title: How to select best kernel parameters for radial basis function of SVM with fingerprint recognition?

text: Here, I am using RBF function of SVM for fingerprint verification and matching. Please tell me What is the approximate range of sigma and gamma values in RFB for fingerprint recognition.


title: Two-sample one-sided Kolmogorov-Smirnov test vs one-sided Wilcoxon-Mann-Whitney test

text: I've read in the library manual for test `ks.test` that > The possible values "two.sided", "less" and "greater" of alternative specify > the null hypothesis that the true distribution function of x is equal to, > not less than or not greater than the hypothesized distribution function > (one-sample case) or the distribution function of y (two-sample case), > respectively. This is a comparison of cumulative distribution functions, and > the test statistic is the maximum difference in value, with the statistic in > the "greater" alternative being D^+ = max[F_x(u) - F_y(u)]. Thus in the two- > sample case alternative = "greater" includes distributions for which x is > stochastically smaller than y (the CDF of x lies above and hence to the left > of that for y), in contrast to t.test or wilcox.test. Unfortunately I failed to understand this difference between (I guess one sided) `wilcox.test` and `ks.test`. It seems, that both test for displacement of one distribution versus another. Does anyone can shed some light on it, please?


title: Calculating R-squared values from a semivariogram

text: I have some spatially autocorrelated vegetation data, and would like to know the how well tree size measured at one location can predict tree size in plots 100m away. I've made a semivariogram of the data, but am wondering if I can calculate an R squared value from this? From my (basic) understanding, I'm thinking the semivariance at the sill gives the total variance in the data, whereas the semivariance at 100m measures the variance at that distance; so $\small 1-\frac{\text{semivariance at 100m}}{\text{semivariance at sill}}$ might be equivalent to an $R^2$ value for the relationship between plots separated by 100m? For example, if I have a semivariance of 10 at the sill and 2 at 100m, would this suggest that 1-2/10 = 80% of the variation in tree size across all plots can be explained by tree size in plots 100m away?


title: Collaborative filtering and implicit ratings; normalization?

text: I want to use the time a user spends viewing an article as an implicit rating of how much the user likes the article. My question is how do I normalize this information across all users. At the moment, I'm subtracting the time spent by the user-specific mean, and dividing by the standard deviation. Is this the right way to go about it? It doesn't seem so, as the ratings can still take any values. Maybe I should scale the ratings into some interval (like [1-10]) after? Thanks


title: Creating point clusters of a specified size using R

text: Given a set of points `(x,y)`, I would like to create K clusters of these points similar to the k-means algorithm. The difference is that I would like to define the number of points contained in a cluster, in my situation I would like them to all be equal. For example, take 30 points and put them into 3 sets of 10 points so that the total distance between points in all clusters is minimized. Ideally I would like to implement the solution using R. Thanks for your help!


title: Logistic regression in R with large amounts of data

text: I am currently carrying out an investigation to find if certain factors such as playing home or away or position of a footballer affects overall pass completion using logistic regression. I am using R to compute my data. In my current section in which I am trying to analyse uses the data of every player to convey a general conclusion to whether or not the position of a player affects the successfulness of pass completion. so far I have computed:               test.logit <- glm( cbind(Total.Successful.Passes.All,Total.Unsuccessful.Passes.All) ~                        as.factor(Position.Id), data=passes.data, family = "binomial")          summary(test.logit)      and my output was:               Coefficients:                             Estimate Std. Error z value Pr(>|z|)              (Intercept)              0.28482    0.01256   22.67   <2e-16           as.factor(Position.Id)2  0.99768    0.01438   69.38   <2e-16           as.factor(Position.Id)4  1.06679    0.01398   76.29   <2e-16           as.factor(Position.Id)6  0.68090    0.01652   41.23   <2e-16           Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 32638  on 10269  degrees of freedom     Residual deviance: 26499  on 10266  degrees of freedom          AIC: 60422          Number of Fisher Scoring iterations: 4      the intercept is goalkeepers,position.Id 2 is for a defender, 4 = midfielder and 6 = striker Is this a good set of results to come to a conclusion? and with the large deviances?


title: Minimizing variance of an estimator under sampling cost penalty

text: I have an estimator $t$, whose variance depends on the dimension of my sample $x_{1:n}$: $$ \text{Var}(t(x_{1:n})) = f(n). $$ Suppose that the form of $f(n)$ is known. I would like to determine what is the optimal sample size $n$, given that every sample has a cost $c$. I was considering: a) To minimize a loss function of the type: $$ \text{Loss}(n) = f(n) + n \times c. $$ b) To plot $f(n)$ against $n \times c$ to try determine a good sample size $n$ visually. In my case $c$ is the time needed to simulate $x_i$, for example 0.1sec. I don't think that minimizing a function were different units are summed together makes much sense, so I was wondering whether there is any way to translate CPU times into something more sensible. Some additional information about my application: $t(x_{1:n})$ is the estimate of a gradient, which is calculated at each step of a Stochastic Steepest Descent algorithm. So if $n$ is small the algorithm is fast but unstable, while if $n$ is large it is slow but stable. So far I have just experimented with many values of $n$, but maybe somebody knows a better way of determining $n$. Thanks a lot!


title: Explain ridge in the log-likelihood for Logistic Regression classifier

text: What does the ridge parameter change in a Logistic Regression classifier as for example implemented in Weka Logistic classifier "Parameter -R ridge". The paper describing the underlying theory: Ridge Estimators in Logistic Regression. Any help is much appreciated!


title: Draw multiple plots on one graph in R?

text: Using the following code, I tried to draw four plots on a graph in `R`. I am not happy with the figure because there is a lot of space between the plots, therefore the width of the plots are not sufficient enough to analyze the plots.   1. Could someone help me to produce a nice graph having four plots?    2. How can I keep x-axis labels from 1 to 10 instead of the default 5 labels?  Data: a1: 11.013 13.814 13.831 13.714 13.787 13.734 13.778 13.771 13.823 13.659 a2: 5.181 7.747 8.314 8.061 7.920 8.153 8.540 8.845 7.881 8.301 I have used a1 data for b1, c1, and d1; a2 data for b2, c2, and d2 just for here. Figure: ![enter image description here](http://i.stack.imgur.com/0FgvA.png) Code:                       op=par(mfrow=c(4,1), mar=c(5.5,5.1,4.1,2.1))             plot(a1, type="b", ylim=c(0,14.5), xlab="Time (secs)", ylab="", cex.axis=1.4, cex.lab=1.3,cex=1.2,lwd=2.5,col="red1",lty=2,pch=1, main="A")             lines(a2,type="b",pch=3,lty=3,col="darkblue",lwd=2.5,cex=1.2)             par(xpd=T)             legend(1,26.5,c("X","Y"),bty="n",horiz=T,cex=1.5,col=c("red1","darkblue"),text.col=c("red1","darkblue"),pch=c(1,3),lty=c(2,3),x.intersp=0.4,adj=0.2)             plot(b1, type="b", ylim=c(0,14.5), xlab="Time (secs)", ylab="", cex.axis=1.4, cex.lab=1.3,cex=1.2,lwd=2.5,col="red1",lty=2,pch=1, main="B")             lines(b2,type="b",pch=3,lty=3,col="darkblue",lwd=2.5,cex=1.2)             plot(c1, type="b", ylim=c(0,14.5), xlab="Time (secs)", ylab="", cex.axis=1.4, cex.lab=1.3,cex=1.2,lwd=2.5,col="red1",lty=2,pch=1, main="C")             lines(c2,type="b",pch=3,lty=3,col="darkblue",lwd=2.5,cex=1.2)             plot(d1/1000, type="b", ylim=c(0,14.5), xlab="Time (secs)", ylab="", cex.axis=1.4, cex.lab=1.3,cex=1.2,lwd=2.5,col="red1",lty=2,pch=1, main="D")             lines(d2,type="b",pch=3,lty=3,col="darkblue",lwd=2.5,cex=1.2)        > mtext("Price", side=2, at=100,line=3,cex=1.1)


title: Confusion related to multicollinearity, FA and regression of heterogeneous data

text: I am currently working with a data set that contains about 26 IVs of almost all sorts of scale of measurement (binary, nominal, ordinal and interval scale variables). There are strong reasons to suspect that some variables are probably highly correlated, while some may not be related to any other IVs to a great extent. I came across with great suggestions to resolve this problem in this site (which was an useful advice to use optimally scaled variables in the FA procedure and to use derived factor score as the IVs). But due to my inexperience in this field I am in need of some expert advice on the following issues: **How should I check if Multicollinearity really exists?** I am not sure how to check Multicollinearity with such a heterogeneous data. I may calculate the Heterogeneous Correlation Matrix (or Spearman's Rank Correlation) by somehow forcing me to consider the nominal variables as ordinal but even if I do it what should be the value of the correlation coefficient at which Multicollinearity can be ignored? I am also not sure if it is going to give any insight at all, as I am missing something like a VIF measure! **Should I take only those variables for a FA which are highly correlated?** Say, if I can find two sets of variables (one set containing 8 IVs and another containing 4 IVs) quite highly correlated to each other within each set, then should I use only those 12 variables for FA and derive FA scores for those two factors to use them as IVs? Clearly my intention is to use the other 14 variables separately as IVs along with the two derived scores. I am confused if I should actually use not the 12, but all 26 variables in the FA in this scenario. Remember in that case my FA scores are weighted by the other 14 unrelated variables too! **Is there any problem to categorize a proportion type DV for an ordinal logistic regression?** I've actually found people using logistic regression instead. But I want to mention here that unfortunately I don't know the number of cases (or trials) out of which each proportion was calculated. So I cannot use the number of trials as the weights in the logistic regression. In that case a logistic regression may not be accurate enough. So, as I only know the proportions, won't it be good to categorize the proportions by median split or by quartile split? So that I can use it as a DV in a logistic or in an ordinal logistic regression? I am thankful for reading this thread patiently and hoping some expert advice. Regards.


title: Correcting biased polling

text: Let's say I'm polling for a binary election in different states with known biases. Furthermore, let's say I only manage to poll only a small sample of people in each of these states. How would you combine the results to come up with the best estimate for the election? My best guess is the normalized vote is: $$ \sum_i \frac{\text{# polled in state}_i}{\text{# total polled}} * \text{normalized vote} $$ where $$ \text{normalized vote} = \frac{p(d) - \text{average}(p(d))}{\text{stdev}( p(d))} $$ Assuming $d$ is one of the binary choices, and I have a history of $p(d)$ in each state.


title: Sample space of pmf

text: I have several exercises to solve that deal with sample space of PMF. One is: Let $X_1,\dots , X_n$ be independent random variables with pmf $p(x;\pi) = (1-\pi)^x \pi$ What is the sample space of $X_1?$ Try to give a probabalistic interpretation of such a sample space. Hint: for example. a Bernoulli random variable can be used to model a coin with probability of success $p \in \;\rbrack0,1\lbrack$. Any hint? I do not know how to "see" from the distribution what the sample space is without knowing more about this distribution... ## Possible Solution? I have the following idea: If I plug in all values in the sample space and calculate the sum I have to get $1$ (because it is a probability). But obviously here we can have just one value for the sample space which is $\log({\frac{1}{\pi}})/{\log(1-\pi)}$ (I set the formula for pmf to one and solve the equation for x). Is that correct?


title: Recall and AUC of a binary classifier

text: Is it possible for a binary classifier to have a recall of `0.0` for one of the classes and at the same time an area under the ROC curve (AUC) of `1.0` for the same class?


title: "Bayesglm", p-values and degrees of freedom?

text: I am trying to perform some logistic regressions (and I am a neophyte user of R). Initially I used "glm" to compute coefficients, AIC and p-values; this worked great until I ran across a data set suffering from complete separation. In [1], Gelman et alia suggest using an (informative) prior to address this problem; the corresponding algorithm is implemented in R as "bayesglm" (in the ARM package). Here is my problem. Previously, with "glm", I would compute p-values as follows:               mylogit <- bayesglm(a ~ b+c+d+e+f+g+h, data = mydata, family="binomial")     with(mylogit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))      There are 53-48=5 degrees of freedom:               Null deviance: 71.188  on 53  degrees of freedom     Residual deviance: 37.862  on 48  degrees of freedom      However, if I use "bayesglm" instead of "glm", the resulting degrees of freedom are a bit surprising to me:               Null deviance: 22.279  on 53  degrees of freedom     Residual deviance: 39.030  on 54  degrees of freedom      If I plug in the preceding formula for a p-value, I have -1 degrees of freedom! Can someone help me get a more sensible answer (or help me interpret this)? By the way, the documentation on the "bayesglm" command includes the following ominous comment: > We include all the glm() arguments but we haven’t tested that all the > options (e.g., offests, contrasts, deviance for the null model) all work. [1] Gelman, Andrew, et al. "A weakly informative default prior distribution for logistic and other regression models." The Annals of Applied Statistics (2008): 1360-1383.


title: How to derive the conditional variance and CDF of a bivariate GARCH process?

text: **Question:** Assume a bivariate GARCH process as follows: \begin{align} r_{mt} &= \sigma_{mt}\epsilon_{mt} \ \ \ \cdots \ \ \ \text{(1)} \\ r_{it}&=\sigma_{it}\rho_{it}\epsilon_{mt}+\sigma_{it}\sqrt{1-\rho_{it}^2}\xi_{it} \ \ \ \cdots \ \ \ \text{(2)} \\ (\epsilon_{mt}, \xi_{it}) & \sim S \end{align} where: $\sigma_{mt}$ is the conditional standard deviation of $r_{mt}$ $\sigma_{it}$ is the conditional standard deviation of $r_{it}$ $\rho_{it}$ is the conditional correlation between $r_{it}$ and $r_{mt}$ $(\epsilon_{mt}, \xi_{it})$ are the shocks that drive the system. The shocks $(\epsilon_{mt}, \xi_{it})$ are independent and identically distributed and have zero mean, unit variance and zero covariance. However they are not necessarily independent of each other. The distribution $S$ is unspecified. From Eqns. $(1)$ and $(2)$ we have: $$r_{mt} = \frac{\sigma_{mt}}{\sigma_{it}\rho_{it}}r_{it} - \frac{\sigma_{mt}\sqrt{1-\rho_{it}^2}}{\rho_{it}}\xi_{it}$$ Define $CoVaR$ as: $$Pr\left(r_{mt} \le CoVaR \mid r_{it} = C \right)=\alpha$$ or equivalently, $$1-Pr\left(r_{mt} \le CoVaR \mid r_{it} = C \right)=1-\alpha \\ Pr\left(\xi_{it} \le \frac{\rho_{it}}{\sigma_{mt}\sqrt{1-\rho_{it}^2}} \left(\frac{\sigma_{mt}}{\sigma_{it}\rho_{it}}C - CoVaR \right)\mid r_{it} = C \right)=1-\alpha \ \ \ \cdots \ \ \ (3)$$ * * * **Need clarification #1** : The question then goes on to state: In the special case where the conditional mean function of $\xi_{it}$ is linear in $r_{it}$, the first two conditional moments of $\xi_{it}$, given $r_{it} = C$, can be expressed as: $$\displaystyle \mathbb{E}\left(\xi_{it} \mid r_{it} = C\right) = \frac{cov(\xi_{it}, r_{it})}{\sigma_{it}^2} \times C \\ \mathbb{V}(\xi_{it} \mid r_{it} ) = \mathbb{V}(\xi_{it}) - \mathbb{V}_{r_{it}}\left[\mathbb{E}\left(\xi_{it} \mid r_{it} \right) \right] = \mathbb{V}(\xi_{it})\times\left[1-\left( \frac{cov(\xi_{it}, r_{it})}{\sigma_{it}^2}\right)^2 \sigma_{it}^2 \right]$$ Where does the expression $\displaystyle \mathbb{V}(\xi_{it})\times\left[1-\left( \frac{cov(\xi_{it}, r_{it})}{\sigma_{it}^2}\right)^2 \sigma_{it}^2 \right]$ come from? How do you derive this expression? **Need clarification #2** : The question continues to state: Consider $G(\cdot)$, the conditional (location-scale) demeaned and standardized cdf of $\xi_{it}$ such that: $$ \mathbb{E}\left[\frac{1}{\rho_{it}}\left(\xi_{it} - \frac{\sqrt{1-\rho_{it}^2}}{\sigma_{it}} \times C \right) \mid r_{it} = C\right] =0 \\ \mathbb{V}\left[\frac{1}{\rho_{it}}\left(\xi_{it} - \frac{\sqrt{1-\rho_{it}^2}}{\sigma_{it}} \times C \right) \mid r_{it} = C \right] =1$$ Then Eqn. $(3)$ can be expressed as: $$\frac{1}{\rho_{it}} \left[ \frac{\rho_{it}}{\sigma_{mt}\sqrt{1-\rho_{it}^2}} \left(\frac{\sigma_{mt}}{\sigma_{it}\rho_{it}}C - CoVaR \right) - \frac{\sqrt{1-\rho_{it}^2}}{\sigma_{it}} \times C \right] = G^{-1} \left(1-\alpha\right) \ \ \ \cdots (4)$$ Where does Eqn. $(4)$ come from? How do you derive it? Thanks in advance!


title: Simple probability question

text: A person from group A has 20% chance of having some characteristic A person from group B has 30% chance of having the same characteristic How can I calculate the probability of a person belonging to both groups having the given characteristic? **edit** Since a number of people pointed out that it's impossible to know, I'll change my question. Let's assume that there are 10 different groups. Do I need to know the probabilities for each possible combination, or can I infere at least some probabilities? **edit** I've added a solution witch I think is plausible for some cases.


title: What is the appropriate test for comparing means?

text: I have a case and control group, both of which had a blood sugar level tested multiple times over 48 hours. I took the mean of the levels for each person and then calculated the mean of the means (with an SD etc) to compare the two groups. I am not sure if this is a valid method for comparison and whether another method is more appropriate? Thanks


title: Common parameters for conditional likelihood

text: I am trying to understand the concept of conditional likelihood in the context of logistic regression. One paper I am reading defines $L(\theta; y|x) = f(y|x; \theta)$, then goes on to point out that there may be different distributions $f(y|x; \theta)$ for every different $x$, and also that these different distributions share the same parameters $\theta$. The first part of this makes sense that there may be different distributions for each $x$, however, I am failing to understand why $\theta$ must be the same for each of these distributions. Can anyone shed some light on this for me?


title: Design fitness function for polynomial approximation

text: I'm trying to apply a `polynomial approximation` for a given function (via Genetic Algorithms), and so far the results are not so good:                # or any other GA package       require(gaoptim)            # for polyval       require(pracma)            # polynomial of degree 10       ndeg = 10         ndim = ndeg + 1            # search limits       search.low = rep(-1, ndim)       search.up = rep(1, ndim)            # no. of data points       m = 101        xi = seq(-1, 1, length = m)       yi = 1 / (1 + (5*xi)^2)            # fitness function       pfn = function(p) max(abs(polyval(c(p), xi) - yi))            # gaoptim perform maximization, so transform the fitness function       pfninv = function(p){ 1/(pfn(p) + 1) }            ## set up the ga       ga = GAReal(pfninv, search.low, search.up, popSize = 500)       ga$evolve(100)       y2 = polyval(ga$bestIndividual(), xi)            plot(xi, yi, ylim = range(c(yi, y2)), type = 'l', main = 'Runge function')       lines(xi, y2, col = 'red')      ![enter image description here](http://i.stack.imgur.com/cJ5kX.png) Is there any strategy i can apply here, or this is a No-no approach? Maybe a better fitness function, or expand the search limits? Higher values of _popSize_ doesn't seem to help too much. Thanks for any insight!


title: How to report significance of factor levels shown in summary of glm?

text: The summary of my GLM shows day 12 of factor Date to be significant, but anova(model, test="Chisq") shows Date to be not significant overall. I know how to report the statistics from the Chisq table, but as I have z values in the summary table I am unsure how to report, or if I should report, that day 12 is significant. Similarly when finding that Date is significant, how to report what specific dates seem to be important. Thanks in advance Lara I have Fertility and Fecundity of female flies measured over 12 days. I want to check if there is a decline (or otherwise) in fertility/fecundity over this time. Additionally, at day 13 females are mated (with a male from one of two different groups), and fertility and fecundity are measured until day 20. I want to use Date as a factor to identify significant peaks in fertility/fecundity i.e. after mating, and potential difference in peaks between groups.               Call:     glm(formula = dda$Fertility.Absolute ~ sqrt(dda$Fecundity) +          dda$Group + dda$Date + sqrt(dda$Fecundity):dda$Group + sqrt(dda$Fecundity):dda$Date,          family = poisson)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.1397  -0.6786  -0.4797   0.3596   3.7588            Coefficients:                                        Estimate Std. Error z value Pr(>|z|)             (Intercept)                    -1.91501    0.51539  -3.716 0.000203 ***         sqrt(dda$Fecundity)             0.72372    0.12441   5.817 5.99e-09 ***         dda$Group2                      0.19540    0.19230   1.016 0.309585             dda$Date4                       0.18117    0.62648   0.289 0.772439             dda$Date6                      -0.28952    0.68983  -0.420 0.674706             dda$Date8                       0.07111    0.60531   0.117 0.906480             dda$Date10                      0.19557    0.62232   0.314 0.753325             dda$Date12                      0.79619    0.60710   1.311 0.189696             dda$Date14                      1.93702    0.53938   3.591 0.000329 ***         dda$Date16                      0.75623    0.58296   1.297 0.194554             dda$Date18                     -0.05392    0.67805  -0.080 0.936618             dda$Date20                     -0.26291    0.68841  -0.382 0.702530             sqrt(dda$Fecundity):dda$Group2 -0.07309    0.04822  -1.516 0.129583             sqrt(dda$Fecundity):dda$Date4   0.27388    0.17555   1.560 0.118734             sqrt(dda$Fecundity):dda$Date6   0.37684    0.22832   1.651 0.098836 .           sqrt(dda$Fecundity):dda$Date8   0.13017    0.13861   0.939 0.347674             sqrt(dda$Fecundity):dda$Date10  0.04552    0.15345   0.297 0.766722             sqrt(dda$Fecundity):dda$Date12 -0.16593    0.14861  -1.117 0.264186             sqrt(dda$Fecundity):dda$Date14 -0.24864    0.12754  -1.949 0.051240 .           sqrt(dda$Fecundity):dda$Date16  0.05496    0.14578   0.377 0.706170             sqrt(dda$Fecundity):dda$Date18  0.15439    0.19341   0.798 0.424715             sqrt(dda$Fecundity):dda$Date20 -0.02006    0.16314  -0.123 0.902161         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for poisson family taken to be 1)              Null deviance: 1327.36  on 359  degrees of freedom     Residual deviance:  298.26  on 338  degrees of freedom     AIC: 873.65          Number of Fisher Scoring iterations: 5          ###          Analysis of Deviance Table         Model: poisson, link: log         Response: dda$Fertility.Absolute         Terms added sequentially (first to last)                                                Df Deviance Resid. Df Resid. Dev  Pr(>Chi)         NULL                                            359    1327.36                       sqrt(dda$Fecundity)            1   893.88       358     433.48 < 2.2e-16 ***         dda$Group                      1     0.03       357     433.45    0.8699             dda$Date                       9    82.16       348     351.29  6.01e-14 ***         sqrt(dda$Fecundity):dda$Group  1     0.07       347     351.22    0.7859             sqrt(dda$Fecundity):dda$Date   9    52.96       338     298.26  2.97e-08 ***


title: Parallel logistic regression

text: I need to perform stepwise binary logistic regression (The horror! The horror!) on 1.5 million observations. This takes far too long in SAS, so I'm wondering if I can use R to process it in a multicore environment. Apparently package gmulti (http://www.jstatsoft.org/v34/i12/paper) will do the trick, but it's not clear to me if it will do that outside of its genetic algorithm. That still might work for me, but I don't have a large number of variables (about 30) so it's not necessary. As long as the results of the brute force and ga approach could be assured to be similar, then I might try it. However, I see others have had problems getting the parallel feature to run: https://stat.ethz.ch/pipermail/r-help/2013-April/351820.html. Any other suggestions on how to parallelize logistic regression in R? A web search turned up a couple of papers, but not much that seemed specific to R. And please spare me a lecture about stepwise regression-I'm very well aware of the pitfalls. I'm replicating someone else's analysis. I'm using a Windows 64 bit system.


title: What is an acceptable value of the Calinski & Harabasz (CH) criterion?

text: I have done a data analysis trying to cluster longitudinal data using R and the kml package. My data contains of around 400 individual trajectories (as it is called in the paper). You can see my results in the following picture: ![enter image description here](http://i.stack.imgur.com/8cB8m.png) After reading chapter 2.2 "Choosing an optimal number of clusters" in the corresponding paper I didn't get any answers. I would prefer having 3 clusters but will the result still being Ok with a CH of 80. Actually I even don't know what the CH value represents. So my question, what is an acceptable value of the Calinski & Harabasz (CH) criterion?


title: Obtaining predictions from linear model

text: This is probably a pretty simple question, but I have been having some trouble interpreting the documentation for the `predict` function. I am generating a simple linear model from a data frame containing (X, Y) pairs, which I would then like to use to predict Y given new X. My code looks something like this:               my_lm = lm(Y ~ X, data=my_data)     new_Y = predict(my_lm, new_data, interval="confidence", type="response")     print(new_Y[,'fit'])      This does seem to give me values already, but I am a bit confused about the `type` parameter. The documentation only mentions that `type` specifies "response or model term", but it doesn't say what that actually means. Searching the web for examples I have seen both methods mentioned on different sites, with no clear explanation why they were used. Can someone tell me what the right way to solve my problem is and where I can read up on the actual meaning of the type?


title: Multivariate random walks in BUGS

text: I want to jointly estimate a very simple MV-Normal two-dimensional AR[1] process, $[x_t,y_t]=[x_{t-1},y_{t-1}]+\text{[Bivariate Gaussian error]}$, in BUGS. But the syntax has been impossible to figure out. Here's the problem part of the code:                ## transition model (aka random walk prior)      for(i in 2:NPERIODS1){                  mu.vector[i,1:2]<-vector[i-1,1:2]         vector[i,1:2]~dmnorm(mu.vector[i,1:2], omega[1:2,1:2])         }      The compiler throws up a "Expected a multivariate node" error. Looking through some examples, there doesn't seem to be any easy way to introduce a structured mean or covariance variables for the multivariate normal function. How should I proceed? **Edit:** Changed omega[,] to omega[1:2,1:2] for clarity.


title: Inspect generators and defining relations of a fractional factorial design

text: Is there a way to give (in R or Minitab or Statgraphics) a fractional factorial design like that and inspect the generators and the complete defining relation ($2^4 - 1$ relations)?               A   B   C   D   E   F   G   H     -1  1   1   1   -1  -1  1   -1     -1  -1  -1  1   1   1   1   -1     -1  1   1   -1  1   1   -1  -1     -1  -1  -1  -1  -1  -1  -1  -1     -1  1   -1  -1  1   -1  1   1     -1  -1  1   -1  -1  1   1   1     -1  1   -1  1   -1  1   -1  1     -1  -1  1   1   1   -1  -1  1     1   1   -1  -1  -1  1   1   -1     1   -1  1   -1  1   -1  1   -1     1   1   -1  1   1   -1  -1  -1     1   -1  1   1   -1  1   -1  -1     1   1   1   1   1   1   1   1     1   -1  -1  1   -1  -1  1   1     1   1   1   -1  -1  -1  -1  1     1   -1  -1  -1  1   1   -1  1      Thank you **EDIT** (This is not an answer. It's just a workaround)               df<- read.table(textConnection("     A B C D E F G H     -1  1   1   1   -1  -1  1   -1     -1  -1  -1  1   1   1   1   -1     -1  1   1   -1  1   1   -1  -1     -1  -1  -1  -1  -1  -1  -1  -1     -1  1   -1  -1  1   -1  1   1     -1  -1  1   -1  -1  1   1   1     -1  1   -1  1   -1  1   -1  1     -1  -1  1   1   1   -1  -1  1     1   1   -1  -1  -1  1   1   -1     1   -1  1   -1  1   -1  1   -1     1   1   -1  1   1   -1  -1  -1     1   -1  1   1   -1  1   -1  -1     1   1   1   1   1   1   1   1     1   -1  -1  1   -1  -1  1   1     1   1   1   -1  -1  -1  -1  1     1   -1  -1  -1  1   1   -1  1     ")->con,header=T)     close(con)      ABCDEFGH gives a column of ones, so I=ABCDEFGH. We need the 14 remaining defining relations.               # start with four way interactions     four.way <- combn(c("A","B","C","D","E","F","G","H"),4)          res <- apply(four.way,2,function(x) {     apply(df[,x],1,prod)})     colnames(res) <- apply(four.way,2,paste,collapse="")      The res matrix has 70 columns of the products of the 4-combinations. Get the column names with a colSum of either 16 or -16 (they are defining relations)               def <- colnames(res[,colSums(res) == 16 | colSums(res) == -16])     [1] "ABCH" "ABDE" "ABFG" "ACDF" "ACEG" "ADGH" "AEFH" "BCDG" "BCEF" "BDFH" "BEGH" "CDEH" "CFGH" "DEFG"      So I got another 14 defining relations (including the generalized interactions). No need to look for other "words" since I already got 15. The resolution is IV (min word length). It's easy to observe that the BGHA columns form the typical $2^4$ design (with opposite signs for B and G). Using the above defining relations you can get the generators: C=ABH, D=AGH, E=BGH and F=ABG. To get the alias structure I did the following (I applied it only for the main effects and 2-way interactions, but you can get whatever you ask for)               # two way interactions     two.way <- apply(combn(c("A","B","C","D","E","F","G","H"),2),2,paste,collapse="")          test <- c("A","B","C","D","E","F","G","H",two.way)     gen <- c(def,"ABCDEFGH")          mat <- character(length(test)*(length(gen)+1))     dim(mat) <- c(length(test),length(gen)+1)     colnames(mat) <- c("Effect",paste("I=",gen,sep=""))     for (j in 1:length(test)){     mat[j,1] <- test[j]     for (i in 1:length(gen)) {     res <- paste(sort(c(unlist(strsplit(gen[i],"")),unlist(strsplit(test[j],"")))),collapse="")             le <- rle(unlist(strsplit(res,"")))$lengths             va <- rle(unlist(strsplit(res,"")))$values             mat[j,i+1] <- paste(sort(va[le %% 2 == 1]),collapse="")     }}     noquote(mat)      And here is the result                     Effect I=ABCH I=ABDE I=ABFG I=ACDF I=ACEG I=ADGH I=AEFH I=BCDG I=BCEF I=BDFH I=BEGH I=CDEH I=CFGH I=DEFG I=ABCDEFGH      [1,] A      BCH    BDE    BFG    CDF    CEG    DGH    EFH    ABCDG  ABCEF  ABDFH  ABEGH  ACDEH  ACFGH  ADEFG  BCDEFGH         [2,] B      ACH    ADE    AFG    ABCDF  ABCEG  ABDGH  ABEFH  CDG    CEF    DFH    EGH    BCDEH  BCFGH  BDEFG  ACDEFGH         [3,] C      ABH    ABCDE  ABCFG  ADF    AEG    ACDGH  ACEFH  BDG    BEF    BCDFH  BCEGH  DEH    FGH    CDEFG  ABDEFGH         [4,] D      ABCDH  ABE    ABDFG  ACF    ACDEG  AGH    ADEFH  BCG    BCDEF  BFH    BDEGH  CEH    CDFGH  EFG    ABCEFGH         [5,] E      ABCEH  ABD    ABEFG  ACDEF  ACG    ADEGH  AFH    BCDEG  BCF    BDEFH  BGH    CDH    CEFGH  DFG    ABCDFGH         [6,] F      ABCFH  ABDEF  ABG    ACD    ACEFG  ADFGH  AEH    BCDFG  BCE    BDH    BEFGH  CDEFH  CGH    DEG    ABCDEGH         [7,] G      ABCGH  ABDEG  ABF    ACDFG  ACE    ADH    AEFGH  BCD    BCEFG  BDFGH  BEH    CDEGH  CFH    DEF    ABCDEFH         [8,] H      ABC    ABDEH  ABFGH  ACDFH  ACEGH  ADG    AEF    BCDGH  BCEFH  BDF    BEG    CDE    CFG    DEFGH  ABCDEFG         [9,] AB     CH     DE     FG     BCDF   BCEG   BDGH   BEFH   ACDG   ACEF   ADFH   AEGH   ABCDEH ABCFGH ABDEFG CDEFGH         [10,] AC     BH     BCDE   BCFG   DF     EG     CDGH   CEFH   ABDG   ABEF   ABCDFH ABCEGH ADEH   AFGH   ACDEFG BDEFGH         [11,] AD     BCDH   BE     BDFG   CF     CDEG   GH     DEFH   ABCG   ABCDEF ABFH   ABDEGH ACEH   ACDFGH AEFG   BCEFGH         [12,] AE     BCEH   BD     BEFG   CDEF   CG     DEGH   FH     ABCDEG ABCF   ABDEFH ABGH   ACDH   ACEFGH ADFG   BCDFGH         [13,] AF     BCFH   BDEF   BG     CD     CEFG   DFGH   EH     ABCDFG ABCE   ABDH   ABEFGH ACDEFH ACGH   ADEG   BCDEGH         [14,] AG     BCGH   BDEG   BF     CDFG   CE     DH     EFGH   ABCD   ABCEFG ABDFGH ABEH   ACDEGH ACFH   ADEF   BCDEFH         [15,] AH     BC     BDEH   BFGH   CDFH   CEGH   DG     EF     ABCDGH ABCEFH ABDF   ABEG   ACDE   ACFG   ADEFGH BCDEFG         [16,] BC     AH     ACDE   ACFG   ABDF   ABEG   ABCDGH ABCEFH DG     EF     CDFH   CEGH   BDEH   BFGH   BCDEFG ADEFGH         [17,] BD     ACDH   AE     ADFG   ABCF   ABCDEG ABGH   ABDEFH CG     CDEF   FH     DEGH   BCEH   BCDFGH BEFG   ACEFGH         [18,] BE     ACEH   AD     AEFG   ABCDEF ABCG   ABDEGH ABFH   CDEG   CF     DEFH   GH     BCDH   BCEFGH BDFG   ACDFGH         [19,] BF     ACFH   ADEF   AG     ABCD   ABCEFG ABDFGH ABEH   CDFG   CE     DH     EFGH   BCDEFH BCGH   BDEG   ACDEGH         [20,] BG     ACGH   ADEG   AF     ABCDFG ABCE   ABDH   ABEFGH CD     CEFG   DFGH   EH     BCDEGH BCFH   BDEF   ACDEFH         [21,] BH     AC     ADEH   AFGH   ABCDFH ABCEGH ABDG   ABEF   CDGH   CEFH   DF     EG     BCDE   BCFG   BDEFGH ACDEFG         [22,] CD     ABDH   ABCE   ABCDFG AF     ADEG   ACGH   ACDEFH BG     BDEF   BCFH   BCDEGH EH     DFGH   CEFG   ABEFGH         [23,] CE     ABEH   ABCD   ABCEFG ADEF   AG     ACDEGH ACFH   BDEG   BF     BCDEFH BCGH   DH     EFGH   CDFG   ABDFGH         [24,] CF     ABFH   ABCDEF ABCG   AD     AEFG   ACDFGH ACEH   BDFG   BE     BCDH   BCEFGH DEFH   GH     CDEG   ABDEGH         [25,] CG     ABGH   ABCDEG ABCF   ADFG   AE     ACDH   ACEFGH BD     BEFG   BCDFGH BCEH   DEGH   FH     CDEF   ABDEFH         [26,] CH     AB     ABCDEH ABCFGH ADFH   AEGH   ACDG   ACEF   BDGH   BEFH   BCDF   BCEG   DE     FG     CDEFGH ABDEFG         [27,] DE     ABCDEH AB     ABDEFG ACEF   ACDG   AEGH   ADFH   BCEG   BCDF   BEFH   BDGH   CH     CDEFGH FG     ABCFGH         [28,] DF     ABCDFH ABEF   ABDG   AC     ACDEFG AFGH   ADEH   BCFG   BCDE   BH     BDEFGH CEFH   CDGH   EG     ABCEGH         [29,] DG     ABCDGH ABEG   ABDF   ACFG   ACDE   AH     ADEFGH BC     BCDEFG BFGH   BDEH   CEGH   CDFH   EF     ABCEFH         [30,] DH     ABCD   ABEH   ABDFGH ACFH   ACDEGH AG     ADEF   BCGH   BCDEFH BF     BDEG   CE     CDFG   EFGH   ABCEFG         [31,] EF     ABCEFH ABDF   ABEG   ACDE   ACFG   ADEFGH AH     BCDEFG BC     BDEH   BFGH   CDFH   CEGH   DG     ABCDGH         [32,] EG     ABCEGH ABDG   ABEF   ACDEFG AC     ADEH   AFGH   BCDE   BCFG   BDEFGH BH     CDGH   CEFH   DF     ABCDFH         [33,] EH     ABCE   ABDH   ABEFGH ACDEFH ACGH   ADEG   AF     BCDEGH BCFH   BDEF   BG     CD     CEFG   DFGH   ABCDFG         [34,] FG     ABCFGH ABDEFG AB     ACDG   ACEF   ADFH   AEGH   BCDF   BCEG   BDGH   BEFH   CDEFGH CH     DE     ABCDEH         [35,] FH     ABCF   ABDEFH ABGH   ACDH   ACEFGH ADFG   AE     BCDFGH BCEH   BD     BEFG   CDEF   CG     DEGH   ABCDEG         [36,] GH     ABCG   ABDEGH ABFH   ACDFGH ACEH   AD     AEFG   BCDH   BCEFGH BDFG   BE     CDEG   CF     DEFH   ABCDEF


title: Exchangeable Processes over the Simplex

text: You are likely all familiar with Polya Urn process. I initially start with an urn containing $b$ black balls and $w$ white balls. At each step, I sample a black ball with probability $\frac{b}{b+w}$ (alternatively, a white ball with probability $\frac{w}{b+w}$), then put back the same ball and another of the same color. Equivalently, if $x_i$ is the $i$-ith ball drawn, I update $w \leftarrow w + \delta(x_i = w)$ and $b \leftarrow b +\delta(x_i = b)$ The Polya Urn process has some wonderful properties, first that the distribution of $\frac{b}{b+w}$ converges to a sample from a Beta distribution, and second that the probability of any sequence of draws is exchangeable. I would like to use a similar process wherein each draw isn't a discrete random variable, but rather a continuous variable on $(0,1)$. My initial intuition was to replace each draw from the urn with a draw $p_i \sim Beta(b, w)$, then to update $b \leftarrow b + p_i$ and $w \leftarrow w + (1-p_i)$. While $\frac{b}{b+w}$ converges much the same way it does in the Polya urn process, the pdf of any sequence of draws is **not** exchangeable. My question is, do there exist any urn-like processes over the probability simplex that are exchangeable?


title: Random Investment

text: $\DeclareMathOperator{\E}{\mathrm{E}}$ An investor has $\$100,000$. When the current interest rate is $i$% (compounded continously so that the growth per year is $e^{i/100}$) she invests her money in an $i$ year CD, takes the profits, and then immediately reinvests the $\$100,000$. Suppose that the $k$th investment has an investment rate of $X_k$% where the $X_k \overset{\mathrm{iid}}{\sim} \operatorname{Uniform}\{1,2,3,4,5\}$. In the long run, how much money does she make per year? **My Thoughts** Modeling this as a renewal process, if I set the reward in the $i$th period $\E r_i$ and the length of the $i$th period, $\E t_i$ equal to $$\E r_i = \sum_{i = 1}^5 100,000(e^{i/100} - 1) \cdot i \approx 56,141.5 \qquad \E t_i = {1 + 2 + 3 + 4 + 5 \over 5} = 3$$ Then by one of our theorem, I should get the long-run rate at which rewards are earned as $${\E r_i \over \E t_i} = 18,713.83 $$ This should be my answer. Is this a valid approach, or have I overlooked something?


title: Joint distribution of a random variable and the sample maximum

text: This is one necessary part of a slightly larger problem, but this part has me stumped. We have that $X_1, X_2, ..., X_n\stackrel{iid}{\sim} U(0,\theta)$. What is the joint density of the first observation and the maximum order statistic: $f_{X_1,X_{(n)}}(s,t)$? Is this doable, any ideas? Many thanks for any help.


title: Casting multidimensional data in R into a data frame

text: I'm trying to do some percentage based comparisons across different groups in a criminal sentencing data set (http://dl.dropbox.com/u/1156404/wightCrimRecords.csv) I have a data in an array of the form:                   $Female         x                                  Burglary                 Criminal Damage and Arson                                0.004950495                               0.017326733                           Driving Offences                                    Murder                                0.371287129                               0.000000000         $Male         x                                  Burglary                 Criminal Damage and Arson                                0.013001083                               0.058504875                           Driving Offences                                    Murder                                0.303358613                               0.000000000              $`Not Stated`         x                                  Burglary                 Criminal Damage and Arson                                  0.0000000                                 0.0000000                           Driving Offences                                    Murder                                  0.1111111                                 0.0000000       This was derived from the original data as follows:               iw=read.csv("~/data/recordlevel.csv")     iwp=tapply(iw$Offence_type,iw$AGE,function(x){prop.table(table(x))})      What I would like to do is generate a single data frame that contains a gender column, a frequency column, and rows corresponding to Burglary, Murder etc. I can extract a single datatable from the multidimensional array, eg using:               iwpF =data.frame(iwp['Female'])      which generates a separate row for each offence and columns referring to offence type and frequency, but can't see how to generate a single datatable. PS I was also wondering whether it's possible to pull out even more structured data, that for example counts the percentages of offence type sex and age group, so for example I could lookup up what percentage of convictions for males in the 35+ age range are related to murder.


title: Fixed, random effect and nested factor in lme

text: I have a dataset with the following variables:   * Treatment : (fixed,3 levels)    * Location : (fixed, 4 levels)    * Sample (random, 5 levels): 5 samples are taken in each location (randomly)    * Subsample (random - 3 levels): 3 subsamples are taken in each sample. One treatment is assigned (randomly) to one of the subsample. The study was conducted every Wednesday during 15 weeks in 2013, from week 32 to week 47 (no observations for week 33). For the first location, samples are taken week 32, 36, 39, 42 and week 49.   For the second location, samples are taken week 46   For the third location, samples are taken week 35, 38, 41, 44   For the 4th location, samples are taken at week 34, 37, 40, 43, 47. For each location _i_ , week _j(i)_ , sample _k(ij)_ , and treatment _l_ , the weight of the subsample has been recorded. The aim of the study is to see if there is an effect of the factor treatment on the weight. I am new in mixed effect model so I am reading the book "Mixed effects Models in S and S-plus" from Pinheiro & Bates, but I still have some doubts on how to treat the different variables (random, fixed, nested). I would like to use the lme function from the nlme package. I think the formula will look like this : `weight ~ treatment , random = ~ 1| location/week/sample` So I consider that   * sample is nested in week    * week is nested in location Is my model correct? Thank you for you help


title: How do I sort an ordinal list of user-generated ratings data?

text: The Ubuntu Software Center uses a 1 to 5 star rating system for its App Reviews. However, it's current ratings sorting algorithm looks very fishy. I believe this is a different question from this one, which seems to assume the mean of the ratings is a meaningful number. If I assume that user supplied ratings are ordinal data, then doing things like adding the ratings together or taking their mean are not correct. The primary method of sorting must be the sample median. Unfortunately, this leaves me with a lot of duplicates since 30 applications are being pigeon-holed into 5 stars of ratings, so it must be possible to further subsort the apps with identical median stars. I believe I want:   * Rating at the median should be better than rating below it.      * {2,3,3,3,4} > {2,2,3,3,4}   * Similarly, rating at the median should be worse than rating above it.      * {2,2,3,3,4} < {2,2,3,4,4}   * Ratings above and below the median should be equivalent.      * {1,1,3,4,4} = {2,2,3,4,4} = {2,2,3,5,5}   * Among two apps with the same median, higher confidence that the median is at least that large should rank higher.      * {2,3,3,3,4} > {2,3,4} Are these reasonable desires? What algorithm can get me there? My intuition tells me I want something like sample median + {lower bound probability estimate someone ranks that app higher than median} - {upper bound probability estimate someone ranks that app lower than median}. So:   * a large data set composed of 20 % 1's, 40% 3's, and 40% 4's would approach (3)+(2/5)-(1/5) = 3.2   * a set composed of equal parts 2,3, and 5 would approach (3)+(1/3)-(1/3) = 3.0   * and a set composed of 40% 1's and 60% 3's would approach (3)+(0)-(2/5) = 2.6. Is this reasonable?


title: Notation for nested factors that are random variables in R

text: I need to compute the estimate the variance component for my data. I have the following model, y = gene, cell line, gene*cell line, DNA extract[cell line] How do I write this for the lme function? I am confused on how to indicate that "DNA extract" is nested within "cell line" and also have it marked as a random variable.


title: Calculating the Standard Error of the Mean (SEM) of a sample from the SEMs of two subsets of that sample

text: Can you calculate the Standard Error of the Mean (SEM) for a set of data given the SEM from two subsets of the data? How would you do this? (The two subsets are mutually exclusive and exhaustive.) Let's say I have the mean and SEM for both men and women on some test, but I do not have access to the original data. I know there are _n1_ men and _n2_ women, and each person took the test _m_ times. It is trivial to calculate the mean of all test takers, but can I calculate the SEM of the entire sample the same way? Wikipedia tells me that "If the standard error of several individual quantities is known then the standard error of some function of the quantities can be easily calculated in many cases" but does not say anything about which quantities or functions.


title: What are the important conditions in ANOVA fixed effects?

text: I am working with an ANOVA model. I want to run a fixed effects ANOVA in which I have a ratio dependent variable and three independent variables with two and three levels. Obviously, before analyzing the results, I want to check the assumptions of factorial ANOVA. Reviewing some handbooks, I found some divergences in their explanations about assumptions of ANOVA. Moreover, I have doubts about the underlying assumptions of factorial ANOVA. The major ones are:   1. All of the handbooks that I check point out that the dependent variable in ANOVA models should be at least an interval variable. I work with a count variable, in fact converted to a ratio variable (i.e. a percentage). So, is ANOVA appropriate in this case?   2. All of the handbooks stress the importance of checking the assumptions of the ANOVA model for inference: mainly a) independence, b) normality, and c) homogeneity of variances. However, they examine this aspects in different ways. Some of them check the data, i.e. the independence of cases, normality of each group and homogeneity of variances between groups. But others examine only the residuals (error) derived from the analysis (i.e. the independence, normality and homoscedasticity of the residuals). So, I am confused about the appropriateness of my approach, but also about the assumptions that I should review. What does the ANOVA model require? Parametric assumptions for variables, only for residuals or both? References are welcome.


title: comparing means between different studies

text: In a questionnaire, I have asked some Likert scale questions. These questions were asked to the same population 6 years ago in another study. I only have the mean scores and standard deviation of the old study. Is there a way to compare the answers of my study to that of the old one statistically? I have limited statistical knowledge and am using SPSS. Below are the means (standard deviation) of each. Scale is from 0-4 The first study has n=170 vs the second study having n=184 Q1 2.11 (1.352) vs 1.83 (1.112) If you could possibly help me do one and I could work out the rest. Thanks!


title: Extrapolating the amount of data missing from the amount of data partially missing

text: Suppose I have a fleet of two-car trains riding around, and that each car is equipped with a data recording device. Unfortunately, some of the recording devices aren't working. I don't know either the exact size of the fleet or the percentage of failed recording devices. I'd like to make a reasonable guess about how many cars I'm missing data from. In particular:   * Total fleet size $S$ is unknown.   * Failure rate (failed units / total cars) $F$ is unknown.   * I have data from both cars from $S_2$ trains and from exactly one car from $S_1$ trains.   * Therefore, I know that there are at least $S_1$ failed units. What I don't know is $S_0$, the number of trains with failed units on both cars.   * Let's assume that the distribution of cars with failed units is random, and whether one car has a failed unit is independent of whether it's mate does. Does the following make sense for a first-order approximation?   * Guess that the failure rate $F'$ is equal to the proportion of missing cars that I know about to total cars that I know about: $F' = S_1 / (2*(S_2 + S_1))$   * Assume that the likelihood of a train having two failed units is $F'^2$.   * Therefore $S_0 = F'^2 * (S_0 + S_1 + S_2)$   * Therefore $S_0 = (S_2 + S_1) * F'^2 / (1 - F'^2)$


title: Finding the change point in data from a piecewise linear function

text: Greetings, I'm performing research that will help determine the size of observed space and the time elapsed since the big bang. Hopefully you can help! I have data conforming to a piecewise linear function on which I want to perform two linear regressions. There is a point at which the slope and intercept change, and I need to (write a program to) find this point. Thoughts?


title: Keywords to find academic lectures on multivariate time series analysis on youtube

text: If I search for "HOTT homotopy type theory" on youtube, I find numerous (advanced/state-of-the art) academic lectures on the topic. For instance the following lectures are found on youtube:   http://www.youtube.com/watch?v=O45LaFsaqMA   http://www.youtube.com/watch?v=wUJ83KnsQnw I would like to find academic lectures on multivariate time series analysis. Unfortunately, my searches usually lead to industry lectures or introductory academic lectures. What keywords can I try to use to find (advanced) academic lectures on multivariate time series analysis?


title: Are there bounds on the Spearman correlation of a sum of two variables?

text: Given $n$-vectors $x, y_1, y_2$ such that the Spearman correlation coefficient of $x$ and $y_i$ is $\rho_i = \rho(x,y_i)$, are there known bounds on the Spearman coefficient of $x$ with $y_1 + y_2$, in terms of the $\rho_i$ (and $n$, presumably)? That is, can one find (non-trivial) functions $l(\rho_1,\rho_2,n), u(\rho_1,\rho_2,n)$ such that $$l(\rho_1,\rho_2,n) \le \rho(x,y_1+y_2) \le u(\rho_1,\rho_2,n)$$ **edit** : per @whuber's example in the comment, it appears that in the general case, only the trivial bounds $l = -1, u = 1$ can be made. Thus, I would like to further impose the constraint:   * $y_1, y_2$ are permutations of the integers $1 \ldots n$.


title: Time series modeling with R on weekly data

text: I am trying to do time series modeling and forecasting using R based on weekly data like below -               biz week     Amount        Count     2006-12-27   973710.7     816570     2007-01-03  4503493.2    3223259     2007-01-10  2593355.9    1659136     2007-01-17  2897670.9    2127792     2007-01-24  3590427.5    2919482     2007-01-31  3761025.7    2981363     2007-02-07  3550213.1    2773988     2007-02-14  3978005.1    3219907     2007-02-21  4020536.0    3027837     2007-02-28  4038007.9    3191570     2007-03-07  3504142.2    2816720     2007-03-14  3427323.1    2703761     ...     2014-02-26  99999999.9   1234567      Regarding my data, as seen above, each week is labeled by first day for the week (my weeks start on Wednesdays and end on Tuesdays). When I construct my `ts` object, I tried:               ts <- ts(df, frequency=52, start=c(2007,1))      The problems I have are:   1. Some years may have 53 weeks, so `frequency=52` will not work for those years.   2. My starting week / date is `2006-12-27`, how should I set the start parameter? Should I use: `start=c(2006,52) or start=c(2007,1)`, since week of `2006-12-27` really crosses the year boundary? Also, for modeling, is it better to have complete year worth of data (say for 2007 my start year if I only have partial year worth of data, is it better I should not use 2007, instead to start with 2008. What about 2014 since it is not complete year yet, shall I use what I have for model or not? Either way, I still have issue of whether or not to include those weeks in the year boundary like 2006-12-27, shall I include it as week 1 for 2007 or last week of 2006?   3. When I use `ts <\- ts(df, frequency=52, start=c(2007,1))` and then print it, I got results shown below, so instead of `2007.01, 2007.02, 2007.52, ...`, I got `2007.000, 2007.019, ...` which it gets from `1/52=0.019`, which is mathematically correct but not really easy to interpret. Is there a way to label it as the date itself just like data frame or at least `2007 wk1, 2007 wk2, ...`? =========              Time Series:     Start = c(2007, 1)      End = c(2014, 11)      Frequency = 52                Amount        Count     2007.000   645575.4     493717     2007.019  2185193.2    1659577         2007.038  1016711.8     860777     2007.058  1894056.4    1450101     2007.077  2317517.6    1757219     2007.096  2522955.8    1794512     2007.115  2266107.3    1723002         4. My goal is to model this weekly data, then try to decompose it to see seasonal component, it seems like I have to use `ts()` function to convert to `ts` object then I can use `decompose()` function, I tried `xts()` function, and I got error stating "time series has no or less than 2 periods" I guess reason is because `xts()` won't let me specify the frequency?               xts <- xts(df,order.by=businessWeekDate)        5. I looked for the answer in this forum and other places as well, most of the examples are monthly, there are some weekly time series question, but none of the answers are straight forward.


title: Help with JAGS model using matrix multiplication and Dirichlet prior

text: I want to set up a simple JAGS model using matrix multiplication, but I am not sure how to set up the matrices. My overall equation is $c=yx$ where $c$ is a 2x1 matrix defined in R of observed data values, $y$ is a 2x2 matrix of prior distributions, and $x$ is a 2x1 matrix of values to be estimated with the prior of the Dirichlet distribution. Here is my model file:               model{              c[] ~ dnorm(mu[],tau)         mu[] <- y[,] %*% x[,]              y[1,1] <- mN         y[1,2] <- fN         y[2,1] <- mO         y[2,2] <- fO              # priors         tau ~ dunif(0,1)         x ~ ddirch(alpha[1:2])         # d15N manure and septic source          mN ~ dunif(0,25)         # d15N fertilizer source         fN ~ dunif(-10,5)         # d18O manure and septic source         mO ~ dunif(-20,25)         # d18O fertilizer source         fO ~ dunif(-20,25)     }      When I try to compile this model I get the error: > Dimension mismatch in subset expression of c


title: statistical test for 3-preference choice design

text: We collected data from 40 animals in which we recorded the total time the animal spent in each of the 3 possible chambers. Thus, for each animal, we have three times (one for each chamber) and those 3 times add up to 10 minutes. Does anyone have any recommendations on statistical test that incorporates the magnitude of the differences in times? Much thanks.


title: Understanding (1-mean(L) - c*(stderr(L))

text: I inherited some old programming code for a scoring mechanism and I am trying to understand it. The premise is that every user $u$ has receives many ratings from other users, all in the range $[0, 1]$. We store all of user's ratings in a list $L$. The code then computes this "score" for the user $u$: $(1 - \mathrm{mean}(L)) - 1.96 \times \mathrm{stderr}{(L)}$ Can anyone give me any insight into what this technique is called, or where I can learn more about it? Any idea what purpose does the 1.96 serve?


title: Least-squares, pivoting/rotating around a reference point?

text: Suppose I have two datasets $y1$ and $y2$ sampled equally on $x$. I would like to bring $y2$ in register with $y1$ in a least-squares sense, and referenced to one point at $x_{ref}$ in $y1$ (i.e. the point $y2$ at $x_{ref}$ should be the same as $y1$ at $x_{ref}$. How would go I about doing that?


title: Probability exercise

text: I was doing one of the Papoulis problems and got a bit confused if I got it right, because as you will see in the end, the conditions given in the problem can be relaxed and still be obtained. The problem is to calculate $P(A\overline{B}+B\overline{A})$ given that $P(A) = P(B) = P(AB)$. The Papoulis solutions that I have available states:   1. $P(A) = P(AB+A\overline{B}) = P(AB) + P(A\overline{B})$ (given that $AB$ and $A\overline{B}$ are mutually exclusive)    2. $P(B) = P(BA+A\overline{A}) = P(BA) + P(B\overline{A})$ (given that $BA$ and $B\overline{A}$ are mutually exclusive) therefore, if $P(A) = P(B) = P(AB)$, from eq. (1) we have $P(A\overline{B})=0$ and from eq. (2) $P(\overline{A}B) = 0$, and we conclude that $P(A\overline{B}+B\overline{A})=0$ The issue is that when I developed my answer, I found a way that P(AB) was not needed for this effect to happen. I found that $P(A\overline{B}+B\overline{A})=0$ if only P(A) = P(B), and I would like to know if it is ok what I did: 3.$P(A\overline{B}+B\overline{A})=P(\overline{A}B+A\overline{B})=P(A\overline{B})+P(\overline{A}B) $ (it is quite easy to see that $A\overline{B}$ and $\overline{A}B$ are mutually exclusive, but it is also provable: Let $\overline{A}B=C$ and $A\overline{B}=D$, they are mutually exclusive if, and only if, CD=$\emptyset$. We have $CD = \overline{A}B\overline{B}A = \overline{A}\emptyset A = \emptyset$, thus $\overline{A}B$ and $A\overline{B}$ are mutually exclusive) We can develop that: 4.$P(A\overline{B})=P(\overline{\overline{A\overline{B}}})=P(\overline{\overline{A}+B})$, and we know that 5.$P(\overline{K})=1-P(K)$, then we can rewrite (4) as: 6.$P(\overline{\overline{A}+B})=1-P(\overline{A}+B)=1-(P(\overline{A})+P(B)+P(\overline{A}B))$. ( **edit** : as stated, everything below is wrong due to the fact that this equation should be $P(\overline{\overline{A}+B})=1-(P(\overline{A})+P(B)-P(\overline{A}B))$ leading (7) to $P(A) - P(B) + 2\times P(\overline{A}B)$) Using (6) in (3) we have: $P(A\overline{B}+B\overline{A})=\overbrace{1-(P(\overline{A})+P(B)+P(\overline{A}B))}^{P(A\overline{B}) \text{ rewritten}}+\overbrace{P(\overline{A}B)}^{P(B\overline{A}) \text{ rewritten}}=1-P(\overline{A})-P(B)$ and we can rewrite that using (5) to obtain: 7.$P(A\overline{B}+B\overline{A})=\overbrace{1-P(\overline{A})}^{P(A)}-P(B)=P(A)-P(B)$. Please observe that the equation (7) could be obtained as: 8.$P(A\overline{B}+B\overline{A})=P(B)-P(A)$, if we had used $P(\overline{B}A)$ in (4). The correct equation to be used is the one that obey the first axiom ($P(A\overline{B}+B\overline{A})>0$), that is eq. (7) when $P(A)>P(B)$ and (8) when $P(B)>P(A)$. So, if $P(A) = P(B)$, we have that both (7) and (8) lead to 0. Therefore, the condition for $P(A)=P(B)=P(AB)$ can be relaxed to $P(A)=P(B)$. The case when $P(A) = P(B) = P(AB)$ leads to $P(A\overline{B}+B\overline{A})=0$, but also to $P(A\overline{B})=0$ and $P(B\overline{A})=0$, which also means that $A=B$ (that is, A and B are the same subset from $S$)!


title: How to test group differences on a five point variable?

text: I have a series of observations that fall into bins (or "scores"); that is, the data can be 0, 1, 2, 3 or 4. There are two groups of such data, control and treated. I know the number of individuals with each score for each group. **What is the best way to determine whether these groups are different or not?** A colleague suggested just arranging the data as individual data points with the given score, and doing the analysis on those two columns of data. Since there are ten individuals per group, this is not difficult, but I do not believe that I am getting a valid answer.


title: Normal Data Distribution (Visual tests accept and statistical test ks reject normal distribution) Need Help

text: I have four variables want to run regression while i check for data distribution i found that histogram and qq plot provide evidance of normal data distribution where as ks test is significant for all variables. one more thing i see there are no negative values in data set all are positive can it be the reason of ks fail. please look at graph and ks test results below and give me suggestions how can i proceed. Can i run regression or not thanks **Histogram and QQ plot for PDM** ![enter image description here](http://i.stack.imgur.com/fmbhl.jpg) ![enter image description here](http://i.stack.imgur.com/7EjJF.jpg) **Histogram and QQ plot for MEOC** ![enter image description here](http://i.stack.imgur.com/cBTst.jpg) ![enter image description here](http://i.stack.imgur.com/Yavfq.jpg) **Histogram and QQ plot for CLT** ![enter image description here](http://i.stack.imgur.com/PJs1j.jpg) ![enter image description here](http://i.stack.imgur.com/ZuI9C.jpg) **Histogram and QQ plot for CTY** ![enter image description here](http://i.stack.imgur.com/5tOOC.jpg) ![enter image description here](http://i.stack.imgur.com/OflCz.jpg) **KS test** ![enter image description here](http://i.stack.imgur.com/gujPy.jpg)


title: Why use G-test and the likes for AB testing at all?

text: I am learning about AB testing using the G-test. In my example, I have a 2x2 contingency table.               >print(T)            response     AB    no yes Sum       A   29   7  36       B   23  16  39       Sum 52  23  75      Event A is the red background of a website. Event B is the blue background of a website. I showed the website to the total of 75 people. Yes is like, and no is the opposite. After running the G-test I get               >likelihood.test(T)          Log likelihood ratio (G-test) test of independence without         correction          data:  T     Log likelihood ratio statistic (G) = 4.1914, X-squared df = 1,     p-value = 0.04063      The p-value is pretty small (significant at the 5% level), so I reject the null that the samples A and B have the same performance. Now, I have two questions:   1. How do I know what background color is better and has the higher performance?   2. Why do I need the G test at all? I can just compare the percentage of likes for each background color. For red (A), 7/36=19.4%, and, for blue background (B), 16/39=41%. So, clearly B has the higher like percentage, hence, is better. So, why use G test at all? Remark: I also use the Fisher exact test since one of the measured values, A_yes, is smaller than 10. The output is               > fisher.test(T)              Fisher's Exact Test for Count Data          data:  T     p-value = 0.04948     alternative hypothesis: true odds ratio is not equal to 1     95 percent confidence interval:      0.9167091 9.6380957     sample estimates:     odds ratio        2.841057       The p-value is nearly identical to the G-test's. So again why use the G-test or Fisher test at all when one can just compare the yes percentage for each event. Thanks.


title: Dropping predictor variables, based on variable of importance, effect of Random Forest Accuracy

text: I am trying to use Random Forest to accurately predict forested land cover classes using Landsat 7, climatic and geographical data. I have 23 predictor variables and 1 response variable. When I drop the least important predictor variables as reported by the Random Forest variable of importance list, I get what seems to be an overly important influence of the remaining predictor variables. That is, I can see obvious misclassification between forest classes based on, for example, the climate data, when viewed in GIS. This results in lower accuracy. My thinking is that the model should improve when the least important predictor variables are removed. **Why does my Random Forest model accuracy decrease when low importance predictor variables are dropped from the model as compared to leaving them in the model?**


title: Regression coefficients with longitudinal data yields different results in R and SAS

text: I have a question about SAS and R. For a research, I used a longitudinal data and I initially used SAS (`GLIMMIX`) and then I analyzed the data with R (`glmer`) programming. There are differences between p-values of SAS and R. I expected that regression coefficient and standard error could be different for R and SAS. But there are differences for p value for some variables, which are significant in R, are not significant in SAS. My R model and SAS model are respectively :               #R     m3.glmm <- glmer(y ~ timebefore + timeafter + x1 + x2 +...+ x11 +                                (1+timebefore+timeafter|id),                       data=data, family=binomial(link="logit"), nAGQ=3)          #SAS     proc glimmix data=data METHOD=QUAD(QPOINTS=3) NOCLPRINT ;       class id x2 x3 x4 x5;       model y(event='1')=timebefore timeafter x1 x2 x3 x4 x5              x6 x7  x8 x9 x10 x11 /solution CL link = logit dist = binary;       random intercept timebefore timeafter/subject = id GCORR SOLUTION;     run;      Eg: variable "x1"(defined as age) was significant (p val= 0.04) in SAS but not in R (p val=0.1). But others were similar. It means that significant variables in SAS are found significant in R, or insignificant variables in SAS are insignificant in R. Does anybody know about the differences?


title: Kernel density estimation (KDE)

text: Due to an assignment I need to implement a algorithm based on _KDE_ to schedule an input data in different servers. So far, I studied statistics in my bachelor but we did not go that far and they did not explain me about estimation functions and Kernel. So my questions are the next: Assuming this unidimensional population:               1 2 5 5 3 4 5 6 7 8 1 3 9 6 7 5 3 1 2 3 5 7 8 2 1 3 4 5 1 7 6 3        1. What is a kernel and why is useful to estimate the trend in a population?   2. As far as I know, _KDE_ can use different kernel functions, such as triangle, uniform, etc... How do I know which one should I choose?   3. What is the _bandwidth_ smoothing parameter?


title: Data mining, domain knowledge and visualization

text: In my book ( _Introduction to data mining_ by Tan, Steinbach and Kumar - chapter 3), in the section about visualization, it's written: > Another general motivation for visualization is to make use of the domain > knowledge that is "locked up in people's heads." While the use of domain > knowledge is an important task in data mining, it is often difficult or > impossible to fully utilize such knowledge in statistical or algorithmic > tools. In some cases, an analysis can be performed using non-visual tools, > and then the results presented visually for evaluation by the domain expert. > In other cases, having a domain specialist examine visualizations of the > data may be the best way of finding patterns of interest since, by using > domain knowledge, a person can often quickly eliminate many uninteresting > patterns and direct the focus to the patterns that are important. I don't fully understand it. What's meant by _domain knowledge_? Who is a _domain expert_ or specialist? Could someone explain the paragraph in simple terms?


title: How to find important features in this problem

text: I was just thinking how ML techniques can be applied in the retail industry. Suppose we have data from a retailer who deals with apparel and cloth in this format and for each item there are pre-defined features, for example a shirt will have features like color, half/full sleeves etc.. I want to understand how to extract important features of "most sold items" from this data. For example if I find that shirts are the most commonly sold item, then which kind of shirt : Black in color, Half Sleeves, etc... I was thinking to use decision trees here. I am not sure if this is a good approach.


title: When can we look at the significance of pairwise comparisons rather than an interaction term?

text: I will defend my Master's thesis on Thursday, and I have a doubt about an analysis that I have to present. In my experiment, I had two independent variables:   1. Age category (SU = senior unemployed subjects and YU = young unemployed subjects)    2. Experimental condition (ST = stereotype threat; NST = No threat condition). My dependent variable is the performance achieved by subjects on a memory test. My hypothesis was that there would be an effect of Experimental condition only on SU subjects. Consequently, I performed a two-way 2x2 ANCOVA and expected the interaction between Age category and Experimental condition to be significant. Unfortunately, the interaction turned out to be non significant (p=0.615). However, when I looked at pairwise comparisons, only SU subjects' performance varied (almost) significantly from one experimental condition to another (p=0.058), whereas YU did not vary significantly across experimental conditions (p=0.213). This helped me confirm my hypothesis and conclude that there was an effect of the experimental condition only on SU participants. However, I'm not quite sure about the conditions under which one is allowed to look directly at pairwise comparisons (as I did) and overlook the global interaction is not significant. I don't have any textbook at hand, so if anyone could indicate me a stat published article arguing in favor of this method and indicating the conditions under which it is doable, I'd be very grateful!


title: Probability statistics exercise, with not enough data?

text: I have a question connected with this exercise. In question a) I am supposed to calculate the probability using this formula $P(A \text{ and } B)=P(A) \times P(B|A)$. I think there is not enough data. Can anybody please help me? In March 2005, the Australian government announced that an extra 450 troops would be sent to Iraq. A special Morgan poll was conducted to evaluate public opinion on Australia’s military presence in Iraq in the middle of March 2005 using an Australia-wide representative sample. It was found that 0,85 % of the respondents did not approve the decision to send 450 extra troops to Iraq. Of the respondents who are Liberal party voters, 0,18% approved the decision to send an extra 450 troops to Iraq. Assume that, of the people surveyed, 0,84 % were not Liberal party voters. One Australian is selected randomly. a) What is the probability that the person is not a Liberal party voter and approves the decision to send an extra 450 troops to Iraq? b) What is the probability that the person is not a Liberal party voter or disapproves the decision to send an extra 450 troops to Iraq?


title: Residualize a binary variable to remedy multicollinearity?

text: Imagine a regression model where there is a continuous-valued response variable and three continuous-valued explanatory variables. For concreteness, imagine that we are interested in the effects of "Predictability," "Length," and "Frequenty," on the reading times "RT" of words. Suppose further that we know that there is collinearity between "Frequency" and the other two explanatory variables. A way of dealing with this in the model is to residualize "Frequency," along the lines of the R code given (where the second line is the form of the model inputted in to a suitable regression algorithm):               r. <- function (formula, ...) rstandard(lm(formula, ...))     RT~Predictability+Length+r.(Frequency~Length+Predictability)      Now suppose that one of the the variables of interest was binary-valued. Let's say that we wanted to run a model with the form:               RT~Education+ReadsDaily      but we found that neither variable had a significant effect when both were included, but when a saturated systematic component, `Education*ReadsDaily` was used, the model found highly significant results for both coefficients as well as the interaction term. When the two models are compared, the inclusion of the interaction term decreases the deviance by `16.7`. When a VIF analysis is run on the model with an interaction terms, very high values are reported:               ReadDailyY            Education          ReadDailyY:Education      8.693957             4.266084            15.607665       The model without an interaction term, which fits poorly, has low values in the VIF analysis:               ReadDailyY  Education      1.227842    1.227842      I _think_ this means that there is multicollinearity between the two predictors (cf. Agresti 2007:138): > ...models with several predictors often suffer from multicollinearity – > correlations among predictors making it seem that no one variable is > important when all the others are in the model. A variable may seem to have > little effect because it overlaps considerably with other predictors in the > model, itself being predicted well by the other predictors. Deleting such a > redundant predictor can be helpful, for instance to reduce standard errors > of other estimated effects. We would, however, like to try to separate their effects on "RT" rather than combining the variables or discarding one, if at all possible, since the distinction between the two variables happens to be of theoretical importance. It would be possible to do something like the following in R with no error messages:               r. <- function(formula) rstandard(glm(family=binomial(),formula))     RT~Education+r.(ReadsDaily~Education)      The residualized "ReadsDaily" parameter in the model would then be continuous- valued. I suspect that this strategy for dealing with collinearity is methodologically suspect, but I do not know enough to give a reason for why or why not the strategy is to be dispreferred. Agresti, A. (2007) An Introduction to Categorical Data Analysis. Wiley.


title: How to refer to variables which lie beyond causal pathway

text: Causal diagrams are extremely good tools for discussing research plans for multivariate modeling between statisticians and non-statisticians. It's easy after some deliberation to decide which variables should go into a model based on their classification as confounders, precision variables. At least this applies insofar as identifying additive models, I'm unsure what language/visualizations applied to identifying interactions graphically. However, it's generally accepted that adjusting for variables which are _caused_ by the outcome will bias the estimated causal relationship between the main effect and outcome of interest. For instance, if you're estimating links between smoking and obesity, researchers might consider adjusting for forced expiratory volume (don't ask me why... maybe because of it's relationship to exposure). However obese people might have less FEV because of the compressing effect of fat on the chest cavity and also lack of physical exercise. So FEV is causally determined by obesity. We show that FEV lies beyond the causal pathway. Is there a term that refers to such variables?


title: Applying Spearman's rank correlation on average values, does it make sense?

text: I have a sample of 2000 projects randomly extracted from a population of 10,000,000 following an exponential distribution. From this sample I would like to check if a statistical dependency exists between the number of lines of code and the number of users. Instead of comparing these two features with Spearman's rank correlation directly on the 2000 projects, I've grouped the projects in the sample according to the number of users (obtaining 10 groups), and for each group I've calculated the average of lines of code. The corresponding rho and p-values show that there is a significant dependency between such variables. On the contrary, if I run the Spearman's rank correlation on the 2000 projects, this dependency does not appear. Can I say that:   1. the 2000 projects are representative enough with respect to the population?   2. there exists a dependency between number of users and lines of code, only if the number of lines of code for a given project is above the average?


title: Summarizing confidence intervals when there are many levels

text: How can I summarize the confidence intervals for several levels of factors? This is what I got:               interest rate (1.280393e+01 16.1980844628);      fico range 645-649 (-3.581231e+00  1.9591171919);      fico range 650-654 (-4.533451e+00  3.7764777115);          fico range 655-659 (-2.797205e+00  2.2914947096);          fico range 660-664 (-2.822852e-01  3.1803914727);          fico range 665-669 (-9.659140e-01  2.4864556261);          fico range 670-674 (-1.954018e+00  1.4892769161);          fico range 675-679 (-2.694203e+00  0.7515725574);          fico range 680-684 (-3.457226e+00 -0.0085546990);          fico range 685-689 (-3.965619e+00 -0.5093765443);          fico range 690-694 (-4.283243e+00 -0.8263499212);          fico range 695-699 (-4.808411e+00 -1.3570095978);          fico range 700-704 (-5.564438e+00 -2.1040755402);          fico range 705-709 (-6.290848e+00 -2.8315079720);          fico range 710-714 (-6.422048e+00 -2.9508840053);          fico range 715-719 (-7.563818e+00 -4.0779259228);          fico range 720-724 (-7.744044e+00 -4.2745503539);          fico range 725-729 (-8.117765e+00 -4.6328158788);          fico range 730-734 (-8.919313e+00 -5.4341114666);          fico range 735-739 (-9.447742e+00 -5.9220440264);          fico range 740-744 (-9.497308e+00 -5.9439118677);          fico range 745-749 (-9.437903e+00 -5.8873782684);          fico range 750-754 (-1.001172e+01 -6.4804655589);          fico range 755-759 (-1.017971e+01 -6.6033187178);          fico range 760-764 (-1.000456e+01 -6.4294990859);          fico range 765-769 (-1.063397e+01 -7.0112845774);          fico range 770-774 (-1.171413e+01 -7.8514910920);          fico range 775-779 (-1.051257e+01 -6.7492730448);          fico range 780-784 (-1.121741e+01 -7.5319951164);          fico range 785-789 (-1.081962e+01 -7.0030665863);          fico range 790-794 (-1.121984e+01 -7.4234794495);          fico range 795-799 (-1.138416e+01 -7.3877810394);          fico range 800-804 (-1.103957e+01 -7.0461654289);          fico range 805-809 (-1.204510e+01 -7.9475491304);          fico range 810-814 (-1.105703e+01 -6.7223043818);          fico range 815-819 (-1.204715e+01 -7.4519985436);          fico range 820-824 (-1.183545e+01 -3.5254524137);          fico range 830-834 (-1.393567e+01 -5.6193498315);          loan length60 months  (3.066167e+00  3.4715640126);          amount requested (1.330939e-04  0.0001547571);      As you can see fico range has a lot of levels. What is the appropriate way to report them?


title: How to report multiple experimental observations each with their own error.

text: Hi I have experimental observations that look like this. ![enter image description here](http://i.stack.imgur.com/XF2dM.png) The output is not a function of concentration, so I thought within each experiment they could just be considered repeated observations. However the machine figures out this "output" by fitting a line to the curve and gives you a error of that fit (referred to here as output error). My question is, how do I report one value of output with some type of error including all these values? Can I just pool X1 - X9 and take a SEM? Or will it be more complicated? Please let me know if I'm being to vague or you need more information from me. Here are similar questions here and here but each assumes that you have one observation and one error from independent observations. I have three and am wondering how to incorporate them. Thanks!


title: McNemar's test for $3 \times 2$ data

text: I adopt a within-group experiment, such that the same participant uses 3 devices to perform a task, which has a result "pass" or fail. My data looks like this:                       Pass  Fail     deviceA   21    13     deviceB   15    20     deviceC   9     25      My questions are:   1. Is McNemar's test appropriate for this data?   2. In R, when I run `mcnemar.test(data)`, I got the following error:              "'x' must be square with at least two rows and columns"      How can I run this $3 \times 2$ data for R's `mcnemar.test`?


title: How to find Influence function?

text: > Derive $IF(x;T,F)$ when $$\displaystyle > T(F)=\int_{F^{-1}(\alpha)}^{F^{-1}(1-\alpha)}x ~dF(x)$$ Here $IF$ stands for > Influence function. Trial: Here $$\begin{align}IF(x;T,F) &=\lim_{t\to 0}\frac{T((1-t)F+t\Delta_x)-T(F)}t \\ &=\lim_{t\to 0}\frac{g(t)-g(0)}t=\frac{d}{dt}g(t)|_{t=0} \end{align}$$ Then I try to simplify $T(F)$ as $$\int_{F^{-1}(\alpha)}^{F^{-1}(1-\alpha)}x ~dF(x) \\ =\int_{\alpha}^{1-\alpha}F^{-1}(y) ~dy$$ Then I am stuck. Please help.


title: Comparisons of two rankings

text: There are two methods which produce different ranking of a group of hypotheses. I want to compare the two methods showing that one method generally produces lower ranking than the other method for a specific subset of the considered group of hypotheses. I was thinking of comparing the quantiles of the two rankings on this subset, but not sure whether that was the right thing to do. Can anyone familiar with this give some hint? Thanks. Hanna


title: What is a quadravariate distribution?

text: **Question:** What is a quadravariate distribution? **Motivation:** I found a reference to quadravariate distributions in Priestley (1981) _Spectral Analysis and Time Series_ on p325, but no definition (the author seems to think it is common knowledge). 10 minutes of googling revealed that the term gets used in quite a few papers, but I am still unable to find a definition anywhere.


title: Relevance of overall absolute values in covariance analysis of two variables

text: I am performing K means clustering on a gene expression dataset. I am aware of the fact that the Pearson correlation metric allows to group trends or patterns irrespective of their overall level of expression. I was wondering if the same concept stands for Covariance metric (I believe that the only difference between the two metrics is the fact that covariance returns unbounded values, while Pearson maps value in interval [-1,1])


title: What is a good reference book for the matrix algebra linear algebra used a lot in applied statistics?

text: I have been working in R for a bit and have been faced with things like PCA, SVD, QR decompositions and many such linear algebra results (when inspecting estimating weighted regressions and such) so I wanted to know if anyone has a recommendation on a good comprehensive linear algebra book which is not too theoretical but is mathematically rigorous and covers all of these such topics...


title: RMSE, MAD of vectors

text: I have a set of $N$ high-dimensional vectors. I use some approximation routine to make my output faster. Now I would like to evaluate the error of the approximation. Typically I use the RMSE to compute errors, but I'm getting stuck here on how to do it right. For each vector $\vec x_i \in \mathbb R^D$ I have approximation $\vec y_i$. I would compute $RMSE = \sqrt{\frac{\sum_i^N (\|\vec x_i - \vec y_i\|_2)^2}{N}} = \sqrt{\frac{\sum_i^N \sum_j^D (x_{i_j} - y_{i_j})^2}{N}}$. Or is it more reasonable to scale like this: $\sqrt{\frac{\sum_i^N \sum_j^D (x_{i_j} - y_{i_j})^2}{N\cdot D}}$? Or would it make more sense to just do $\frac{\sum_i^N \|\vec x_i - \vec y_i\|_2}{N}$ and how would that be called? How would I compute the maximum absolute deviation (MAD)? $MAD = \arg\max_i \|\vec x_i - \vec y_i\|_2$ or $MAD = \arg\max_i \|\vec x_i - \vec y_i\|_\infty$? I know that MAD typically stands for median absolute deviation, is there a better acronym?


title: Equal frequency vs. equal weight discretization

text: How does one choose between the two? How do you decide which one is more appropriate for your data? I have a numeric data and I would like to discretize it using one of these. I know what they do but I have no idea which one I should pick and why.


title: Classification algorithm robust to variance?

text: tldr version: What's an algorithm I can use to classify different series with very different standard deviations across a single variable? I'm working with data that follows changes of stock prices through non- sequential, uncorrelated, instances of time. I would like to group stocks together that tend to change price at the same time. My first thought was K-nearest neighbors, but I would need it to be radial and the different variances make it difficult to set a threshold radius. Maybe I just need a good way to normalize the data? Classifying based on a correlation matrix seems like a good idea to me. Most of the data could be sorted by simply setting a threshold correlation, but then there are some stocks that belong to more than one group and others that don't belong to any groups. Rather than try to design something from scratch, I'm guessing there's a widely used algorithm that would work well for my purpose. Any suggestions are welcome. Edit: In case it wasn't clear, I'm working with a single variable, which is price.


title: Stationarity of $\nabla x_{t}$

text: How should I prove that for a given series $x_{t}=\frac{3}{2}x_{t-1}-\frac{1}{2}x_{t-2}+w_{t}$ that is non-stationary, $y_{t}=\nabla x_{t}$ is stationary? I attempted to get an equation based on the lag operator $B$, defined as $\nabla x_{t} = (1-B)x_{t}$ but I can't obtain roots greater than 1 (which is the condition for stationarity). However, since I have only seen examples dealing with $x_{t}$ instead of $\nabla x_{t}$, it's very likely I'm missing something. What is the correct procedure in this case? **UPDATE** : whuber helped me a lot with this question. Now I have another question as I posted in a comment below: Is there a general procedure when it is not possible to express $x_{t}$ in terms of $y_{t−i}$'s for some $i$?


title: References for weighted linear regression with linear constraints on the coefficients?

text: I would like to estimate a multiple linear regression model $N$ observations (with $\beta$ of length $k$): $$Y = X \beta + \epsilon$$ subject, however, to some linear constraints on the coefficients. I.e a constraint of the form $$M \beta = 0$$ where (at least in my case) $M$ is of dimension $2 \times k$ but more generally can be expressed as $$M \beta = c$$ where $M$ is of dimension $p \times k$ with $p < N$ and $c$ of dimension $p \times 1$. I have not been able to find much about this online and would appreciate references to read up on this.. Thanks


title: What is the definition of covariance under the fixed-effects assumption for a meta-analysis of correlation?

text: There is a doubt about the meaning of covariance and its definition. Does it change under different conditions such as primary data or statistic such as r in case of meta-analysis?


title: ARIMA with difficult seasonality in R

text: I have non-stationary time series. It has evident trend in means and seasonality. These raw data are measured every second. On the plot of original series I see trend and seasonality about 80,81 (after each 80 second we can see new growth). It is possible that I will never obtain model with white noise? I have tried first difference and first seasonal difference on 80 and second seasonal difference with 81 pediod, and twice of 81. And stiil I have regular data in residuals! It is possible?? I have used R. I have tried different ways: step-by-step, ARIMA with different settings for differences, auto.arima. I do not know why I cannot receive white noise in the residuals. I look on ACF of residuals. Has the task got solution? Below I present my data.               3785     4542     4542     4542     4542     4542     4542     4542     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     4731     5772     8043     8516     8516     8516     8516     8516     8516     8516     8516     8516     8516     8516     8516     8516     8516     8516     8705     8705     8705     8705     8705     8705     8705     8705     8705     8705     8800     8800     8800     8800     8800     8800     8800     8800     8800     8800     8800     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     8989     10314     12585     13058     13058     13058     13058     13058     13058     13058     13058     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     13247     14288     16843     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17032     17316     17316     17316     17316     17316     17316     17316     17316     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     17505     18073     21101     21290     21290     21290     21290     21290     21290     21290     21290     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21574     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     21763     25359     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25548     25832     25832     25832     25832     25832     25832     25832     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     26021     29333     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30090     30279     30279     30279


title: Simulating financial strategies with given value-at-risk (VaR) and mean

text: I would like to generate thousands of distributions with respect to: VaR@5%=-7% and Median=0%. A first idea:   1. Choose a family, say the Gaussian one.   2. Estimate its parameters. This is straightforward for the Gaussian but let's write the general algorithm (see mpiktas' answer on Estimating a distribution based on three percentiles):              library(nleqslv)     VaR_Threshold <- c(0.05,0.5)     VaR_Value <- c(-0.07,0)          ufn <- function(x,q) q-pnorm(VaR_Value, x[1],x[2])     usol <- nleqslv(c(0,1), ufn,q=VaR_Threshold)     usol$x          [1] -1.262016e-11  4.255698e-02 (estimated parameters)          plot(x<-seq(-1,1,by=0.01),dnorm(x,usol$x[1],usol$x[2]),type="l",col=2)     points(p<-VaR_Value,dnorm(p,usol$x[1],usol$x[2]))          pnorm(VaR_Value,usol$x[1],usol$x[2])          [1] 0.05 0.50 (check)      So this is for the Gaussian family: one unique solution. My question: is there a parametric family flexible enough to generate thousands of distributions with a given VaR and median? I tried with the stable distributions (4 parameters : 2 manually fixed + 2 optimized), in vain, probably because I don't master it...


title: What is the best software for time series analysis and forecasting?

text: Is MATLAB better than R for time series analysis and forecasting or vice versa? What other software is considered best for time series analysis?


title: What do the residuals in a logistic regression mean?

text: In answering this question John Christie suggested that the fit of logistic regression models should be assessed by evaluating the residuals. I'm familiar with how to interpret residuals in OLS, they are in the same scale as the DV and very clearly the difference between y and the y predicted by the model. However for logistic regression, in the past I've typically just examined estimates of model fit, e.g. AIC, because I wasn't sure what a residual would mean for a logistic regression. After looking into R's help files a little bit I see that in R there are five types of glm residuals available, c("deviance", "pearson", "working","response", "partial"). The help file refers to Davison, A. C. and Snell, E. J. (1991) Residuals and diagnostics. In: Statistical Theory and Modelling. In Honour of Sir David Cox, FRS, eds. Hinkley, D. V., Reid, N. and Snell, E. J., Chapman & Hall, of which I do not have a copy. Is there a short way to describe how to interpret each of these types? In a logistic context will sum of squared residuals provide a meaningful measure of model fit or is one better off with an Information Criterion?


title: Use of expectation propagation for inference

text: I have a PGM as described by the attached diagram. ![Graphical model](http://i.stack.imgur.com/t8LSK.png) $y$ is observed and I want to infer the joint posterior distribution as given by $P(w,\lambda, \phi |y)$. Distribution on $w$ is modelled using a multivariate 0 mean with a covariance structure. $\lambda$ needs to be inferred from the data and the prior on it is modelled using a gamma distribution with scale and shape parameters $s_0$ and $C_0$. The noise variance $\phi$ also needs to be estimated and its prior is also modelled using a Gamma distribution with parameters $a_0$ and $b_0$ I have been looking at using Expectation Propagation to perform this approximate inference on this model. As I understand, EP works by minimising the KL divergence between the true posterior and the approximated distributions. So, my questions are: 1: Can I use EP in such a case? 2: Can someone help me understand what objective function I need to design in this model? As far as I can tell, the EP will approximate the try posterior by independent distributions from the exponential family. So, $P(w,\lambda, \phi) \approx q(w) q(\lambda) q(\phi)$ And now do we minimise the KL divergence between true p and q's? I am really having trouble understanding this setup and would be really grateful if someone can help me understand where and how to start! Any help/pointer/comments will be appreciated. Thanks, Luca


title: Can complex Kruskal-Wallis results be reported in tables for APA-style publication?

text: Advice? I ran a Kruskal-Wallis H test in SPSS and have a lot of variables to report for a paper for publication in APA format with limited word count. Can I show all the paired comparisons in a table? (the independent variable has 4 groups, the dependent variable has 5 groups; that`s a lot of statistics). If so, would it be sufficient to report the Chi-Square and p-value in one table, and post- hoc median and p-values for paired comparisons in a second, more complex table? Or would this be frowned upon, and why? Thanks!