title: Predicting the Success of a tweet

text: I want to predict the success of a tweet. In my case a tweet is successful when the sum of the number of favorites and the number of retweets is greater than 5. So my outcome value y is: **y= if((retweets + favorites) >5 ==1** else == 0 Therefore I created the following dataset: (This sample just shows the data for one twitter user and so friendCount etc. is always the same. In my real dataset there are thousand of users and their last tweets)               head(dataTrain)          +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | friendsCount | followersCount | nbTweets     | at | rt | link | nchar           | nwords           | r_nachar_words  | hash | no_hash | hash_mean    | y |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | 0.7105428815 | 0.7105428815   | 0.7105428815 | 0  | 0  | 1    | -20.471.177.338 | -125.525.459.396 | -22.192.436.207 | 0    | 0       | -0,777283372 | 1 |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | 0.7105428815 | 0.7105428815   | 0.7105428815 | 1  | 0  | 1    | -0.5838923895   | -0.68260381719   | 0.2862346486    | 1    | 1       | -0,330839739 | 0 |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | 0.7105428815 | 0.7105428815   | 0.7105428815 | 1  | 1  | 1    | -0.6884084855   | -125.525.459.396 | 19.565.534.949  | 1    | 5       | 0,562047528  | 1 |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | 0.7105428815 | 0.7105428815   | 0.7105428815 | 0  | 0  | 1    | -0.2006667041   | -0.10995304042   | -0.2562620706   | 1    | 3       | 0,6364548    | 0 |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | 0.7105428815 | 0.7105428815   | 0.7105428815 | 0  | 0  | 1    | 0.9490103521    | 122.623.210.539  | -0.6388650200   | 0    | 0       | -0,777283372 | 1 |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+     | 0.7105428815 | 0.7105428815   | 0.7105428815 | 0  | 0  | 1    | 0.3567524747    | 0.08093055184    | 0.4350079872    | 1    | 2       | 0,227214803  | 1 |     +--------------+----------------+--------------+----+----+------+-----------------+------------------+-----------------+------+---------+--------------+---+      So the dataset contains the features:   * friendsCount (how many people does the user follow?)    * followersCount (how many people do follow the user?)   * nbTweets (how many tweets did the user send in the last weeks?)   * at (does the tweet contain a @? )   * link (does the message contain a link? )   * nchar (number of characters in the tweet)   * nwords (number of words in the tweet)   * r_nachar_words (ratio betweent chars and words)   * hash (does the messsage contain a hashtag? )   * no_hash (how many hashtags does the message contain? )   * hash_mean (what is the mean length of the hashtags? )   * y (was the sum of number of favorites and number of retweets greater than 5? ) At the moment I use a Support Vector Machine for creating the model and the SVMlight algorithm. But it does not scale really well and is very slow. I get a accuracy of about 80% but very poor values for Sensitivity and Kappe. My Questions:   * What would be a better approach for creating the model? Would you prefer for example a deep learning network instead of a SVM?   * Do I use enough features? Could you think of other features which could increase the performance of my model?


title: Mann-Kendall test STATA

text: I am new in this forum. I am beginning to work with time series, I have a daily (25,000+ observations) temperature dataset (01/01/1946 - 07/01/2014) I want to test for the following:   1. **Trends:** So far I used OLS, but I have heard about using Mann-Kendall test may be useful.   2. **seasonality:** In the same way that in the previous point so far I tried to analyze them with OLS, any recommendation is very welcome.   3. **structural breaks:** the problem with this aspect is that graphically it cannot be observed any breaks, therefore I have been trying with CUSUM and CUSUMSQ to identify this structural breaks, some additional recommendations? As additional point I need to perform this on STATA. I saw in the Web about Kendall package for R, but I cannot use R. Currently I use Gretl and STATA, I hope someone might help me. Thanks in advance for your help. Jorge


title: Paired Differences T-tests for Odd number of Observations

text: Suppose we have $m$ baseline data points and $n$ post-baseline data points where $m < n$. This is for one subject. Is it possible to perform a paired t-test on this data even though $m+n$ is odd? Can we thrown away the data that is not paired and perform the test? For example, suppose the baseline data is $(1,2,3)$ and the post-baseline data is $(6,7,8,9,10)$ for one subject. Then we can form the pairs $(1,6)$, $(2,7)$ and $(3,8)$. The $9$ and $10$ would be thrown out. I guess if $m$ is not much greater than $n$, then we can still perform a paired t-test without a loss of power? Or perhaps do some data imputation?


title: How do I define generators in Minitab

text: How do I define generators in Minitab 16? I have an experiment design (DOE) with 5 factors and 2 generators. I entered D=ABC E=AC in the box of generators but it says only F through G ?


title: ANOVA uncertainty

text: If I have 10 classification accuracies for 3 different parameter, and I perform an ANOVA test which yields $F=1.19$ and $\text{Prob>F}$ $0.3201$. The parameter value indicates the number of features. Does that mean none of the 3 parameter values yields a significantly different result than the other? More importantly, how can I justify choosing one parameter value over the other? The thing is with 3 features I get an accuracy of 96%, while with 2 and 4 features I get 92% and 88% respectively. Hence, I want to justify using 3 features instead of 2.


title: quantile transformation with t distribution

text: I'm wondering about the monotonic quantile transformation in case of a t distribution. First, for explanation, lets consider a normal distribution: \begin{align*} f(l | \mu , \sigma ^2)=\frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{(l-\mu)^2}{2 \sigma ^2}} \end{align*} A quantile at $\alpha$ can be calculated by $Quantile_\alpha=\mu + \sigma \Phi^{-1}(\alpha)$ so e.g. lets say $\mu = -0.0000615$ and $\sigma=0,0145$, the 0,95 quantile of the normal distribution is 1,644854 so $Quantile_{0.95}=-0.0000615+0,0145*1,644854=0.02379$ So now for the t distribution, I use the following notation: \begin{align*} f(l|\nu ,\mu ,\beta) = \frac{\Gamma (\frac{\nu+1}{2})}{\Gamma (\frac{\nu}{2}) \sqrt{\pi \nu} \beta} \left(1+\frac{1}{\nu}\left(\frac{l - \mu}{\beta}\right)^2 \right)^{\text{$-\frac{1+\nu}{2}$}} \end{align*} now, I have this from a book: ![quantile](http://i.stack.imgur.com/dnrep.png) ignore the "VaR", it only sais, we focus on the right side of the distribution. And now I am not really sure about the quantile computation anymore. I think my $\beta$ is the $\sigma^2$ right? So is this correct: $Quantile_{\alpha}^{t}=\mu + \beta t_\nu^{-1}(\alpha) $ or this: $Quantile_{\alpha}^{t}=\mu + \sqrt{\beta} t_\nu^{-1}(\alpha) $ As an example, I have the following values: $\mu=-0.0012827$ $\beta=0.009891722$ $\nu = 3.34331$ and therefore $t_{3.34331}^{-1}(0.95)=2.258039$ now, is $Quantile_{0.95}^{t}=-0.0012827+0.009891722*2.258039=0.02105319$ or $Quantile_{0.95}^{t}=-0.0012827+\sqrt{0.009891722}*2.258039=0.2232954$ correct? Mh, I am not sure, since I am confused about the $\sigma^2$ they use, but 0.2232954 CANNOT be empirical correct, so it MUST be the 0.021? Where is my mistake? the book: Alexander J. McNeil, Rüdiger Frey, Paul Embrechts: Quantitative Risk Management 2005 page 40.


title: Interpreting linear regression with NA for estimates

text: I have 3 variables. Income is the dependent variable(continuous), sex(binomial), workhrs(categorical) i fitted a linear model and obtained the result: My question is: Is this model usable? The NA values come from the fact that workhrs is missing 4 categories. The fifth one is i think due to coliniarity. The adjusted R squared is pretty high, but is the model Usefull to predict? or the interaction is not useful and should be dropped?               Call:     lm(formula = income ~ factor(sex) + factor(workhrs) + factor(sex):factor(workhrs))          Residuals:        Min     1Q Median     3Q    Max      -6.757 -1.261  0.000  1.531  5.690           Coefficients: (5 not defined because of singularities)                                    Estimate Std. Error t value Pr(>|t|)         (Intercept)                      31.330      4.645   6.744 8.42e-05 ***     factor(sex)1                      0.160      6.569   0.024 0.981101         factor(workhrs)1                  2.930      5.689   0.515 0.618960         factor(workhrs)2                 13.240      8.691   1.523 0.161968         factor(workhrs)4                 10.935      5.689   1.922 0.086773 .       factor(workhrs)5                  6.960      6.569   1.059 0.316997         factor(workhrs)6                 30.045      8.691   3.457 0.007192 **      factor(workhrs)7                 18.315      5.689   3.219 0.010503 *       factor(workhrs)8                 26.947      5.364   5.024 0.000715 ***     factor(workhrs)9                 25.340      6.569   3.857 0.003863 **      factor(workhrs)10                38.850      9.291   4.182 0.002370 **      factor(sex)1:factor(workhrs)1     9.160      8.691   1.054 0.319347         factor(sex)1:factor(workhrs)2        NA         NA      NA       NA         factor(sex)1:factor(workhrs)4     9.060      8.046   1.126 0.289276         factor(sex)1:factor(workhrs)5    23.035      8.691   2.651 0.026451 *       factor(sex)1:factor(workhrs)6        NA         NA      NA       NA         factor(sex)1:factor(workhrs)7     3.325      8.691   0.383 0.710895         factor(sex)1:factor(workhrs)8        NA         NA      NA       NA         factor(sex)1:factor(workhrs)9        NA         NA      NA       NA         factor(sex)1:factor(workhrs)10       NA         NA      NA       NA         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 4.645 on 9 degrees of freedom     Multiple R-squared:  0.9263,    Adjusted R-squared:  0.8117      F-statistic: 8.083 on 14 and 9 DF,  p-value: 0.001742


title: AUC is equivalent to a Mann-Whitney U-score, is the basic multiclass AUC related to the Kruskal-Wallis test statistic?

text: I've read that the area under the ROC curve is equivalent to a Mann-Whitney U-score. Is a multiclass AUC score (which averages the AUC scores for pairs of classes) related to the Kruskal-Wallis test statistic?


title: Non parametric measure of strength of association between a discrete and a continuous random variable

text: I'm throwing here the problem as I received it. I have two random variables. One of which is **continuous** (Y) and the other one which is discrete **ordinal** (X). I put below the plot I received together with the query. ![enter image description here](http://i.stack.imgur.com/bidGa.png) The person who send me the data wants to **measure the strength of the association** between X and Y. I'm looking for ideas that would not come front loaded with assumptions about what process generated the data. Note that this is not about finding a non parametric _test_ the strength of the relationship (as in bootstrap) but about finding a non parametric way to _measure_ it. On the other hand, efficiency is not an issue since there is a lot of data points.


title: R mvpart - splitting index

text: I need to understand what `mvpart` is doing. Which index does it use as a splitting criterion (in my case, method=class)? Does it use simultaneous partitioning along multiple axis (multivariate)? Are there any references I can read on the subject?


title: Calculate normal distribution of a calculation, function of several normally distributed parameters

text: I usually make calculations of the estimated profitability of some businesses. These calculations are function of some parameters: estimated prices, estimated demand, estimated exchange rates, ... I often assume that all of these input parameters are normally distributed, and I assess their means and standard deviations. My questions are: Which is the best way to predict the distribution of the estimated profitability? Can I determine mean + standard deviation of the profitability? I happened to create, by using `R` or `NumPy`, 1,000 random values of prices (based on the normal distribution), other 1,000 random values for the demand, etc., make $1,000^N$ calculations, and obtain the mean + standard deviation of the $1,000^N$ values, always assuming that the result is a normal distribution as well.


title: Use of a Poisson distribution from a distance matrix to determine dbscan parameters

text: I´ve been researching about automatic determination of parameters for DBSCAN (a density-based clustering algorithm -- http://en.wikipedia.org/wiki/DBSCAN), especially _eps_ , and have found the following paper, which proposes a method for determining the best _eps_ from the data themselves: http://www.joics.com/publishedpapers/2012_9_7_1967_1973.pdf However, the paper claims, in page 1969, section 3.1, about a distance matrix DIST$_{n \times n}$: > We calculate the value of each element in DIST$_{n\times n}$ successively, > and then sort them in an ascending order line by line. We use DIST$_{n > \times i}$ to express the $i$th column value in DIST$_{n \times n}$. As all > data in DIST$_{n \times i}$ obeys the Poisson distribution[6], we use the > maximum likelihood estimation in mathematics to estimate the value of > parameter $\lambda$. That is to say, $\lambda$ can be obtained by means of > the geometrical mean of the value of DIST$_{n \times i}$. The paper goes to claim that this $\lambda$ can be used as _eps_. Is there a reason for this claim? I have read about the Poisson distribution, but nowhere I have found found support for the claim that distances should follow this distribution (I admit to not being very knowledgeable about probability distributions, however). If this claim is not spurious, it would enormously simplify working with DBSCAN, as it would diminish the guesswork necessary to determine parameters.


title: Correcting data for heteroscedasticity in a regression model

text: I applied OLS on a regression model that looks as follows: $$ y = b_0 + b_1x_1 + b_2x_2 $$ and found that signs of heteroscedasticity. In an econometrics text book, I found that I can divide each variable by the square root of an independent variable. So the new model looks as follows: $$ \frac{y}{\sqrt{x_1}} = \frac{b_0}{\sqrt{x_1}} + \frac{b_1x_1}{\sqrt{x_1}} + \frac{b_2x_2}{\sqrt{x_1}} $$ This is only for $\sqrt{x_1}$ but I will have to do it for $\sqrt{x_1}$ as well. The problem I have is that some of the $x_1$ or $x_2$ values is zero. So what do I do if I want to divide by the $\sqrt{0}$?


title: Steps followed when Binary logistic regression when both dependent and independent variables are binary

text: I had set of binary variables. To apply logistic regression, I have checked association between dependent and independent variables and considered only those independent variables in the model which came to be associated with dependent variable. My query is whether it is an appropriate way of fitting logistic regression model.


title: Looking for measurement for nonlinear model

text: I use KNN regression to train out a model. The model estimates running time of a program based on different inputs, and the output is a single variable, which is time (double type). I want to evaluate this model. I have a list of observation times and a list of prediction times. I would like to know how I can evaluate the model based on these two lists. What kind of measurements can I use? Two lists for example: Observation List: {0.01, 0.03, 0.02, 0, 0.02} Prediction List: {0.05, 0.02, 0.04, 0.37, 0.05}


title: Recall and precision in classification

text: I read some definitions of recall and precision, though it is every time in the context of information retrieval. I was wondering if someone could explain this a bit more in a classification context and maybe illustrate some examples. Say for example I have a binary classifier which gives me a precision of 60% and a recall of 95%, is this a good classifier ? Maybe to help my goal a bit more, what is the best classifier according to you ? (dataset is imbalanced. Majority class has twice the amount of examples of the minority class) I'd personally say **5** because of the area under the receiver operator curve. (as you can see here model 8 has a low precision, very high recall, but one of the lowest AUC_ROC, does that make it a good model ? or a bad one ?) ![enter image description here](http://i.stack.imgur.com/LZAHw.jpg) * * * **edit:** I have an excel file with more information: https://www.dropbox.com/s/6hq7ew5qpztwbo8/comparissoninbalance.xlsx In this document the area under the reciever operator curve can be found and the area under the precision recall curve. Together with the plots.


title: How to generate a pair of ellipses shaped as Saturn in $\mathbb{R}^p$

text: I'm trying to generate multivariate data shaped like Saturn (long story). More formally;   1. *cluster 1 is a rank-p Gaussian with correlation matrix R (where every out of diagonal entry of R is the same number in $(0,1)$).   2. cluster 2 is a bunch of points distributed on a hyper-plan of rank p-1, around the equator of cluster 1, all located beyond a certain mahalanobis distance wrt to cluster1 of say $\zeta$. Rejection sampling is ok, so long as the rejection rate is small: this configuration has to be generated many times for $n$ and $p$ large reliably. **Update:** So i implemented Bob Durrant's solution:               x0<-matrix(rnorm(n*p),n,p)     b0<-matrix(rnorm(n*p),n,p)     b1<-qchisq(0.99,df=p);     b2<-sqrt(c(b1,b1*1.25))     b0<-b0/sqrt(rowSums(b0*b0))*runif(n,b2[1],b2[2])     plot(rbind(x0,x1))     mahalanobis(cbind(0,b0),colMeans(x0),var(x0))/(qchisq(0.99,df=p))      yielding the attached picture, which look indeed like Saturn ![Saturn](http://i.stack.imgur.com/Rv6VE.jpg). with the rings located at least @ qchisq(0.99,df=p) away from the center. Now however, i want my saturn to be ellipse shaped, i.e. to have correlation structure R. The problem is that is i pre-multiply $\verb+x0+$ and $\verb+cbind(0,x1)+$ by $R^{1/2}$, the rings are no longer on the equator :( **Update2:** For example, things already go awry when i replace the diagonal variance structure above with something else. For example:               library(MASS)     x0<-mvrnorm(n,rep(0,p),diag(rchisq(p,p),p))     b0<-mvrnorm(n,rep(0,p-1),diag(rchisq(p-1,p-1),p-1))     b1<-qchisq(0.99,df=p);     b2<-sqrt(c(b1,b1*1.25))     b0<-b0/sqrt(rowSums(b0*b0))*runif(n,b2[1],b2[2])          mahalanobis(cbind(0,b0),colMeans(x0),var(x0))/(qchisq(0.99,df=p))      most of the distances are way too large (compare with the output to the same call after using diagonal co-variance matrix above)! **Update3**               delta<-0.9     p<-3     R<-matrix(runif(p^2,delta*0.99,delta),p,p)      #to avoid repeating eigen-values, i jitter a bit.     diag(R)<-1


title: Accurate prediction of binomial group level outcome with lme4

text: I have a dataset in which individuals have rated a number of subjects. I am trying to predict a subject-level outcome (binary) on the basis of these ratings. While analyzing these data, I noticed that whenever I include a random intercept for subject, this random effect accounts for almost all variance in the subject-level outcome, so there is no prediction. Trying to find out whether this is caused by no relation between predictor and outcome, or by a problem in a way that I modelled the data, I wrote the following very simple syntax that demonstrates the problem. This generates the group-level data.               set.seed(1234)     pred <- rnorm(n=200, mean=5, sd=2)     outcome <- rbinom(length(pred), size=1, prob=.5+(pred / 25))     summary(glm(outcome ~ pred, binomial))      The relevant B estimate and p-value:               pred         0.21782    0.07951   2.740  0.00615 **      Now both predictors and outcomes are all repeated 5 times, to resemble that there are five (identical) observations of the predictor in a single group, and the group-level outcome does not vary:               pred <- rep(pred, each=5)     outcome <- rep(onoff, each=5)     grp <- factor(rep(1:(length(pred) / 5), each=5))     summary(glmer(outcome ~ pred + (1 | grp), family=binomial))      The output:               Random effects:      Groups Name        Variance Std.Dev.      grp    (Intercept) 4038.8   63.552       Number of obs: 1000, groups: grp, 200          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)     (Intercept)   12.389     22.094   0.561    0.575     pred           0.204      4.793   0.043    0.966      The same model using aov (I know, not suitable for binary outcomes) gives results that look much better:               > summary(aov(outcome ~ pred + Error(grp)))          Error: grp                Df Sum Sq Mean Sq F value  Pr(>F)        pred        1   8.86   8.858   8.078 0.00495 **     Residuals 198 217.12   1.097                         OK, so what about adding some intra-group variance? Maybe some intragroup variance would help to get this to work:               indiv <- factor(rep(1:5, times=length(unique(grp))))     pred <- pred + rnorm(length(pred), 0, .2)     summary(glmer(outcome ~ pred + (1 | grp/indiv), family=binomial))           indiv:grp (Intercept) 2.6142e-10 1.6168e-05      grp       (Intercept) 4.0416e+03 6.3574e+01     Number of obs: 1000, groups: indiv:grp, 1000; grp, 200          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)     (Intercept)  12.4160    20.1808   0.615    0.538     pred          0.1986     4.3130   0.046    0.963      And again, the aov output is better:               summary(aov(outcome ~ pred + Error(grp/indiv)))          Error: grp                Df Sum Sq Mean Sq F value  Pr(>F)        pred        1   8.73   8.731   7.958 0.00528 **     Residuals 198 217.24   1.097                         Inspecting the random effects in the two glmer models shows that the group intercept is high when the outcome is 1, and low when the outcome is 0 - thus explaining all group-level variance. It might be the case that I am doing something totally wrong here, but I was wondering whether this is indeed way that lme4 was intended to behave. If you have any suggestions on how to improve the model specification, I would be happy. PS: I know I can just aggregate the data to the group level, but the actual dataset is not complete; i.e., not all judges have judged all subjects, and not all subjects were judged by all judges. Thus, I would like to use a model that I can use random intercepts for both subjects and judges in; and aggregating the data to the group level would not take the inter-judge variance into account.


title: Clustering of points based on vector feature similarities in R

text: I have as an input a number of points that I need to partition into clusters. Each point has a number of features that are ideally to be used to find the similarity between each point and the others. Some of these features are scalar values (one number) and others are vectors. For example, assume that each point has the following features:   1. **S1** : scalar value   2. **V1** : 48 $\times$ 1 vector   3. **V2** : 48 $\times$ 1 vector For example one point may have (S1,V1, V2) as (100, {0, 100, 20, 30}, {75,0,10, 5}) My hypothesis is to use cosine similarity to find how similar the vector V1 or V2 of one point is to the vector V1 or V2 of another point. I have already computed the similarity matrices between all points in terms of V1 and V2 similarities. By exploring the standard clustering algorithms in R, I have found that k-means turns to use the Euclidean distance, which might be suitable for clustering points according to their scalar values, because [subject unclear] doesn't work for the situation where I have hybrid types of features (scalars and vectors). Also the K-medoid clustering seems to be supporting only the Euclidean and the Manhattan distances. I think what should be done is to generate one more distance/similarity matrix between all points based on the scalar value, so that we end with three similarity matrices that show the similarity between each point and the other points according to each feature regardless of it being a scalar or a vector, and use those matrices for finding the neighbourhood of points while clustering. I wonder if there is an implementation for a clustering algorithm that accepts as an input the similarity matrices (or alternatively the dissimilarity/distance matrices) between vector features of multiple points and uses them for clustering?


title: Change in binomial proportion confidence interval

text: I'm having trouble calculating 95% confidence intervals for a change in binomial proportion. For example, in group $A$, there are $4$ successes out of $n =20$. In group $B$, there are $12$ successes out of $n =20$. I am able to calculate the CIs for the proportions separately in R using:               binom.test(4,20,(1/2),alternative="two.sided")     binom.test(12,20,(1/2),alternative="two.sided")      However, I'm interested in the change in proportion for these two groups, group $A = 0.2$ and group $B = 0.6$ with the change in proportion being $0.4$. What I'm really lost on is how to calculate the CIs for this change ($0.4$)?


title: How to look for valleys in a graph?

text: I'm examining some genomic coverage data which is basically a long list (a few million values) of integers, each saying how well (or "deep") this position in the genome is covered. I would like to look for "valleys" in this data, that is, regions which are significantly "lower" than their surrounding environment. Note that the size of the valleys I'm looking for may range from 50 bases to a few thousands. What kind of paradigms would you recommend using to find those valleys? **UPDATE** Some graphical examples for the data: ![alt text](http://i.stack.imgur.com/5UEOg.png) ![alt text](http://i.stack.imgur.com/tlaV1.png) **UPDATE 2** Defining what is a valley is of course one of the question I'm struggling with. These are obvious ones for me: ![alt text](http://i.stack.imgur.com/MaR6O.jpg) ![alt text](http://i.stack.imgur.com/GP9UZ.jpg) but there some more complex situations. In general, there are 3 criteria I consider: 1\. The (average? maximal?) coverage in the window with respect to the global average. 2\. The (...) coverage in the window with respect to its immediate surrounding. 3\. How large is the window: if I see very low coverage for a short span it is interesting, if I see very low coverage for a long span it's also interesting, if I see mildly low coverage for a short span it's **not** really interesting, but if I see mildly low coverage for a long span - it is.. So it's a combination of the length of the sapn and it's coverage. The longer it is, the higher I let the coverage be and still consider it a valley. Thanks, Dave


title: Regression with correlated errors

text: In a standard mean regression setting $y_i=x_i^T \beta + \epsilon_i$, if we assume the errors $\epsilon_i$ are stationary with some covariance $\Sigma$, what can be said about $x_i$? Is $x_i$ also stationary? Also, is $\epsilon_i|x_i$ stationary? What if $x_i$ includes nonstationary covariates, what can be said about $\epsilon_i$.


title: Mean of a sliding window in R

text: I have a vector of values that I would like to report the average in windows along a smaller slide. For example, for a vector of the following values:               4, 5, 7, 3, 9, 8      A window size of 3 and a slide of 2 would do the following:               (4+5+7)/3 = 5.33     (7+3+9)/3 = 6.33     (9+8)/3 = 5.67      And return a vector of these values:               5.33, 6.33, 5.67      Is there a simple function that will do this for me? If it also returned the indices of the window starts that would be an added bonus. In this example that would be 1,3,5


title: Should anova(model, test="Chisq") not produce a test statistic as well as a P value?

text: Doing a poisson regression like this: `model<-glm(y~x*z,family=poisson)` with one predictor being a factor, I would use `anova(model,test="Chisq")` to test the overall effect of the interaction rather than `summary(model)`, which would give me two interaction terms (one for each of two levels in the factor variable that is not the intercept). When I do this, I get an output with five columns: `Df, Deviance, Resid. Df, Resid.Dev, Pr(>Chi)` If I do the same with `test="F"` instead I get an F value, so I assume there is a reason R doesn't report a chi-squared or LRT value. Is this right? I would have thought that I should report a test statistic with my P value. Can I calculate it from the residual deviance and residual degrees of freedom? Should I report the residual deviance and df with my P value and leave it at that?


title: Comparison of subpopulations: do I need normalisation?

text: I have a population of people. Each person has one of three characteristics, say X, Y or Z. I want to compare other characteristics of these people, using the characteristics X, Y and Z as a independent variable. My hypothesis would be something like: "The mean of the salary for people with characteristic X is bigger than the mean of the salary for people with characteristic Y and Z" or "The mean of the number of children of people with characteristic X equals the mean of the number of children of people with characteristic Y and Z". What I know beforehand is that these 3 subpopulations differ greatly in size and that they're independent. My question is: do I need to normalize/standardize/use z-score of the numerical variables, for instance the salary/number of children, so I can do a fair comparison? Would making a random sample of the entire population make anything easier?


title: Draw biplot on nonclassical MDS interpration?

text: How can I draw similar biplot like you can on PCA, so I could compare these two visualization techniques / methods? I got dissimilarities and distances of my data in simple plot, but I would like it to visualize in biplot with arrows that are used in PCA? Is it possible to visualize like that? Or it's not logical to do that? I was following these examples: For PCA: http://www.mathworks.se/help/stats/princomp.html ![enter image description here](http://i.stack.imgur.com/I3z6f.gif) FOR MDS: http://www.mathworks.se/help/stats/mdscale.html ![enter image description here](http://i.stack.imgur.com/1aYOT.png) You can see for PCA example it shows biplot with arrows / lines where you can see direction of all attributes, but in MDS, I only see approximation of distances for dissimilarities (with simple plot). Am I missing there something? I read that these techniques are quite similar, so you should get similar results, but in MDS example I don't get such data that I could use biplot as in PCA (or maybe I don't know how).


title: Good reference on sample autocorrelation?

text: I'm not a statistician but I'm writing my thesis on mathematical finance and I think it would be neat to have a short section about independence of stock returns. I need to get better understanding about some assumptions (see below) and have a good book to cite. I have a model for stock prices $S$ in which the daily ($t_i - t_{i-1}=1$) log-returns $$X_n = \ln\left(\frac{S(t_n)}{S(t_{n-1})}\right), \ \ n=1,...,N$$ are normally distributed with mean $\mu-\sigma^2/2$ and variance $\sigma^2$. The autocorrelation function with lag 1 is $$r = \frac{\mathrm{Cov}(X_1,X_2)}{\mathrm{Var}(X_1)}$$ which I estimate by $$\hat{r} = \frac{(n+1)\sum_{i=1}^{n-1} \bigl(X_i - \bar{X} \bigr)\bigl(X_{i+1} - \bar{X} \bigr)}{n \sum_{i=1}^{n}\bigl(X_i - \bar{X} \bigr)^2} $$ where $$\bar{X} = \frac{1}{n}\sum_{i=1}^N X_i$$ Now I understand that under some **some assumptions** it holds that $$\lim_{n \rightarrow \infty} \sqrt{n}\hat{r} \in N(0,1)$$ I would be very glad if someone could point me towards a good book which I can cite in my thesis and read about these assumptions (I guess it has something to do with the central limit theorem). Thank you in advance! Crossposting at: Mathematics: http://math.stackexchange.com/questions/139408/good-reference-on- sample-autocorrelation Quantative Finance: http://quant.stackexchange.com/questions/3390/good- reference-on-sample-autocorrelation


title: R: generate binomial sample with (pretty) exact probability

text: When I simulate normal data in R, I make sure that the sample have the exact mean and sd of the sampling distribution: `x = scale(rnorm(n))*sd + mean`. I want to do the same for binomial data, making the sample express the near- exact probability that they were generated from. Of course it can't be exact when the probability is continuous and the sample is discrete but something that gets pretty close would be nice `x = rbinom(n, 18, 0.5)` can potentially give samples where an MLE estimate would indicate a probability of p=0.2 or p=0.8 which is pretty far from p=0.5. Purpose: I'm building a bayesian model where I infer a binomial rate from a small sample. To test that the model works, I'd like to simulate well- specified data, in order to diagnose whether a strange inferential result is due to chance in the simulation or in the model.


title: What is probability to get through to the value

text: I have data set of ints with values 1..10. Before any sample is added to the dataset it can be observed for a period of time. During this period X can grow from 1 to 10 ( **Note : it only grows in one directio through all values 1,2,3...10** ). **Once added it cant change**. For example if you are observing next X currently being 3 that means it went throug 1 and 2 and it can stay as 3 till the moment it added to the dataset or it can go higher to 4, 5 .. up to 10. I have count of each X value from the past. eg               X  Y (count of X in the past)     1  12     2  21     3  11     ..     10 100      _I need to calculate probability of X which **hasnt been added to the dataset yet** to grow up to 10 (knowing that X cant go lower than the value it has at this particular moment)._ Obviously if the obeserved X is 10 at the moment it has 100% probability to grow to 10. **EDIT:** (I am trying to explain the process). This is discrete process. Initial value of X is 1 always. Then it can move to 2 then 3 ... up to 10. When value goes up it guaranteed makes it through each consequtive value (ie 10 goes though each - 1,2,3..9 then 10) Process can stop at any time. There is no additional information on why value goes from X1 to X2. The process is completely random in that direction.


title: Machine learning techniques for spam detection, and in general for text classification

text: I am going to configure a system for spam detection. What I have is a dataset of labeled (spam/not-spam) strings containing, mostly, sentences. I have a background in machine learning techniques, but no background in machine learning applied to text. One approach could be creating vectors of extremely high dimensions, in which the features are boolean and each feature represent one possible word (present/not present). Of course such approach is unsatisfactory due the high dimensionality, but also for the extreme sparsity of one feature across the instances and extreme sparsity across the features. What I am asking is just a few pointers (to tutorials, for example) to simple entry-level techniques that may address the aforementioned shortcomings of boolean per-word encoding. Any ideas on which tools may be more suitable for this task? Maybe RapidMiner?


title: How do you find the best feature combination to achieve maximum accuracy?

text: I have a set of 9 possible features that I can extract from an image. Does anyone know of a method that I can employ to utilize the best feature combination to get the maximum possible classification accuracy?


title: Logistic regression using penalized likelihood (lasso?) in Matlab/R

text: I am trying to use `logistic regression` in a scenario where there are very few positives. I'm aware that maximum likelihood suffers from small sample bias. So MATLAB's `glmfit` doesn't work for me. I tried using `firth` regression in R but it simply hangs up my powerful PC (150,000 observations, 9 dummy variables, 1500 positives). MATLAB also has the function `lassoglm` but I'm not sure if it can be used for logistic regression with few positives. Can you please confirm/suggest alternative? Some useful links - here and here.


title: Generate two categorical variables with a certain association?

text: Say, I wish to simulate two categorical variables who are associated with each other like in the table:               prop.table(table(mtcars$gear, mtcars$vs),2)           vs     gear          0          1        3 0.66666667 0.21428571        4 0.11111111 0.71428571        5 0.22222222 0.07142857               # I can generate two new variables which have the same probabilities as the original variables:          gearSim <- sample(x=0:2, size=1000, replace=TRUE, prob = table(mtcars$gear)   )          vsSim <- sample(x=0:1, size=1000, replace=TRUE, prob = table( mtcars$vs)   )               # But they are not (off course) associated with each other:          prop.table(table(gearSim, vsSim),2)                 vsSim     gearSim         0         1           0 0.4766355 0.4365591           1 0.3495327 0.3956989           2 0.1738318 0.1677419      How do I get them to associate with each other, like in the table? (mtcars is in the base dataset)


title: Regression with sequence as predictor

text: **Background:** Obviously, some musical compositions are better than others (measured by charts, saleability, last.fm site, ...). Each composition is, in fact, ordered sequence of discrete values (notes). Putting all together, it looks like supervised problem - sequence of notes (independent "variable") can predict popularity of composition (dependent variable). Hovewer, each sequence is ordered and of different length, what is not usual in classical regression setting, where we have for each object the same number of predictors and set of predictors can be permuted arbitrarily. **Questions:** Is this also supervised problem? If yes, how to predict popularity of composition from sequence of notes. Are there any special techniques, how to solve this problem? If yes, how to efficiently generate set of new sequences of notes with high popularity?


title: Item ranking algorithm

text: I am developing a website with a lot of different items that can be reviewed and rated, although I need a solid ranking algorithm which I can sort after. The ratings will be performed by the user being able to give a score from 1 to 10 where 10 is the best. Now I need the ranking algorithm to be able to calculate a ranking score, so that I can sort the items accordingly. I might introduce a believability / creditability factor of the user giving a score, but I am afraid it would make the ranking calculation more complex. I would appreciate if someone could provide me with both algorithm formulas to calculate the ranking score of an item and with some calculation examples.


title: Unexpected correlation, Method to account for the shared variance

text: **Summary:** While looking at the correlations I noticed a wired, relation. Variables: User comments - A; User community activity- B; Comments on user's wall - C (control) I expected that A and B will be highly correlated, and later on will load into the same latent factor. However, it turns out that B is strongly correlated to C (significant), and moderately correlated to A (significant). That assumes that most of the user's comments A are indeed on his own wall (C). But obviously there is some shared variance between A and B which will be in line with the latent factor. So that user's comments outside his own wall. **Question:** In order to construct the factor how I can find that shared variance between A and B and to account only for it?


title: Significance test for randomly generated task

text: I have conducted an empirical study as part of my master thesis. Each test subject was asked to find 10 cards (one at a time, without knowing what the next card was) among 50 cards divided into 10 piles. Because of different reasons did I chose to shuffle all cards for each subject, and my hypothesis is that this had **_no significant_** affect on the completion time. How I can test if the randomly shuffled piles affected the result? This is how each subjects cards were randomly put, and what completion time it resulted in:                         |--------------------------- PILE ----------------------------|               |   1     2     3     4     5     6     7     8     9    10   |                 |-SUBJECT-|-------------------------------------------------------------|--- SEC --|     |    1    |   1     1     1     0     2     1     2     0     1     1   |    171   |     |    2    |   0     1     1     2     1     2     0     2     1     0   |     78   |     |    3    |   1     1     0     3     0     1     0     1     2     1   |    141   |     |    4    |   1     1     2     1     2     0     0     1     1     1   |    122   |     |    5    |   2     1     2     0     1     0     0     0     4     0   |    227   |     |    6    |   2     0     1     1     0     0     2     2     1     1   |    100   |     |    7    |   1     2     2     0     0     0     2     1     2     0   |    148   |     |    8    |   1     1     1     1     0     1     1     1     1     2   |    113   |     |    9    |   0     2     0     0     2     2     0     1     1     2   |    156   |     |   10    |   1     0     2     1     1     1     0     1     1     2   |    130   |     |   11    |   0     2     3     0     1     2     0     1     1     0   |    112   |     |   12    |   1     1     2     1     1     0     1     1     2     0   |    208   |     |   13    |   1     0     2     0     0     1     1     1     2     2   |    105   |     |   14    |   1     0     1     0     1     2     2     1     2     0   |    216   |     |---------|-------------------------------------------------------------|----------|      _Each number in the table represents how many cards were "hidden" in that pile_ I believe the randomness _could_ affect the result, if for example all 10 cards were put in the first two piles a test subject tried. However, my gut feeling tells me it had no significant affect, but I would like to write that I actually have verified this.


title: IRLS working weights proportional odds

text: According to Reduced-rank vector generalized linear models the parameter estimates are obtained by Fisher-scoring algorithm. In the `VGAM` package the IRLS algorithm is used: $\beta^{(a+1)} = (\sum X_i^\top W_i ^{(a)} X_i)^{-1} (\sum X_i^\top W_i ^{(a)} z_i^{(a)})$ My aim is to reproduce the IRLS steps by extracting the working weights $ W_i ^{(a)}$ for each step $a$ and then fit the parameters for the $(a+1)$-th step. I am failing with the first step because the `vglm` function does not accept the working weights.               pneumo <- transform(pneumo, let = log(exposure.time))     fit1 <- vglm(cbind(normal, mild, severe) ~ let,                  cumulative(), data=pneumo, maxit=1)     weight1 <- weights(fit, type = "working") # working weights          # the weights can't be used to fit the second step     fit2 <- vglm(cbind(normal, mild, severe) ~ let,                  cumulative(), data=pneumo, maxit=1, weights=weight1)     weight2 <- weights(fit, type = "working")      `?vglm` tells me: If the VGAM family function handles multiple responses (q > 1 of them, say) then weights can be a matrix with q columns. Each column matches the respective response. `?weightsvglm` tells me: Working weights are used by the IRLS algorithm. They correspond to the second derivatives of the log-likelihood function with respect to the linear predictors. The working weights correspond to positive- definite weight matrices and are returned in matrix-band form, e.g., the first M columns correspond to the diagonals, etc.


title: Interpretation of Bayesian 95% prediction interval

text: Assume the following bivariate regression model: $$ y_i = \beta x_i + u_i, $$ where $u_i$ is i.i.d $N(0, \sigma^2 = 9)$ for $i = 1,\ldots, n$. Assume a noninformative prior $p(\beta) \propto \text{constant}$, then it can be shown that the posterior pdf for $\beta$ is $$ p(\beta|\mathbf{y}) = (18\pi)^{-\frac{1}{2}}\left(\sum_{i=1}^n x_i^2\right)^{\frac{1}{2}} \exp\left[-\frac{1}{18}\sum_{i=1}^n x_i^2 (\beta - \hat{\beta})^2\right], $$ where $\hat{\beta} = (\sum_{i=1}^n y_ix_i)/(\sum_{i=1}^n x_i^2).$ Now consider the value of $y$ with a given future value of $x$, $x_{n+1}$: $$ y_{n+1} = \beta x_{n+1} + u_{n+1}, $$ where $u_{n+1}$ is i.i.d $N(0, \sigma^2 = 9)$, then we can show that $$ p(y_{n+1}|x_{n+1},\mathbf{y}) = \int_{\beta} p(y_{n+1}|x_{n+1}, \beta, \mathbf{y}) p(\beta|\mathbf{y})d\beta $$ is a normal density with expectation and variance $$ E[y_{n+1}|x_{n+1},\mathbf{y}] = \hat{\beta}x_{n+1},\quad {\rm var}[y_{n+1}|x_{n+1},\mathbf{y}] = \frac{9[x_{n+1}^2 + \sum_{i=1}^n x_i^2]}{\sum_{i=1}^n x_i^2}. $$ Thus the posterior probability density function for $y_{n+1}$, conditional on $x_{n+1}$, is \begin{align} p(y_{n+1}|x_{n+1},\mathbf{y}) = &\left(\frac{18\pi\left[x_{n+1}^2 + \displaystyle{\sum_{i=1}^n x_i^2}\right]}{ \displaystyle{\sum_{i=1}^n x_i^2}}\right)^{-\frac{1}{2}} \\\ & \times \exp\left\\{-\frac{\displaystyle{\sum_{i=1}^n x_i^2}}{18\left(x_{n+1}^2 + \displaystyle{\sum_{i=1}^n x_i^2}\right)}\left(y_{n+1}-\hat{\beta}x_{n+1}\right)^2\right\\} \end{align} * * * Now the question is: Specify a 95% prediction interval for $y_{n+1}$ and carefully interpret it. What aspect(s) of the data generating process does the interval fail to accommodate our uncertainty about? * * * I'm not exactly too sure on how to answer the question but here's my attempt: So essentially we need to find some $a$ and $b$ such that $P(a < y_{n+1} < b) = \int_{a}^b p(y_{n+1}|x_{n+1}, \mathbf{y}) dy_{n+1} = 95\%$ Now we know that $y_{n+1}|x_{n+1}, \mathbf{y} \sim N(m, v^2)$ where $m = E[y_{n+1}|x_{n+1},\mathbf{y}]$ and $v^2 = var[y_{n+1}|x_{n+1},\mathbf{y}]$, Hence: $$\frac{y_{n+1}-m}{v} \sim N(0,1)$$ $$P(-1.96 < \frac{y_{n+1}-m}{v} < 1.96) = 95\%$$ $$P(-1.96v+m < y_{n+1} < 1.96v+m) = 95\%$$ Now because we are conditioning on $x_{n+1}$ and looking at the expression for $v$ and $m$, we see that both $v$ and $m$ are known values. So we can take $a = -1.96v+m$ and $b = 1.96v+m$. i.e., we can select many other possibilities of $a$ and $b$ which yield a probability of $95\%$... but how does this relate to answering the part of the question that asks what aspects of the data generating process this interval fails to accommodate?


title: Estimation and comparison of flowering curves

text: I am a research scholar and monitoring phenological events of timber line at Himalayan region from past 4 years. During data analysis I found a research paper Estimation and comparison of flowering curve similar to my work. In this paper bbmle package was used and five parameters ($β_0, ... β_4$) describe: (i) the height; (ii) the peak date; (iii) the range; (iv) the symmetry; and (v) the peakedness of the regression curve was calculated. I also read the appendix table and followed the code to estimate these parameter but as a newbie I failed to calculate. Appendix               library(bbmle) ## depends R(≥ 2.0.0)     ## Collect data together as a data frame,     fdat <- data.frame(Days=c(212:238,250:271),     Count=c(0,2,2,6,10,18,29,39,59,75,104,130,     145,169,193,209,216,227,231,214,212,226,242,     225,214,202,211,104,90,70,55,52,45,38,29,22,     14,15,14,13,11,6,5,4,3,2,1,1,0))     ## Define function     GESN <- function(x,b0,b1,b2,b3,b4) {     exp(b0-abs(((x-b1)/(b2∗(1+b3∗(sign(x-b1))))))ˆ(b4))     }     ## Get reasonable starting values     startvals <- list(b0=log(250),b1=230,b2=15,b3=0,b4=2)      Please help me to find these parameters as they are mentioned in figure number 5 (link 1).


title: Counterpart of Chow test for Least absolute deviations (regression analysis)

text: Does is exist counterpart of Chow test for Least absolute deviations method in analytical or not version ? In least absolute deviations we minimize sum of absolute values of errors.


title: R function which uses innovations algorithm?

text: I can't seem to find much info on the following: I have a dataset _D_ at time _t_ which I use to fit an ARIMA model. I forecast the value of the time series at time _t_ \+ 1. Now, when I'm in _t_ \+ 1, I would like to predict the value of my time series at _t_ \+ 2 using data up until time _t_ \+ 1. However, I don't want to refit the whole model. Is there a function in R which does this? I would have thought this is a common thing to do with time series analysis.


title: What does a post hoc power analysis say about the significance of the study results?

text: I orginially posted this over at AskMetafilter, and a commenter suggested I ask it here. I work for a dietary supplement company that also makes skin care products, and some of those skin care products are tested clinically. Now they are talking about repeating some of the clinical tests in another region of the world in which the products will be sold because the marketing department thinks that would be good. The question came up: How do you decide how many subjects to include? The answer, I learn, is statistical power. You include at least as many subjects as are needed to give you an appropriate power, say 80%, given your chosen significance level alpha and expected effect size, for the type of statistical test you are performing. What if p is small, lower than alpha, but power is also low? Does that invalidate the results? For example, one of the studies they want to repeat used a two-tailed paired t test to compare before and after treatment means for a measurement. Alpha was 0.05 and population size was 30. The pooled SD for the two data sets was 0.688. After the fact, I calculated a Cohen's d of 0.494. All this gives a power of less than 50%, which means the study was underpowered. At the same time, p was 0.000004. I can tell the people at work that, when the study is repeated, we are going to need more subjects or we risk missing the effect that we saw in the first trial, but what can we conclude about that first trial? Power was low, but p was much lower than alpha. Are the results no good? Or can we still trust that p? Or, what is also possible, am I completely confused and all this doesn't work the way I think? Thanks for any help you can provide!


title: How to choose the right distance matrix for clustering?

text: I am attempting simple Ward type clustering. However, the R package is proving several choices to use for the distance matrix. I am wondering how I am supposed to determine the right distance matrix method. Are there any generally acceptable criteria for specific sets of problems?


title: JavaScript/PHP libraries for calculating common

text: I am looking for a `JavaScript` or a `PHP` library to calculate common business/statistical metrics (e.g. `Total`, `Avg`, `YTD`, `YoY`, `MoM`, `CAGR` growth, standard deviation, regression line coordinates) from simple arrays/objects with a structure such as: ![enter image description here](http://i.stack.imgur.com/hhSyz.png) The idea is that I would just provide the arrays with some basic parameters (e.g. time period I'm interested in) and would get all the numbers crunched automatically in return. **Can anybody point me to a library they've used and trust from experience?**


title: Statistical testing against empirical density function

text: I have a feeling this is a somewhat common problem but reading about tests for normal distributions hasn't really helped me as their assumptions didn't seem to fit mine. Consider the following case:   * I have data (a lot) that is somewhat normally distributed but quite leptokurtic.    * I collect new data samples and want to test if the new data is from the same distribution. What test should I use here? To give some context the testing would be done online (e.g. repeated after 10, 20 etc. samples) and I am interested to find out how many samples I need to reject the null hypothesis (where appropriate). Edit: Alternatively, is there some way to forgo the distribution/goodness-of- fit testing and evaluate new incoming data according to its likelihood? I.e. data point 1 is in the 90-percentile of the distribution, point 2 in the 95% and so on. Eventually, I want a test to tell me that the behaviour is likely according to the previously collected data or not. Thanks!


title: Compare skewness of many distributions with few observations

text: I have a dataset with page view data for about 500,000 users, divided into two groups. Each user can visit up to 5 pages, each as many or as few times as they want. So for each user, I have the distribution of number of visits to each page. I would like to compare the 'average skew' in distribution between the two groups. Roughly, users in one group are more likely to have distributions that look like this {3,0,0,0,0} and users in the other group are more likely to have distributions like this {1,1,0,1,0}. How can I compare the degree of skewness between the two groups? I thought of using the average Gini coefficient or entropy for each group, but each user has so few observations. How can I do this?


title: How to change data between wide and long formats in R?

text: You can have data in wide format or in long format. This is quite an important thing, as the useable methods are different, depending on the format. I know you have to work with `melt()` and `cast()` from the reshape package, but there seems some things that I don't get. Can someone give me a short overview how you do this?


title: Anova for a factor with only 2 levels

text: I run an experiment with some control condition and some test conditions, in a repeated measures design, and I want to know whether for my dependent variable there is a significant difference between the two conditions. Is it enough to test for the main effect of the repeated measures one-way ANOVA? Or do I need a post-hoc test? What about if the main effect is not significant? Can I conclude that there is not significance in the differences between the two conditions? Is ANOVA the best method in this case? Or do you suggest a different analysis method?


title: How can I use Park test for heteroskedasticity?

text: For example, I have a model: $$Y= \beta_0 + \beta_1*X_1 + \beta_2*X_2 + \beta_3*X_3 + \beta_4*X_4 + \epsilon_1$$ I suppose that $X_1$ may be the factor causing heteroskedasticity, and I regress Y on $X_1$ $$Y= \alpha_0 + \alpha_1*X_1 + \epsilon_2 $$ Now, which residual, $\epsilon_1$ or $\epsilon_2$ could be used to regress on $X_1$ to test for heteroskedasticity ?


title: Stability of cross-validation in Bayesian models

text: I'm fitting a Bayesian HLM in JAGS using k-fold cross-validation (k=5). I'd like to know whether estimates of parameter $\beta$ are stable across all folds. What's the best way to do this? One idea is to find the differences of the posteriors of $\beta$ and to see if 0 is in the 95% CI of the difference. In other words, is 0 in the 95% interval of $\beta_{k=1}-\beta_{k=2}$ (and then repeat for all pairs of folds). Another idea is to treat the posteriors from each fold as different MCMC chains, and to compute Gelman's $\hat{R}$ (Potential Scale Reduction Factor) across these pseudo-chains. Is one of these preferable, and are there alternatives?


title: Limit of a convolution and sum of distribution functions

text: I need to prove an induction step. $X_i$ are independently distributed with the distribution function $1-F_i=x^{-\alpha}L_{i}(x)$ where $\alpha \geq 0$ and $L_{i}(x)$ is regularly varying (If the limit $g(a)=\lim\limits_{x\rightarrow\infty}\frac{L(ax)}{L(x)}$ is finite and nonzero for $a >0$, then L is regularly varying). $\lim\limits_{x\rightarrow \infty} \frac{P(X_1+...+X_n>x) }{P(X_1 > x)+...+P(X_n>x)} = 1$ is true. Now we have to show: $\lim\limits_{x\rightarrow \infty} \frac{P(X_1+...+X_{n+1}>x) }{P(X_1 > x)+...+P(X_{n+1}>x)} = 1.$ How do we show this?


title: What is a method to calculate precisely $P(Y \geq X, Y\leq Z)$, given three independent random variables $X, Y$, and $Z$

text: For three independent normally distributed continuous random variables X, Y, and Z (each with its own mean and standard deviation), I need a way to calculate $P(Y \geq X, Y \leq Z)$ I know that I can do this by the following: $P(Y \geq X, Y \leq Z) = P(Y \geq X) \cdot P(Y \leq Z | Y \geq X)$ I am able to calculate $P(Y \geq X)$ using the following relation: $P(Y \geq X) = P(Y - X \geq 0)$ However, I'm having trouble calculating precisely $P(Y \leq Z | Y \geq X)$


title: Why t-stat is normally reported but not p-value?

text: I just wonder why t-stat is normally reported with an asterisk to indicate significance but not p-value.


title: logistic regression. How to get dual function?

text: Given pairs $(x_i, y_i), x_i \in R^n , y_i \in R$ we want to solve minimization problem (logistic regression):$\min \frac{1}{2} ||w||^2 + \sum_i^{i=m}\log(1+\exp(-y w\cdot x_i))$. How to do that? I know the dual form is: $ \min_{\alpha} D(\alpha)= \frac{1}{2}\sum_{i, j}\alpha_i \alpha_j y_i y_j x_i \cdot x_j + \sum_i\alpha_i \log(\alpha_i) + (C-\alpha_i)\log(C-\alpha_i)$, subject to: $0\le \alpha_i\le C$. How to get dual?


title: Normalize a periodic parameter

text: I am using inverse modelling software (PEST) to estimate a periodic parameter for the direction of anisotropy, $\hat{\theta}$, which is somewhere in $[0^{\circ}, 180^{\circ})$ (i.e., has a wavelength of 180°). With PEST, parameters must be defined with lower/upper bounds and an initial guess. The inverse model will then seek to minimize the objective function by adjusting the parameter up/down within the bounds. The problem is that my parameter is periodic, and this is difficult for inverse modelling software to solve. If the ideal $\theta$ is near 90°, then there isn't any issue solving for this parameter, since it is far from 0° or 180°, and it behaves somewhat linearly. However, if the ideal $\theta$ is near 0° or 180°, then the program will have difficulty approaching this answer to either the lower or upper bounds, since it does not understand periodicity. Is there a technique to transform a periodic parameter, such that it behaves linear? I already have both pre- and post-processing routines in Python that can be adapted to transform this parameter.


title: Is it possible to measure the independent variable with part of the dependent variable

text: I have Beta as my independent variable and Economic value added (EVA) as my dependent variable. To calculate EVA I need to use Cost of capital and to calculate that I have to use Beta, so is it possible to use EVA as the dependent variable. thank you > EVA= Net Operating Profit After Taxes (NOPAT) - (Capital * Cost of Capital)


title: Flat, conjugate, and hyper- priors. What are they?

text: I am currently reading about Bayesian Methods in Computation Molecular Evolution by Yang. In section 5.2 it talks about priors, and specifically Non- informative/flat/vague/diffuse, conjugate, and hyper- priors. This might be asking for an oversimplification but, could someone explain simply the difference between these types of priors and how that affects the outcome of an analysis/decisions I would make during the process of a Bayesian analysis? (I'm not a statistician and I am just starting out on the road to learning Bayesian analyses so the more it is in layman terms the better)


title: How to extract long run and short run coefficients from ARDL (UECM) estimates?

text: I have estimated ARDL(UECM) in eviews but I dont know how to specify or extract the long run an short run estimates/coefficienst? what is the standard procedure to do so?


title: Unimodality of continuous univariate distributions

text: What are the different methods to prove that a continuous univariate distribution is unimodal? Some of them are available at Wikipedia but not much detail is given. Solution or suggestion of any continuous univariate distribution is welcome.


title: Friedman test to report on occurrence of categories in multiple articles

text: To investigate reasons for certain behavior, we have defined three categories for these reasons, and have determined how often these were mentioned across 15 (very different) scientific articles. These articles were selected by systematic review of titles and abstracts. They span the whole (broad) medical field and are very heterogeneous in methodology and setting. I want to report on which category of reasons is the most 'important' category for the investigated behavior. **I.e. if one category is mentioned significantly more often than the others.** Our data is as follows, we have a table with the 15 articles for rows, and the categories as columns. Each cell shows how often an article has mentioned a category in absolute numbers. Like so:               Art C1  C2  C3     1   10  3   2     2    7  18  1     3    5  2   7     etc...      After much consideration my idea would be to use the Friedman test to determine if one category is mentioned significantly more often than the others, or to show the mean rank of one category is significantly higher than the others. First question, assuming the data, would this be the right test? Second question, should I use the actual occurrences in the articles as data in the table? EG:               Art C1  C2  C3     1   10  3   2     2    7  18  1     3    5  2   7     etc...      Or ranked?               Art C1 C2  C3     1    3  2   1     2    2  3   1     3    2  1   3      And then I would report the mean rank and the result of the Friedman test e.g. the probability that the difference in mean rank is based on chance" So two questions:   1. Friedman, good idea?    2. List the data with the absolute numbers, or ranked?


title: Use of Weibull model or Weibull distribution

text: I have a basic question regarding use of the terms "Weibull model" and "Weibull distribution". (I take Weibull as an example, but the question could apply to any distribution.) I read in books or research papers frequent use of "Weibull model" or "Weibull distribution", so my questions are 1) Are the meanings of Weibull model and Weibull distribution the same? 2) If not so, in what situations do we use Weibull model or Weibull distribution?


How can I estimate 95% confidence intervals using profiling for parameters estimated by maximising a log-likelihood function using optim in R? I know I can asymptotically estimate the covariance matrix by inverting the hessian, but I am concerned that my data do not meet the assumptions required for this method to be valid. I'd prefer to estimate confidence intervals using some other method. Is the profile likelihood method appropriate, as discussed in Stryhn and Christensen, and in Venables and Ripley's MASS book, §8.4, pp. 220-221? If so, are there any packages that can help me do this in R? If not, what would the pseudo code for such a method look like?


title: Stratified classification with random forests (or another classifier)

text: So, I've got a matrix of about 60 x 1000. I'm looking at it as 60 objects with 1000 features; the 60 objects are grouped into 3 classes (a,b,c). 20 objects in each class, and we know the true classification. I'd like to do supervised learning on this set of 60 training examples, and I'm interested both in classifier accuracy (and related metrics) as well as feature selection on the 1000 features. First, how's my nomenclature? Now the real question: I could throw random forests on it as stated, or any number of other classifiers. But there's a subtlety - I really only care about differentiating class c from classes a and b. I could pool classes a and b, but is there a good way to use the _a priori_ knowledge that all non-c objects likely form two distinct clusters? I'd prefer to use random forests or a variant thereof, since it's been shown to be effective on data similar to mine. But I could be convinced to try some other approaches.


title: Cross validation for multivariate imputation

text: I am currently using the mice: Multivariate Imputation by Chained Equations in R (JSS 2011 45(3)) package. Consider the following example. I am using Sites B to Z and mice() to help infill missing data for Site A. The observations for Site A span the period 2002 to 2005 (there are some missing data in this period). I have cross validated the infilled data over the period 2002 to 2005. Because Sites B to Z cover the period 1995 to 2005, I would like to extend the record of Site A to also cover a similar period. My question is: how I can check that the imputted values for Site A over the period 1995 to 2002 are reasonable? Obviously I cannot cross validate any data over this period. Thanks for your help in advance. Matt.


title: Missing cells with Type III SS

text: I have a linear model of the form y ~ x + z + x:z. I have unbalanced data and further have a few missing cells (data would be unbalanced even without those missing cells). My understanding of the different Types of SS comes mostly from how you can manually obtain the SS using different model comparisons in R. For example, to obtain the Type III SS associated with z, you would compare a model including x and x:z to a model including x, z, and x:z. More details on this procedure can be found in this answer. I understand this model comparison method (for Type III SS) as throwing away any error that cannot be unambiguously assigned to only one term. I can perform this method on data with missing cells, and the Anova function in the car package in R will compute the anova table for missing cells, but I have read some places (e.g. here) that Type III SS should not be used for missing cells data.   1. Is it really wrong to use Type III SS with missing cells, or just not normally appropriate, and why?   2. How is the interpretation of an anova table using Type III SS with missing cells different from one without missing cells (but still having unbalanced data)? ## EDIT I performed the following simulation in R comparing type I error rates for different patterns of missing observations in synthetic data (balanced, observations missing completely at random, cells missing completely at random, and both cells and observation missing at random, all with the same number of observations). From the results of this I don't see a problem with using Type III SS on data with missing cells. The type I error rates don't seem to be violated for any of the data sets. For one run I got the following observed Type I error rates:                                      intercept       effect A        effect B      balanced               .056            .058            .039     missing obs            .052            .052            .039      missing cells          .050            .055            .044        missing obs and cells  .056            .050            .049      In a different setup I added a main effect for B and the neither the power nor the observed alpha level seem to be affected by missing cells. Am I not understanding something important here? Here is the code I used:               library(car) # for Type III SS          # Data properties (simulation may not work if changing these)     nLevelsA <- 5     nLevelsB <- 5     nreps <- 5     factA <- rep(1:5, times = nLevelsA*nreps)     factB <- rep(1:5, each = nLevelsB*nreps)     cell <- factor(paste(factA, factB, sep = ""))          nSims <- 1000 # Number of repetitions of the simulation     # Empty output objects     mod1P <- matrix(rep(NA, times = nSims * 3), ncol = 3)     mod2P <- matrix(rep(NA, times = nSims * 3), ncol = 3)     mod3P <- matrix(rep(NA, times = nSims * 3), ncol = 3)     mod4P <- matrix(rep(NA, times = nSims * 3), ncol = 3)          # Create four datasets, all based on the same values, but different patterns of cell counts     for(i in 1:nSims){         myData0 <- data.frame(y = rnorm(nLevelsA*nLevelsB*nreps, 0, 1), A = factA, B = factB, cell = cell)         # Use this to add a main effect for B         # for(m in 1:nLevelsB){         #       myData0[which(myData0$B == m), 1] <- myData0[which(myData0$B == m), 1] + m         # }              # Randomly remove one observation from each cell         myData1 <- myData0         for(j in 1:length(levels(cell))){             # randomly pick one row from each cell             rowNums <- 1:nrow(myData1)             rowToDrop <- sample(rowNums[which(myData1$cell == levels(cell)[j])], 1)             myData1 <- myData1[-rowToDrop, ]         }              # Randomly remove observations without regard to cell         myData2 <- myData1[sample(1:nrow(myData0), nrow(myData1), replace = FALSE), ]               # Randomly empty cells from balanced data         myData3 <- myData0         cellToDrop <- sample(levels(myData0$cell), 5)          for(k in 1:length(cellToDrop)){           myData3 <- myData3[which(myData3$cell != cellToDrop[k]), ]         }              # Randomly empty cells from unbalanced data         myData4 <- myData0         cellToDrop <- sample(levels(myData0$cell), 3)          for(l in 1:length(cellToDrop)){           myData4 <- myData4[which(myData4$cell != cellToDrop[l]), ]         }         myData4 <- myData4[sample(1:nrow(myData4), nrow(myData1), replace = FALSE), ]                    # nrow(myData0  # basis to start with to get the same number observations in each other dataset         # nrow(myData1) # balanced         # nrow(myData2) # observations missing at random         # nrow(myData3) # cells missing at random         # nrow(myData4) # cells and observations missing at random         #          #          # xtabs(~ A + B, data = myData0)         # xtabs(~ A + B, data = myData1)         # xtabs(~ A + B, data = myData2)         # xtabs(~ A + B, data = myData3)         # xtabs(~ A + B, data = myData4)                   mod1 <- lm(y ~ A + B, data = myData1)         mod2 <- lm(y ~ A + B, data = myData2)         mod3 <- lm(y ~ A + B, data = myData3)         mod4 <- lm(y ~ A + B, data = myData4)              # P-values for intercept, factor A, and factor B         mod1P[i, ] <- Anova(mod1, type = 3)$'Pr(>F)'[1:3]          mod2P[i, ] <- Anova(mod2, type = 3)$'Pr(>F)'[1:3]         mod3P[i, ] <- Anova(mod3, type = 3)$'Pr(>F)'[1:3]          mod4P[i, ] <- Anova(mod4, type = 3)$'Pr(>F)'[1:3]          }          # Count how many times a significant test result is found     p1 <- mod1P <= .05     p2 <- mod2P <= .05     p3 <- mod3P <= .05     p4 <- mod4P <= .05          # Get proportion of times a significant test result is found     pobs1 <- apply(p1, MARGIN = 2, FUN = sum) / nSims     pobs2 <- apply(p2, MARGIN = 2, FUN = sum) / nSims     pobs3 <- apply(p3, MARGIN = 2, FUN = sum) / nSims     pobs4 <- apply(p4, MARGIN = 2, FUN = sum) / nSims          # Proportion of experiments observed as significant at 0.05 level.     # Intercept, then factor A, then factor B     pobs1 # balanced     pobs2 # observations missing at random     pobs3 # cells missing at random     pobs4 # cells and observations missing at random


title: Either quadratic or interaction term is significant in isolation, but neither are together

text: As part of an assignment, I had to fit a model with two predictor variables. I then had to draw a plot of the models' residuals against one of the included predictors and make changes based on that. The plot showed a curvilinear trend and so I included a quadratic term for that predictor. The new model showed the quadratic term to be significant. All good so far. However, the data suggest that an interaction makes sense, too. Adding an interaction term to the original model also 'fixed' the curvilinear trend and was also significant when added to the model (without the quadratic term). The problem is, when both the quadratic and the interaction term are added to the model, one of them is not significant. Which term (the quadratic or the interaction) should I include in the model and why?


title: How to calculate the scale parameter of a Cauchy random variable

text: Let $(X_n)$ be iid random variables and suppose they have mean 0 and follow Cauchy distribution. I know I can set the location parameter to 0. My question is how to find the corresponding scale parameter. Thank you!


title: What are some illustrative applications of empirical likelihood?

text: I have heard of Owen's empirical likelihood, but until recently paid it no heed until I came across it in a paper of interest (Mengersen et al. 2012). In my efforts to understand it, I have gleaned that the likelihood of the observed data is represented as $$L = \prod_i p_i = \prod_i P(X_i=x) = \prod_i P(X_i \le x) - P(X_i \lt x)$$ , where $\sum_i p_i = 1$ and $p_i > 0$. However, I have been unable to make the mental leap connecting this representation with how it can be used to make inferences about observations. Perhaps I am too rooted in thinking of a likelihood w.r.t parameters of a model? Regardless, I have been searching Google Scholar for some paper employing empirical likelihood that would help me internalize the concept... to no avail. Obviously, there is Art Owen's book on Empirical Likelihood, but Google Books leaves out all the yummy bits and I'm still in the slow process of getting an inter-library loan. In the meantime, can somebody kindly point me to papers and documents that clearly illustrate the premise of empirical likelihood and how it is employed? An illustrative description of EL itself would also be welcome!


title: Computing effect for each firm for one of the explanatory variable in the panel data (Stata)

text: Suppose I have panel data with firms $i=1,...,N$ and year $t=1,...,T$. Also, assume that I have the model like below: $y_{it}=\beta x_{1it}+\alpha_{i}x_{2it}+v_{i}+u_{it}$ My simple question is how to estimate $\alpha_{i}$. Note that $x_1$ and $x_2$ are continuous variables. In Stata, I did like below:               xtset firm year     reg y x1 i.firm#c.x2 i.firm          Is this the correct way to estimate $\alpha_{i}$?


title: Understanding ACF - PCF Plots of ARIMA model

text: ![enter image description here](http://i.stack.imgur.com/eG6l1.jpg) I used `auto.arima` function of R software to get following order of arima model: ARIMA(2,0,2). Then I used the following commands on my data to generate ACF and PACF plots.               par(mfrow=c(1,2))     acf(mydata[,1],main="")     pacf(mydata[,1],main="")      My question is how we can observe the values of p=2, and q=2 from the given plots. I tried to understand from the book page: https://www.otexts.org/fpp/8/5 , but it mentions about arima(p,d,0) or arima(0,d,q), and i need to understand arima(p,0,q) case.


title: Algorithm for scoring co-varying traits

text: I am sure this has been done, but I can't find quite the right approach. **EDIT:** Trying to explain better. The rows of colored boxes below are _columns_ of molecular sequence data -- positions in a protein. The circles represent individuals (species, etc) which have that amino acid at that position. An actual data set consists of many of these positions for the same set of taxa. I want to quantify how well each position (set of traits, A1,A2, etc) is correlated with the trait shared by the species (black or white color of dots, e.g., indicating habitat or other external factor). There is an underlying relationship between the species which is going to drive the major pattern of similarity, and I am trying to find where things differ from that pattern, being instead correlated to the trait (black/white). In the diagrams and question below tried to reduce this question down to what seems to be the underlying statistical question: generate a score for how well correlated the colored squares are with the black and white circles, without regard to the actual color of the square. (The correlation could be the opposite one position over, as in B1 and B2, but that is still strongly linked to the black/white trait.) * * * **Original description below:** I am trying to find a way to score a suite of traits (colored boxes in figure) by how well they coincide with another categorical trait (black and gray circles below). The categorical traits are shown as binary, but would ideally scale to 3 or more. One approach would be to generate expected frequencies based on the number of traits and the number of categories, and compare the actual frequencies for each. To get expectations, for example in scenario A1 one state is present in 1/3rd of the cases, its expectation would be the number of each category times that probability. ![enter image description here](http://i.stack.imgur.com/MneCU.png) My questions are (1) whether this is an accepted way to generate the expected probabilities and (2) what the metric should be for assessing them. I want to have a value for each black/white trait state indicating how strongly it is linked to the colored state. There is no way to get a zero observed in this situation, because one or the other trait will be found in each category. (For example, in B1 and B2, they should both score as high as possible because the traits sort perfectly to the categories, even though they are in opposite order.) Scenario A4 and B4 should score as low as possible and A1 should also have the highest score. To put it another way, I am not testing direct correlations ("Do females prefer Democrats?") but rather for each potential preference, which ones separate by category ("What topics of conversation show male/female differences?"). The expectation has to take into account that the associations should score high with either association. Another thing that might/might not be relevant for summing the trait scores is that, they are not independent: If A and B are traits, then there is an expectation that individual samples (all the 0 and 1 circles) will have an underlying pattern of similarity that is not driven by one- or zero-ness. The point of this metric is to sort out which sets of traits show the most influence of category, even if the majority are driven by some other relationship. Currently I am assuming I would compare the expected total and the observed for each category using a Chi-squared calculation of `(O - E)^2/E` and then sum them for each trait? Does this look like a standard test? I would greatly appreciate any links or suggestions.


title: Calculate the correlation between a principal component and a feature of the original data set in R?

text: I am interested in seeing the correlation between a particular principal component and a particular independent variable in my 'original' data set, that is, I'd like to calculate $\rho_{Y_{i},X_{k}} = \frac{e_{i,k}\sqrt{\lambda}}{\sqrt{\sigma_{k,k}}}$ where $Y_{i}$ is the $i$th principal component and $X_{k}$ is the $k$th independent variable. I have been looking around online and can't find a package in $R$ that has this function. Has anyone seen this function in any package?


title: Dynamic Time Warping and normalization

text: I'm using Dynamic Time Warping to match a "query" and a "template" curve and having reasonable success thus far, but I have some basic questions:   1. I'm assessing a "match" by assessing whether the DTW result is less than some threshold value that I come up with heuristically. Is this the general approach to determining a "match" using DTW? If not, please explain... Assuming the answer to (1) is "yes", then I'm confused, since the DTW result is quite sensitive to a) the difference in amplitudes of the curves and b) the length of the query vector and the length of the "template" vector. I am using a symmetric step function, so for (b) I'm normalizing my DTW result by dividing by M+N (width + height of DTW matrix). This seems to be _somewhat_ effective, but it seems that it would penalize DTW matches that are further from the diagonal (i.e., which have a longer path through the DTW matrix). Which seems kind of arbitrary for a "normalization" approach. Dividing by the number of steps through the matrix seems to make intuitive sense, but that doesn't appear to be the way to do it according to the literature.   2. So is there a better way to adjust the DTW result for the size of the query and template vectors?   3. Finally, how do I normalize the DTW result for the difference in amplitudes between the query and the template vectors? As it is, given the lack of reliable normalization techniques (or my lack of understanding), there seems to be a lot of manual effort involved in working with the sample data to identify the best threshold level for defining a "match". Am I missing something?


title: Degrees of Freedom for Ridge regression without knowing the Ridge Parameter?

text: There is a very nice post here that gives a neat solution to the problem of finding the ridge parameter when the degrees of freedom are known: How to calculate regularization parameter in ridge regression given degrees of freedom and input matrix? My question is how can you know what the degrees of freedom are before knowing the ridge parameter value is? ( I have added a comment to the other thread but its quite old so thought it best to create a new topic).


title: Muthén's robust weighted least squares factoring method for binary items...in R?

text: I am working on an exploratory factor analysis of dichotomous items. I've found this post extremely helpful. However, Mplus seems to be a standard that journals in my field are expecting, and I would love to replicate Mplus results in R. I've gotten my results to come very close to Mplus, using the fa function on a polychoric correlation matrix (created with the hetcor function) and the geominQ rotation. It seems that the major difference between the fa function and Mplus is that the latter uses a robust weighted least squares factoring method (WLSMV - a diagonal weight matrix), whereas the former uses a regular weighted least squares (WLS) factoring method. Has anyone managed to use R to replicate Mplus factor analysis for binary items? Has anyone figured out how to implement WLSMV in R? Thanks!


title: Deriving posterior distribution with a normal likelihood and a normal-gamma prior

text: How can I derive the posterior distribution of a likelihood function $N({\bf{t}}|\phi(x)^{T}{\bf{w}}, \beta^{-1})$ and a prior $N({\bf{w}}|{\bf{m_{0}}}, \beta^{-1}{\bf{S_{0}}})Gam(\beta, a_{0},b_{0})$? I'd like to use something similar to the general relations for conditional distributions in the normal case: $$p(x|y) \propto p(y|x)p(x) = N(y|Ax+b, L^{-1})N(x|\mu, \Lambda^{-1})$$ which gives the result of: $$p(x|y) \propto N(x|\Sigma\\{A^{T}L(y-b)+\Lambda \mu\\}, \Sigma)$$ However, this doesn't work with a normal-gamma distribution as a prior, and I'm not sure if there is such equivalent relation. Any thoughts about how to approach this?


title: Drop a predefined percentage of outliers

text: I have a dataset that follows something between a power and exponential law. I'm not happy with the IQR method of detection of outliers because on small sets of data (<50-100), it does not give you an idea of the percentage of outliers that have been droped.   I thought of an iterative method that would drop on each iteration the most outlying number until the the desired percentage is reached. I could run this method only when the IQR method drops more than the desired percentage.   Does it make sense and is there a standard way of dealing with this problem ?


title: Estimating HMM parameters from a sequence with missing observations

text: The BW algorithm to estimate HMM parameters works on a consecutive sequence of observations. But what should be done if only a partial sequence is available? From another point of view, the observable may change its internal state without producing any output. Can we say that a continuous-time HMM is less sensitive than a discrete-time one to this missing information when estimating the parameters? P.S. I have faced the problem with a discrete-time HMM. I don't know whether we can estimate the parameters of a continuous-time one efficiently, or not.


title: On ARIMA and its implementation

text: I am interested in ARMA/ARIMA models. First off, is there any material online about how to identify such models? It is done via acf and pacf, but I find the material I found online somewhat misleading. I would like to have a set of rule of thumbs I can use for that. Secondly, do you know of any C++ library, equivalent to R for time series? Finally, is it possible to include dependency across different time series? For example, suppose we know that series $\\{x\\}_n$ depends upon series $\\{y\\}_n$. How would I include this when using ARIMA models?


title: Kendall Tau or Spearman's rho?

text: In which cases should one prefer the one over the other? I found someone who claims an advantage for Kendall, for pedagogical reasons, are there other reasons?


title: Random forests for multivariate regression

text: I have a multi-output regression problem with $d_x$ input features and $d_y$ outputs. The outputs have a complex, non-linear correlation structure. I'd like to use random forests to do the regression. As far as I can tell, random forests for regression only work with a single output, so I would have to train $d_y$ random forests - one for each output. This ignores their correlations. Is there an extension to random forests that takes output correlations into account? Maybe something like Gaussian process regression for multi-task learning.


title: Standardizing variables for k-means?

text: I only have two variables and they are on the same scale. However, the variance corresponding to the first variable is approximately 0.609, whereas for the second variable is 0.154. So my question is should I standardize the observations (subtracting off the mean and dividing by the standard deviation) or leave it as is? My professor said if the variables are on the same scale (i.e. measuring the same thing) standardization is not critical, and you could do it if you want to, but it will not have a big effect. However, I think it does have a somewhat significant impact if the variances differ significantly. Though I'm not sure if still is best to standardize or not. Would standardizing cause me somehow lose what the original data was trying to convey or is it the best course of action when applying K-means?


title: What does it mean that random effects are highly correlated?

text: What does it mean when two random effects are highly or perfectly correlated?   That is, in R when you call summary on a mixed model object, under "Random effects" "corr" is 1 or -1.               summary(model.lmer)      Random effects:     Groups   Name                    Variance   Std.Dev.  Corr                      popu     (Intercept)             2.5714e-01 0.5070912                                     amdclipped              4.2505e-04 0.0206167  1.000                              nutrientHigh            7.5078e-02 0.2740042  1.000  1.000                       amdclipped:nutrientHigh 6.5322e-06 0.0025558 -1.000 -1.000 -1.000      I know this is bad and indicates that the random effects part of the model is too complex, but I'm trying to understand   * 1)what is doing on statistically    * 2)what is going on practically with the structure of the response variables. **Example** Here is an example based on "GLMMs in action: gene-by-environment interaction in total fruit production of wild populations of Arabidopsis thaliana" by Bolker et al Download data               download.file(url = "http://glmm.wdfiles.com/local--files/trondheim/Banta_TotalFruits.csv", destfile = "Banta_TotalFruits.csv")     dat.tf <- read.csv("Banta_TotalFruits.csv", header = TRUE)      Set up factors               dat.tf <- transform(dat.tf,X=factor(X),gen=factor(gen),rack=factor(rack),amd=factor(amd,levels=c("unclipped","clipped")),nutrient=factor(nutrient,label=c("Low","High")))      Modeling log(total.fruits+1) with "population" (popu) as random effect               model.lmer <- lmer(log(total.fruits+1) ~ nutrient*amd + (amd*nutrient|popu), data= dat.tf)      Accessing the Correlation matrix of the random effects show that everything is perfectly correlated               attr(VarCorr(model.lmer)$popu,"correlation")                                   (Intercept) amdclipped nutrientHigh amdclipped:nutrientHigh     (Intercept)                       1          1            1                      -1     amdclipped                        1          1            1                      -1     nutrientHigh                      1          1            1                      -1     amdclipped:nutrientHigh          -1         -1           -1                       1      I understand that these are the correlation coefficients of two vectors of random effects coefficients, such as               cor(ranef(model.lmer)$popu$amdclipped, ranef(model.lmer)$popu$nutrientHigh)      Does a high correlation mean that the two random effects contain redundant information? Is this analogous to multicollinearity in multiple regression when a model with highly correlated predictors should be simplified?


title: Subjects:condition interaction random effect in a growth model

text: I'm investigating the effect of 'Condition' (3 levels: Quiet, Intelligible, Unintelligible) on pupil response over time (intercept, linear, cubic, quadratic, quartic and quintic terms). When I use this as my mixed growth curve model;               m.full <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) * ConditionFactor + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber), data=Study2rbind, REML=F)      I get significant effects on the linear, quadratic, cubic and quintic terms when I conduct pairwise comparisons between the 3 levels of my 'Condition' variable. However, when I use this as my model (adding a subject:condition random effect);               m.fullCondition <- lmer(PupilBaseCorrect ~ (ot1+ot2+ot3+ot4+ot5) * ConditionFactor + (ot1+ot2+ot3+ot4+ot5 | ParticipantNumber) + (ot1+ot2+ot3+ot4+ot5 | ConditionFactor:ParticipantNumber), data=Study2rbind, REML=F)      I only get a significant effect on the quintic term in just one of my pairwise comparisons (intelligible vs unintelligible). So, I have two questions: (1) Should I be including a subject:condition random effect in my model? My experiment is a within-participant design so I was just wondering whether, by virtue of my design, I should definitely include this is a random effect or not? (2) If you agree that a subject:condition random effect should be included in the model, is an effect on the quintic term (between 2 levels) even meaningful? From my understanding this reflects differences in the shape and timing of the secondary peaks in the data...


title: Can one compare the results of exploratory and confirmatory factor analysis?

text: I have been asked to compare an exploratory factor analysis (EFA) to one based on a hypothesized model. The answer in Why would the results of PCA differ from a confirmatory factor analysis? seems to suggest that this is not conventionally done. I have very little knowledge of this field, but it seems reasonable to want to compare how a model based on prior assumptions/knowledge compares to the one returned by an EFA. Can anyone expound on this and point me to some specific references? Also can this be done in R (or SAS)?


title: Predict after running the mlogit function in R

text: Here's what I want to do, but there seem to be no `predict` method for the mlogit. Any ideas?               library(mlogit)     data("Fishing", package = "mlogit")     Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode")          Fish_fit<-Fish[-1,]     Fish_test<-Fish[1,]     m <- mlogit(mode ~price+ catch | income, data = Fish_fit)     predict(m,newdata=Fish_test)


title: Choice of weight function in Moran's I

text: I'm doing an autocorrelation analysis for a spatially distributed collection of observations. To perform my analysis, I am using Moran's I statistic. My questions are: **(1)** What are the implications and benefits of using different weighting functions, i.e. $d^{-1}$, $d^{-2}$, $\exp(-d)$, and **(2)** Is there any (perhaps informal) answer to which of the possible weighting functions is used most frequently in the geo-statistics literature (and for what purposes)? As for why I care: I am trying to explore whether there is clustering in my data set at different scales of structure, following some of the methodology of Fauchald 2000. I am plotting Moran's I versus aggregation scale. The interesting thing that the resulting correlation curves show very different qualitative behavior when calculating using $d^{-1}$ and $d^{-2}$ weighting functions ($d^{-1}$ has a discontinuity point, for example). I'm having a hard time understanding why this would be true -- does anyone have experience with this who may be able to point me to some references?


title: A "systematic" part of a random time series component?

text: I've seen at least 3 sources on time series* state that the component of a series that is variously called random, stochastic, or noise (something clearly separate from any deterministic, patterned component) itself consists of 2 parts, a systematic part and an unsystematic part. I can't for the life of me figure out how the random part can in turn have a systematic part. I can see how the _series_ can have unsystematic and systematic parts, but not how the random component of the series could have both of these. *For example, Burns and Grove, The Practice of Nursing Research (2005, p. 475-6).


title: How do I generate a QQ-Plot for data fitted using fitdistr?

text: I am doing the following to fit my data using an exponential function:               # Define the data     data <- c(67, 81, 93, 65, 18, 44, 31, 103, 64, 19, 27, 57, 63, 25, 22, 150,                31, 58, 93, 6, 86, 43, 17, 9, 78, 23, 75, 28, 37, 23, 108, 14, 137,                69, 58, 81, 62, 25, 54, 57, 65, 72, 17, 22, 170, 95, 38, 33, 34, 68,                38, 117, 28, 17, 19, 25, 24, 15, 103, 31, 33, 77, 38, 8, 48, 32, 48,                26, 63, 16, 70, 87, 31, 36, 31, 38, 91, 117, 16, 40, 7, 26, 15, 89,                67, 7, 39, 33, 58)          # Fit the data to a model     params = fitdistr(data, "exponential")      I got the following for the params:               > params           rate           0.019694623       (0.002087626)      I want to draw a QQ-plot to see how good the fit was. I am guessing I need to generate exponentially distributed data using the parameter generated and then use some function to draw the QQ plot but am not sure how to go about doing this. Can someone tell me how to do this in R?


title: Within Subject ANOVA?

text: I want to study effect of three factors (Size, Stress and distortion) on readability of graphs. Each factor has two level of difficulties such as easy and hard. Every time I will fix two factors and change another factor to generate a graph. For example, I will fix Size and Stress and will generate a graph with low distortion (easy) and high distortion (hard). Thus In general I will generate 6 graphs. (3 factors * 2 levels). Then I will ask 12 users to perform three tasks with all graphs. I record time and accuracy. So what is the best technique to analyse the data? 2*3 within subjects ANOVA? Repeated ANOVA? and why?


title: Use of R-squared and Slope for determining time-stabilized Data

text: Can I determine cutoffs for the values of R-squared and its slope that are not arbitrary? Can I do anything to improve my method? Exposition: I work with a type-2 superconducting magnet for the measurement of very small magnetic moments. There is a problem which occurs when the the magnet passes through "critical magnetic fields" in which magnetic flux becomes trapped on the surface, this affects measurements. This trapped flux issue can only be alleviated by waiting. I am trying to construct an algorithm which takes into account the R correlation between the retrieved data points (each individual point takes 14 seconds to measure) and the slope of the R linear-regression line to an acceptable set of data points. ![Datapoints between 0 and 5](http://i.stack.imgur.com/yGl5q.png) Above, for x in minutes between 0 and 5, is the type of situation I am trying to avoid by using the linear regression slope to determine when the data is leveled off. ![Data points between 10 and 15](http://i.stack.imgur.com/wc3Ni.png) Above, for x in minutes between 10 and 15, is "good" data. (Neither of these is actual data, just for example) My current algorithm:   1. Take N points (minimum N = 10)   2. Throw away points which are 1 Standard Deviation away from the mean   3. Obtain R^2 value   4. Obtain Slope of R   5. Check Rslope and R^2 with arbitrary limits  If this fails, the algorithm takes data points to fill those which were taken away in Step 2 and cycles. If it meets the criteria then the points are averaged out and their mean and Std. Dev. are exported as one point. Which is why I need Rslope ~> 0. I've found this algorithm is actually pretty good for the high-noise, small- sample data that I extract. The limits/cutoffs that are placed on R^2 and Rslope are arbitrary and were obtained by "eye-balling" the data. To restate the question: Can I determine cutoffs for the values of R-squared and its slope that are not arbitrary? Can I do anything to improve my method? Can anyone direct me toward other methods to consider? (Algorithm written in Delphi/Pascal and used for lots of varying samples)