title: Extracting the model p-value for a multiple regression in R

text: When fitting multiple variables to one outcome via the `lm()` function in R, `summary(lm)` gives me the p-values for individual regressors but not for the full model in an easily extractable (as in, just accessing fields) kind of way. According to this question, it is possible to extract the p-value via `summary(lm)$fstatistic` by using the command:               pf(x$fstatistic[1],x$fstatistic[2],x$fstatistic[3],lower.tail=FALSE)      However, while in the example linked this provides the same p-value as is printed, I get a different one:               > summary(model)     # ...     Residual standard error: 1.533 on 371 degrees of freedom       (555 observations deleted due to missingness)     Multiple R-squared:  0.3364,    Adjusted R-squared:  0.2864      F-statistic: 6.718 on 28 and 371 DF,  p-value: < 2.2e-16      and:               f = summary(model)$fstatistic     > pf(f[1],f[2],f[3],lower.tail=F)            value      5.948007e-20       What are possible reasons for these values to be different, and which one is the "right" one for the significance of the whole model?


title: How get pvalue with a glmer model (lme4)

text: I am desperately looking a way to obtain the p-values for the fixed effects on my glmer model. Any idea ?? I found a lot of things for lmer model but nothing worked with the glmer ... thank all !


title: Power calculation for a test of trend over time

text: This paper (open access) Babies born at the threshold of viability: changes in survival and workload over 20 years (table 2) presents data for survival of extremely premature babies over the last 20 years in 4-year time periods. They state that the survival of 23-week gestation babies has not improved over this time (p=0.08). Does this test of trend have sufficient power to detect a trend if one existed? How could I determine that?


title: Cubic clustering criterion in R

text: Does anybody know if any package calculates the cubic clustering criterion (CCC) index in R to aid the selection of optimal number of clusters?


title: Repeated Measures Anova in R

text: If I have multiple entries for a certain subject-treatment combination, will the aov() function correctly handle these entries when running a repeated measures ANOVA? For example, suppose the table looks like this, Subject Treatment Time Result   Subj1 Coke Morning 15   Subj1 Coke Morning 14   Subj1 Water Morning 17   Subj1 Water Morning 20   Subj2 Coke Morning 45   Subj2 Coke Morning 46   Subj2 Water Morning 58   Subj2 Water Morning 75   Notice that there are repeated subject-treatment combinations. Will these be correctly handled when the aov(subj1 ~ Treatment*Time, + Error(Subject/(Treatment*Time))); code is run?


title: Max-Diff Using Rank-Order Data?

text: Many discussions of acquiring data for Max-Diff describe the following steps.   * Show all items and ask the respondent to pick the best   * Remove the item identified as the best, and ask respondents to pick the best item from those items remaining   * Repeat until all items but one have been selected Would a standard rank-order question not produce the identical results? Thanks in advance to all for any info. relevant links: http://surveyanalysis.org/wiki/Counts_Analysis_of_Max- Diff_Data http://joelcadwell.blogspot.com/2013/01/if-spss-can-factor-analyze- maxdiff.html


title: Is visual inspection the only way to compare large datasets?

text: I have two large data sets, in fact, one of them is even much larger than the other. Visually, there doesn't seem to be that much difference between them: ![comparing gene targets from different algorithms depending upon whether they have a miRNA inside](http://i.stack.imgur.com/seW9P.jpg) The actual data underlying the box plot isn't normally distributed and doesn't normalise well to transformations. They are roughly the same distribution (i.e. the YES and NO distributions for each algorithm), but the large data size differences make other tests a bit useless. I have applied the Two-sample Kolmogorov-Smirnov test, however this is probably wrong and it gives extremely significant results. My questions are: 1) Do statistical tests on large datasets produce significant results given even slight differences between the two samples? The 'slightness' being magnified given huge data points. 2) Is visual inspection better with large datasets rather than applying non- parametric and parametric tests in which certain underlying assumptions may be violated. 3) For this data, what is the best course of action? Edit My data has structure like : My data is of the form:               Name    Bind    miRNA     a       300     NO     b       500     YES     c       140     YES     d       2345    NO


title: Boosted decision trees using Matlab

text: I would like to experiment with classification problems using boosted decision trees using Matlab. In the paper An Empirical Comparison of Supervised Learning Algorithms this technique ranked #1 with respect to the metrics the authors proposed. My question is, is there a library in Matlab for this type of supervised classification? The function fitensemble(...) has multiple techniques which I'm finding difficult to understand. For example, the 'LSBoost' technique is relevant for regression problems while I'm interested in classification only. Moreover, in a similar question posted this exact question was asked and answered with regards to R, not Matlab.


title: Predicting when & how much is the retailer going to buy

text: I am fairly new to statistical modelling. Currently I have the retailer-wise & month-wise sales data for past 12 months for a particular product. I am trying to predict when the dealer is going to purchase next & approximately how much quantity of the said product. Please suggest an approach/model that will be suitable for achieving the objective. Thanks..


title: What is a second order prior?

text: What is a second order prior distribution in Bayesian statistics? What are some good examples of this concept?


title: Trend in residuals vs dependent - but not in residuals vs fitted

text: I am fitting a linear model to a problem, and a little confused by what is going on. Without the details here are the two plots confusing me: Residuals vs Fitted ![Residuals vs Fitted](http://i.stack.imgur.com/gWcb1.png) Residuals vs Y ![Residuals vs Y](http://i.stack.imgur.com/LsRTw.png) Now the residuals vs fitted looks good to me. Fairly evenly dispersed, no clear pattern. However, the $y$ vs fitted do not look good. I would have expected there to also not be a clear trend in this relationship. This seems basic, but I'm rather lost. I don't want to fit a model that allows large positive errors at high values of $Y$ to make up for large negative errors at low values of $Y$. Is this just an indication of a poorly fitting model? Or is something else going on here? EDIT: Asked for fitted vs y (changed to have labels) ![fitted vs y](http://i.stack.imgur.com/QPl5t.png) EDIT2: Just want to point out that this question has already been asked and answered (albeit in a more abstract sense) What is the expected correlation between residual and the dependent variable? .


title: tight bound of bernoulli sums with unknown dependency

text: Consider n random variables $X_1, \ldots, X_n$ all follow same bernoulli distribution of mean $p$. But the dependency of these variables are unknown (i.e., cannot assume that they are independent). Define $Y=\sum_1^n X_i$, then can I get a lower bound of $\Pr[Y \leq t\cdot n]$ (e.g., I want to get $p_\tau$ such that $\Pr[Y \leq t\cdot n] > p_\tau$) which is tighter than applying the Markov inequality? Thanks.


What do you call an average that does not include outliers? For example if you have a set:               {90,89,92,91,5} avg = 73.4      but excluding the outlier (5) we have               {90,89,92,91(,5)} avg = 90.5      How do you describe this average in statistics?


title: Getting critics to recognize that two similar input patterns refer to the same output-performance relationship

text: The actor-critic model is used within temporal difference learning, which is a method within reinforcement learning, to optimize a process on a state-by- state basis by using the difference between performance and expected performance for each respective state. In other words, there is a set of possible input patterns to the process, each of which constitutes a state (the state in which the agent finds itself). For each state, the actor computes an output, which is evaluated by a performance metric (which can theoretically define a different target output for each state) that returns the reward (performance) signal. Meanwhile, for each state, the critic tracks the expected reward, which is a function of the past reward signals for that state. At each iteration, the difference between the reward signal and the expected reward signal is taken and used to update both the actor's policy (hopefully up a performance gradient) and the critic's performance expectation for that state. Hypothetically, for some problem there exists two or more states (input patterns) that are similar and where the information that differs between the two is irrelevant. In other words, the target output for the states is the same and the input is similar. It would then be potentially useful for the critic to have a method for recognizing that all states that share a certain similarity are effectively the same state, so that the set of those states can share the same predicted outcome. In particular, if you are dealing with a system that features a massive amount of inputs, and if most of those are non- pertinent at any given point in time, then you don't want to relearn the expected performance for each possible combination of non-pertinent inputs. My question, then, is how you do this. How do you get a critic to organize different states into state classes, and in a way that improves performance over the long run?


title: How can I calculate confidence intervals and determine statistical significance for each cell in a contigency table?

text: I'm creating a large contigency table to summarise the difference between a bunch of (categorical) demographics across groups. I would like to know what the confidence intervals are for each cell in the proportional cross-tab, and if the proportion of each demographic group in each sample is statistical different from random. Say I have the following data:               > set.seed(1)     > y <- sample(factor(c("Outcome 1", "Outcome 2", "Outcome 3", "Outcome 4")), 1000, replace = T)     > x1 <- sample(c("Male", "Female"), 1000, replace=T)     > x2 <- sample(c("Employed full-time", "Employed part-time", "Not in labour force", NA), 1000, replace=T)     >      > table(x1, y)             y     x1       Outcome 1 Outcome 2 Outcome 3 Outcome 4       Female       123       126       124       108       Male         121       150       108       140     > round(prop.table(table(x1, y), 2) * 100, 1)             y     x1       Outcome 1 Outcome 2 Outcome 3 Outcome 4       Female      50.4      45.7      53.4      43.5       Male        49.6      54.3      46.6      56.5     >      > table(x2, y)                          y     x2                    Outcome 1 Outcome 2 Outcome 3 Outcome 4       Employed full-time         57        64        65        77       Employed part-time         66        80        57        53       Not in labour force        57        61        51        58     > round(prop.table(table(x2, y), 2) * 100, 1)                          y     x2                    Outcome 1 Outcome 2 Outcome 3 Outcome 4       Employed full-time       31.7      31.2      37.6      41.0       Employed part-time       36.7      39.0      32.9      28.2       Not in labour force      31.7      29.8      29.5      30.9      How would I compute such a table in R? I have read papers where Fisher's exact test is used for this purpose, but it seems to me only to be useful for 2 x 2 tables.


title: Difference between skewed continuous variable and/ or ordinal variable by their binary group allocation

text: I need to select the most appropriate hypothesis test to answer: "whether there is any difference between variable 1 (ordinal) and/or variable 2 (continuous) by group (binary)". This is a randomised study of just over 500 individuals. Variable 1 is an ordinal scale of severity for the first event that occurred to each individual (0 = no event, 1 = mild, 2 = severe, 3= extremely severe). The data is positively skewed with 45 % being 0, 52 % being 1, 2% being 2 & <1% being 3. Variable 2 is the total number of events that each individual experienced over one year (ranges from 0 to 12). 45% of individuals had no event. There is a strong positive skew (0 = 241, 1= 120, 2 = 84, 3 = 33, 4=18, 5= 9...etc until only 2 individuals had 12 events). The individuals are randomised to one of two groups (placebo vs intervention) I have so far done two separate Mann-Whitney U tests (one each for variable 1 and variable 2 - each by group i.e.placebo and intervention). I have also used SPSS to calculate Kendall's tau and Spearman's rho correlation coefficients, which are negative. I have little experience of using these but I understand that Spearman's is ok for continuous data such as variable 2 and and Kendall's tau is ok for ordinal data such as variable 1. However I am not sure that it is appropriate to use this correlation to compare variables 1 or 2 with the binary grouping variable? Does anyone have any suggestions on whether this seems appropriate or is there a better way of doing this?


title: Does anyone know how to make this infograph with R?

text: This picture was posted on a Brazilian newspaper called "Folha de São Paulo". The graph relates each candidate's name with their allied political parties ![enter image description here](http://i.stack.imgur.com/xqV4m.jpg)


title: Multiple comparison on subgroups and overall comparison

text: I am trying to work out is there is any association between occurrence or not of an event in around 500 individuals and their randomised binary grouping (intervention vs placebo). However I also want to see if this holds true when broken down by one of 5 locations in which the individuals were treated. A chi square test on the entire group shows a significant result. However if I repeat the test on each of the 5 subgroups only one subgroup is significant (<0.05) and all of the others are non-signficant. Further still I have then corrected for multiple testing by multiplying the p-values by 6 (5 tests for each centre and then 1 test for entire group). This does not affect the subgroups but does affect the overall group which is no longer significant after multiple testing correction. I know this issue is controversial. Had I just stuck with the original hypothesis and performed one chi square test I would be saying it is significant but once I do subgroup testing I am then having to say this is not the case. Does anyone have any suggestions on how to approach this problem?


title: What is the Probability Distribution of NLTK Naive Bayes?

text: As I know Naïve Bayes has various distributions, as said in Sci-kit learn manual: > The different naive Bayes classifiers differ mainly by the assumptions they > make regarding the distribution of $P(x_i \mid y)$. So, there are various methods of Naïve Bayes like, Gaussian NB, Multinomial Naïve Bayes, and I found Sci-kit learn has different libraries for them. There is another famous Naïve Bayes library in Python named NLTK. I was trying to find out the distribution of NLTK based Naïve Bayes. After searching I felt:   1. It is not Multinomial (see here, and here)   2. But in the documentation, I am finding the following examples, Example 6.4, Example6.6, it seems NLTK’s Naïve Bayes Classifier is supporting, `nltk.FreqDist`.  So I am bit confused. If any one of the esteemed members may kindly guide me.


title: Limits of zero-inflated negative binomial in % of zeroes

text: What are the limits of a zero-inflated regression? Specifically, if more than 80% of the data is zeroes, is the ZINB still valid? What is a good rule of thumb or educated way to understand the right threshold for the number of zeroes?


title: Rule of thumb when drawing N samples from a discrete distribution with N possible values with replacement

text: I'm looking for an explanation and possibly the name of a rule of thumb: When drawing N samples with replacement from a discrete uniform distribution of N values, it is very likely that:   * 1/3 of the values are not drawn   * 1/3 of the values are drawn exactly once   * 1/3 of the values are drawn more than once The way I heard that was with throwing 64 pieces of rice onto a chess board and then observing roughly what I described above. I'm seeing this very often, but I want to read up on the theory and calculations behind it. Why is this rule of thumb valid and does this rule of thumb have a name? And if yes, which one?


title: Hidden Markov model for event prediction

text: **Question** : _Is the set-up below a sensible implementation of a Hidden Markov model?_ I have a data set of `108,000` observations (taken over the course of 100 days) and approximately `2000` events throughout the whole observation time- span. The data looks like the figure below where the observed variable can take 3 discrete values $[1,2,3]$ and the red columns highlight event times, i.e. $t_E$'s: ![enter image description here](http://i.stack.imgur.com/QkIn0.png) As shown with red rectangles in the figure, I have dissected {$t_E$ to $t_{E-5}$} for each event, effectively treating these as "pre-event windows". **HMM Training:** I plan to train a Hidden Markov Model (HMM) based on all "pre-event windows", using the multiple observation sequences methodology as suggested on Pg. 273 of Rabiner's paper. Hopefully, this will allow me to train an HMM that captures the sequence patterns which lead to an event. **HMM Prediction:** Then I plan to use this HMM to predict $log[P(Observations|HMM)]$ on a new day, where $Observations$ will be a sliding window vector, updated in real-time to contain the observations between the current time $t$ and $t-5$ as the day goes on. I expect to see $log[P(Observations|HMM)]$ increase for $Observations$ that resemble the "pre-event windows". This should in effect allow me to predict the events before they happen.


title: Number of samples vs Number of features

text: I've got a set of two classes with 4000 observations total. I've a set of 63 features to construct a predictor. My question is, is there a relation that would prevent overfiting for having too much dimensions? I mean a maximum number of features per number of observations in training?


title: Software to produce confidence interval error bars from summary statistics without access to the original data

text: A bit like a box plot. I mean not necessarily the standard upper confidence interval, lower confidence interval, mean, and data range-showing box plots but I mean like a box plot with just the three pieces of data I have specified in my question. This is a screenshot of a journal article which had exactly what I want: ![enter image description here](http://i.stack.imgur.com/mAuzT.png) I would also like to know how I would use the software the answerer mentions to create such a plot. Something I said in prior versions of this question that has since been edited out is **I ONLY KNOW THE 95% CI INTERVAL AND MEAN; NO OTHER DATA!**


title: Using an odds ratio when data is sparse

text: Suppose I have around 20 exposures that potentially affect an outcome and I want to see which exposures have bigger impacts on the outcome. So I want to calculate each exposures' odds ratios by exponentiating the coefficients obtained from logistic regression. So I have the following input set and output set where 1 means it (exposure or outcome) is present and 0=not present: ![enter image description here](http://i.stack.imgur.com/4GTNy.png) So, for example, the first row represents a sample where exposure 1 wasn't present, exposure 2 was present,...exposure 20 was present and the outcome was present. I fit a logistic regression model to this data and exponentiate the coefficients to get odds ratios. The potential problem is that I am going to be working with a VERY sparse data set with many samples. There are many instances where almost all exposures except one or maybe two is going to be present in a sample. My question is if this sparsity is something to be concerned about and if this will make my method of comparing exposures using odds ratios a bad idea. Page 6 of this paper Greenland 1987 seems to imply that sparsity won't matter too much but I want to see what the statisticians here say. Any links to papers would be appreciated.


title: What precisely is the interpretation of distance correlation?

text: The key feature of distance correlation is that if it is $0$ then two variables are independent. If the pearson correlation is independent, however, this does not imply independence. How would one interpret distance correlation?


title: Understanding the derivation of an equation in LDA modeling

text: When reading the derivation of LDA models, I usually get the following equations. ![enter image description here](http://i.stack.imgur.com/hZYho.jpg) I do not quite understand the second step, where $p(\mathbf{z}_{-i},\mathbf{w}|\alpha,\beta)$ was removed. Is that because it can be thought of as a constant when studying $z_{i}$. I know this is a common practice in Bayesian statistics, but I am not sure whether $p(\mathbf{z}_{-i},\mathbf{w}|\alpha,\beta)$ is a constant.


title: ARIMA forecast with seasonality and trend, strange result

text: as I am stepping into forecasting with ARIMA models, I am trying to understand how I can improve a forecast based on ARIMA fit with seasonality and drift. My data is the following time series ( over 3 years, with clear trend upwards and visible seasonality, which seems to be not supported by autocorrelation at lags 12, 24, 36??).                   > bal2sum3years.ts                  Jan     Feb     Mar     Apr     May     Jun     Jul     Aug                   2010 2540346 2139440 2218652 2176167 2287778 1861061 2000102 2560729          2011 3119573 2704986 2594432 2362869 2509506 2434504 2680088 2689888          2012 3619060 3204588 2800260 2973428 2737696 2744716 3043868 2867416                   Sep     Oct     Nov     Dec         2010 2232261 2394644 2468479 2816287         2011 2480940 2699780 2760268 3206372         2012 2951516 3119176 3032960 3738256      The model that was suggested by `auto.arima(bal2sum3years.ts)` gave me the following model:                   Series: bal2sum3years.ts          ARIMA(0,0,0)(0,1,0)[12] with drift                       Coefficients:                   drift               31725.567         s.e.   2651.693              sigma^2 estimated as 2.43e+10:  log likelihood=-321.02         AIC=646.04   AICc=646.61   BIC=648.39      However, the `acf(bal2sum3years.ts,max.lag=35)` does not show acf coefficients higher than 0.3. The seasonality of the data is, however, pretty obvious - spike at the beginning of every year. This is what the series looks like on the graph: ![Original Time Series](http://i.stack.imgur.com/kQi5N.png) The forecast using `fit=Arima(bal2sum3years.ts,seasonal=list(order=c(0,1,0),period=12),include.drift=TRUE)` , called by function `forecast(fit)`, results in the next 12months's means being equal to the last 12 months of the data plus constant. This can be seen by calling `plot(forecast(fit))`, ![Actual and Forecasted Data](http://i.stack.imgur.com/GJqcG.png) I have also checked the residuals, which are not autocorrelated but have positive mean ( non zero). The fit does not model the original time series precisely, in my opinion ( blue the original time series, red is the `fitted(fit)`: ![Original vs fit](http://i.stack.imgur.com/ux3i7.png) The guestion is, is the model incorrect? Am I missing something? How can I improve the model? It seems that the model literally takes the last 12 months and adds a constant to achieve the next 12 months. I am a relative beginner in time series forecasting models and statistics. Thank you very much for your answers!


title: How to handle negative Y values with exponential regression?

text: I have data points that include one instance of a negative y value. With this data, is it possible to generate a exponential decay function in a way that is mathematically sound? The only variable is X.


title: Regression using dummy variables

text: I am working on a credit scoring modelling project and we decided to use dummy variables for regression. The way we create dummy variables are: For each predicting numeric variable,   1. We create by default 10 equal-size bins, examine the weights of evidence (WOE) for each bin and merge adjacent bins if their WOEs are similar (meaning their risk signals are similar.   2. Then we create dummy variables fomr this numeric variable according to the final bins. Say if we end up with 3 bins for this variable, then we create 3 dummy variables representing the bins. Then we pick one dummy variable as reference to avoid perfect multicollinearity.   3. Then we run logistic regression. The above is the highlevel description of how we do the modelling. My question is not specific to this procedure but more of a general dummy variable regression question: In the case when we have many variables (~20) for regression, the number of dummy variables will be even more and the regression result often says some dummy variables are insignificant. How do you treat insignificant coefficient estimates? What if a variable is significant in variable level, having some significant dummy variables but not all? Thank you.


title: Probability mass function of the sample mean

text: Suppose that $X_1,X_2,\ldots,X_n$ is a random sample from a distribution with probability function: $$p_X(x)=\begin{cases} 1/2 &x=-1,1 \\\ 0 & \text{otherwise} \end{cases}$$ Now if we define the average $\bar{X}= \frac{1}{n} \sum X_i$, I need to show that that for an odd $n$ the probability function for $\bar{X}$ is: $$p_{\bar{X}} (x)= \frac{\binom{n}{\frac{n}{2} (x+1)}}{2^n}$$ for $x=\pm 1/n,\pm 3/n,\ldots,\pm 1$, $0$ otherwise. * * * Since $P[\sum X_i=k] =P [\bar{X}=k/n ]$, I thought it would be easier to first derive the probability function for the sum. Since they are still $2^n$ n-tuples and each one is equilikely, I need to count the ones whose sum is $k$ for: $k=\pm 1,\pm 3,\ldots \pm n-2,\pm n$ (since even numbers are not a feasible combination) The problem here though is that I do not immediately recognize how I could count all these tuples. Could you please help me with this? If it is easier to proceed another way, I am of course all ears.


How do you compute confidence intervals for positive predictive value? The standard error is: $$SE = \sqrt{ \frac{PPV(1-PPV)}{TP+FP}} $$ Is that right? (here my concern is the denominator) Does that formula work for any similar ratio in a 2x2 table. E.g. for sensitivity, it would be $$SE = \sqrt{ \frac{SENS(1-SENS)}{FP+TN}} $$ Is that right? (here my concern is that it is generalizable to other ratios as long as you get the denominator right) And the for the 95% confidence intervals: $$CI_{PPV} = PPV \pm 1.96*SE$$ Is that right? (my concern here is how to go from SE to the confidence interval) (of course with all the cell restrictions like $n\cdot p\cdot (1-p) \ge 5$)


title: Estimating event probability from historical time series with clear seasonality

text: I would like to predict the average number of days in a year for which two conditions are true:   * daily average temperature is below zero celsius   * the day was preceded by at least four days with daily average temperature below zero celsius I've historical daily average temperature data for the location available for about 10 years. My initial approach was to use the one sided Chebyshev inequality which can be used to approximate a probability if the distribution is not known. However in this application I am interested in the probability of a special condition, can I use the Chebyshev inequality for a dummy time series as well? I.e.: 1 if condition is fullfilled, otherwise 0, --> the dataset would therefore look something like 0,0,0,0,0,1,1,0,0,0,0,0,1,1,1,1,1,1,etc. How would you approach a problem like that from a different angle, the data clearly has seasonality, is there any distribution which I could use to have a better estimate than Chebyshev?


title: How to carry out a Likert scale analysis?

text: ### Context: My research aims to assess whether parents who have been involved in a relationship with a social worker are at higher risk of child abuse/neglect. I am trying to establish whether there exists a relationship between social workers' attitude toward the parents that reflect an overall outcome of the case (positive or negative). I want to assess whether social worker attitudes have an effect on parental responses.   * The questionnaire was only given to parents who have a relationship with a social worker.    * There were 20 participants (16 females and 4 males).   * I aggregated my 13 Likert items/questions into a single attitude (word) per each.    * I also reversed the order of the Likert scale when I entered the data, so that a score of 5 stands for "strongly disagree", whereas a score of 1 reads "strongly agree". This was done because of the way my Likert items were asked. ### Question:   * Can I run a test between genders given that there are an uneven number of males and females?   * What analyses can I run to test my research question?


title: Solutions to E. T. Jaynes Probability Theory for self study?

text: I am doing a self study of E.T.Jayne's Probability Theory, the logic of science. I wanted to see that does any one know of a web page providing solutions to the exercises? I have no teacher or anyone experienced in this book around me, I need to make sure that my solutions to them are correct.


title: How can I predict a value like this from a bayesian network model?

text: Let $\theta$ be the set of parameters for a Bayesian Network model. The model observes some data $D$. I want to then predict the values of some new data $D$. Why is the following true? $$ p(D'|D) = \int p(D'|\theta)p(\theta|D) d\theta $$ Does it stem from some manipulation of Bayes rule? Or the formula for expected value?


title: What does "degree of freedom" mean in neural networks?

text: In Bishop's book "Pattern Classification and Machine Learning", it describes a technique for regularization in the context of neural networks. However, I don't understand a paragraph describing that during the training process, the number of degrees of freedom increases along with the model complexity. The relevant quote is the following: > An alternative to regularization as a way of controlling the effective > complexity of a network is the procedure of early stopping. The training of > nonlinear network models corresponds to an iterative reduction of the error > function defined with respect to a set of training data. For many of the > optimization algorithms used for network training, such as conjugate > gradients, the error is a nonincreasing function of the iteration index. > However, the error measured with respect to independent data, generally > called a validation set, often shows a decrease at first, followed by an > increase as the network starts to over-fit. Training can therefore be > stopped at the point of smallest error with respect to the validation data > set, as indicated in Figure 5.12, in order to obtain a network having good > generalization performance. **The behaviour of the network in this case is > sometimes explained qualitatively in terms of the effective number of > degrees of freedom in the network, in which this number starts out small and > then to grows during the training process, corresponding to a steady > increase in the effective complexity of the model.** It also says that the number of parameters grows during the course of training. I was assuming that by "parameters", it refers to the number of weights controlled by the network's hidden units. Maybe I'm wrong because the weights are prevented to increase in magnitude by the regularization process but they don't change in number. Could it be referring to the process of finding a good number of hidden units? What's a degree of freedom in a neural network? What parameters increase during training?


title: Finding the best fit parameters for only a certain segment of data

text: I have some experimental data that are observations of an $n$ dimensional quantity, $x = \\{x_1, x_2, \ldots, x_n\\}$. I also have a complicated parameterized model $f(a,b)$ for the data, which for a given choice of $a$ and $b$, will spit out an $n$ dimensional vector $y = \\{y_1, y_2, \ldots,y_n\\}$. $y$ is supposed to be "similar" to $x$ for certain segment with indices $1< i, i+1,\ldots,j<n$ ($i$ and $j$ are known). The objective is to find the $a$ and $b$ that best describes $x$ (more on how after the list). Since my actual application is complicated and time consuming I won't be able to share the data, but I'll mention some specifics and provide sample data:   * I don't have an analytic form for $f$; it's a complex numerical simulation that just spits out the list given $a$ and $b$.   * The underlying distribution for $x$ is not known.   * All $x_i$ and $y_i$ are non-negative Here's some sample data ( _Mathematica_ code, but feel free to use any language you like):               test = Sort@RandomReal[ChiSquareDistribution[1], {1000}];     f[a_, b_] := a Sort@RandomReal[ChiSquareDistribution[b], {1000}]^2;      Now if we look at the plots of the test data and the model for some values of `a` and `b` that I eyeballed to look good:               ListLinePlot[Log10@{test, f[0.064, 3.3]}]      ![Image 1](http://i.stack.imgur.com/KvuZL.png) you see what I mean by "similar" for a certain segment (here indices between 300 to 900). Compare this with an arbitrary choice of `a` and `b` that are "not similar"               ListLinePlot[Log10@{test, f[0.4, 1.3]}]      ![Image 2](http://i.stack.imgur.com/cXjUY.png) Assuming that this similarity segment 300 – 900 is known, how would I select the right values for a and b? My idea was to discretize the space spanned by `a` and `b` (sensible values for both are known _a priori_ , and for the above, if necessary, you can pick $a\in[0.05, 1], b\in[2.5,4]$) and try to find the combination that would minimize $\Vert \\{x_i,\ldots,x_j\\} - \\{y_i,\ldots,y_j\\}\Vert_p$, where $p$ is either 1 or 2 depending on how I want to weight the outliers. Trying this out,               res = With[{t = Take[test, {300, 900}]},          Table[Norm[t - Take[f[a, b], {300, 900}], 2],              {a, 0.5, 1, 0.01}, {b, 2.5, 4, 0.01}]];      unfortunately gives something like this ![Image 3](http://i.stack.imgur.com/AuOqj.png) which does not indicate a clear minimum. Similar result for the $L_1$ norm as well... I also tried taking the $\log_{10}$ of each and then finding the error (is that even valid?), but that wasn't helpful either. I've seem to have hit a wall on finding a way to identify the right (or close to the right) parameters. Would it be possible to, instead, use the slope and curvature of the segment to find the parameters instead of a point-by-point error? I'd appreciate any pointers towards a possible solution. I can do the programming myself, and can understand logic in most common programming languages, so that's not a barrier. > **NOTE:** My data is _not_ $\chi^2$ distributed, nor is the model. This is > just an example that I picked, since I could illustrate how the shape is the > same in the segment, but not elsewhere.


title: Tests to establish correlation

text: I'm trying to analyze if there exists any correlation between two time series? So far I have not been able to notice any significant correlation. So if I were to assert in my report that these time series are uncorrelated what are the minimum experiments/tests/metrics/results that I must provide to support my claim?


title: Transforming / Evaluating (Probably) Normal Data from Sample with Linear Probability Distribution

text: I have a sample of data for which the count of a given value increases roughly in proportion to the actual value. For example, the value 22 will occur around 22 times and the value 30 will occur around 30 times. If I recall this is what is described by a linear probability model. However, I am almost certain that the underlying distribution is normal, and that the issue is with the sample that has been provided. Unfortunately, there is no easy way to change this sample, and any suggestions that involve a new sample should be considered a last resort. Is there a way to use intuition about the data's underlying distribution to model the sample? Are there known methods to handle this sort of problem?


title: Finding the corresponding bayesian network of a predefined joint probability distribution

text: Given a joint probability distribution over the variables $X_1,X_2,\dots,X_n$. Is there an algorithm for constructing the corresponding Bayesian Network?


title: Algorithm for real-time normalization of time-series data?

text: I'm working on an algorithm that takes in a vector of the most recent data point from a number of sensor streams and compares the euclidean distance to previous vectors. The problem is that the different data streams are from completely different sensors, so taking a simple euclidean distance will dramatically overemphasize some values. Clearly, I need some way to normalize the data. However, since the algorithm is designed to run in real time, I can't use any information about any data-stream as a whole in the normalization. So far I've just been keeping track of the largest value seen for each sensor in the start-up phase (the first 500 data vectors) and then dividing all future data from that sensor by that value. This is working surprisingly well, but feels very inelegant. I haven't had much luck finding a pre-existing algorithm for this, but perhaps I'm just not looking in the right places. Does anyone know of one? Or have any ideas? I saw one suggestion to use a running mean (probably calculated by Wellford's algorithm), but that if I did that then multiple readings of the same value wouldn't show up as being the same, which seems like a pretty big problem, unless I'm missing something. Any thoughts are appreciated! Thanks!


title: A measure to identify discriminating subgraphs of a larger graph?

text: I have genes from four species of bacteria. I constructed an undirected graph from the sequence similarity between all the genes from these bacteria. Each node in the graph is a gene and each edge is the similarity between two genes. I partitioned this graph into subgraph clusters of similar genes using the MCL software. The four species of bacteria I am using are split into two groups (A & B) based on their response to an environmental condition. I would like to identify the clusters of genes from my graph that are producing this response. To put it another way I would like to create a measure to find the clusters which discriminate A from B. Generally each cluster has genes from all four species. What would be the best way to do this? One idea I had was to determine the mean weight between all A-type nodes in a cluster and then divide by the mean edge weight of the cluster. Values greater than 1 would indicate the A-type genes are more closely related to each other.


title: Which test to see if two set of continuous data are comparable

text: Scenario: Two different persons have measured lesions in millimeters before and after medicine, so the data are continuous. I have calculated both of these readers' difference in measurements and want to use the mean difference for further calculations. To do this I think I have to test whether or not these measurements are not too different? How do I do that? I thought about a paired t-test since the measurements have been done on the same patients ... but I have the feeling it is not right. Ideas?


title: Reliability engineering: Measuring the proportion of defective products at several time points using the same test units

text: A product being designed will have a failure-free warranty requirement of 1 year and the Engineers are planning a series of tests to estimate our warranty exposure. Based on the binomial confidence interval analysis procedure described in section 7.2.4.1 of the NIST/SEMATECH e-Handbook of Statistical Methods I have computed the following in advance of the test. I used the exact intervals for small numbers of failures.   * If out of a random sample of 20 units, 4 units fail the 1 year warranty simulation test, I believe we can say that we are 90% confident that the population defect rate will lie between 7.1% and 40.1%    * If out of a random sample of 20 units, 1 unit fails the 1 year warranty simulation test, I believe we can say that we are 90% confident that the population defect rate will lie between 0.26% and 21.6%    * If out of a random sample of 20 units, 0 units fail the 1 year warranty simulation test, I believe we can say that we are 90% confident that the population defect rate will be 13.9% or less.   * If out of a random sample of 1 unit, zero units fail the 1 year warranty simulation test, I believe we can say that we are 90% confident that the population defect rate will be 95% or less. Now assuming no failure occurs at the 1 year point in the test(s), continuing the test(s) beyond the 1 year warranty period with a some corresponding benefit in confidence would be desirable as the test units are costly to build. Can the above analysis procedure be extended to answer the following questions when the simulation test time exceeds the actual warranty period?   1. At the 90% confidence level, what are the upper and lower limits of the population defect rate if a 1.5 year warranty simulation test is run on 1 unit and that unit does not fail?   2. At the 90% confidence level, what are the limits of the population defect rate if 2 units are warranty simulated. One unit fails at 1.5 years and the other fails at 4 years?   3. Can Monte Carlo be used here to help answer/validate these questions? Many thanks for any guidance.


title: is there any difference between taking more samples and a sample with more observations?

text: I'm quite confused with the difference between taking the average of more samples and taking the average of a sample with more observations. Do you get unbiased estimates by taking more sample of more observations?


title: Skewed but bell-shaped still considered as normal distribution for ANOVA?

text: This could be a pretty basic question, I'm a little rusty on my stats knowledge. **Background:** I am monitoring website load time performance. To do so, I have a script running and capturing data points (About 400) on load time through various Agents. Every Agent is located in different geographic locations, but they measure the same steps. I would like to determine if there is statistical difference between the agents. So if one is consistently reporting slower load time performance I would like to know if its because of the Agent or not. I would include images but I need 10 reputation points and I just found out about this website. **Problem:** I have two sets of data from different agents measuring the seconds it takes a website to download, both are bell-shaped but are heavily skewed to the right. Can I still perform ANOVA to determine if there is difference, even though they are skewed? Thanks in advance


title: Performance evaluation of a model: MAPE in R

text: I'm trying to compare two models. I have got the root mean square error (RMSE) but it's not satisfying. This is why I'm trying to work out the mean absolute percent error (MAPE). Is there any package in R where I can compute the mean absolute percent error?


title: Identity in Simple Linear Model

text: I'm working on this identity $$\sum_{i-1}^n (y_i - \hat {\beta_0} - \hat {\beta_1}x_i)^2 = \sum_{i=1}^n y_i^2 - \hat {\beta_0}\sum_{i=1}^n y_i - \hat {\beta_1} \sum_{i=1}^n x_iy_i$$ I have these relationships to work with: $$\hat {\beta_1} = \frac { n\sum_{i_1}^n x_iy_i - \left ( \sum_{i=1}^n x_i \right ) \left ( \sum_{i=1}^n y_i \right )}{ n \left ( \sum_{i=1}^n x_i^2 \right ) -\left (\sum_{i=1}^n x_i \right )^2 }$$ $$ \hat {\beta_0}= \overline {y} - \hat {\beta_1} \overline {x}$$ A little manipulation also shows $$\hat {\beta_1}= \frac { \sum_{i=1}^n ( x_i - \overline {x}) y_i}{ \left ( \sum_{i=1}^n x_i^2 \right ) -n \overline {x}^2 }$$ My strategy has to substitute $\hat {\beta_1}$ out of the equation, but I just keep getting $\overline {y}$ in too many terms and not any terms with $y_i$. The formula for $\hat {\beta_1}$ is too complicated to consider working with. What am I missing?


title: Is my understanding of this “confidence interval” solution correct?

text: I was faced with a problem of automatically detecting regions of continuity in a vector. I have a lot of these vectors, hundreds. One example of such vector is here. Basically looking at the vector, one can see that there are regions of continuity, then it breaks, jumps some value, and then has another region of continuity and so on. I was provided with this solution by some nice fellow on the newsreader:               v = [ 9 18 21 58 59 60 63 66 69 70 72 74 ...     dv = diff([0 v]);                                       % Create difference vector     dvs = [mean(dv)  std(dv)];                              % Determine mean & std     dvd = [0 find(dv > 1.96*dvs(2)+dvs(1)) length(v)];      % Use ‘dvs’ to detect discontinuities & create index reference vector     for k1 = 1:length(dvd)-1         vs{k1} = v(dvd(k1)+1:dvd(k1+1)-1);                  % Create cell array of regions-of-interest         vi{k1} = [dvd(k1)+1 dvd(k1+1)-1];                   % Create reference array of start-end indices for ‘vs’ regions     end      this works wonderfully and uses the 95% confidence interval to find the start and end points of the indexes of the regions. I know it works; I just want to know why and how it works. Statistics is not my strong point, so I have been reading in detail anything I can find related to the provided solution. Statistics is really fascinating me as to how one can make assumptions intelligently based on a sample and find them to be true. I want to fully understand the working of the solution provided. Based on what I have read and understood, this is my understanding of the explanation (please note I have no base in stats at all, so I may be going very basic here). Can anyone strong in statistics shed some light on this? A 95% confidence interval is basically saying that I am 95% sure that given mean of a probability is within these limits [a, b]. And this can be calculated with the formula: $x \pm 1.96 (std/\sqrt n)$ where $x$ is the mean, $std$ is the standard deviation and $n$ is the sample size. And this confidence interval applies to a normal distribution only. Judging by my vector, the samples were observed to be normally distributed samples (between the jumps). Knowing that these are normal distributions, simply finding their confidence intervals will give the start and endpoint of a normally distributed sample. And that is what is going on in this line:               dvd = [0 find(dv > 1.96*dvs(2)+dvs(1)) length(v)];      where `find` finds the indices of the points where the confidence interval equation holds true, but why the difference vector `dv`? and why use it in the condition statement? Why even create a difference vector? And why use it in the line above? This is based on what I read today, so I may be way off. Am I on the right path? Is my understanding correct? Thank you.


title: How do you calculate eigenvector/value of a non-square matrix in R?

text: I have been wanting to calculate eigenvalues and eigenvectors for a non-square matrix and I know that `svd` method is used. But, given my poor background, I don't understand how to figure out eigenvalues and eigenvectors from the `u`, `d` and `v` matrices.


title: Spectral clustering using RBF Kernel function in R

text: I have extracted user-features and item features in my recommender system using a modified SVD approach built on ALSE (loosely based on Yehuda Koren's paper). I now want to cluster items not directly based on distance/cosine similarity measure alone but learning from the extracted features and capturing latent conditions which may turn out to useful in clustering. I am trying out the Kernel based Spectral Clustering method as part of the R package 'kernlab' using the specc method. I have 2 doubts: 1\. How do I choose the 'k' (equal to both the reduced dimensional space as well as the # of clusters)? For a given choice of 'k', sigma is optimized for minimizing the within cluster SOS. But in a broader sense, what should I do for choosing the # of clusters? People suggest eigen gap method (ie based on difference between consecutive eigen values. But for that you need a value of sigma to get the eigen vectors. So what value of sigma I give decides the difference. I am confused on what to do).   1. Is there a way to update the model incrementally rather than re running the clustering algorithm from scratch in case new items are added? Would be grateful if somebody could reply quickly. Thanks.


title: Elastic net package for mixed effects models?

text: I know about `glmmLasso` but would prefer to use elastic net. I wonder if there are any glmm analogues of `glmnet` out there, or if anyone is working on any. I couldn't find any in my searches. Thanks


title: Covariates in regression models

text: Should covariates be included in regression analyses if they are correlated with the dependent variable or if they are correlated with the predictor variable/s. Alternatively, should they be included as a result of past (fairly robust) findings that they are significantly related to the outcome and/or predictor variables?


In the book “OpenIntro Statistics” it is specified that samples of voluntary patients cannot be generalized. Why cannot a sample of voluntary patients be generalized?


title: Nonparametric test for largely skewed count data

text: My research design looks as follows: an experimental game with 4 participants (human subjects), repeated for 20 rounds. During each round, participants are allowed to form bilateral coalitions which are “frozen” until the end of a round. Hence, six different possible coalitions exist, but maximum two coalitions can emerge during a single round. The same setup is tested through a simulation, where simulation agents follow pre-set decision rules to form coalitions, the difference from the experiment being that all simulation rounds are independent. For the experiment and for the simulation, I have the data on the number of times each particular coalition was formed. It can be visually represented on a bar chart, on which each colored bar corresponds to each of the possible coalitions, and the vertical axis shows the number of times it occurred (count variables). ![enter image description here](http://i.stack.imgur.com/nhKTZ.png) My general idea is to understand whether stable coalition relationships emerge between participants in the experiment as they play repeatedly: for example, whether participant 1 tends to form coalitions with participant 2 more often than with other players. To that end, I would like to test:   1. Whether there are statistical differences between the frequency of occurrences of different coalitions;   2. Whether the experimentally observed frequency distribution statistically differs from a theoretical benchmark (each coalition occurs with the equal 1/6=0.15 probability). What would be the appropriate test in this case given that I am working count data which is (a) not normally distributed, (b) is largely skewed (a numeric outcome is frequently zero)?


title: How to deal with Z-score greater than 3?

text: In a standard normal distribution how do I deal with a $Z$ value greater than 3? I know that z-score ranges form -3 to 3 Consider this one ... mean = 70, standard deviation = 4 I need to find $P(65 < X < 85)$. Transforming to standard normal gives $P(-1.25 < Z < 3.75)$ How to deal with the $3.75$? Edit *Actually it's not greater than 3, $z < 3.75$. I meant smaller than, sorry. Should I assume it's just 0.4990 or what?


title: Dividing up Training Data into Test Set

text: Suppose we have a training data set. We want to learn some hypothesis using some algorithm. Would we divide up the training set differently if we used, for example, logistic regression as opposed to support vector machines? So if we divide the training data set into: $70 \%$ in-sample data and $30 \%$ out-of-sample data, would this hold if we used logistic regression or support vector machines?


title: Bias of estimating a Ratio?

text: I would like to know why 1/N * Summation of(y/x) is a worse estimator than average of y divided by average of x. This is in the context of constructing an estimator of a ratio where the ratio would be something like cost per person etc. I'm drawing from sample data as well. Thanks in advance,


title: negative parameter estimate but positive relationship in gamm

text: I am trying to figure out how to interpret results obtained when estimating GAMM using mgcv. Consider this hypothetical model: M <\- gamm(Y ~ offset(nn) + s(X, by = Z) + Z, random = list(id = ~ 1)) where: Y - dependent variable X - predictor nn - variable included as an offest Z - factor covariate included as interaction with X, then as a parametric term, as well id - factor included as an random term Now considering output graphs, Y and predictor X show nice positive relationship. However, when looking at estimated parameters in lme object, the estimate b for X is negative.. I'm confused, because i always thought, that when b is negative for a given variable, it should be rather obvious, that this variable affect Y negatively. Any ideas why graphs show positive relationship between Y and X while estimate b for X is negative? Cheers!


title: Creating semi-random planting designs in R

text: I want to use R to create 'random' planting plans - drawing plants manually (CAD) is neither random nor efficient - my aim is a list of plants, within an irregular region, (up to 10,000 plants) with their respective numbers and positions. This would use points (trees) of different diameters, with a coverage ratio for each tree, and a predefined level of how the points attract or repel as in attached image. ![enter image description here](http://i.stack.imgur.com/ksrTc.jpg)True randomness would not occur as the points would not be truly independent. The maximum plant coverage (range 50% to 100%) would probably be a circle- packing problem, but where packing is random - so maximum number of circles that fit the region while being as random as possible, the result of this would be 100% coverage. My (unstatistical) mind tells me this number could inform the number of circles for subsequent runs of the model for coverages <100%. Before I take this further is this do-able in R, and if so what steps should I take?


title: Testing main effects and interactions with rma from metafor

text: I'm not sure how to test for main effects and interaction effects between multiple moderator variables in metafor. Am I doing this right? My dataset is quite large, so we'll play with the iris dataset. Create dataset               sp <- gsub("setosa", 0, iris$Species)     sp <- gsub("versicolor", 1, sp)     iris$Species <- as.numeric(gsub("virginica", 2, sp))     iris$var <- rep(c(0,1,2), nrow(iris)/3)      I think this is correctly testing for main and interaction effects between the two moderator variables `Species` and `var`               rma(Sepal.Length, Petal.Length, mods = ~ Species * var -1, data = iris)               Mixed-Effects Model (k = 150; tau^2 estimator: REML)          tau^2 (estimated amount of residual heterogeneity):     3.8158 (SE = 0.8043)     tau (square root of estimated tau^2 value):             1.9534     I^2 (residual heterogeneity / unaccounted variability): 58.63%     H^2 (unaccounted variability / sampling variability):   2.42          Test for Residual Heterogeneity:      QE(df = 147) = 387.0300, p-val < .0001          Test of Moderators (coefficient(s) 1,2,3):      QM(df = 3) = 574.7435, p-val < .0001          Model Results:                       estimate      se     zval    pval    ci.lb    ci.ub          Species        3.8987  0.3038  12.8333  <.0001   3.3032   4.4941  ***     var            3.0421  0.2400  12.6769  <.0001   2.5718   3.5124  ***     Species:var   -1.8851  0.2761  -6.8284  <.0001  -2.4262  -1.3440  ***          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       And this by adding `factor()` call to each moderator is simply getting estimates and CI's for each level of each moderator               rma(Sepal.Length, Petal.Length, mods = ~ factor(Species) * factor(var) -1, data = iris)          Mixed-Effects Model (k = 150; tau^2 estimator: REML)          tau^2 (estimated amount of residual heterogeneity):     0 (SE = 0.2714)     tau (square root of estimated tau^2 value):             0     I^2 (residual heterogeneity / unaccounted variability): 0.00%     H^2 (unaccounted variability / sampling variability):   1.00          Test for Residual Heterogeneity:      QE(df = 141) = 10.7649, p-val = 1.0000          Test of Moderators (coefficient(s) 1,2,3,4,5,6,7,8,9):      QM(df = 9) = 1666.4569, p-val < .0001          Model Results:                                         estimate      se     zval    pval    ci.lb   ci.ub          factor(Species)0                 5.0488  0.2947  17.1332  <.0001   4.4712  5.6263  ***     factor(Species)1                 5.7185  0.4935  11.5867  <.0001   4.7512  6.6858  ***     factor(Species)2                 6.7073  0.5870  11.4267  <.0001   5.5568  7.8577  ***     factor(var)1                    -0.0596  0.4106  -0.1452  0.8845  -0.8645  0.7452          factor(var)2                    -0.1082  0.4221  -0.2562  0.7978  -0.9354  0.7191          factor(Species)1:factor(var)1    0.3329  0.8248   0.4036  0.6865  -1.2838  1.9495          factor(Species)2:factor(var)1   -0.2557  0.9126  -0.2802  0.7793  -2.0445  1.5330          factor(Species)1:factor(var)2    0.3598  0.8174   0.4402  0.6598  -1.2423  1.9620          factor(Species)2:factor(var)2   -0.0776  0.9217  -0.0842  0.9329  -1.8841  1.7288               ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       However, the `0` level of the `var` variable is not in the second call to `rma`. Is there a different way to specify the model to get the estimates for that variable?


title: Do we really need a control group when we can compare the treatment group with the status quo?

text: Here I am not talking about medical research, where there is placebo effect need to be controlled with the blind design. In social science studies, we almost never do blind experiment. In this case, if we have data from all the subjects, do we really need to have a randomly selected control group to be compared to? Say, if I have a teaching strategy to improve 5th grade students' math achievement and I want to test its effect. Student math achievement is measured by statewise standard assessments at every end of school year and all data will be extracted from district database directly. Which of the following is better: 1) randomly select 20 schools from 200 schools in a district, and then randomly assign 10 schools to treatment group (use the new teaching strategy) and 10 schools to control group (do business as usual). Compare the control and treatment. 2) randomly select 10 treatment schools from the same district. Compare the treatment with the rest of the district (control on previous achievement and demographics)? 3) randomly select 10 treatment schools from the same district. Compare the treatment with 10 randomly selected schools (post hoc) from the same district (control on previous achievement and demographics)? Is method 1 a better design and why?


title: Persistence in time series

text: Could someone tell me what the term 'persistence' mean in time series analysis? It's regarding econometrics and applied regression.


title: Bayesian model averaging in the case of large number of predictors

text: In the context of a linear factor model, Bayesian Model Averaging (BMA) is used to obtain the posterior probability of all possible combinations of predictors. A final model is obtained as a weighted average of all models, where the weight of each model is its posterior probability. In the case of a large number of factors K, the total number of possible combinations, i.e. models, is $2^K$. Therefore, it is almost impossible to average out over that many models. Is there an approach to get around this difficulty? Thanks,


title: Variation of Skewness and Excess Kurtosis

text: I apologize ahead of time if this a trivial question. I am building probability models for games - and in order to test the game code, I am comparing the results of large iterations of games with the original probability models. Specifically, I am calculating Standard Deviation, Skewness and Excess Kurtosis and using these calculations for verification. The sticking point lies in calculating appropriate confidence intervals since the underlying distributions are almost always not normal. I have been searching for quite some time for a solution - and the only thing which seems apparent is the BCa bootstrapping method, though I am not sure that this is really appropriate or necessary. It seems to me that there ought to be a way, given the original distribution, to properly calculate the expected variation not only of the mean, but also of the standard deviation, skewness and kurtosis. Is this possible?   * Due to some helpful comments, I am going to try to rephrase this (you may ignore the comments that I tried to make in response as they are a bit disjointed): Okay, let me see if I can rephrase this thought taking away the focus of the games since that is only relevant in terms of context. If I have a discrete distribution which is not normal - and then I proceed to recreate this distribution through some sort of Monte Carlo process - what is the best way to check that the Monte Carlo process is accurately recreating the original distribution? My thought was that since a given distribution is characterized by its moments, I could use this characterization. The problem with this approach is how to estimate whether the observed error in the moments (from the Monte Carlo simulation) is significant. It seems to me that there should be some way to do this.


title: D-efficiency value does not vary in D-optimal design of experiment

text: I am designing an experiment using D-Optimal design. With full factorial design I have to do 200 runs. With D-optimal, I can choose the number of runs I want to do and I evaluate the design using D-efficiency value. My question is the D-efficiency value varies in order of 0.001 when I change the number of runs. Can anyone give me an indication why this is happening? I expected larger number when I chose to run many experiment. Thank you in advance!


title: How to collapse groups in data defined by a group variable in R?

text: I have some multilevel data with a group variable (groups 1-100). Some of my other variables are within-group means, hence they have the same value repeated for every individual within group i. I'd like to reshape the data so that each group has only one line (instead of one for every individual within the group) so I can better use the within-group mean variables. For now, I don't care what happens with the individual-level variables. I'm guessing the reshape package must be the way to go, but I haven't had any success fiddling around with it. _Edit:_ Here's an example of what my data look like:               groupmean = c(1.4, 1.4, 1.4, 1.4, 6.2, 6.2, 6.2, 6.2)     group = c(1, 1, 1, 1, 2, 2, 2, 2)     mydata = data.frame(group, groupmean)      The repetition is because I have additional variables representing individual- level data (that I don't care about for now). I would like to reshape the above into a dataframe with only two lines, one for each group. Thanks in advance!


title: Sample size planning with required CI width and guessed propotion, using "Wilson" method

text: I refer to this link for calculating the sample size needed based on the width confidence interval needed, and the guessed proportion. I realize that it is use the "asymptotic" method when calculating the sample. Having know that the "wilson" method is a better way for estimating the CI for proportion[a]. I wonder if there is any "wilson" method for estimating the sample size? (Preferably a package used in R), thanks. [a]: A. Agresti and B.A. Coull, Approximate is better than "exact" for interval estimation of binomial proportions, American Statistician, 52:119–126, 1998.


title: Expected number of coin tosses to get N consecutive, given M consecutive

text: Interviewstreet had their second CodeSprint in January that included the question below. The programmatic answer is posted but doesn't include a statistical explanation. (You can see the original problem and posted solution by signing in to the Interviewstreet website with Google creds and then going to the Coin Tosses problem from this page.) **Coin Tosses** You have an unbiased coin which you want to keep tossing until you get N consecutive heads. You've tossed the coin M times and surprisingly, all tosses resulted in heads. What is the expected number of additional tosses needed until you get N consecutive heads? Input:   The first line contains the number of cases T. Each of the next T lines contains two numbers N and M. Output:   Output T lines containing the answer for the corresponding test case. Print the answer rounded to exactly 2 decimal places. Sample Input:   4   2 0   2 1   3 3   3 2 Sample Output:   6.00   4.00   0.00   8.00 Sample Explanations:   If N = 2 and M = 0, you need to keep tossing the coin until you get 2 consecutive heads. It is not hard to show that on average, 6 coin tosses are needed.   If N = 2 and M = 1, you need 2 consecutive heads and have already have 1. You need to toss once more no matter what. In that first toss, if you get heads, you are done. Otherwise, you need to start over, as the consecutive counter resets, and you need to keep tossing the coin until you get N=2 consecutive heads. The expected number of coin tosses is thus 1 + (0.5 * 0 + 0.5 * 6) = 4.0 If N = 3 and M = 3, you already have got 3 heads, so you do not need any more tosses. All of the mathematical equations I came up with had the right answers for the sample input data listed above, but was wrong for all of their other input sets (which are not known). Their programmatic solution appears to solve the problem far differently from my try-to-come-up-with-an-equation method. Can someone please explain how to come up with an equation that would solve this?


title: How can I draw a boxplot without boxes in R?

text: Using R, I'd like to plot two boxplots without the boxes—just the points. Creating clean boxplots in R is trivial:               business <- runif(50, min = 65, max = 100)     law <- runif(50, min = 60, max = 95)          boxplot(business, law, horizontal=TRUE, names=c("Business", "Law"), col=c('green', 'red'), main="Salary example (boxplot)")      ![Regular boxplots](http://i.stack.imgur.com/MTqxe.png) However, the only way I've found to plot just the points in the two random distributions seems needlessly complicated: I overlay two scatterplots with each variable plotted against either 1 or 2, to make a flat line:               plot(business, rep(1, length(business)), xlim=range(business, law), ylim=c(0, 3), pch=20, col='green', main="Salary example (dots)")     points(law, rep(2, length(law)), col='red', pch=20)      ![Sample dotted boxplot](http://i.stack.imgur.com/ROmNa.png) While this works, it will require a ton more tweaking to get the axes, tickmarks, and labels to match what R does with `boxplot()`. It seems that there has to be a simpler, more R-like way to do this. What's the best way to draw a boxplot without the box and whiskers—just the individual points?


title: Approximating binomial distribution by normal distribution

text: I got a very small probability in my homework so are there any mistakes? Benford's law says that if we have a collection of numbers, the first digit has the distribution $\mathbb{P}(D=k)=\log_{10}(1+1/k).$ Take a collection of 1200 numbers. Approximate by using normal distribution that the probability that those numbers has first digit 1, 2 or 3 at least 680 cases. My approach: One can compute easily that $\sum_{k=1}^9\log_{10}\left (1+\frac{1}{k}\right )=1.$ Benford's law says that the probability that the first digit is 1, 2 or 3 is $\sum_{i=1}^3\log_{10}(1+1/i)\approx 0.778$. Now $1200\cdot 0.778>10$ and $1200(1-0.778)>10$, so it is reasonable to approximate binomial distribution by a normal distribution. Now $\operatorname{Bin}(1200,0.778)$ is approximately $N(1200\cdot 0.778,1200\cdot 0.778(1-0.778))=N(933.6,207.3)$. Therefore, $P(X\leq 680.5)\approx P\left (Z\leq \frac{680.5-1200\cdot 0.778}{\sqrt{1200\cdot 0.778\cdot 0.222}}\right )\approx P(Z<-17.6)$. Now $\Phi(-17.6)\approx 0.$


title: How to interpret F-measure values?

text: I would like to know how to interpret a difference of f-measure values. I know that f-measure is a balanced mean between precision and recall, but I am asking about the practical meaning of a difference in F-measures. For example, if a classifier C1 has an accuracy of 0.4 and another classifier C2 an accuracy of 0.8, then we can say that C2 has correctly classified the double of test examples compared to C1. However, if a classifier C1 has an F-measure of 0.4 for a certain class and another classifier C2 an F-measure of 0.8, what can we state about the difference in performance of the 2 classifiers ? Can we say that C2 has classified X more instances correctly that C1 ? Sorry if it is a stupid question...


title: Non-increasing Log Likelihood with Expectation Maximization

text: I am struggling with my implementation of the expectation maximization (EM) algorithm for a certain model. The sequence of log likelihood values is not increasing, which is contradicting the theory. The measured outcome variable y is binary and is generated by one of two distributions:   * p(y=1|z=0) = xi, a fixed probability   * p(y=1|z=1,x) = 1/(1+exp(theta'*x)) So either there is a fixed probability of obtaining a 1, or the probability depends on the instance's characteristics x, and is modeled with a logistic regression. You can find a more detailed derivation (with the actual LL and the E and M step) in the following document: https://dl.dropbox.com/u/46839575/MMapproach.pdf The setting is customer churn prediction, but I guess no background knowledge is needed to check the modeling approach. I have implemented this in Matlab, and the problem is that the LL does not increase monotonically. I have checked my code several times, and I do not find any mistakes until now, so I am starting to wonder whether my derivation is correct. Maybe someone of you could point to any incorrect assumptions? In particular, I started to doubt whether the assumption that my latent class membership is modeled with the person's characteristics is correct (in most (Gaussian) mixture models this is a fixed probability, not varying over different persons). Furthermore, the same characteristics are used to model latent class membership and churn probability (although my model also gives non increasing log likelihood values when the set of characteristics in both models is mutually exclusive). I hope my explanation is clear. Thanks and best regards, Thomas


title: How to maximize (optimize) a response, given only a dataset of responses and features?

text: Let's say I have an n x p dataset. For each n, I have the response, 'y', and p - 1 features associated with it. What is the best way to determine the values of the features that will maximize 'y'? The only way I can think of is to build a model, regression or otherwise, and then perform multivariable optimization on the resulting equation. Is there a better way to do this?


title: How remove variations in a time series X caused by another time series Y?

text: I have a time series on a monthly basis (a commodity) of which much variation is caused by the weather. I want to adjust this commodity for weather changes. I use Heating degree day as a proxy for the weather conditions (If the average temperature of one day is below a threshold value, the absolute value of the difference is the heating degree for one day. This is summed up for all days in a month.) Important: I already adjusted for the series for seasonality. How should I proceed?


title: How to validate classifier (built by using MLN method)?

text: I have developed a method (let's call it `Method X`) that has a classifier function. The classifier function was built by using MLN (Markov logic network). I need to validate the result of applying the `Method X` to identify the root cause of problem in a textile factory. Since I use MLN, the data I use is non-parametric data. So, depending on the data and MLN method that I use, what method do you suggest to validate the `Method X`?


title: lme / lmer - split-plot with Non orthogonal subdesign

text: I have an agricultural field experiment (testing a plant protection agent): **Split plot design** with:               2 whole plot treatments "infestation": "high" & "low"      8 split-plot treatments ("treat"):           1. Untreated Control ("Ctrl1")     2. Reference Product ("Ctrl2")     3. 1 x Test-Product 1     4. 2 x Test-Product 1     5. 3 x Test-Product 1     6. 1 x Test-Product 2     7. 2 x Test-Product 2     8. 3 x Test-Product 2          and 4 replicates ("block"):      The parameter of interest in this example is grain ( **yield** ): First, I could model this:               lme(yield ~ infestation * treat, random = ~ 1 | block/infestation, data)      or                lmer(yield ~ infestation * treat + (1 | block/infestation), data)      But as can be seen treatments 3-8 can and have to be recoded as 2 products ("prod") being tested 1-3 times ("times"), so I have a 2x3 subdesign. One possibility would be subsetting the data:                 data2 <- subset(data, !data$treat == "Ctrl1" &  !data$treat == "Ctrl2")      and recode the resting treatments to "prod" = 1,2 and "times" = 1:3 then run:               lme(yield ~ infestation * form * times, random = ~ 1 | block/infestation, data)      Afterwards I could still do contrasts to compare the control treatments with the treated ones. But (here my actual problem starts): I read an article of **H.P. Piepho** : " _A Note on the Analysis of Designed Experiments with Complex Treatment Structure_ ", HortScience 41(2):446--452. 2006 The author wants to show " _how a meaningful analysis can be obtained based on a linear model by appropiate coidng of factors. (...) Our main objective is to demonstrate that the introduction of dummy variables can conveniently solve a wide variety of inferential problems that would otherwise either require ... multiple linear contrasts... or not make fully eficient use of the data, e.g when only data from orthogonal subdesigns are analysed._ " A very similar example (Example 1 in the article) is discussed within, and an alternative analysis in SAS is proposed - which I wanted to try to realise in R. The author adds a dummy variable ( **ctrl_vs_trt** ) to the data and codes it: "control", "trt" (in my case **trt** , **Ctrl1** , **Ctrl2** ". The he uses: (in his case **prod** is **form** ulation, and **times** is **conc** entration)               PROC GLM;     CLASS block contr_vs_trt form conc;   ##      MODEL set = block contr_vs_trt             contr_vs_trt * form         contr_vs_trt * conc         contr_vs_trt * form * conc;     RUN.                          I cite a further paragraph: " _Of course, a test for contr_vs_trt is not produced with this model, and one cannot compute simple means or marginal means. Also, the Type I SS for **form** , **conc** , and **form x conc** are not the same as with Type III SS. With Type III SS, the test for form is adjusted for **conc** , as fitting **conc** takes out the control when coding factors as in Table "xy" (as I did here). Similarly, the test for **conc** is adjusted for **form** , because fitting of **form** takes out the control. As a result, the Type III ANOVA for the model **form x conc** turns out to be that for the 3x2 factorial subdesign. (...) It seems much more stringent and transparent to use the nested model **contr_vs_trt/(form x conc)** , as this properly reflects all nesting and crossing features of the design._" Now, how to do that in `lme` or `lmer`? lme does not run at all, even if I simplify to:                lme(yield ~ prod * times, random = ~1|block, data),       I get      Error in MEEM(object, conLin, control$niterEM) :       Singularity in backsolve at level 0, block 1      The term **prod * times** cannot be run ( **prod + times** logically can). Eliminating both controls from the data set resolves this problem. `lmer` runs with **prod * times** , but always given the message:                fixed-effect model matrix is rank deficient so dropping "x" columns / coefficients      I understand that the subdesign is not orthogonal and therefor dropping is occuring, but I cannot say if the analysis after dropping can still be right. Also, I do not know how to specify the full model (leaving out the "infestation" whole plot for a second):               lmer(yield ~ prod * times + (1|block/ctr_vs_trt), data)      **prod * times** is nested inside **ctr_vs_trt** but both are nested inside the same block (or whole plot). Is nesting of fixed effects possible in `lme` or `lmer` \- does it work as I proposed? Does it even make sense to run the full model? With `aov()` I get the model running, even the partitioning of Df's is right. But due to strong non-orthogonality it is not possible to assume that the results are right. I can get meaningful results subsetting and using contrasts, but I found the authors approach interesting and it would help in the analysis of some of my other trials. Thanks in advance for any help; I hope this question is not too long...


title: Instability of one-pass algorithm for correlation coefficient

text: What should I know about the instability of the Pearson product-moment correlation coefficient? When might I experience problems using this calculation? I will quote the following Wikipedia article for some background information: http://en.wikipedia.org/wiki/Pearson_product- moment_correlation_coefficient#Mathematical_properties The Pearson correlation can be expressed in terms of uncentered moments. Since $μX = E(X)$, $$σX2 = E[(X − E(X))2] = E(X2) − E2(X)$$ and likewise for Y, and since $$E[(X-E(X))(Y-E(Y))]=E(XY)-E(X)E(Y)$$ the correlation can also be written as $$\rho_{X,Y}=\frac{E(XY)-E(X)E(Y)}{\sqrt{E(X^2)-(E(X))^2}~\sqrt{E(Y^2)- (E(Y))^2}}$$ Alternative formulae for the sample Pearson correlation coefficient are also available: $$r_{xy}=\frac{\sum x_iy_i-n \bar{x} \bar{y}}{(n-1) s_x s_y}=\frac{n\sum x_iy_i-\sum x_i\sum y_i} {\sqrt{n\sum x_i^2-(\sum x_i)^2}~\sqrt{n\sum y_i^2-(\sum y_i)^2}}$$ The second formula above needs to be corrected for a sample: $$r_{xy}=\frac{\sum x_iy_i-n \bar{x} \bar{y}}{(n-1) s_x s_y}=\frac{n\sum x_iy_i-\sum x_i\sum y_i} {\sqrt{(n-1)\sum x_i^2-(\sum x_i)^2}~\sqrt{(n-1)\sum y_i^2-(\sum y_i)^2}}$$ The above formula suggests a convenient single-pass algorithm for calculating sample correlations, but, depending on the numbers involved, it can sometimes be numerically unstable.


title: algorithms to pick best collection of classifiers

text: I have a bunch of methods that classify a binary outcome. I'm trying to figure out if some combination of those classifiers is better than any others. I'm hoping to run a bunch of methods. I've thought of running random forests/cart and perhaps a SVM model. What are some other common methods typically done for this?


title: How to control intercept of only one IV-DV relationship in multiple linear regression

text: If I have a relationship between three variables (survivorship (outcome), time and dose (continuous) with interaction) how can I control the intercept for time but not dose in MLR? The reasoning for this is that survivorship must be 100% at time 0. If I set the intercept to 100 in the equation survivorship = b1*time + b2*dose + b3*time*dose + intercept it will also constrain the dose which is not necessarily true.


title: Neural network with skip-layer connections

text: I am interested in regression with neural networks. Neural networks with zero hidden nodes + skip-layer connections are linear models. What about the same neural nets but with hidden nodes ? I am wondering what would be the role of the skip-layer connections ? Intuitively, i would say that if you include the skip-layer connections, then the final model will a sum of a linear model + some non-linear parts. Is there any advantage or disadvantage in adding skip-layer connections to neural nets ?


title: in nonlinear binary classification problems, which is the optimal dimension for make it lineary separable?

text: My question pertains to linear separability with hyperplanes in a support vector machine. Is posible to determinate the optimal dimension in which i have to transform a training data set for make it lineary separable? How can i determinate the dimension of the transformed space?


title: Most similar combination of partitions to the rest of group

text: I have a number of sample sequences of data that I want to compare with another large number of sequences of data. These samples are partitioned in between and I want to see which combination of partition for these sequences are the best match. I have some criteria for comparing these sequences, but I am not sure what is the most efficient way of finding the best grouping of partition for each sample. For example we have a list of sequences like:               1,3,4,4,4,4,4,5     2,4,4,4,4,4,4,5,3,3,4     5,6,3,5,3,3,3,3     2,2,2      and a sample sequence of:               3,4,4,4,4,4,5||2,2,2      Now I want to check whether the best combination of these partitions [(3,4,4,4,4,4,5) and (2,2,2)] is [(3,4,4,4,4,4,5) and (2,2,2)] or (3,4,4,4,4,4,5,2,2,2) as a sequence of data. I might be able to do it easily for a combination of two partition, but the number of alternatives for the combination will rise very rapidly. What is the best way to do this? I might be able to set an ending the number of partitions like 6 or something, but I want to see if there is a better way to do this.


title: large variables and low sample (p > n) problem: ridge , LASSO, PLS, PCR which is most suitable for predictions

text: I am trying see whether to go for **ridge regression** , **LASSO** or **principal component regression** (PCR) or **Partial Least Squares** (PLS) in a situation where there are large number of variables / features (p) over samples (n) and my objective is prediction. This my understanding: (1) **Ridge regression** shrinks the regression coefficients, use all coefficients without making them 0. (2) **LASSO** also shrinks the coefficients, but also makes them 0 means that can do variable selection too. (3) **Principal component regression** truncates the components so that p < n, will discard p-n components. (4) **Partial least square** \- also constructs a set of linear combination of the inputs for regression, but unlike PCR it uses y (in addition to X) for construction. The main practical difference between PCR and PLSR is that PCR often needs more components than PLSR to achieve the same prediction error[see]. Let's see the following dummy data - that I am trying to work in similar one:               #random population of 200 subjects with 1000 variables           M <- matrix(rep(0,200*100),200,1000)     for (i in 1:200) {     set.seed(i)       M[i,] <- ifelse(runif(1000)<0.5,-1,1)     }     rownames(M) <- 1:200          #random yvars      set.seed(1234)     u <- rnorm(1000)     g <- as.vector(crossprod(t(M),u))     h2 <- 0.5      set.seed(234)     y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))          myd <- data.frame(y=y, M)      Implementation of four methods:                require(glmnet)             # LASSO         fit1=glmnet(M,y, family="gaussian", alpha=1)             # Ridge           fit1=glmnet(M,y, family="gaussian", alpha=0)            # PLS       require(pls)       fit3 <- plsr(y ~ ., ncomp = 198, data = myd, validation = "LOO")       # taking 198 components and using leave-one-out cross validation        summary(fit3)       plot(RMSEP(fit3), legendpos = "topright")           # PCR       fit4 <- pcr(y ~ ., ncomp = 198, data = myd, validation = "LOO")      The best description of the data is: (1) $P > n$, most of time > 10 times (2) Variables (X and Y) are **correlated** with each other with different degree My question is which strategy may be best for this situation ? Why ?


title: Comparing percentages from one sample and between samples

text: I've conducted a large survey with about 600 responders. Most of the questions were multiple choice questions. I'm interested in comparing several kinds of percentages. First, I'm simply interested in determining whether the percentage of people choosing answer A is the same or different from those choosing answer B, C, D. etc. Most of the online stuff I've read suggests this is a one sample test of proportions but those seem to require a theoretical or expected value to compare the percentage to and that's not the case here. My second broad category of analyses would be to compare BETWEEN populations. So, do males and females answer the questions differently. There I would be interested in seeing whether the array of answer choices differ (e.g. males chose A 40% of the time, B 50% of the time, and C 10% of the time while females chose A 30% of the time, B 10% of the time, and C 60% of the time) and whether specific percentages differ (e.g. is the percentage of males that chose A [40%] different from the percentage of females that chose A [30%]). Thanks for your advice.


title: Simple SVM Question

text: For a linear SVM, the documentation tells me the formula is: $$ \frac{1}{2}w^Tw+C\sum\limits_{i=1}^l\xi_i$$ **Please explain to me in layman's terms what _w_ (and _ξ_ ) represent.** Is _w_ the distance to the hyperplane, the weight of different variables,... ? - PS: Is there a difference with the following formula or is it just a different notation: $$f(x)=w^T x+b$$


title: How would you merge different hierarchies over the same domain?

text: Say that you have $K$ different tree hierarchies. Each tree $k_i$ has associated a subset of elements of a given domain, $S_{k_i} \in D$ and finally on the leafs of each tree you have subsets of elements of the elements attached to that tree, $S_{k_i,1},S_{k_i,2},...,S_{k_i,n_k}$. As an example, think of the modern phylogenetic tree, $k_i$, classifying animals from Asia, $S_{k_i}$ where $D$ is the animal kingdom. Another tree could be the Aristotelian classification of animals known to the Greeks. Say that we decide that the correct hierarchy for the elements in the domain is the phylogenetic tree, what would be the best procedure for classifying European animals that don't exist in Asia? Intuitively I'm assuming that one has to use the co-occurrence animals existing in Asia and known to the Greeks to influence the assignment of "known-to-the-Greek" animals. I also think that the Aristotelian hierarchy might also guide probabilities. How would you work out this problem? Keep in mind that my problem hasn't anything to do with the animals example, but I though that it was a good example to flesh out the abstract question of how to merge hierarchically grouped instances of a domain into a reference hierarchy.


title: Should the 'TPR of adaboost' be better than base classifiers'?

text: There are two different base classifiers which produce true positive rate (TPR) values 99.46% and 91.79%. When I use these base classifiers in adaboost, what should the new TPR be? Better than two of them? Better than one of them?


title: How to combine the chi-square distributions (with one model parameter) of many items in a sample?

text: I have a sample of items, for each of which I have fitted models to obtain the best-fitting ($\chi^2$-minimising) value of a parameter $\alpha$. So for each, I have the values of $\chi^2_i(\alpha)$ for every possible value of $\alpha$ in the grid. I want to combine all these distributions to obtain some sort of stacked probability distribution for the sample. If I just take the mean of $\chi^2$ at each $\alpha$, I get something that looks sensible out - but is this meaningful? I initially expected that looking at the "error" on $\alpha$ from this distribution (by looking at the values of $\alpha$ enclosed by $\Delta\chi^2 = 1$ away from the minimum) would yield an error similar to the standard deviation of the original best-fitting $\alpha$ values, but it doesn't. That is perhaps reasonable since the width of the $\chi^2$ minima for each item are quite varied (i.e. the items have quite different errors). How can I get a meaningful distribution out?


title: Feature selection with k-fold cross-validated least angle regression

text: I am using the least angle regression (LARS) to extract the most important predictors ($x_1, x_2,...,x_p$) for my response variable ($y$). I have seven predictors ($x_1,x_2,...,x_7$) for each response variable. I did 10-fold cross validation by using R package `lars`. I am making a cut point at 0.62 (approx.). Am I doing right? Is there any criteria to select important predictors? Any help in this regard would be highly appreciated! Thanks!


title: Why does ridge regression classifier work quite well for text classification?

text: During an experiment for text classification, I found ridge classifier generating results that constantly top the tests among those classifiers that are more commonly mentioned and applied for text mining tasks, such as SVM, NB, kNN, etc. Though, I haven't elaborated on optimizing each classifier on this specific text classification task except some simple tweaks about parameters. Such result was also mentioned Dikran Marsupial. Not coming from statistics background, after read through some materials online, I still cannot figure out the main reasons for this. Could anyone provide some insights on such outcome?


title: Weighted mean of two 2D Gaussian random variables

text: I have two 2D Gaussian random variables. I'd like to find the weighted average of the two (based on their covariance matrix around the mean, meaning that the mean of the final Gaussian should be closer to the one with the smallest variance). What is the proper way of doing that? Any pointer to web urls around this topic would be highly appreciated. Clarification: I have a robot that has two sensors. The robot reads two bi- variate normal observations of an object. I want to find the location of the object by finding the weighted mean of the two samples. So I want my final location that I calculate to be closer to the observation that has a smaller spread. What is the proper way of doing that?


title: Is it valid to model discrete numerical test scores as coming from a continuous random variable?

text: I'm working with a sample of test scores which range from 0 to 100. These scores are generated from a set of 100 binary responses (0 or 1), so the higher the resulting sum, the better the performance. The distribution of these scores is definitely not normal: it is highly skewed. I'm interested in carrying out a Box-Cox transformation to make this data more normal. Since the data is very discrete, but ordinal, I'm considering doing a Kernel Density Estimation (which would model a continuous distribution), and apply the transformation on the resulting (continuous-valued) model. Is this a valid approach? If not, I'd appreciate any pointers on modeling this sort of probability distribution.


title: Parallel solving Ax=b?

text: > Cross posted on StackOverflow. I have some extremely large sparse matrices created using spMatrix function from the matrix package. Using the solve() function works for my Ax=b issue, but it takes a very long time. Several days. I noticed that http://cran.r-project.org/web/packages/RScaLAPACK/RScaLAPACK.pdf appears to have a function that can parallelize the solve function, however, it can take several weeks to get new packages installed on this particular server. The server already has the snow package installed it. So   * Is there a way of using snow to parallelize this operation?    * If not, are there other ways to speed up this type of operation?   * Are there other packages like RScaLAPACK? My search on RScaLAPACK seemed to suggest people had a lot of issues with it. Thanks. Additional details   * The matrices are about 370,000 x 370,000. I'm using it to solve for alpha centrality, http://en.wikipedia.org/wiki/Alpha_centrality.    * I was originally using the alpha centrality function in the igraph package, but it would crash R.   * This is on a single machine with 12 cores and 96 gigs of memory (I believe)    * It's a directed graph along the lines of paper citation relationships.   * Calculating condition number and density will take awhile. Will post as it comes available.


title: Binomial probability function

text: What is the probability of rolling exactly two sixes in 6 rolls of a die? Solution by the Binomial Probability formula is $$\binom{6}{2} \left(\frac{1}{6}\right)^2 \left(\frac{5}{6}\right)^4 = \frac{15 \times 5^4}{6^6} = \frac{3125}{15552} \approx 0.200939.$$ But by basic probability understanding > Probability = successful outcome/ total possible outcomes So that way the probability should be $15/ (6^6) \approx 0.0003215.$ I know that the second calculation is not giving the right answer, but somehow am not convinced on why it is wrong as I don’t see the approach being wrong. Can someone help me understand why the second approach is wrong?


title: Multivariate EWMA

text: Is there any package in R which computes the Multivariate EWMA? I have a data frame of 4 series and I do not want to use a simple rectangular method to compute the covariance estimator. So is there any preimplemented function for this? Or could you give me help to implement this? My programming skills in R are very limited.


title: borrow strength accross studies (FDR)

text: Good evening everyone I have following problem. I hope someone can help me out a bit :) Consider a study were you perform 1000's of tests. You correct for multiple testing with FDR. Obviously you will miss a lot of true significant hits (false negatives). There are also 1000's of studies were similar tests are conducted. Some tests will test non-significant across all studies because of the multiple correction in all studies while the nulhypothesis is actually not true (false negatives). But it's not hard to see that if a test consistently has a relativly low p-value across different studies, the probability that the alternative hypothesis is true increases. It should be possible to use this interstudy information to increases the power of each test and ultimatly identify test that has relativly little proof for the alternative hypothesis to be truly positive/significant. I tried googling it but I'm not sure what to look for and can't find a clear solution. Does someone has a clue for me? best regards Thanks in advance


title: Can a growth model of diseases incorporate vital signs data?

text: I intend to study the growth rate of diseases among patients of different age groups (cohorts) and comment on the growth of diseases in the next 5 years. Can I incorporate vital sign readings such as body temperature, systolic blood pressure, etc. in the growth model? Is any other approach more suitable for this kind of analysis?


title: How to calculate Kullback-Leibler divergence/distance?

text: I have three data sets X, Y and Z. Each data set defines the frequency of an event occurring. For example: Data Set X: E1:4, E2:0, E3:10, E4:5, E5:0, E6:0 and so on..   Data Set Y: E1:2, E2:3, E3:7, E4:6, E5:0, E6:0 and so on..   Data Set Z: E1:0, E2:4, E3:8, E4:4, E5:1, E6:0 and so on.. I have to find KL-divergence between X and Y; and between X and Z. As you can see for some of the events there will be 0 and non-zero values. For some of the events all three data sets are 0. I would appreciated if someone can help me find the KL divergence for this. I am not much of a statistician, so I don't have much idea. The tutorials I was looking at online were a bit too complex for my understanding.


title: Sampling from conditional distribution in general case

text: I'm dealing with Gibbs Sampling now. Let's consider the example: I know the distribution of X|Y and the distribution of Y. They are some known - Binomial or Beta or other but particular. Thus I have in analytical view f(X|Y), f(Y) and I can calculate joint distribution f(X,Y). To provide Gibbs Sampling I need to calculate many times X∼X|Y and Y∼Y|X. The question is: which technique can I use to sample Y∼Y|X in general case, for any given f(X|Y), f(Y)?


title: Removing the difference between treatments

text: Is there any statistical method to remove the difference between treatments observed before starting the treatments? I have measured the length of shoots on trees of three different treatments but unlucky they were different even before treatment onset!


title: How to find residuals and plot them

text: I have been given data               x = c(21,34,6,47,10,49,23,32,12,16,29,49,28,8,57,9,31,10,21,26,31,52,21,8,18,5,18,26,27,26,32,2,59,58,19,14,16,9,23,28,34,70,69,54,39,9,21,54,26)     y = c(47,76,33,78,62,78,33,64,83,67,61,85,46,53,55,71,59,41,82,56,39,89,31,43,29,55,           81,82,82,85,59,74,80,88,29,58,71,60,86,91,72,89,80,84,54,71,75,84,79)      How can I obtain the residuals and plot them versus $x$? And how can I test if the residuals appear to be approximately normal? I'm not sure if I do the original linear fit correctly as I got the equation $y=6.9x-5.5$ but the lecture notes says that the linear regression line should be of the form $y_i=\beta_0+\beta_1x+\epsilon$.


title: How to Use Neural Networks to Forecast Time Series Data with Predictor Variables?

text: I have browsed a lot of topics here, but the ones I see were all about forecasting a single variable, depending on its historical values. Whereas I want to predict a variable, by estimating a relationship between multiple predictor variables as well. So I'm trying to find a function $f(x_1,x_2, ... , x_p) = y$ -How- Can I apply a neural network approach for such a task? If not, what alternatives can I use for this? ## Edit : I am trying to implement this in MATLAB, so I would really appreciate some MATLAB implementation of such methods.