title: SAS - REG vs GENMOD; OLS vs MLE

text: I'm using a very simple data set from an article in trying to further my understanding of GLMs. I've input the data using SAS, and I've run both the PROC REG and PROC GENMOD procedures on the data. In the PROC GENMOD procedure, I used a log link with a normal distribution; in the PROC REG procedure, I used the log of the response variable in the model. My question is, why don't the parameter estimates of the two procedures match? My understanding is that PROC REG uses OLS to estimate the parameters, whereas PROC GENMOD uses MLE with a Newton-Raphson iterative process for estimation. But I had thought that, when the assumed distribution is normal and the relationship is linear (which, after the log transformation, it is in the GLM, right?), MLE is equal to OLS. Here are the resulting parameters from the run:                    REG     GENMOD     A1   4.623    4.579      A2   4.688    4.730      A3   4.654    4.654      B1   (0.735)    (0.741)     B2   (0.487)    (0.436)      And here is my code:               data GLM;      input Y A1 A2 A3 B1 B2;      lnY = LOG(Y);      datalines;     95 1 0 0 0 0     115 0 1 0 0 0     105 0 0 1 0 0     55 1 0 0 1 0     45 0 1 0 1 0     30 1 0 0 1 1     ;           proc genmod data=GLM;      model Y = A1 A2 A3 B1 B2 / dist=normal link=log scale=deviance noint ;      weight Y;     run;          proc reg data=GLM;      model lnY = A1 A2 A3 B1 B2 / noint;      weight Y;     run;      Any insight that anyone can contribute is greatly appreciated! Bonus question - in my data I have 6 equations and 5 variables. Why is an iterative process needed to solve that?


title: correlation analysis with collapsed subjects

text: I conducted an experiment with 21 subjects. Each subjects performed 80 trials for 3 different categories of test (total 240 trials). What I need to do now is look for a correlation between category 1/2/3 for each of the 80 trials and I need to do that by collapsing all subjects data as if they were produced by one single subject. Does anyone have any idea on how to do this? Also I have a categorical dependent variable for each of the 80 trials, how can I do a regression which again collapsed all the subjects into one and uses category 1/2/3 as independent variables? Finally, in that case do I need to use a regressor which take into account the fact that data where produced by 21 different individuals? Thanks Best wishes, Andrea


title: State-dependent independence in graphical models

text: I have a graphical model with binary variables Y, X1, X2, and observed data D. D depends on Y X1 and X2 depend on Y When Y is false, X1 and X2 are independent:               (D) <-- (Y) --> (X1)              |              V             (X2)      However, when Y is true, X1 and X2 are not independent and must be modeled in a joint manner:               (D) <-- (Y) --> (X1,X2)      I've tried putting in an "edge" variable E between X1 and X2 that is influenced by Y; I can get the correct behavior when Y is false, but not when Y is true:               (D) <-- (Y) --> (X1)              |       ^              V       E            (X2) <- /      (that is Y --> E and E --> X1 and E --> X2) If I reverse the directed edges so that they point from X1 and X2 to E, then the moral graph once again always treats X1 and X2 as joint. So I can't seem to encode state-dependent independence... Is there established theory about such networks (if so what is this relationship called)? Is there a way to encode or transform the graph with temporary variables or something so that, if you know given D, Y must be false, that you can use the first graph (where X1 and X2 are independent given Y), and have faster inference? Of course, I should point out that this is a toy example; the real one is huge. Thank you very much for your help. And sorry for the ascii graphs, I was raised by wolves.


title: Ideal learning sample in machine learning

text: I am constructing a model for the prediction of a binary (Yes/No) outcome. I have a learning sample that gives the machine 1500 examples of the "Yes" group and 500 example of the "No" group. Should I be using all the data I have for input to learn the machine? Would this be biased towards the "Yes"? I had the thought of giving 500 "Yes" and 500 "No" examples, but I am not sure if this is going to positively or negatively my future predictions. Thanks.


title: Non-Parametric Test controlling for third variable

text: There are two laboratories measuring identical samples. One of the lab's data is non-normal, therefore I am looking for a non-parametric test. After comparing overall differences in the dependent variable (continuous) using Wilcoxon Mann-Whitney, I wish to test for differences between labs based on a third variable. There are a number of third variables that I wish to test for, and some of them are categorical (2 levels, and 3 levels), and another is ordinal (by month). Which test would be appropriate? I read about the van Elteren test, but the SAS output did not make sense. Thanks in advance!


title: Which statistical test to test if counts from different populations are actually different? t-test versus chi-squared

text: I'm having a hard to deciding which of a two-sample t-test and a chi-squared test is more appropriate for these sorts of questions. Examples would be: 70 out of 120 americans like hotdogs. 50 out of 95 europeans like hotdogs. Do the two groups differ in their hotdog opinions? 1000 each of caucasians, african-americans, and asians were asked if they preferred red, white, or yellow onions (counts provided in a table). Do respondents from different ethnicities differ in their responses? I realize that those questions are worded slightly differently - I think it may be the wording that differentiates them? In either case, i'm confused about which of the two tests is more appropriate.


title: In R, why are results different for 1-tailed wilcoxsign_test, wilcox_test, and wilcox.exact with paired data?

text: I am trying to determine why I am getting different results with these different Wilcoxson test implementations in R version 2.15.2. I have paired data with some ties. I have read the response for What is the difference between wilcox.test and wilcox_test in R, but some of the p-values are extremely different in my data. Also, I do not understand why the coin package would have two functions (wilcoxsign_test and wilcox_test) to produce the same result. What am I missing? Here is my example data:               library('exactRankTests') # for wilcox.exact     library('coin')          data <- data.frame(         y = c(2770.00, 3160.00, 4120.00, 4510.00, 3320.00, 3170.00, 3340.00, 3810.00, 3760.00, 6350.00, 2720.00, 3740.00, 5210.00, 3330.00, 4230.00, 3490.00, 3138.07, 2630.88, 4058.70, 3521.11, 4941.09, 5762.31, 3565.89, 3517.91, 3413.32, 3415.98, 3439.96, 2602.11, 2659.36, 3099.79, 2820.00, 2830.00, 4310.00, 4010.00, 2780.00, 2730.00, 3130.00, 2700.00, 3510.00, 3460.00, 2470.00, 2920.00, 3230.00, 3370.00, 3290.00, 2380.00, 2845.69, 2137.58, 3477.96, 3128.84, 3117.77, 4949.78, 3061.60, 2942.57, 3149.46, 3067.10, 3164.60, 2135.01, 2275.32, 3154.66),         condition = factor(c(rep('A', 30), rep('B', 30))),         participant = factor(rep(1:30, times=2))     )      When trying to get exact results I compared wilcoxsign_test, wilcox_test, and the depreciated wilcox.exact.               wilcoxsign_test(y ~ condition | participant, data=data, alternative='less',          distribution = "exact")     ##  Exact Wilcoxon-Signed-Rank Test     ##      ## data:  y by x (neg, pos)      ##   stratified by block      ## Z = -4.577, p-value = 3.818e-08     ## alternative hypothesis: true mu is less than 0          wilcox_test(y ~ condition | participant, data=data, alternative='greater',          distribution = "exact")     ##  Exact Wilcoxon Mann-Whitney Rank Sum Test     ##      ## data:  y by     ##   condition (A, B)      ##   stratified by participant      ## Z = 4.15, p-value = 0.05913     ## alternative hypothesis: true mu is greater than 0          wilcox.exact(y ~ condition, data=data, alternative='greater', exact=TRUE,          paired = TRUE)     ##  Exact Wilcoxon signed rank test     ##      ## data:  y by condition      ## V = 455, p-value = 3.818e-08     ## alternative hypothesis: true mu is greater than 0      Even the approximate methods sometimes differ as seen below.               wilcox_test(y ~ condition | participant, data = data, alternative = "greater",          distribution = "approximate")     ##      ##  Approximative Wilcoxon Mann-Whitney Rank Sum Test     ##      ## data:  y by     ##   condition (A, B)      ##   stratified by participant      ## Z = 4.15, p-value < 2.2e-16     ## alternative hypothesis: true mu is greater than 0     wilcoxsign_test(y ~ condition | participant, data = data, alternative = "less",          distribution = "approximate")     ##  Approximative Wilcoxon-Signed-Rank Test     ##      ## data:  y by x (neg, pos)      ##   stratified by block      ## Z = -4.577, p-value < 2.2e-16     ## alternative hypothesis: true mu is less than 0     wilcox.exact(y ~ condition, data = data, alternative = "greater", exact = FALSE,          paired = TRUE)     ##  Asymptotic Wilcoxon signed rank test     ##      ## data:  y by condition      ## V = 455, p-value = 2.362e-06     ## alternative hypothesis: true mu is greater than 0     wilcox.test(y ~ condition, data = data, alternative = "greater", paired = TRUE)     ## Warning: cannot compute exact p-value with ties     ##      ##  Wilcoxon signed rank test with continuity correction     ##      ## data:  y by condition      ## V = 455, p-value = 2.481e-06     ## alternative hypothesis: true location shift is greater than 0     pairwise.wilcox.test(data$y, data$condition, p.adj = "none", paired = TRUE)     ## Warning: cannot compute exact p-value with ties     ##      ##  Pairwise comparisons using Wilcoxon signed rank test      ##      ## data:  data$y and data$condition      ##      ##   A         ## B 5e-06     ##      ## P value adjustment method: none      As a related question. How can I get the extract observed Wilcoxon statistic for wilcoxsign_test? I know you can do the following for wilcox_test, but it does not seem to work for wilcoxsign_test.               w_t <- wilcox_test(y ~ condition | participant, data = data, alternative = "greater",      distribution = "exact")     statistic(w_t, "linear")     ##            ## A 1128


title: Decision Tree - Splitting Factor Variables

text: I'm new to decision trees and I have some confusion about how factor variables and non-ordered character/string variables get handled in a split. Suppose I have a factor such as "tiny, small, medium, large, huge" where the levels are important. How does a decision tree try to find the best split? Will it only check the 4 obvious splits, or will it check splits for weird combinations like, "tiny or huge but not small medium or large"? Similarly, how does a decision tree check for a split for an unordered character variable such as "New Orleans, Birmingham, Jackson, Miami, Atlanta"? I'm using the rpart package in R as I try to learn this stuff, so any references to rpart's implementation would be helpful. Thanks


title: Dimensionality reduction for classifying textures with MATLAB

text: I am trying to apply dimensionality reduction on an a set of images (~3000 pixels) using Matlab's Dimensionality Reduction toolbox. However, I know very little about dimension reduction. So I tried several functions by trial and error. PCA returned a matrix with complex numbers, and the others froze MATLAB. Can I get some advice as to which method works good on images? Here are some of the images: ![enter image description here](http://i.stack.imgur.com/fi7kJ.jpg)


I know this is the general specification for a single group: $Y_i$ = $\beta_0$ + $\beta_1$*time + $\beta_2$*post + $\beta_3$*timepost where `time` is a continuous variable indicating the time (e.g., month) of every observation in the study; `post` is an indicator for whether the observation is after the introduction of the intervention; and `timepost` is time since the introduction of the intervention (0 if pre-intervention). What is the correct specification for a segmented regression of interrupted time series data with a non-equivalent comparison group?


title: Initial Probabilities of an HMM

text: I have a learned Hidden Markov Model (HMM) from a certain sequential data using Gibbs sampling. I have managed to obtain the transition probabilities (transition matrix) of the Markov chain and the params of the probability distributions of the hidden states. The only set of params that I have yet been able to obtain are the initial probabilities $\pi$ that specify where the Markov chain would start. So given the parameters that I have estimated, is there a way to figure out $\pi$? Or how can I specify the Gibbs sampler to find out $\pi$?


title: Why would concatenating feature vectors lead to better estimates?

text: I wish to estimate the state of a system from two separate and disparate observations. A simple approach that I have seen in some literature is to combine the feature vectors (observations) by simply concatenating them. So if the first observation is denoted $x \in \Re^{10}$ and the second is denoted $y \in \Re^5$, we obtain a single (composite) feature $z \in \Re^{15}$. I have only seen this technique being applied (seemingly blindly) without any discussion as to what is actually going on "under the hood". **Questions**   1. When does this approach make sense and when does it not make sense? It would be great to see an example of what is taking place and why the approach works or does not work.   2. Is there anything that must be done to the two features ($x$ and $y$) before concatenating them? (e.g. scaling/normalisation?)   3. Is there anything that might assist this approach? (e.g. dimensionality reduction)


title: Experiment design

text: We have x1 and x2 representing variations in temperature and pressure and the response is yield of penicillin. What is the practical reason for using the points (-1,-1), (-1,1), (1,-1), (1,1) rather than using (0,0) in place (-1,-1)? Is it because using (0,0) would force it through the origin? Also, is replicating those four points three times a sensible design?


title: AR1 and Law Of Iterated Expectations : No serial correlation

text: In the AR(1) model $y_{t}=\beta_{0}+\beta_{1}y_{t-1}+u_{t}$, assuming $E(u_{t-1}|y_{t-1},y_{t-2}...)=0$, how does the law of iterated expectations ensure that the errors must be uncorrelated: $E(u_{t}u_{s}|x_{t},x_{s})=0$?


title: Learning a value of a parameter u given "true" or "false" prediction for each data-point x

text: We have a data-point x and many classes. Let $P(c|x)$ the probability that $x$ is of class $c$. We note $c_1$ the most probable class for $x$ (i.e. $P_1=P(c_1|x)$ is the highest probability), $c_2$ the second most probable class for $x$ (i.e. $P_2=P(c_2|x)$ is the second highest probability ($P_1> P_2$)). Obviously, when $P_1$is close to P2 (i.e. $P_1P_2$ is small or closer to 0) we are not sure whether $c_1$ or $c_2$is the true class of $x$, then it is useful to ask for the true label of $x$. We define $u$ to be a threshold value for the difference between the highest and the second highest probability (a threshold for $P_1$-P2); thus, if $P_1-P_2 > u$ then we can (with some cost) ask for the true class of $x$ (let's note it $c_x$). $c_1$(the predicted class for $x$) is usually equals to cx, but may sometimes not be equal. Given this problem, I want to learn a good value for the parameter $u$. To do that, currently, I just set u at an initial value (e.g. $u=0.2$) and then adjust this value according to whether or not $c_1$ equals $c_x$: if $c_1 = c_x$then we get more confident and thus decrease the value of $u$ (e.g. $u=u-\epsilon$), otherwise (when $c_1$!= cx) we get less confident and thus increase the value of u (e.g. $=u+\epsilon$), where $\epsilon=0.01$ for example. **Question:** Is there any better why to "learn" a value for $u$ ? (Assuming that we can start with a hight initial value of u in order to get labelled data at the beginning, or assuming that I have a subset of labelled data that I can to use to learn the value of $u$).


How to estimate the out crosses rate in crops? in have experiment included; main plot (1- seeds from open field 2- seeds from free pollinators cage 3-seeds from pollinator cage), sub plot (5 genotypes) and 50 plants grown for each genotype in each plot. This in three replications. It easy to estimate variance between main plots and sub plots but my question about the genetic variance between plants in same plot. and this will be due to out cross.


title: PCA when the dimensionality is greater than the number of samples

text: I've come across a scenario where I have 10 signals/person for 10 people (so 100 samples) containing 14000 data points (dimensions) that I need to pass to a classifier. I would like to reduce the dimensionality of this data and PCA seems to be the way to do so. However, I've only been able to find examples of PCA where the number of samples is greater than the number of dimensions. I'm using a PCA application that finds the PCs using SVD. When I pass it my 100x14000 dataset there are 101 PCs returned so the vast majority of dimensions are obviously ignored. The program indicates the first 6 PCs contain 90% of the variance. Is it a reasonable assumption that these 101 PCs contain essentially all the variance and the remaining dimensions are neglectable? One of the papers I've read claims that, using a similar (though slightly lower quality) dataset than my own, they were able to reduce 4500 dimensions down to 80 retaining 96% of the original information. The paper hand-waves over the details of the PCA technique used, only 3100 samples were available, and I have reason to believe less samples than that were used to actually perform PCA (to remove bias from the classification phase). Am I missing something or is this really the way that PCA is used with high dimensionality-low sample size dataset? Any feedback would be greatly appreciated.


title: 2x2x2 plus control group: nested ANOVA?

text: I recently ran an experiment assessing the impact of three system factors (factor A; factor B; factor C) on user attitudes (4-item scale). Each system factor had two levels (high vs. low). As a result, my experiment had nine conditions: eight treatment conditions (reflecting the 2x2x2 design) plus one control group. Subjects in the control group used a "baseline" system in which none of the three system factors were present. To analyze my data, I am considering to run two ANOVAs: (1) a one-way ANOVA comparing the control group with the overall treatment average and (2) a 2x2x2 ANOVA comparing the three factors and their interactions. However, this would not allow me to simultaneously compare the three factors with the control group. How could this be done? A nested ANOVA? Any suggestions are greatly appreciated.


title: How to cluster LDA/LSI topics generated by gensim?

text: I'm an enthusiastic single developer working on a small start-up idea. I reduced a corpus of mine to an LSA/LDA vector space using `gensim`. Now I have a bunch of topics hanging around and I am not sure how to cluster the corpus documents. I see that some people use k-means to cluster the topics. Can someone please elaborate? * * * So i had some to properly read up LDA/LSA and took a look at the gensim source. i did not realize that the Similarity Matrix was actually an MXM matrix where M is the number of documents in my corpus, i thought it was MXN where N is the number of topics. Using the matrix as an input to scikit's linkage function, i was able to create a heirachical cluster. Well, i was able to do it for a 1000 docs in my corpus (need more ram to handle a large corpus).


title: How to quantify statistical insignificance?

text: I'm relatively new to statistics, and understand that my question may be completely misworded. I am testing my own algorithm versus another. While the outputs are not identical, I want to show that the differences are "statistically insignificant." How can I quantify this, to make my point?


title: How to prove a particular identity involving the CDF of a continuous, strictly positive random variable?

text: Can someone show me why for a distribution function $F$ on $(0,\infty)$, $0 < y \leq x$, the following holds: $\frac{F(x)-F^{*2}(x)}{\overline{F}(x)}=\int\limits_{0}^{x}\frac{\overline{F}(x-t)}{\overline{F}(x)}dF(t)$ Edit: $F^{*2}$ stands for the convolution of $F$ with itself: $F^{*2}= F \star F = \int\limits_0^x F(x-y)dF(y)$ $\overline{F}=1-F$


title: How to perform recursive clustering in R?

text: Please let me know how to perform recursive clustering in R? How to give the input of the 1st iteration of cluster to the next one and so on. Any tutorial in this would be of great help.


title: What does "case-control" and "cross-sectional" mean in the context of logistic modeling?

text: While studying logistic modeling, I read the following statement > The fact that only odds ratios, not individual risks, can be estimated from > logistic   >  modeling in case-control or cross-sectional studies is not surprising. I do not know what do the "case-control" and "cross-sectional studies" stand for in the statistical analysis? Moreover, I do not quite understand what does the above statement mean from the viewpoint of statistical analysis. Any explanations will be appreciated.


title: Is "every blue t-shirted person" a systematic sample?

text: I'm teaching an intro stats class and was reviewing the types of sampling, including systematic sampling where you sample every kth individual or object. A student asked if sampling every person with a particular characteristic would accomplish the same thing. For example, would sampling every person with a blue t-shirt be random enough and provide enough of a representation of the whole population? At least, if you're asking a question other than "What color t-shirt do you prefer wearing?" My sense is no, but I wondered if anyone here had any thoughts on this.


title: Whether to delete cases that are flagged as outliers by statistical software when performing multiple regression?

text: I am performing multiple regression analyses and I am not sure whether outliers in my data should be deleted. The data I am concerned about appear as "circles" on the SPSS boxplots, however there are no asterisks (which makes me think they are not 'that bad'). The cases I am concerned about do appear under the "casewise diagnostics" table in the output - therefore should I delete these cases?


title: Using ksvm in R with differents constraints

text: I would like to use **ksvm** from the **kernlab** package in **R** to train and test a support vector machine. I'm able to use it but I would like to use different values of _C_ for each observations. To be more precise, I want to solve the following minimization problem (using standard SVM notation) $ \frac{1}{2}w\cdot w +C \sum_{i=1}^l \xi_i + C\_ \sum_{i=l+1}^n (\xi_i \cdot 1_{\{y_i=-1\}}) + C_{+} \sum_{i=l+1}^n (\xi_i \cdot 1_{\{y_i=1\}}) $ instead of the more standard $ \frac{1}{2}w\cdot w +C \sum_{i=1}^l \xi_i $ (I'm indeed working on a transductive support vector machine). Question: do you know if it is possible to solve this problem with ksvm (or another R svm package)? Thank you


title: How to deal with non-normality in MANOVA?

text: The dependent variables (DV) have to be normally distributed. I have a problem because some of them aren't. I have one independent variable (IV), namely _type of education_. The DV's are _Externalizing problems_ , _Internalizing problems_ , _Self-image_ , _Motivation_ , _Neuroticism_ , _Perseverance_ , _Social anxiety_ , _Visciousness_ and _Dominance_. The research question is: for which DV's do the levels of the IV differ. I assessed univariate normality with Kolmogorov-Smirnov, and the tests indicated that _Perseverance_ and _Social anxiety_ are not normally distributed. I also ran multivariate analysis based on Mahalanobis distance in order to find multivariate outliers. There are none. There are, however, several univariate outliers. I tried to transform the data but when I did it, I transformed both the levels of the IV so the Kolmorgorov-Smirnov test remained significant. My question is whether I now need to take the _Perseverance_ and _Social anxiety_ out of the MANOVA and run a non-parametric test on it. If that is so, which one might be the most appopriate? And should I do something like Bonferroni correction because I conduct several analysis now? Thank you very much in advance!


title: How to compute prediction bands for non-linear regression?

text: The help page for Prism gives the following explanation for how it computes the prediction bands for non-linear regression. Please excuse the long quote, but I am not following the second paragraph (that explains how $G|x$ is defined and $dY/dP$ is computed). Any help would be greatly appreciated. > The calculation of the confidence and prediction bands are fairly standard. > Read on for the details of how Prism computes prediction and confidence > bands of nonlinear regression. > > First, let's define G|x, which is the gradient of the parameters at a > particular value of X and using all the best-fit values of the parameters. > The result is a vector, with one element per parameter. For each parameter, > it is defined as dY/dP, where Y is the Y value of the curve given the > particular value of X and all the best-fit parameter values, and P is one of > the parameters.) > > G'|x is that gradient vector transposed, so it is a column rather than a row > of values. > > Cov is the covariance matrix (inversed Hessian from last iteration). It is a > square matrix with the number of rows and columns equal to the number of > parameters. Each item in the matrix is the covariance between two > parameters. > > Now compute c = G'|x * Cov * G|x. The result is a single number for any > value of X. > > The confidence and prediction bands are centered on the best fit curve, and > extend above and below the curve an equal amount. > > The confidence bands extend above and below the curve by: = > sqrt(c)*sqrt(SS/DF)*CriticalT(Confidence%, DF) > > The prediction bands extend a further distance above and below the curve, > equal to: = sqrt(c+1)*sqrt(SS/DF)*CriticalT(Confidence%, DF)


title: Ordinal variable with 0.3% of observations in one category - delete, ignore?

text: The response variable in my ordered logit model has 5 categories ranging from '1 = strongly disagree' to '5 = strongly agree'. However, only 0.3% of observations fall into category 1 (4 observations out of 1,320). This has created some problems, such as the inability to test the assumption of parallel lines (Brant test) in Stata. Also, treating these 4 observations as missing data changes the coefficients/significance levels of the model significantly, which worries me insofar as 4 observations should not have such a large effect on the model. Is it appropriate to run the model with the 4 observations set to missing? Or, should I leave the dependent variable unaltered? Any other suggestions are welcome. Thanks.


title: Matrix image regression or regression of subscales

text: A person is given 10 sample images (based on the same image) with corresponding techniques on how the image was achieved. I want to regress the matrix bitmap of the image of the 10 images on the image created by the person as a gross way of knowing which images the person copied from. I am assuming that if the person copied one particular person, they will have more corresponding pixels than other people. Is there such a thing as matrix regression? Can anyone point me to other possible techniques? I also have a rating scale of the features present in each image. With a similar intent, I would like to regress the subscores of the rating scale of the sample images to the subscore of the particular person's image. I imagine the same regression technique could be used. Hoping someone out there can point me to a technique that would let me do this.


title: Autocorrelation and heteroskedasticity in panel data

text: In the research, both autocorrelation and heteroskedasticity are detected in panel data analysis. I can solve them separately in stata with command "xtregar" and "robust", respectly. However, I cannot find a way to solve both problems at the same time.   If possible, please show me how to repair autocorrelation and heteroskedasticity problem for panel data in one estimaion. It will be great by using Stata, but SPSS is also fine.


title: Is there a reason to prefer a specific measure of multicollinearity?

text: When working with many input variables, we are often concerned about multicollinearity. There are a number of measures of multicollinearity that are used to detect, think about, and / or communicate multicollinearity. Some common recommendations are:   1. The multiple $R^2_j$ for a particular variable    2. The tolerance, $1-R^2_j$, for a particular variable    3. The variance inflation factor, $\text{VIF}=\frac{1}{\text{tolerance}}$, for a particular variable    4. The condition number, $\sqrt{\text{max(eigenvalue(X'X))}/\text{min(eigenvalue(X'X))}}$, of the design matrix as a whole  (There are some other options discussed in the Wikipedia article, and here on SO in the context of R.) The fact that the first three are a perfect function of each other suggests that the only possible net advantage between them would be psychological. On the other hand, the first three allow you to examine variables individually, which might be an advantage, but I have heard that the condition number method is considered best.   * Is this true? Best for what?    * Is the condition number a perfect function of the $R^2_j$'s? (I would think it would be.)    * Do people find that one of them is easiest to explain? (I've never tried to explain these numbers outside of class, I just give a loose, qualitative description of multicollinearity.)


title: Manual computation of Gower's similarity coefficient

text: There is an example on the computation of Gower's similarity coefficient on the page, Gower's similarity coefficient I am trying to work out the similarity manually between patient 1 and 2, however my computation is giving me completely different results. Here is the formula I use, in the order of the variables and their types $$ \frac{\left(1*\left(1 - \left(\frac{150-120}{150-110}\right)\right)+ 0 + 1 + 1 + 0\right)}{5} $$ I still cannot get the 0.0625 as indicated in the example. What am I doing wrong?


title: question about price-elasticity and endogenity

text: In _Davidson & McKinnon - Estimation and Inference_ I read that in a competitive market which is always in equilibrium we observe: $Q^d_t = Q^s_t = Q_t$ where $Q^d_t$ is the quantitiy of demand, $Q^s_t$ the quantitiy supplied. If one would estimate the equation (1) by OLS in a model given by (1) $Q_t^d = \alpha P_t + Z^d_t\beta + u^d_t$ (2) $Q_t^s = \gamma P_t + Z^s_t\delta + u^s_t$ with $P_t$ the price in period t, $Z^d_t$ correspods to a bunch of exogenous variables which determine the quantity, we would get a problem caused by the endogenity of $P_t$. This can be seen if we rewrite the equations (1) and (2) in terms of the observable variables $Q_t$ and $P_t$. **My question is:** Is this true in general should accounted for if one does a regression with real market data?


title: Spherical K-means Clustering in R

text: I have a large data set that I would like to cluster using spherical K means algorithm. However, I am relatively new to this subject and R in general. Most of my knowledge is self taught and I am still in the beginning stages- I have read all about K means clustering in the past few days and I would like to apply it to my own project (184 rows of 4000+ columns containing measurements in decimal values). How and where do I start? I am trying to teach myself how to go about this in R but I can't seem to find examples online on how to do this, or at least they don't apply to my case. I was wondering if anyone here would know about a step by step tutorial or has access to a script in which a spherical k means clustering has been conducted in R, and if you could share it with me. I have found a few papers but they are very advanced, and most are for text clustering and not numeric values. I hope this question is not too vague. Does anyone have experience in this subject and could guide me on which steps I have to take to get started? Absolute beginner here, so I apologize if I am completely wrong to ask such a question on this platform. Thank you for taking time out of your day to read this!


title: Stata error: "panels are not nested within clusters"

text: I have an unbalanced panel data set with more than 400,000 observations over 20 years. My panel variable is a person id and my time series variable is the year. The persons are from all over Germany which means that they are from different regions. To account for possible correlations between the persons within the same regions, I would like use clustered standard errors in my fixed effects regression. The Stata code looks like this:               xtreg depvar var1 var2 ... varN if year > 1990, fe vce(cluster clustvar)      First, I used the variable "region" which is a numeric variable indicating in which region the observation (person) lives as the cluster variable (clustvar). Of course, I got the error "panels are not nested within clusters" because I forgot that people move across regions such that they appear in more than a single region. To fix this problem, I wanted to use a region-year variable which is a numeric variable indicating the year as well as the region. To visualize the variable, the follwing extract from my dataset including the variable "region-year" may help:               region year region-year     1      1991 2000     1      1992 3000     2      1991 2500     2      1992 3500     2      1993 4500      Note that this extract does not contain the person id. The values of "region- year" are chosen arbitrarily in this example. In my opinion, it is just important that all observations belonging to the same region in the same year have the same number no matter how this number looks like. Am I right? Running the regression               xtreg depvar var1 var2 ... varN if year > 1990, fe vce(cluster region-year)      also causes the error "panels are not nested within clusters". But I don't know why.


title: Validity of maximising log-likelihood for maximum likelihood estimation

text: For reasons owing to mathematical convenience, when finding MLEs (maximum likelihood estimates), it is often the log-likelihood function---as opposed to the standard likelihood function---which is maximised. From what I've gathered, this approach is deemed valid as a result of the monotonically increasing nature of the (natural) logarithm function. My understanding of a monotonically increasing function is: for all $x$ and $y$ (defined on a subset of the reals), if $x \leq y$ then $f(x) \leq f(y)$. This, however, does not appear to be the case for all log-likelihood functions; for example: for the log-likelihood $\text{Gamma}(3, 5)$ function, if $x = 0.15$ and $y = 0.46$, then $f(x) = -0.59$ and $f(y) = -1.19$. Clearly, I've misunderstood this concept. Fundamentally, I guess I'm asking can somebody (preferably mathematically) demonstrate why: $$ \hat{\theta} = \text{max} \text{ } L(\theta) = \text{max} \text{ } \text{log} \text{ } L(\theta) $$ where $\hat{\theta}$ is the MLE for a given likelihood function.


title: Can variance in error term be estimated without fitting a line?

text: In this exercise, I was asked to prove that the simple linear regression - ( $Y_i=\beta_0+\beta_1X_i+\epsilon_i$ with all the usual basic conditions) - of $\{Y_{ji}\}$ for each $X_i$ is the same as the linear regression of ${\bar{Y}_i}$ for each $X_i$, where $\bar{Y}_i=\sum_j^nY_{ji}$ . This I managed to prove. Afterwards, they ask me to see if the variance of the error term, in this setting, can be estimated without fitting a line. To this last question, I think the answer is no, but I'm not sure. The error term variance unbiased estimator is $s^2=\frac{\sum(Y_i-\hat{Y}_i)^2}{n-2}$, where $\hat{Y}_i$ is the usual fitted value. Adapting to this setting the usual formula, we know that $s^2=\frac{\sum(\bar{Y}_i-\hat{Y}_i)^2}{n-2}=\frac{\sum_i\sum_j(Y_{ij}-\hat{Y}_i)^2}{n-2}$. I have no clue on how to tackle this exercise. Any help would be appreciated.


title: Covariance of binary and continuous variable

text: Suppose $y$ is a continuous random variable and $d$ is a binary random variable that takes the value $1$ with probability $p$ and $0$ with probability $1-p$. How do I show that $\text{Cov}(y,d)=(E[y|d=1]-E[y|d=0])p(1-p)$?


title: I have two sets of data (regular time intervals) is there any way to find out when they correlated the most and when they don't?

text: Sorry if the title is a bit vague however, i'm not sure exactly how to make my sentence concise. I have two times series: Amount invested into Iraq across time (in months) Price of a stock across time (averaged each month) There area few things i want to find out however, i'm only at an intro level stats and econometrics in Uni. Firstly how would i go about correlating these two time series i've read a few posts on here i'd like some clarification on the best approach to finding a correlation between them. so far i believe i have to use cross correlation which i hope to do using minitab. I want to do this to see if there is an over correlation between the two sets of data. Secondly is there anyway i can find out when the data correlates the most and when it does not i have over 150 data points and i really like to see during which months the correlation grew or was at its strongest. I want to do this so i can analyse why correlation might have decreased and so on (perhaps attributing it to other macro factors). I hope this makes sense and if you need any more clarification i'll answer straight away. Thanks for any points and i look forward to learning a little more today.


title: Pivotal quantities, test statistics and hypothesis tests

text: We are learning pivot functions, test statistics, and hypothesis testing at university but it makes no sense. I've tried reading my text book/notes, going through examples, etc., but the concepts seem like a random guess and I'm clueless about how to even start guessing what the answer could be. ### 1st part Can you please explain how to calculate a pivot function? E.g $X_{1},\dots,X_{n} \sim N(\mu, \sigma^2)$. Pivot function for $\sigma^2$ when $\mu$ is known and when $\mu$ is unknown. Why does $\mu$ being un/known matter? Also how would you calculate the pivot function for the ratio of two variances ($\sigma_{x}^2$ and $\sigma_{y}^2$)? Is it an F distribution? Assume $\mu_x$ and $\mu_y$ are known and $X_1,\dots,X_n \sim N(\mu_x,\sigma_x^2)$ and $Y_1,\dots,Y_n \sim N(\mu_y, \sigma_y^2)$. ### 2nd part Can you please explain how to calculate a test statistic (I get how to show it's a test statistic but don't know how to form one from scratch). Last, I have a few questions about hypothesis testing. I don't really understand how to calculate the power of a test or even what it means, to be honest. There is a whole bunch of theory and many definitions out there but they are rather abstract so I don't get it at all... I don't understand the notation or how to calculate the size/power of a test (generic form - not just with numbers). Example: $X \sim N(\mu, \sigma^2)$. $H_0: \mu >= \mu_0$ and $H_1: \mu < \mu_0$. $\mu$ and $\sigma$ unknown. Calculate the power and size of this test. How do I even start? I'm so confused. :( I'm really stuck with all of this and I hope you can help me! :) if there's a better resource out there to help please do let me know. EDIT1: Thanks for your reply. I did ask my lecturer to clarify… but ended up even more confused. He agrees the notes are unclear but will not rectify them because everyone else seems to get them! :( I also went to my class teacher and read through the Statistical Inference chapter several times – I get the basics but still don’t really understand most of it. I have googled around – and read Wikipedia – but it’s just more and more theory with no step by step examples explaining what to do. Everything seems randomly chosen and guesswork and hence my massive confusion. One thing though – I can’t read this: $$T_{X}=\sum_{i=1}^{N}\Big(\frac{X_{i}-\mu_{X}}{\sigma{X}}\Big)^{2} \sim \chi^{2}(N)$$ Is it supposed to be in mathematical notation? How do I view it properly? Yes- you are correct – the pivot function is used to calculate the confidence interval. The thing is – once I have the pivot it’s quite straightforward to calculate the CI. But it’s the pivot that’s causing the problems. I still don’t get the following: Pivot function for sigma^2 when mu is known and when mu is unknown. Why does mu being un/known matter? How would you calculate the pivot function for the ratio of 2 variances (sigmax^2 and sigmay^2)? Is it an F distribution? Assume mu x and mu y are known and X1....Xn - N(mu x,sigmax^2) and Y1...Yn -N(mu y, sigmay^2). Is it 1/Fn-1, m-1 = Fm-1,n-1? And the hypothesis testing questions above please…. Can you shed some light on this please? EDIT2: I did ask my lecturer to clarify… but ended up even more confused. He agrees the notes are unclear but will not rectify them because everyone else seems to get them! :( I also went to my class teacher and read through the Statistical Inference chapter several times – I get the basics but still don’t really understand most of it. I have googled around – and read Wikipedia – but it’s just more and more theory with no step by step examples explaining what to do. Everything seems randomly chosen and guesswork and hence my massive confusion. One thing though – I can’t read this: $$T_{X}=\sum_{i=1}^{N}\Big(\frac{X_{i}-\mu_{X}}{\sigma{X}}\Big)^{2} \sim \chi^{2}(N)$$ Is it supposed to be in mathematical notation? How do I view it properly? Yes- you are correct – the pivot function is used to calculate the confidence interval. The thing is – once I have the pivot it’s quite straightforward to calculate the CI. But it’s the pivot that’s causing the problems. I still don’t get the following: Pivot function for sigma^2 when mu is known and when mu is unknown. Why does mu being un/known matter? How would you calculate the pivot function for the ratio of 2 variances (sigmax^2 and sigmay^2)? Is it an F distribution? Assume mu x and mu y are known and X1....Xn - N(mu x,sigmax^2) and Y1...Yn -N(mu y, sigmay^2). Is it 1/Fn-1, m-1 = Fm-1,n-1? And the hypothesis testing questions above please…. Can you shed some light on this please?


title: How to use an initial posterior for recursive / sequential updating in WinBUGS

text: I am using WinBUGS to estimate / update the parameters of a model. The model is: $$ \begin{aligned} D(T,B,a)&= B*(a_0+a_1T+a_2T^2+a_3T^3)+error(B,T,a) \\ error &= \mathcal N(0, B^{0.5}a_4(a_0+a_1T+a_2T^2+a_3T^3)) \end{aligned} $$ where $D$ is an observed response, $B$ and $T$ are observed input parameters and $a_i$s are non-observable model parameters. To estimate the model parameters with the first set of data, I used non- informative normal priors for $a_i$'s. My question is: **If I want to do sequential / recursive updating,** that is to use a set of observations $(D,B,T)$ to compute posterior distributions for the model parameters and then use these posteriors as priors to update model parameters given a new set of data $(D,B,T)$ (and so on), **can I directly use the posterior distributions that WinBUGS finds or should I instead care about the joint posterior distribution?**


title: Variance explained for averaged values

text: I am trying to work out how good my particular model is at explaining some observed data. The problem here is that the observed data takes the form of averaged (mean) values for each of my predictive scores. When performing a simple correlation, I get a really high R-squared value (and this is repeated for independent data sets), which I am assuming means that there is a consistent relationship between my predictive value and the observed data. However, if I want to estimate how much of the variance in observed data I explain with my predictive score, how would I do this with averaged observed values? My underlying concern is that although I have a good correlation in all cases, there may be many underlying factors that also drive my observed data that are simply being 'averaged out' within each category, and as such, my r-squared value is meaningless. Just as an example, if my underlying observed data was: x<-c(3,4,3,2,1,6,5,7,5,4,9,7,8,10,11) And my predictive score for these values is: y<-c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3) I get r-squared of 0.911. However, if I average my observed data (which is all I can measure in my case), I would get: x1<-c((13/5),(27/5),9) y1<-c(1,2,3) R-squared = 0.997 Thus telling me that my predictive score explains almost all of the variation in the observed data, when a better representation above (x vs y) tells me I am explaining 91% of the variation. Thanks.


title: How to correlate two time series with gaps and different time bases?

text: I asked **this question** over on StackOverflow, and was recommended to ask it here. * * * I have two time series of 3D accelerometer data that have different time bases (clocks started at different times, with some very slight creep during the sampling time), as well as containing many gaps of different size (due to delays associated with writing to separate flash devices). The accelerometers I'm using are the inexpensive GCDC X250-2. I'm running the accelerometers at their highest gain, so the data has a significant noise floor. The time series each have about 2 million data points (over an hour at 512 samples/sec), and contain about 500 events of interest, where a typical event spans 100-150 samples (200-300 ms each). Many of these events are affected by data outages during flash writes. So, the data isn't pristine, and isn't even very pretty. But my eyeball inspection shows it clearly contains the information I'm interested in. (I can post plots, if needed.) The accelerometers are in similar environments but are only moderately coupled, meaning that I can tell by eye which events match from each accelerometer, but I have been unsuccessful so far doing so in software. Due to physical limitations, the devices are also mounted in different orientations, where the axes don't match, but they are as close to orthogonal as I could make them. So, for example, for 3-axis accelerometers A & B, +Ax maps to -By (up-down), +Az maps to -Bx (left-right), and +Ay maps to -Bz (front-back). My initial goal is to correlate shock events on the vertical axis, though I would eventually like to a) automatically discover the axis mapping, b) correlate activity on the mapped aces, and c) extract behavior differences between the two accelerometers (such as twisting or flexing). The nature of the time series data makes Python's numpy.correlate() unusable. I've also looked at R's Zoo package, but have made no headway with it. I've looked to different fields of signal analysis for help, but I've made no progress. Anyone have any clues for what I can do, or approaches I should research? Update 28 Feb 2011: Added some plots **here** showing examples of the data.


title: High dimensional model estimation with outliers

text: I have a set **H** of `k` `m`-dimensional hyperplanes in `n` dimensional space, where `1<<m<<n<<k`. I want to find an `m`-dimensional hyperplane **p** that's closest to **H**. However, a significant portion of **H** are considered outliers. Any inlier's distance to the true **p** is bounded, while the outliers' aren't predictable in advance. What's the canonical way of approaching this problem? RANSAC and other EM- based model estimation methods immediately come to mind, but sampling in hope of finding an all-inlier subset of **H** becomes exponentially more inefficient with the growth of dimension. Pointing me to an existing research paper is all I need.


title: Can we use log-likelihood to cluster classes?

text: I have an SVM classifier for m classes and n data points (somewhat evenly distributed across each class). Could I use the resulting MxN log likelihood matrix to merge classes that are similar?


title: When is MCMC required?

text: I am trying to understand when MCMC algorithms are utilized. I can do density estimation with MLE for non-linear cases just like linear cases, and EM for hidden variables and missing values. What are the scenarios when these tools are not sufficient and we use MCMC? Thanks


title: Entire data considered as support vector

text: I am currently learning to use support vector machine as classification. I have a data set with 161 observation and 18 dimension. I get 160 support vectors using svm function form R package, e1071. The obtained accuracy is 1. Can someone tell whether the result is accurate or not? actually I've found svm example from internet and I followed their strategy which is using svm on untrained data set (entire data) and second on trained data set(splitting data to train and test set). For the former and latter, model accuracy is 1 and 0.51. What should I do to increase the model accuracy?


title: Standard deviation in regression trees

text: In a regression tree, it is often assumed that each leaf is a Gaussian distribution $\mathcal{N}(\mu_i, \sigma)$, where $i$ is the index the leaf. Is $\sigma$ calculated as the standard deviation _within a leaf_ or the standard deviation for the dataset? If a tree is grown such that each leaf contains one instance (as is the case for bagged trees) then the _within leaf_ calculation seems ill-posed. However, using the whole dataset seems like it could induce a lot of bias. Is there something I am missing here?


title: calculating probability or filtering that certain subject is not in the particular cluster

text: I have a situation where there are n individuals and p features (variables). I do have their cluster information. Here is an example:               myd <- data.frame (      sub1 = c(1, "AB", "AB", "BB", "BB", "AA", "BB", "AB", "AA", "BB", "AB", "AA", "AB"),     sub2 = c(1, "AB", "BB", "BB", "BB", "AA", "BB", "AB", "AA", "BB", "AB", "AA", "AB"),     sub3 = c(1, "AB", "BB", "AB", "BB", "AA", "BB", "AB", "AA", "BB", "AB", "AA", "AB"),     sub4 = c(1, "AB", "BB", "AB", "BB", "AA", "BB", "AB", "AA", "BB", "AB", "AA", "AB"),     sub5 = c(1, "AB", "BB", "BB", "AB", "AB", "BB", "AB", "AA", "BB", "AB", "AA", "AB"),     sub6 = c(1, "AB", "BB", "BB", "BB", "AA", "AB", "AB", "AA", "BB", "AB", "AA", "AB"),     sub7 = c(2, "AA", "AA", "BB", "AB", "AA", "AB", "AB", "AB", "AB", "AB", "AB", "BB"),     sub8 = c(2, "AB", "AA", "AB", "AB", "AA", "AB", "AB", "AB", "AB", "AB", "AB", "BB"),     sub9 = c(2, "AA", "AA", "BB", "AB", "AA", "AB", "AB", "AB", "AB", "AB", "AB", "BB"),     sub10 = c(2, "AB", "AA", "BB", "AB", "AA", "BB", "AB", "AB", "AB", "AB", "AB", "BB"),     sub11 = c(2, "AA", "AA", "BB", "AB", "AB", "AB", "AB", "AB", "AB", "AB", "AB", "AB"),     sub12 = c(2, "AA", "AB", "BB", "BB", "AA", "AB", "AB", "AB", "AB", "AB", "AB", "BB"),     sub13 = c(3, "AB", "AB", "AB", "AB", "AB", "BB", "AB", "AA", "AB", "AA", "AB", "BB"),     sub14 = c(3, "AB", "BB", "BB", "AB", "AB", "BB", "AB", "AA", "AB", "AA", "AB", "BB"),     sub15 = c(3, "AA", "AB", "BB", "BB", "AA", "AB", "AB", "AB", "AB", "AB", "AB", "BB"))          rownames(myd) <- c("cluster", paste("var", 1:12, sep=""))      Thus this data has 15 subjects `sub1 : sub15` classified by categorical variables `var1: var12`. The cluster id is in the first row - with row name - `cluster`.               myd              sub1 sub2 sub3 sub4 sub5 sub6 sub7 sub8 sub9 sub10 sub11 sub12 sub13 sub14 sub15     cluster    1    1    1    1    1    1    2    2    2     2     2     2     3     3     3     var1      AB   AB   AB   AB   AB   AB   AA   AB   AA    AB    AA    AA    AB    AB    AA     var2      AB   BB   BB   BB   BB   BB   AA   AA   AA    AA    AA    AB    AB    BB    AB     var3      BB   BB   AB   AB   BB   BB   BB   AB   BB    BB    BB    BB    AB    BB    BB     var4      BB   BB   BB   BB   AB   BB   AB   AB   AB    AB    AB    BB    AB    AB    BB     var5      AA   AA   AA   AA   AB   AA   AA   AA   AA    AA    AB    AA    AB    AB    AA     var6      BB   BB   BB   BB   BB   AB   AB   AB   AB    BB    AB    AB    BB    BB    AB     var7      AB   AB   AB   AB   AB   AB   AB   AB   AB    AB    AB    AB    AB    AB    AB     var8      AA   AA   AA   AA   AA   AA   AB   AB   AB    AB    AB    AB    AA    AA    AB     var9      BB   BB   BB   BB   BB   BB   AB   AB   AB    AB    AB    AB    AB    AB    AB     var10     AB   AB   AB   AB   AB   AB   AB   AB   AB    AB    AB    AB    AA    AA    AB     var11     AA   AA   AA   AA   AA   AA   AB   AB   AB    AB    AB    AB    AB    AB    AB     var12     AB   AB   AB   AB   AB   AB   BB   BB   BB    BB    AB    BB    BB    BB    BB      I would like to see if the **cluster is proper** by calculating some of the measure (such as p-value) that can provide **the probability that particular subject belong to certain cluster**. _How can I achieve this_ ? _I do have a large number of variables / subjects, so I would need quicker solution as well_.


title: How to calculate the robust standard error of predicted y from a linear regression model in R?

text: How can I calculate the robust standard error of the predicted _y_ from a linear regression model in R? Any suggestion is appreciated.


title: Binary event probability optimization

text: I have a relatively small sample of binary events (50-100 events) that occurred during a time of day (the success rate is closely related to the time of day). I'm grouping these events into hour intervals, for example:   * ..   * [08:00-09:00) success: 5; failure: 7; 42% success   * [09:00-10:00) success: 1; failure: 0; 100% success   * [10:00-11:00) success: 30; failure: 20; 60% success   * [12:00-13:00) success: 11; failure: 19; 37% success   * ... After that I need to select n time ranges in order to schedule an event within the selected time ranges and have something like a minimum success probability. My question is if there is a way to take in consideration that some time intervals have very small data available (like 2nd interval in the example) and weight somehow this characteristic in order to avoid selecting only these intervals assuming a very high probability of success. To give some context: when I schedule an event, the event is executed more or less at the beginning of the time interval and in case of failure it is rescheduled 2 or 3 more times randomly within the selected time interval, the goal is to get a successful result using the less possible attempts. Any advice would be greatly appreciated.


title: Converting SAS EWMA code to Python

text: I have the following SAS code that uses PROC FORECAST that I would like to replicated with the Python Pandas pandas.stats.moments.ewma               proc forecast data=ewma method=expo        interval=weekday weight=(0.05) /*i.e. lambda=0.95*/       nstart=20 lead=0       out=out_ewma outfull;       id date;       var sq_gspc cross_returns;     run;      My DataFrame in Pandas is set up to mimic the `ewma` dataset in SAS -- same column names and everything, including same starting values. However, I can't seem to get the same values calculated with Pandas, probably due to a lack of understanding of the SAS options and how to set up the same calculation with Pandas. For example, to add the `sq_gspc` and `cross_returns` columns to my Pandas DataFrame, I'm doing the following:               ewmadf['f_sq_gspc']       = pd.ewma(ewmadf["sq_gspc"],       span=20, freq="D")     ewmadf['f_cross_returns'] = pd.ewma(ewmadf["cross_returns"], span=20, freq="D")      How do I replicate the PROC FORECAST here with the Python parameters? Edit: SAS documentation provides the following: > For METHOD=EXPO, n beginning values of the series are used in forming the > exponentially smoothed values S1, S2, and S3, where n is the value of the > NSTART= option. The parameters are initialized by fitting a time trend > regression to the first n nonmissing values of the series.


title: Generate separate plots for each group of records in dataset

text: I have a large csv file with around 45,000 rows and about 6 columns. The columns are set-up like:               Module, File-path, Status, Version-Old, Version-New, Change-diff      _Change-diff_ is simply the result of doing: _New - Old_. There are around 30 different modules and I want to perform analysis and create different graphs for each module based on the change-diff. I want to know the number of "rows" for each 'module' - I guess I can just count the number of each module name for the Module column to work that out. I'm wondering if there is any program that I can use to make my job easier, I'm currently trying to plug it all into Excel but I'm not sure how I can create graphs for each module instead of one big graph for all modules. Would R, SPSS or Matlab be a better option? Or anything else? EDIT: Just expanding on what the data looks like.               Module,  Path,      Status,  Old,  New,  Change     modA,    modA/a/a,  1,       11,   19,   8     modA,    modA/a/b,  1,       2,    4,    2     modA,    modA/a/c,  1,       0,    1,    1     modA,    modA/b/c,  0,       0,    1,    1     modA,    modA/b/a,  1,       7,    12,   5     modB,    modB/a/a,  1,       8,    9,    1     modB,    modB/a/b,  0,       0,    1,    1     modB,    modB/a/c,  1,       4,    10,   6     modC,    modC/a/a,  0,       0,    4,    4     modC,    modC/a/b,  1,       0,    3,    3      So I want to analyse all of modA based on Change (and possibly based on Old vs New etc.) Then analyse all of modB, modC etc. all in the same way but in separate graphs, automatically because of the amount of data I'm dealing with. Thanks in advance


title: How to add symbols on smooth lines in a R's lattice xyplot?

text: In a lattice xyplot with multiple smooth curves, I plotted curves categorized by color (y1, y2, y3 in blue; z1, z2, z3 in red). If we consider the following simple data :               x<-1:50     y1<-rnorm(50,100,20)     y2<-rnorm(50,80,20)     y3<-rnorm(50,90,20)     z1<-rnorm(50,60,20)     z2<-rnorm(50,40,20)     z3<-rnorm(50,65,20)      we can make such a plot with smooth lines and specify colors with this code:               xyplot(y1+y2+y3+z1+z2+z3~x,type='smooth', span=0.4,     par.settings=list(superpose.line=list(col=c('blue','blue','blue','red','red','red'))))      and we get this result: ![multiple smooth curves](http://i.stack.imgur.com/fWU66.jpg) My question : **Do you know how to add symbol (markers) on each curves to identify it ?** For instance, y1 with a triangle marker in blue, y2 with a square marker in blue, etc. It is easy to do this when raw data is plotted, but I can't find how to do it with smooth curves. I tried to use `superpose.symbol` but it seems that it's works only when `type=c('p','smooth')`. But in this case, symbols are not on the curves but in the raw positions. NB. I tried to identify the curves with different line type (cf. `lty`) but but the result does not allow a good visualization with my data.


title: Comparing arrival rate distributions

text: Given two samples of arrival rate of a process (they should end up something like a Poisson distribution):   * One over a long period (something like several days)   * Another from a shorter period time (maybe 30 minutes) How can I compare the shorter distribution to see if it varies significantly from the longer distribution? Ideally using an existing function in R? I believe a Q-Q Plot is showing this visually, but I am looking for a single number that would represent how far the shorter sample is from the longer sample. Also, perhaps better suited for another question. I am a bit curious about how I might say what amount of time would be a representative sample for the longer period of time. For reference the longer period plot looks something like the following, but it might very somewhat as I decide which attributes I might want to filter out: ![enter image description here](http://i.stack.imgur.com/ethnm.png)


title: What are the uses and pitfalls of regression through the origin?

text: Spuriously high R-squared is one of the pitfalls of regression through the origin (i.e. zero-intercept models). If the predictors do not contain zeroes, then is it an extrapolation? What are the uses and other pitfalls of regression through the origin? Are there any peer-reviewed articles?


title: How to implement k-means cluster analysis algorithm correctly?

text: I am trying to implement the K-mean analysis with the Standard algorithm. My implementation seems to work, but I noticed some strange behavior. If the _k_ is close to half of the length of the list to be analyzed, I will get a set that is empty. I am not sure if it is the correct behavior. I think the worst case is _k_ equal to the length of the list, and each result sets has only 1 element. Empty result sets will happen if _k_ is greater than the length of the list, but it is an invalid situation.


title: Is it valid to apply large scale statistics to a single case?

text: I have a small black box that displays a 0 or a 1 every time I press its button. I can tell you, truthfully, that over the past 20,000,000,000 times it has been pressed, 19,999,999,990 times it has shown a 1. I ask you to make a prediction of what number will appear when I press the button next. For the sake of simplicity, I won't attach any form of consequences to the results. You make a guess, and you are either right or wrong. How do you make your choice, and in what way do you find that choice to be valid? * * * Here's some other constraints that you can feel free to apply or discard if you think the result would be interesting:               1. The temporal distribution of the recorded results is unknown.     2. The display is blank until the button is pressed, and the digit will appear for a short amount of time afterwards. It is blank when you see it.      3. The button being pressed causes the display of the digit.


title: Two-sample t-test / ANOVA on functions, with unequal variances

text: Suppose $N$ experiments can be made in varying conditions. Each of them yields an estimate $f_i$ of a continuous (and, if necessary, positive) function of x over some interval. Experiment $i$ is repeated a number of times $n_i$, and each time it yields a function $f_i^j$. I am given the following data $$\forall x \quad \forall i \quad {n_i,f_i(x),\sigma_i(x)}$$ where $$f_i(x) = \frac{1}{n_i}\sum_{j=1}^{n_i} f_i^j(x) $$ $$\sigma_i^2(x) = \frac{1}{n_i}\sum_{j=1}^{n_i} (f_i^j(x)-f_i(x))^2 $$ Suppose experiment $1$ is the reference experiment. I would like to know whether some or all of the $f_i$ agree with function $f_1$. If they do not agree, I would like to know the region in which there is disagreement. Note however that since the experiments were made in different conditions, it is expected that the variance functions do not agree, and have to be treated as different for each experiment. For now, I was planning on doing a two-way ANOVA on each $x$ value followed by a two-sample $t$-test if necessary to determine the disagreeing experiments. However, I was looking for a method that treats the function as a whole, heteroscedastic, and that is bayesian if possible. I could not find anything. Literature references appreciated.


title: Finding weights that remove residual dot product

text: I have three time-series $m_i$, $k_i$, $s_i$, for $i \in [1,N]$. The problem is to find constants $A_m$, $A_k$, $A_s$, such that for $p_i = ( A_m \cdot m_i ) + ( A_k \cdot k_i ) + ( A_s \cdot s_i )$ , and a given integer constant $t$ \[ \sum_{i=1}^{N-t} { (m_i - p_i) (p_{t+i} - p_i) } = 0 \] \[\sum_{i=1}^{N-t} { (k_i - p_i) (p_{t+i} - p_i) } = 0 \] \[\sum_{i=1}^{N-t} { (s_i - p_i) (p_{t+i} - p_i) } = 0 \] Not sure if such constants would always exist. Any pointers ?


title: ICC in a linear mixed model - longnitudial data - why is this the correlation between time periods?

text: I'm confused by the description of the intraclass correlation (ICC) for a linear mixed model with longitudinal data from this material. A screen capture is shown below. What I'm confused about is that the ICC describes the amount of variation in $Y$ (the dependent variable) BETWEEN individuals, in this context (longitudinal data) but then also is the correlation between observations (which I think of as the correlation between an observation for a given individual at time $t$ and that same individual at time $t+k$). This correlation matrix seems to describe correlation within an individual. Why does this make sense? ![enter image description here](http://i.stack.imgur.com/tEdTz.jpg)


title: Using AIC to select between models that use nested and non-nested variables

text: I'm using SPSS to try and find a mixed model that adequate explains the data that I have. Two of the explanatory variables are closely related ('Sample group' and 'individual'), as an individual is only ever part of one sample group, so I've been nesting them if they are in the same model. I've been using the models AIC score to rank the models in order of explanatory power. Some of the models use the nested variables, and some of the models only use either 'Sample group' or 'individual'. My question is: Is it valid to use the AIC to compare between models that use nested variables and those that don't? To clarify by nested variables, I mean that some of the potential variables used in a model are: 1) sample site(individual) 2) sample site 3) individual


title: time series decomposition/dtrending using splines

text: Is there a way/method/approach to decompose a time series data using regression splines:   1. Seasonal time series into trend+seasonal+random component ?   2. A non seasonal time series into trend+random component ? I'm familiar with STL, Census and classical decomposition in R. All these techniques require time series data with seasonal component. We cannot extract trend if the time series is non seasonal (i.e., Frequency = 1). I recently came across this interesting article which is data driven in the recent 2013 ISF. Any insights on methods like these that are data driven decomposition using regression splines and that can be readily programmed in software packages such R would be greatly helpful. Thanks so much **Detrending time series with cycle and seasonal components** _Tatyana Krivobokova and Francisco Rosales_ In this work we discuss a nonparametric and completely data-driven approach to the decomposition of time series into a trend (cycle), seasonal and random components. Two former are modeled with penalized splines, while the latter is assumed to follow an ARMA structure. Empirical Bayesian approach allows to estimate both smoothing parameters and the orders of the ARMA process simultaneously resulting in an efficient, fast and data-driven decomposition procedure. The practical relevance of the approach is illustrated by real-data examples. The work is the extension of Kauermann, G., Krivobokova, T., Semmler, W. (2011) Filtering time series with penalized splines. _Studies in Nonlinear Dynamics & Econometrics._


title: Variable reduction by means of ANOVA?

text: I have a typical problem with several variables and a large amount of data which are not important right now. The goal of the study is to relate variable $Y$ with variables $X_1,X_2,...,X_n$. I have performed an ANOVA so that I see which variables are able to explain the variability on $Y$. However, my question is, to what extent is this procedure (ANOVA) valid to select variables? My selection criteria is to choose those more significant. Is it correct?


title: Poor predictions from lm.ridge?

text: I don't understand why predictions from `lm.ridge()` are so far out, when using the "best" lambda, based upon GCV. Can anyone help me to obtain better predictions? Or, at least, does anyone have a good ridge example with a simple explanation of the results? Below is my R code for wine quality data (from http://archive.ics.uci.edu/ml/machine-learning-databases/wine- quality/winequality-red.csv):               library(MASS)     wine_all <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";", header = T)     #wine_all <- read.table("winequality-red.csv", sep=";", header = T)     wine_train <- wine_all[1:1400,]     wine_test <- wine_all[-(1:1400),]          train.lm <- lm.ridge(quality~., wine_train, lambda = seq(0, 100, 0.1))     plot(x=train.lm$lambda, y=train.lm$GCV)          pred.test <- scale(wine_test[,1:11], center = F, scale = train.lm$scales) %*% train.lm$coef[, which.min(train.lm$GCV)] + train.lm$ym          pred.all <-  scale(wine_all[,1:11], center = F, scale = train.lm$scales) %*% train.lm$coef[, which.min(train.lm$GCV)] + train.lm$ym     cor(wine_test[, 12], pred.test)^2     cor(wine_all[, 12], pred.all)^2


title: Is there an upper limit on the number of intervals in a histogram?

text: I've read several articles and excerpts from books that explain how to choose a _good_ number of intervals (bins) for histogram of a data set, but I'm wondering if there's a hard _maximum_ number of intervals based on the number of points in a data set, or some other criterion. **Background:** The reason I'm asking is that I'm trying to write software based on a procedure from a research paper. One step in the procedure is to create several histograms from a data set, then choose the optimal resolution based on a characteristic function (defined by the authors of the paper). My problem is that the authors don't mention an upper bound for the number of intervals to test. (I have hundreds of datasets to analyze, and each one can have a different "optimal" number of bins. Also, it's important that the _optimal_ number of bins is selected, so manually looking at the results and picking a _good_ one won't work.) Would simply setting the maximum number of intervals to be the number of points in the data set be a good guideline, or is there some other criterion that's typically used in statistics?


title: How to communicate "very negative" parameter estimate in paper?

text: I'm currently writing a thesis and I have to interpret a coefficient estimate that's extremely negative compared to what I hypothesized. How do I best communicate this? Describing it as `"very negative"` sounds a bit strange, but I can't think of anything better.


title: data normalisation problems

text: I am pretty new in machine learning and hence facing a lot of confusion in data normalisation concepts. Someone pls clarify the following doubts : 1) While normalising a data matrix of m-samples x n-features, the normalisation should be done along m-samples or along the n-features i.e we normalise each row or each column here ? 2) We should normalise all the data together or normalisation should be done separately for each class. 3) Suppose I have a vector of 2-samples x n-features as shown               class_1   0.4432       27.19        0.6733   0.0828  0.4134  -0.6662  0.3381  0.0552   0.0    -0.01167    0.0    0.0    0.0    0.00056    0.0    0.0    0.03111    0.05778    0.01333    0.0    -0.02722    -0.02333    0.0    0.00222    -0.02722    0.0    0.0    0.0    10100   10010   11000   10100   10100   11000   11000   10001   11000   10001   10001   10100   10001   10001   10001           class_2   0.7647       16.1073        0.7867   -0.2414  -0.2264  -0.8025  0.8649  -0.3524   -0.01944    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.04611    0.0    -0.01944    0.03556    -0.02722    0.03111    0.0    0.0    0.0    0.0    0.0    10010   10100   11000   10010   10001   10100   10010   10001   10001   11000   10001   10001   10001   11000   10001       How to normalise in this case as you can see that the first few features have float values but there are some which are encoded binary features ( the 00110 types are some position informations) . If we normalise what will happend to such information ? 4) The training and testing data should be separated after normalisation or before normalisation and then normalised ? 5) What about the user sample ? once I train and prepare my classifier the user give a SINGLE input of 1-sample x n-features ? It is quite obvious that the values for each feature will differ in magnitude greatly as compared to normalised data used in training and testing. How to handle this thing ?


title: Computing a circular-linear partial correlation

text: The `CircStats` toolbox for MATLAB (http://bit.ly/18C1SCF) implements a procedure to compute a correlation between a linear and a circular variable. Specifically, the correlation between a linear variable $x$ and a circular variable $\alpha$ is given by $$\rho_{\textrm{cl}} = \sqrt{\frac{r^2_{\textrm{cx}}+r^2_{\textrm{sx}}- 2r_{\textrm{cx}}r_{\textrm{sx}}r_{\textrm{cs}}}{1-r^2_{\textrm{cs}}}},$$ where $r_\textrm{sx}$ is the Pearson correlation coefficient between $\sin\alpha$ and $x$, $r_\textrm{cx}$ is the coefficient between $\cos\alpha$ and $x$, and $r_\textrm{cs}$ is the coefficient between $\cos\alpha$ and $\sin\alpha$. How do I compute the partial correlation of $x$ and $\alpha$, given a third variable, $y$ (linear), as a control variable?


title: wavelet missing data

text: I have multiple time series datasets with missing data (the data is missing due to it not being recorded on these dates), and wish to perform a wavelet transformation on the data (probably a MODWT, due to its properties defined by Percival and Walden [2000]). Given the wavelet details occur at different frequencies - and there is also the wavelet smooth to avail of - should I address the missing data before the wavelet filtering stage (say, interpolation etc.)? Or, due to the properties of the wavelet filtering at different frequencies, will the bias incurred by missing data be reduced at the different filtering frequencies? Thanks


title: Sufficient Statistic for Uniform

text: Let $X$ have the pdf $f_X (x;\theta)=1/(2\theta)$ for $-\theta<x<\theta$, zero elsewhere where $\theta>0$ Is the statistic $Y=|X|$ a sufficient statistic for $\theta$? * * * My book claims it is but I do not see it. A transformation argument shows that the distribution of $|X|$ is uniform with $f_X (x;\theta)=1/\theta$ for $0<x<\theta$. Thus the conditional probability $$\frac{f(x_1;\theta)f(x_2;\theta)\ldots f(x_n;\theta)}{f_{|X|}(|X|;\theta)}$$ _does_ depend on $\theta$ and therefore, $|X|$ cannot be a sufficient statistic for $\theta$. Is there something that I am doing wrong here? Thank you.


title: Determine minimum acceptable step in cut()

text: I'm using `cut()` to transform data like this:               N <- 10     table(cut(iris$Sepal.Length,quantile(iris$Sepal.Length,probs=seq(0,1,1/N))))      However, I do have many variables and I would like to determine what is the maximum acceptable `N` for variable. I thought this will be linked to unique values, but it is not that easy:               length(unique(iris$Sepal.Length)) # 35 unique values         N <- 19         ## Fails if N > 19         table(cut(iris$Sepal.Length,quantile(iris$Sepal.Length,probs=seq(0,1,1/N))))      Is there a way how to calculate `N` or do I have to find this number iteratively? Maybe this will be simpler case to work with:               N <- 4     x <- c(1,2,2,2,2,2,2,2,3,3,3,3,4,5,6,7,8)     table(cut(x,quantile(x,probs=seq(0,1,1/N))))


title: How to create a barplot diagram where bars are side-by-side in R

text: I want to create a bardiagram for these data in R (read from a CVS file): Experiment_Name MetricA MetricB Just_X 2 10 Just_X_and_Y 3 20 to have the following diagram: ![alt text](http://i.stack.imgur.com/kGwFA.png) I am beginner and I do not know even how to start.


title: R packages that work with biased samples

text: I'm working with a biased sample of web users. I'm only able to track responses of users who have navigated my site in a certain way, and I'd like to run an analysis to determine how certain factors (which products they are into, etc) influence how much they buy from my online store. I've searched around online and found some information on the Heckman Correction, which somehow both won a Nobel Prize and also has extremely little information online, including no R support or youtube videos that I can find. What are some ways you work with building models that try to correct for biased data? Are there easy to use software libraries in R? I have a ton of unlabeled data (products that every user has seen), so it's possible that I could create a naive model for determining probability that a user ends up in my sample if that helps, and I'd be willing to work with more advanced learners like SVM (which I've read can work better for these types of problems) in case linear regressions are a bad choice.


title: Accounting for minimum dependent measure in data when fitting a distribution

text: I have what is possible a naive question. I am current comparing various models (i.e. distributions). And the comparisons do not involve different distributions but rather how the model is fed the data. For example, two models might be 1) a mixture of a gamma + exponential and 2) a mixture of a gamma + exponential wherein the minimum data point is subtracted. With respect to model 2). If I subtract the minimum data then I am left with a data point = 0. The model's fit fine and I get my MLE values. However, If I am working with just a gamma and in one model I remove the minimum data point, then that data point = 0 causes errors in the MLE estimator. To get past this issue I simply remove the 0 data point and fit the distribution. My question is, is there a problem with simply removing the data point (my worry is decreasing the size of my data set)? And should I do that for the gamma + exponential as well? Thanks


title: How to integrate observational errors in goodness of fit tests?

text: I have an astrophysical non-linear curve, specifically a power spectrum. I need to fit this curve with a model and obtain the goodness-of-fit (GOF). This gives me expected and observed values. The data also have observational errors related to instrumental uncertainties. Is there any way I can integrate observational errors with a chi-square test (or any other GOF test)?


title: Cross validation and parameter optimization

text: I have a question about the parameter optimization when I use the 10-fold cross validation. I want to ask that whether the parameters should fix or not during every fold's model training , i.e. (1) select one set of optimized parameters for every fold's average accuracy. or (2) I should find the optimized parameter for every fold and then every fold uses different optimized parameters to train its model then test on the fold's test data respectively, and finally average every fold's accuracy as result? Which one is the correct method for cross validation? Thanks a lot.


title: Conditional probability problem

text: Could someone please help explain how to calculate $P(A)$ given: $ P(B) = 0.2 $ $P(A|B) = 0.6$ $P(A|{\rm not} \ B) = 0.6$ Any explaination/help would be greatly appreciated.


title: How should I analyze a data set with 1 IV and 2DVs (all scale variables)?

text: Hello wonderful statisticians, I'm writing a proposal for a dissertation where I'd like to look at whether heart rate can predict 1) narrative coherence (as judged on a scale - but could also be categorized), and 2) self-reported relationship satisfaction (also scalar). Ideally, I'd also like to see if DV1 (narrative coherence) predicts DV2 (relationship satisfaction), as I suspect they're all interrelated. What's the best way to analyse this data? Can I analyze DV1 as an IV without compromising the integrity of my analysis? My thought is to do a MANOVA or MANCOVA. Any advice, however, would be appreciated!


What are some of the best ranking algorithms with inputs as up and down votes?


title: Software enabling GLS estimation with both heteroskedasticity and serial correlation correction

text: Does anyone know such? I tried to find such procedures in Gretl, but there you can use either hsk procedure for heteroskedasticity correction or ar1 procedure for serial correlation correction. I need the GLS procedure that will deal with both. Thanks


title: Statistical models and methods to evaluate and forecast "fine price" increases

text: I want to evaluate the implications of increasing fine prices. I will have a few different scenarios ranging from business as usual, minor increase, proportional increase, categorical increase, to extreme-increases. Each scenario will have different levels of monetary increase depending on the fine details and subprograms. Also, I'd like to do comparative analysis of implementing a strike base system (x times you get a fine, you're out) versus a cooperative non-extreme user lead system. Are there any models that predict such increases, especially considering that there is a likelihood of people to eitehr decrease their or increase there frequency of fines after a new system is in place.


title: What is the difference between empirical variance and variance?

text: As far as I know variance is calculated as $$\text{variance} = \frac{(x-\text{mean})^2}{n}$$ while $$\text{Empirical Variance} = \frac{(x-\text{mean})^2}{n(n-1)} $$ Is it correct? Or is there some other definition? Kindly explain with example or any refence for reading on this topic


title: Low accuracy in out of time validation

text: I prepared a model which had very good accuracy (80.5%) on my out of sample data. However when I ran that model on a population which is some 6mths old the accuracy went down to abysmal 33%. I am talking here about percentage detection of event (say defaulters). So right now my model is only detecting 33 out of 100 defaulters in the out of time dataset. Please suggest what can be the possible reasons behind this. How do I go about improving this? It'd be tough to defend this accuracy before the client. God forbid, if required to defend this accuracy before the client, how could this be justified?


title: How to compute df for ARIMA models?

text: Running this example from Hyndman's textbook in Chapter 9.1:               auto.arima(austa,d=0,xreg=1:length(austa))     ARIMA(2,0,0) with non-zero mean           Coefficients:              ar1      ar2  intercept  1:length(austa)           1.0371  -0.3379     0.4173           0.1715     s.e.  0.1675   0.1797     0.1866           0.0102          sigma^2 estimated as 0.02486:  log likelihood=12.7        I wish to test if the trend is significant. When using `lm`, I would need the degrees-of-freedom to get the p-values, e.g.:               tstats <- coef(test) / sqrt(diag(vcov(test)))     2 * pt(abs(tstats), df = df.residual(test), lower.tail = FALSE)      How do I compute d.o.f for this example (and in general for ARIMA models) to get the p-values?


title: Category selection for text classification

text: It is said that to achieve a good result (many different metrics) for text classification, it is not always a business of selecting the algorithm/classifier. Sometimes, it is even more important to find a good feature selection method given the specific task. However, even before feature selection, selection of the category set actually will determine the final result a lot. Let's say, in many cases, the categories are actually semi-flexible or totally flexible. I can imagine, to get a better set of categories, both the purpose of text classification and the statistical feature of the text play important roles. I want to know, is there any tutorial/discussion/suggestion on how to select the categories? And what methods can be applied to analyse the statistical features of the text collection? (I can think of using clustering.)


title: What are the assumptions for checking the stationarity of a time series?

text: I am checking stationarity or non-stationarity of a time series with R and I am using `adf.test` and `kpss.test` in `tseries` package. What are the assumptions for these tests?   Is following a Gaussian distribution by the data set one of the assumption?   If yes, what can we do for non-Gaussian time series? Also, is using the default options for lag order in ADF which is $k = \text{trunc}((\text{length}(x)-1)^{1/3})$ a safe assumption? Thank you


title: The multiple regression model - The minimum sum of squares

text: I am struggeling with this definition the minimum sum of squares for a Matrix $M$: The minimum value of $f(β)$ where $β = (X′X)^{−1}X′y$ from the system of normal equations $X′Xβ = X′y$ is $$y′(I − X(X′X)^{−1}X′)y = y′My$$, with $M = (I − X(X′X)^{−1}X′)$ a symmetric, idempotent (T × T)-Matrix. The trace of $M$ is $\operatorname{tr(M)} = \operatorname{tr}(I)−\operatorname{tr}(X(X′X)^{−1}X′) = T −K $and is equal to the rank $r(M)$. What is meant by the rank $M$(can this be the "real" rank of the Matrix $M$) and how to come to this conclusion?


title: categorical independent variable with three levels and binary logistic regression

text: I am not very experienced in statistics that's why you may find my question stupid. Anyway, I want to learn which level of a categorical independent variable should I look to interpret the odd ratios in binary logistic regression. For example, I have one independent categorical variable (education ) with three levels. The results of binary logistic regression show that "edu" (reference level) is not significant but one of the levels is significant "edu(2)". Is it possible to ignore the overall significance of "edu" and interpret the "edu(2)"?


title: Multivariate Multiple Classification

text: Is it possible to have multiple (at least 3) dependent variables in a single classification model. I know this can be accomplished in a regression model but I need to perform this with a classification model. My three DVs are correlated so I do not want to simply write three seperate classification models. Does anyone know of a way to do this in either R or python's scikit?


title: Conditional match and combine data in R

text: What would be the most efficient way to achieve the following task in R: combine the values of one variable for those observations which match on one index variable and have specific values on a second index variable. For example:               id1<-c("a","a","b","b")     id2<-c ("x","y","x","y")     variable<-rep(1,4)     d<-as.data.frame(cbind(id1,id2,variable))     d$variable<-as.numeric(as.character(d$variable))      The task is, for the observations (rows) that have the same 'id1', to add the value of 'variable' for the cases with 'id2' equal to 'x' to the value of 'variable' for the cases with 'id2' equal to 'y'. The best I could come up with is double looping which does the job but is very slow and ugly:               for (i in 1:nrow(d)){       if (d$id2[i]=='x') {             for (j in 1:nrow(d)){               if (d$id1[i]==d$id1[j] & d$id2[j]=='y'){                 d$variable[j]<-d$variable[i]+d$variable[j]           } else next  } } else next }       So what would be the most efficient approach?


title: Box's test for equal covariances or Mauchly's test of sphericity?

text: I am writing a report on a 2 x 3 mixed study. There is ONE within IV with TWO levels and there is ONE between subjects IV THREE levels. There are two DVs. Sample sizes for each group is 25, 25, 27. I tried Mauchly's test but it does not produce any results. I read somewhere that you can only use that for studies with 3 or more within IVs. So my question is do I use the Box's test and refer to the Greenhouse Geisser output?


title: How to determine sample size prior to designing a prospective cohort study research of 2 groups?

text: How do I use the a priori test to determine a sample size for my prospective cohort research that compares vaccinated males versus non vaccinated males. I am looking at the effect of vaccination on disease prevalence in the future.


title: Variance of product of k correlated random variables

text: What is the variance of the product of $k$ correlated random variables?


title: What is the best criterion for performance evaluation in classification algorithms?

text: Here are some results for ANN and KNN on abalone data set using Weka:               Result for ANN      Correctly Classified Instances  3183 76.203 %      Incorrectly Classified Instances 994 23.797 %      Mean absolute error      0.214      Root mean squared error  0.3349      Relative absolute error 58.6486 %          Result for KNN      Correctly Classified Instances  3211 76.8734 %      Incorrectly Classified Instances 966 23.1266 %      Mean absolute error      0.2142      Root mean squared error  0.3361      Relative absolute error 58.7113 %      KNN has high accuracy but ANN has low errors. So which of the two algorithms should I say is better? Which is the more preferable criterion, accuracy or error? What I understood was that error should decrease with high accuracy, but the results here are opposite. Why is this so?


title: What are analysts looking for when they plot a differenced, logged time series?

text: So in R, for instance, this would be:               my_ts_logged_diffed = diff(log(some_ts_object))     plot(my_ts_logged_diffed)      This seems to be part of every experienced analyst/forecaster analytical workflow--in particular, a _visual examination of the plotted data_. What are they looking for--i.e., what useful information does this transformation help reveal? Similarly, I have a pretty good selection of time series textbooks, tutorials, and the like; nearly all of them mention this analytical step, but none of them say why it's done (i am sure there's a good reason, and one that's apparently too obvious to even mention). (i do indeed routinely rely on this transformation but only for the limited purpose of testing for a normal distribution (i think the test is called _Shapiro-Wilk_ ). The application of the test just involves (assuming i am applying it correctly) comparing a couple of parameters (a 'W' parameter and the p-value) against a baseline--the Test doesn't appear to require plotting the data).


title: Find the standard deviation

text: **We must find the standard deviation out of the text below:** _The raw material needed for the manufacture of medicine has to be at least $97\%$ pure. A buyer analyzes the nullhypothesis, that the proportion is $\mu_0=97\%$, with the alternative hypothesis that the proportion is higher than $97\%$. He decides to buy the raw material if the nulhypothesis gets rejected with $\alpha = 0.05$. So if the calculated critical value is equal to $t_{\alpha} = 98 \%$, he'll only buy if he finds a proportion of $98\%$ or higher with his analysis. The risk that he buys a raw material with a proportion of $97\%$ (nullhypothesis is true) is $100 \times \alpha = 5 \%$_ This text confuses me to the extreme. From the first part it seems to me that the chance that he gets raw material with a proportion $\geq 98 \%$ is $0.05$ (as it is the critical value). From the final sentence it seems to me that the chance he gets of buying material with a proportion of $\leq$ $98 \%$ is $0.05$ These of course are contradicting and I hope someone can clear up my confusion! ps - A little side question. Let's say you want to reject a hypothesis ($\alpha = 0.05$) and you find that the chance of the $H_0$ being true miraculously is $0.05$. Do you reject the hypothesis or not? In other words, does the boundary count?


title: Statistical test for n x m contingency tables

text: I have a dataset made up of elements from three groups, let's call them G1, G2, and G3. I analysed certain characteristics of these elements and divided them into 3 types of "behaviour" T1, T2, and T3 (I used cluster analysis to do that). So, now I have a 3 x 3 contingency table like this with the counts of elements in the three groups divided by type:                     |    T1   |    T2   |    T3   |     ------+---------+---------+---------+---       G1  |   18    |   15    |   65    |      ------+---------+---------+---------+---       G2  |   20    |   10    |   70    |     ------+---------+---------+---------+---       G3  |   15    |   55    |   30    |      Now, I can run a Fisher test on these data in R               data <- matrix(c(18, 20, 15, 15, 10, 55, 65, 70, 30), nrow=3)     fisher.test(data)      and I get                  Fisher's Exact Test for Count Data          data:  data      p-value = 9.028e-13     alternative hypothesis: two.sided           So my questions are:   * is it correct to use Fisher test this way?   * how do I know who is different from who? Is there a post-hoc test I can use? Looking at the data I would say the 3rd group has a different behaviour from the first two, how do I show that statistically?    * someone pointed me to logit models: are they a viable option for this type of analysis?   * any other option to analyse this type of data? Thank you a lot nico


title: Latent variables in Bayes nets with no physical interpretation

text: In _Pattern Recognition and Machine Learning_ Bishop writes about Bayes networks: > For practical applications of probabilistic models, it will typically be the > highernumbered variables corresponding to terminal nodes of the graph that > represent the observations, with lower-numbered nodes corresponding to > latent variables. The primary role of the latent variables is to allow a > complicated distribution over the observed variables to be represented in > terms of a model constructed from simpler (typically exponential family) > conditional distributions. And after a few lines: > The hidden variables in a probabilistic model need not, however, have any > explicit physical interpretation but may be introduced simply to allow a > more complex joint distribution to be constructed from simpler components. What do you think he means by this type of hidden variables (with no physical interpretation)?   What can be a simple example of this?   I thought about mixture of gaussians, but they don’t correspond to a situation where the variables we are interested are highernumbered.


title: Comparing observed and predicted values across several measurements

text: As a neuropsychology graduate student with some experience in statistics (I'm usually the guy other psychologists come to with statistics problems _after_ trying it themselves but _before_ seeing a statistician), I am investigating whether a medical treatment has an effect on cognitive measures (aside from curing the medical problem). I have been looking for a statistically sound method to approach the following problem (preferably using R) which will also pass muster with (medical) reviewers. I will try to explain the problem, which appears to be somewhat unusual - perhaps because of an error in my thinking. I have data from two populations: a large group of controls, and a smaller group of patients (n ~ 20). In both groups, I have scores on a normally distributed cognitive scale at four time points (T0-T3), which in patients corresponds to a pre-treatment score (T0) and three post-treatment scores (T1-T3). As part of an accepted approach in neuropsychology (SRB-analysis; e.g. Duff, 2012), I have used multiple regression on the control data to model scores at T1, T2 and T3 from demographics and baseline (T0) scores. This regression model is then used to predict scores at T1-T3 in patients, and change in individual patients is evaluated based on the difference between predicted and observed score at T1-T3 (as well as the prediction error/variability in controls). For a single patient that looks something like this:                                  T0      T1      T2      T3     Observed Score     100     105     106     107     Predicted Score    -       110     117     120      However, the analysis is intended for an audience that is used to group (rather than individual) analysis, and to (partially) oblige this audience I would like to compare the predictions in patients **as a group** with the observed scores. This group analysis presents with me with two problems:   1. Due to phenomena such as regression to the mean, practice effects, etc., the distribution of the predicted scores is different from the distribution of the observed scores. For example, regression to the mean decreases the spread in the predicted scores (as compared to observed scores).   2. Although data exploration (histograms, qq-plots, etc) suggest that difference scores (observed-predicted) are normally distributed, the group size (n ~20) is smaller than what is usually recommended as a minimum group size for approximate normality (n ~30). I would like to answer two questions in the group analysis:   1. Are observed scores different from (/lower than) predicted scores?   2. At which time points (T1-T3) are observed scores lower than predicted scores? As a first approximation, I have used one-sample t-tests on the difference scores (observed-predicted) with Bonferroni correction for multiple (3) comparisons. But that is clearly a non-optimal solution: it does not directly provide the omnibus test required to answer the first question and it assumes normality where that is at least somewhat doubtful. Hence my question: what would be the optimal way to perform the group analysis?