title: Choosing number of PCA components when multiple samples for each data point are available

text: **Update** : I posted an answer below describing my current attempts to approach this problem. I am facing a perennial problem of identifying significant principal components. This question has been discussed here a couple of times (one, two and a very impressive thread three), and the usual advice is to use either so called "parallel analysis" (a very misleading name in my opinion; it is essentially Monte Carlo estimation of eigenvalue spectrum under null) or Laplace/BIC approaches from Minka, 2000, Automatic Choice of Dimensionality for PCA. Apart from that, I am aware of some ways to perform cross-validation over data points (see Bro et al., 2008, Cross-validation of component models), plus there is of course a good old approach of staring at the eigenvalue spectrum trying to locate an "elbow". All of these methods assume that you have one set of data points, let us say $N$ points from $\mathbb{R}^D$, and nothing apart from that. In my case, however, each point is an average of $M$ points, which are repeated measurements of the same value. This makes me think that I can somehow infer which PCA components are significant, by using the data that goes into averages; maybe by some cross-validation or bootstrapping approach. My question is: how? **Update**. Following discussion in the comments, I should clarify what I mean by "significant". For this I will introduce some notation. There are $N$ points $\bar{\mathbf{x}}_i \in \mathbb{R}^D$, and each point is an average over $M$ repeated measurements $\bar{\mathbf{x}}_i=\frac{1}{M}\sum_{j=1}^M \mathbf{x}_i^{(j)}$. Each measurement is a noisy observation of a true position of the corresponding point, i.e. $\mathbf{x}_i^{(j)}=\mathbf{a}_i + \mathbf{\xi}_i^{(j)}$ (note that $\xi$ is a $D$-dimensional vector as well, but I cannot make Greek letters bold here). I am interested in principal components of $\\{\mathbf{a}_i\\}$, but can only perform PCA on $\\{\bar{\mathbf{x}}_i\\}$, and the small components can be completely distorted by the noise. The larger the noise, the less principal components of $\\{\mathbf{a}_i\\}$ I will be able to recover. The challenge is to select only those leading PCs that _definitely_ "come from" $\\{\mathbf{a}_i\\}$. In my answer below I describe two approaches that I currently use. I am not fully happy with either.


title: Data transformations of scales

text: Could really do with some help from someone with more stats expertise than I have. I have some data relating to a customer survey along a range of measures which are consolidated into a number of scales. Participants rated the company and its competitor along all measures. I transformed the raw data by squaring it to reduce skew and kurtosis. I now want to analyze the differences between the measures by the company and its competitor using the delta. I used the delats from the transformed variables, but as I have skew in the delta scores do I transform the data again or should I have taken the difference on the raw, untransformed scores, and then transformed the deltas to reduce skew and kurtosis? Also, given that I’ve got negative and positive deltas, could you recommend an appropriate method to transform the data?   Thanks in advance for any help you can give.


title: Fake uniform random numbers: More evenly distributed than true uniform data

text: I'm looking for a way to generate random numbers that _appear_ to be uniform distributed -- and every test will show them to be uniform -- except that they are **more evenly distributed than true uniform data**. The problem I have with the "true" uniform randoms is that they will occasionally cluster. This effect is stronger at a low sample size. Roughly said: when I draw two Uniform randoms in U[0;1], chances are around 10% that they are within a range of 0.1, and 1% that they are within 0.01. So I'm looking for a good way to generate random numbers that are **more evenly distributed than uniform randoms**. Use case example: say I'm doing a computer game, and I want to place treasure randomly on a map (not caring about any other thing). I don't want the treasure to be all in one place, it should be all over the map. With uniform randoms, if I place, say, 10 objects, the chances are not that low that there are 5 or so really close to each other. This may give one player an advantage over another. Think of minesweeper, chances (albeit low, if there are enough mines) are that you are really lucky and win with a single click. A very naive approach for my problem is to divide the data into a grid. As long as the number is large enough (and has factors), one can enforce extra uniformness this way. So instead of drawing 12 random variables from U[0;1], I can draw 6 from U[0;.5] and 6 from U[0.5;1], or 4 from U[0;1/3] + 4 from U[1/3;2/3] + 4 from U[2/3; 1]. Is there any better way to get this extra evenness into the uniform? It probably only works for batch randoms (when drawing a single random, I obviously have to consider the whole range). In particular, I can shuffle the records again afterwards (so it's not the first four from the first third). How about doing it incrementally? So the first is on U[0;1], then two from each halves, one from each third, one from each fourth? Has this been investigated, and how good is it? I might have to be careful to use different generators for x and y to not get them correlated (the first xy would always be in the bottom half, the second in the left half and bottom third, the third in center third and top third... so at least some random bin permutation is also needed. and in the long run, it will be too even, I guess. As a side node, is there a well-known test whether some distribution is _too_ evenly distributed to be truly uniform? So testing "true uniform" vs. "someone messed with the data and distributed the items more evenly". If I recall correctly, Hopkins Statistic can measure this, but can it be used for testing, too? Also somewhat an inverse K-S-Test: if the largest deviation is below a certain expected threshold, the data is too evenly distributed?


title: Poisson regression assumptions and how to test them in R

text: I would like to test in what regression fits my data best. My dependent variable is a count, and has a lot of zeros. And I would need some help to determine what model and family to use (poisson or quasipoisson, or zero-inflated poisson regression), and how to test the assumptions.   1. Poisson Regression: as far as I understand, the strong assumption is that dependent variable mean = variance. How do you test this? How close together do they have to be? Are unconditional or conditional mean and variance used for this? What do I do if this assumption does not hold?    2. I read that if variance is greater than mean we have overdispersion, and a potential way to deal with this is including more independent variables, or family=quasipoisson. Does this distribution have any other requirements or assumptions? What test do I use to see whether (1) or (2) fits better - simply `anova(m1,m2)`?   3. I also read that negative-binomial distribution can be used when overdispersion appears. How do I do this in R? What is the difference to quasipoisson?   4. Zero-inflated Poisson Regression: I read that using the vuong test checks what models fits better.  `> vuong (model.poisson, model.zero.poisson)` Is that correct? What assumptions does a zero-inflated regression have?   5. UCLA's Academic Technology Services, Statistical Consulting Group has a section about zero-inflated Poisson Regressions, and test the zeroinflated model (a) against the standard poisson model (b):  `> m.a <- zeroinfl(count ~ child + camper | persons, data = zinb)`   `> m.b <- glm(count ~ child + camper, family = poisson, data = zinb)`   `> vuong(m.a, m.b)` I don't understand what the `| persons` part of the first model does, and why you can compare these models. I had expected the regression to be the same and just use a different family.


title: Comparing multiple classifiers over a single dataset

text: Im doing a comparison between 4 classifiers : A linear SVM, a random forest classifier, a multivariate gaussian and a neural network, for a dataset with 50 features and 10000 examples, the observations belong to 2 classes. What is the best way to compare them ? I mean, should I use F-score, Precision/Recall or maybe logloss to see which one performs better on the test subset?


title: Significance of average correlation coefficient

text: Disclaimer: if you find this question to be too similar to another one, I happy for it to be merged. However, I did not find a satisfactory answer anywhere else (and do not yet have the "reputation" to comment or upvote), so I thought it would be best to ask a new question myself. My question is this. For each of 12 human subjects, I have computed a correlation coefficient (Spearman's rho) between 6 levels of an independent variable X, and corresponding observations of a dependent variable Y. (Note: the levels of X are not equal across subjects.) My null hypothesis is that in the general population, this correlation is equal to zero. I have tested this hypothesis in two ways:   1. Using a one-sample t-test on the correlation coefficients obtained from my 12 subjects.   2. By centering my levels of X and observations of Y such that for each participant, mean(X) = 0 and mean(Y) = 0, and then computing a correlation over the aggregate data (72 levels of X and 72 observations of Y).  Now, from reading about working with correlation coefficients (here and elsewhere) I have started to doubt whether the first approach is valid. Particularly, I have seen the following equation pop up in several places, presented (apparently) as a t-test for average corelation coefficients: $$t = \frac{r}{SE_{r}} = \frac{\sqrt{n-2}}{\sqrt{1-r^{2}}}$$ where $r$ would be the average correlation coefficient (and let's assume we've obtained this using Fisher's transformation on the per-subject coefficients first) and $n$ the number of observations. Intuitively, this seems wrong to me as it does not include any measure of the between-subject variability. In other words, if I had 3 correlation coefficients, I would get the same t-statistic whether they were [0.1, 0.5, 0.9] or [0.45 0.5 0.55] or any range of values with the same mean (and $n=3$) I suspect, therefore, that the above equation does not in fact apply when testing the significance of an average of correlation coefficients, but when testing the significance of a single correlation coefficient based on $n$ observations of 2 variables. Could anyone here please confirm this intuition or explain why it is wrong? Also, if this formula doesn't apply to my case, does anyone know a/the correct approach? Or perhaps my own test number 2 is already valid? Any help is greatly appreciated (including pointers to previous answers that I may have missed or misinterpreted).


title: How to compare two different predictors

text: I have developed a new predictor based on neural networks for a specific problem in bioinformatics. This predictor takes as inputs several features and returns a boolean target value. Additionally i have runned my dataset through other already published predictors (none of which based on neural networks). From all these results i have generated 9 contingency tables (one per predictor) based on the target value and the predictor response like the one below.                                 Target                   |    0     1 |                ----------------|     Prediction  0 |  123    10 |                 1 |    5   171 |                 ----------------      Is there any way to compare these statistical tables in such a manner that i can state that my predictor is better or worse than any of the other predictors supported by a significant p-value? Note that i have the results table for all cases (Ei) in my dataset for all the predictors (Pj), like:                  P1 P2 P3 P4 P5 P6 P7 P8 P9     E1  0  0  1  0  1  0  0  1  1     E2  0  1  1  0  0  0  0  1  1     E3  0  1  0  0  0  0  0  0  0


title: What is an example of a chi-squared distribution?

text: I'm trying to figure out the usage of the $\chi^2$ distribution. In other words, in what kind of situations it occurs and how is it useful in that situations.   I read the wikipedia definition of $\chi^2$. Can someone give an example of a $\chi^2$ distribution so that I can better understand it?


title: Testing if short sequence of small values drawn from binomial distribution

text: An experiment was run 10 times in configuration a and there were A failures, 10 times in configuration b giving B failures, 10 times ... (in all 9 different configurations). I have a list of 9 numbers giving the number of failures for each configuration (out of 10). 4 4 4 3 3 2 0 2 10 I suspect the configuration makes no difference and that the failure counts might have a binomial distribution (it is tempting to classify that 10 as an outlier and remove it). Feeding the above list of values into R's prop.test function, along with the probabilities for a binomial distribution with p=0.3556 (the mean of failure count/10) generates a warning that the Chi-squared approximation may be incorrect; not surprising giving the small sample size. A qq-plot would also suffer from the same sample size issues. Is there a "Fisher's exact test" equivalent for testing whether a sequence of values might be drawn from a binomial distribution? A related but different experiment produced a failure count sequence of (no outliers here :-): 10 8 7 10 9 9 10 8 7


write program bayesian test of linear regression hypothesis In R or winbugs :               h (θi)= µi     βi~exp-family      except Normal distribution     H0 : µi=α+β1x1+ β2x2     H1 : µi=α+β1x1      (B.F12 = (m(x|μ_2))/(m(x|μ_1)) )=???


title: How to find relationships between different types of events (defined by their 2D location)?

text: I have a dataset of events that happened during the same period of time. Each event has a type (there are few different types, less then ten) and a location, represented as a 2D point. I would like to check if there is any correlation between types of events, or between the type and the location. For example, maybe events of type A usually don't occur where events of type B do. Or maybe in some area, there are mostly events of type C. What kind of tools could I use to perform this ? Being a novice in statistical analysis, my first idea was to use some kind of PCA (Principal Component Analysis) on this dataset to see if each type of event had its own component, or maybe some shared the same (ie were correlated) ? I have to mention that my dataset is of the order of 500'000 points $(x, y, type)$, thus making things a bit harder to deal with. EDIT: As noted in the answers below and the comments, the way to go is to model this as a marked point process, and then use R to do all the heavy- lifting, as explained in details in this workshop report : http://www.csiro.edu.au/resources/Spatial-Point-Patterns-in-R.html


title: Linear regression with shot noise

text: I'm looking for the right statistical terminology to describe the following problem. I want to characterize an electronics device that has a linear response $Y = \beta_0 + \beta_1 X + \epsilon$ where $\epsilon \sim N(0,\sigma^2_{ro})$ is a term due to the read-out noise of the device. In order to determine $\beta_0, \beta_1, \sigma^2_{ro}$ I would measure a series of responses $\\{X_i,Y_i\\}$ and apply the standard linear regression toolbox. But I don't know what the $X_i$ are exactly, because I use a source that is affected by shot noise. That is I know that if I set the dial on the source to a certain value $J_i$ then $X_i \sim N(\mu, \mu)$ (a Gaussian with average $\mu$ and variance $\mu$). This looks like an errors-in-variables model of linear regression (http://en.wikipedia.org/wiki/Errors-in-variables_models), where not for the fact that in order to characterize my device over its entire input range, during the measurements I have to change the value of $J_i$, and now the variance of the $X_i$ is not fixed, but it depends on $X_i$ (through J_i), although because of the shot noise if $X_i=X_j$ this does not mean that the variance of $X_i$ is the same as the variance of $X_j$. What is this model called, and are there articles where I can find out such a problem is approached? Or am I formulating in the wrong manner?


title: Example use of ARMA forecasting method

text: I have no experience in forecasting, so can anyone give me a step-by-step example or link to example-real values with ARMA forecasting method application?


title: Why is chi square used when creating a confidence interval for the variance?

text: This is a very basic question. Why do we use a chi square distribution? What is the meaning of this distribution? Why is this the distribution used for creating a confidence interval for the variance? Every place I google for an explanation just presents this a fact, explaining when to use chi, but not explaining _why_ to use chi, and why it look the way it does. Many thanks to anyone that can point me towards the right direction and that is - really understanding why I am using chi when I am creating a confidence interval for the variance.


title: Introductory statistics video courses for social scientists

text: I have just discovered, through this answer, a nice video resource to share with students in introductory statistics courses. The lectures, produced by the Khan Academy, are fairly didactic and well illustrated. I could not find the author of the videos. Is there another set of video courses which I could recommend to my students, who are social scientists taking an introductory course in statistics? The aforementioned set is software-independent and free to watch, two conditions that I need to maintain.


title: Autocorrelation $p^2(1)<\frac{p(2)+1}{2}$

text: I'm struggling with a problem how to show that stationary AR(2) process autocorrelation: $p^2(1)<\frac{p(2)+1}{2}$ * * * My solution attempt: Because $p(2)=\phi_1p(1)+\phi_2p(0)$ and $p(0)=\frac{Var(x)}{Var(x)}=1$ we can also write that (can we?) $p(1)=\phi_1 + \phi_2p(1)$ and $p(1)=\phi_2 + \phi_1p(1)$ if all above are correct, then we can write $p(\phi_1+\phi_2p(1))<\frac{\phi_1p(1)+\phi_2+1}{2}$ which leads to $p(p(1))<\frac{p(1)+1}{2}$ which cannot be correct, since this is not true: $p^2(1)<\frac{p(2)+1}{2}$ <=> $p(p(1))<\frac{p(1)+1}{2}$ That cannot be correct, since $p(p(1))$ must be larger than $\frac{p(1)+1}{2}$... ie. $p(1)=0,1$ then $p(0,1)=0,9>\frac{0,1}{2}+\frac{1}{2}$


title: Is it acceptable to use Cronbach's alpha to assess reliability of questionnaire composed of categorical and conditional items?

text: **Background:** I have data on a questionnaire made up of categorical factual items. Some of them are binary and others have more than two categories. For example:               1. Can you do painting (Yes- 1 point, No- 0 Point)     2. Do you have qualification in painting         (No - 0 point; Bachelor - 1 point; Postgraduate - 2 points)     3. Is there any painter(s) in your family (Yes- 1 point, No- 0 Point)          Etc.      The questionnaire also contains conditional items. E.g., if someone chooses "No" for question 1, he or she will skip question 2 and jump to question 3. ### Question **Can I use Cronbach's alpha on a survey that contains categorical factual items?** I assume questions like these cannot be measured with Cronbach's Alpha. I can't seem to find anyone using cronbach's alpha for such questionnaire.


title: interpretation from loess graph

text: I have measured events X and Y in two strains (A and B) of same species. After calculating the values of these X and Y, I compared them with wild type values (X and Y) to calculate fold change. I have plotted these fold changes of X and Y in A and B, and curve was fitted with LOESS. ![enter image description here](http://i.stack.imgur.com/0Ntfb.png) Now for event X (top plot) LOESS approximates to straight line, but for event Y (bottom plot), the smaller values (Chromosome size < 600 kb) seem to be elevated. My question is by just showing this plot can I claim that smaller chromosomes get more event-Y than larger chromosomes (chromosome size > 600 kb) or do I need to do statistical tests to say this ?? If so what tests do you suggest ? Sorry I am poor when it comes to advanced statistics. Hope I have clearly stated my problem.


title: What is the intuition behind exchangeable samples under the null hypothesis?

text: Permutation tests (also called a randomization test, re-randomization test, or an exact test) are very useful and come in handy when the assumption of normal distribution required by for instance, `t-test` is not met and when transformation of the values by ranking of the non-parametric test like `Mann- Whitney-U-test` would lead to more information being lost. However, one and only one assumption should not be overlooked when using this kind of test is the assumption of exchangeability of the samples under the null hypothesis. It is also noteworthy that this kind of approach can also be applied when there are more than two samples like what implemented in `coin` R package. Can you please use some figurative language or conceptual intuition in plain English to illustrate this assumption? This would be very useful to clarify this overlooked issue among non-statisticians like me. **Note:**   It would be very helpful to mention a case where applying a permutation test doesn't hold or invalid under the same assumption. **Update:**   Supppose that I have 50 subjects collected from the local clinic in my district at random. They were randomly assigned to received drug or a placebo at 1:1 ratio. They were all measured for paramerter 1 `Par1` at V1 (baseline), V2 (3 months later), and V3 (1 year later). All 50 subjects can be subgrouped into 2 groups based on feature A; A positive = 20 and A negative = 30. They can also be subgrouped into another 2 groups based on feature B; B positive = 15 and B negative = 35.   Now, I have values of `Par1` from all subjects at all visits. Under the assumption of exchangeability, can I do comparison between levels of `Par1` using permutation test if I would:   \- Compare subjects with drug with those received placebo at V2?   \- Compare subjects with feature A with those having feature B at V2?   \- Compare subjects having feature A at V2 with those having feature A but at V3?   \- By which situation this comparison would be invalid and would violate the assumption of exchangeability?


title: Assessing associations between many categorical variable pairs

text: I want to check for associations between categorical variables before using some of them as predictors in a logistic regression. Normally I would simply use chi-square, however I have about 15 variables. What is the best way to test for associations between pairs of all 15 variables?


title: What are tools or library to visualize real-estate price on map?

text: What are appropriate tool or library to plot real-estate price on map. Let's say I have list of real-estate listings that consist of the price and coordinate (latitude and longitude). **UPDATED** I understand this can be very controversial. So here's some background. I have millions of real-estate data that consist of price and coordinate. And I have a proficiency in Java and Python programming. But I don't really have a background on data visualization. Basically I want to visualize the data to the users so that they can understand the distribution of the price in their neighborhood. It would be nice if there is a library (for web, i.e. HTML5/CSS/Javascript) to do visualization.


title: Calculating predicted values from categorical predictors in logistic regression

text: ### Context: I am working with an ordinal logistic model and trying to interpret/present the results. The model has two continuous predictors of interests, and a mix of continuous and categorical controls. I was hoping to graph the predicted likelihood of the top outcome (being accepted into a school) across multiple levels of my IVs of interest. I am using R's predict() function to generate predicted likelihoods. For my IVs of interest, I chose a range of reasonable values (i.e. mean +- 1 SD). For the continuous predictors, I can use sensible baseline values (usually 0) because they are mean-centered or standardized. I am trying to work out how to approach the categorical predictors. I've explored my options by plugging in different values, and in most cases the result is just a small shift in the output curve. For one variable however, the differences are huge, so I need to find a way to present results that are general to the different levels of that variable. Perhaps an example would help clarify. In these two graphs, the two IVs of interest are plotted on the x-axis and as the 3 lines. Each graph shows the output given a single level of my troublesome categorical control, "Admitting School" (which has 4 levels total) ![enter image description here](http://i.stack.imgur.com/5rIXW.png) ![enter image description here](http://i.stack.imgur.com/gu1X8.png) Other graphs and R syntax here if you're curious ### Question:   * How should I represent the model across all levels of the categorical variables in a single graph? ### Initial Thoughts:   * Aggregate predicted values across each level of Admitting School with some sort of weighted average.   * This post suggests using the proportion of cases of each type as the input for each variable. As in, if 32% of my cases came from School 1, I would use .32*B-school1 in the prediction formula. I don't know how to do that in R, since those variables are factors, but if it's an appropriate approach, I'm sure I could figure it out. Sorry for the verbosity and thanks in advance for any help.


title: Doing statistical analysis and charting with large data?

text: I have 4301 lines of data from my science project. I have to do unorthodox statistical analysis for my results because they are 4-dimensional; I have three independent variables. With Excel I can do a lot of analysis, even multivariable linear regression, but what I need is a step up: Multivariable nonlinear regression. I have R. Is there a way I can do it with R? Also, is there a way to generate multiple 3-dimensional graphs with slices of the data automatically? EDIT: I've searched with no luck. Is there _any_ way of doing multivariable nonlinear regression automatically?


title: prove that this expression converges in probability to zero

text: Apparently $T^{-3/2} \sum\limits_{t=1}^T{y_{t-1}}u_t$ converges by law to $0.5\times T^{-1/2}\sigma^2 (X-1)$, where $X$ is a $\chi^2(1)$ random variable. $u_t$ is white noise and $y_t$ is an **AR(1)** process, that is $y_t=y_{t-1}+u_t$ with $u_t\sim iid N(0,\sigma^2)$. I want to prove that the above expression also **converges in probability to zero**. Now it appears the variance of this expression is of order $1/T$. I deduct this by noting that the variance of a a $\chi^2$ variable is finite and we have something finite divided by $\sqrt{T}$. This implies that the variance of the whole expression converges to zero. If I am not mistaken, this then also implies that the above expression converges in probability to its mean (?) Now am a little unsure now how to obtain the mean of this assymptotic expression. To prove $0.5\times T^{-1/2}\sigma^2 (X-1)$ converges in probability to zero, can I simply take the expectation this expression by assumging that the expecation of a $\chi^2(n)$ variable is $n$? That would be then $0.5\times T^{-1/2}\sigma^2 (1-1)=0$. Therefore we have $0.5\times T^{-1/2}\sigma^2 (X-1)\rightarrow_p0$? Is this rigorous enough?


title: How to convert epidemiological Relative Risks (RR) to carcinogenic Inhalation Unit Risk Factor (IUR)?

text: Measures of effect in epidemiological studies on air pollution are normally presented as Relative Risk (RR). Classical Risk Assessment procedures (e.g. US-EPA) require instead a Inhalation Unit Risk Factor (IUR), which is multiplied by the exposure concentration to obtain an excess risk estimation. Is it possible to convert a RR to a IUR to be used in risk assessment? Any guideline/publication?


title: Reporting from a likelihood ratio test

text: I don't know if this belongs here or in StackExchange, it is a mixed but probably pretty simple question. How do I normally report a Likelihood Ratio Test? I would love a good reference in your answer, I have searched but could not find any good answers.               > glmm0 <- glmer(yngel ~ (1|nest), data, family=poisson(link="log"))     > glmm <- glmer(yngel ~ age.level + (1|nest) + 0, data, family=poisson(link="log"))     > anova(glmm0, glmm)     Data: data     Models:     glmm0: yngel ~ (1 | nest)     glmm: yngel ~ age.level + (1 | nest) + 0           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)         glmm0  2 682.33 689.38 -339.16   678.33                                  glmm   3 672.37 682.95 -333.18   666.37 11.959      1   0.000544 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      My best guess so far is: I used likelihood ratio test to compare the model with the fixed effect to a model without it. The model including the fixed effect (age-level) was a better fit ($\chi^2(df=?)=11.96$, $p=0.00054$). And I cannot actually figure out how many df to report from that. The there is 2 for one model and 3 for the other, and 1 between them. Thank you for your help.


title: Clustering: Should I use the Jensen-Shannon Divergence or its square?

text: I am clustering probability distributions using the Affinity Propagation algorithm, and I plan to use Jensen-Shannon Divergence as my distance metric. Is it correct to use JSD itself as the distance, or JSD squared? Why? What differences would result from choosing one or the other?


title: Problem calculating joint and marginal distribution of two uniform distributions

text: Suppose we have random variable $X_1$ distributed as $U[0,1]$ and $X_2$ distributed as $U[0,X_1]$, where $U[a,b]$ means uniform distribution in interval $[a,b]$. I was able to compute joint pdf of $(X_1,X_2)$ and marginal pdf of $X_1$. $$ p(x_1,x_2) = \frac{1}{x_1}, \text{ for }\quad 0\le x_1\le 1, \quad 0\le x_2 \le x_1,$$ $$ p(x_1)= 1, \text{ for } \quad 0\le x_1\le 1.$$ However while computing marginal pdf of $X_2$ I am encountering limits problem. The resultant of integral through marginal of $X_2$ is $\log(X_1)$ and the limits are from 0 to 1. As $\log(X_1)$ is not defined for $X_1=0$, I am facing a difficulty. Am I wrong somwhere? Thanks.


title: T-tests for power-law distributed data?

text: I have empirical data consisting of a variable that is arranged in 4 groups. In each of the groups the variable roughly follows a power-law distribution (Pareto distribution). I want to test for differences in means across the group, as I would with a t-test of the mean for normally distributed data. What is the best way of going about doing this?


title: Multinomial distribution for 4 side dice roll

text: Suppose I roll a 4-side dice 25 time to get a final sum between 25 and 100. How I calculated the distribution of probability for each sum between 25 and 100? Thank you


title: isotropic correlation functions

text: Given $l$: characteristic length-scale, and $r$ is the distance of each points. Could you propose some popular correlation functions that are periodic and characterized by $l$ and $r$?


title: Equivalent events

text: Define the term "equivalent events". If $M$ is the event that the number rolled from a die is a prime number, which event can be equivalent to $M$?


title: Is it possible to append training data to existing SVM models?

text: I'm using libsvm and I noticed that everytime I call svmtrain(), I create a new model and that there seems to be no option to put data in an existing model. Is this possible to do however? Am I just not seeing this aspect in libsvm?


title: Adaboost and factor variables

text: My dataset contains both numerical and categorical features like education level, region etc (i use factor variable type for them). I think these variables are important for predicting the outcome of my experiment and training the final classifier. I was recently learning about Adaboost and I want to train boosting classifier on this data. What I don't understand is how should I treat categorical data within the boosting algorithm. Most of R boosting package do not support factor variables as features (except for R's package gbm which supports up to 1024 levels). I think that the biggest issue here is how factor variables are coded, because most of the packages do not support automatic conversion, is that right? What coding is used gbm package? Is it dummy coding?


title: Why do Laplacian Eigenmaps ignore the second KKT condition in the solution calculation?

text: From the original Laplacian Eigenmaps paper, we find the new manifold embedding $Y$ as follows, given the Laplacian and degree matrices $L$ and $D$: $\arg\min_Ytr(Y^TLY)$ subject to ${Y^TDY=I}$. Regarding the constraint ${Y^TDY=I}$: > For one-dimensional embedding problem, the constraint prevents collapse onto > a point. For the $m$-dimensional embedding problem, the constraint presented > above prevents collapse onto a subspace of dimension less than $m-1$ ($m$ > if, as in the one-dimensional case, we require orthogonality to the constant > vector). Using a Lagrange function the above objective is formulated as: $\arg\min_YC(Y)=\arg\min_Y{tr(Y^TLY)+\lambda(I-Y^TDY)}$ The main KKT conditions defining the solutions are: Condition 1: $\frac{\delta C(Y)}{\delta Y}=0:\quad LY=\lambda DY$ Condition 2: $\frac{\delta C(Y)}{\delta\lambda}=0:\quad Y^TDY=I$ > Standard methods show that the solution is provided by the matrix of > eigenvectors corresponding to the lowest eigenvalues of the generalized > eigenvalue problem $Ly=\lambda Dy$. However, the imposed constraint ${Y^TDY=I}$ will not necessarily be met by the eigenvectors corresponding to the smallest positive non-zero eigenvalues of $Ly=\lambda Dy$, especially since condition 2 ($Y^TDY=I$) is ignored according to the text. **Why is the second KKT condition ignored in the calculation of the solution?**


title: Does Akaike information criterion penalize model complexity any more than is necessary to avoid overfitting

text: The AIC penalizes complex models. Obviously a certain penalty for complex models is necessary to avoid overfitting of statistical models: otherwise we would favour a model which is simply a copy of the data itself, and that would tell us nothing. Does the AIC penalize complexity any more than is strictly _necessary_ to avoid this occurance, however? Edit: As an example of what I mean, wiki page for Occam's Razor references a textbook showing that a presumption in favour of simpler models is not required in Bayesian stats: > There have also been other attempts to derive Occam's Razor from probability > theory, including notable attempts made by Harold Jeffreys and E. T. Jaynes. > The probabilistic (Bayesian) basis for Occam's Razor is elaborated by David > J. C. MacKay in chapter 28 of his book Information Theory, Inference, and > Learning Algorithms,[32] where he emphasises that a prior bias in favour of > simpler models is not required. Is the AIC merely reflecting the frequentist equivalent of the same effect, or does it favour simplicity over and above that?


title: Dataset and papers for baseline

text: I'm doing a project about Topic Detection and Tracking in text. I need to perform a baseline so I can compare existing results with mine. I read some papers where they use datasets that are not so easy to get, for various reasons. For example:   1. They use some sort of API to get text from a website at a certain time;   2. They use a dataset from many years ago and are nowhere to be found;   3. They use a dataset of which they ask hundreds or even thousands of dollars for. Of course, to see if my recreation of their baseline was successful, I need the exact same dataset so I can compare mine with their results. I was wondering if someone can point me to a research paper (more is even better) where they use a dataset which is **free** and **still available** for download. And where they use it with different classifiers (Naïve Bayes, kNN, SVM, Decision Trees, etc.) showing their corresponding results.


title: Paired t-test or Independent t-test

text: I am new to this forum so here goes. I have an assignment question I need help with. A scientist has performed some experiments to determine the relationship between the concentration of phosphorus (P) in water and the concentration of the phosphorus in the water after 200 g of sand is added to the water. The idea is that some of the phosphorus will be absorbed into the sand, and we wish to know how much. As part of this experiment, the scientist added approximately 6000 units of phosphorus to a jug of water, mixed the phosphorus so that it was completely dissolved, and then poured all of this water into 3 smaller jugs so that there was an equal amount of water in each jug. This resulted in each jug containing approximately 2000 units of phosphorus. Initial measurements were then taken on the amount of phosphorus in each jug. Sand was then added to each jug, mixed, and the amount of phosphorus in the water in each jug was recorded again. Thus for each of the 3 jugs, there is an initial reading giving the amount of phosphorus in the water before adding the sand, and a final reading giving the amount of phosphorus in the water after adding the sand. I need to find a 95% CI for the amount of phosphorus absorbed into the sand. I haven't included figures from the example, but I don't know whether this is a paired t-test or not. I am thinking it is as each jug will be its own baseline, but I am only familiar with paired t-tests where we test if the difference is zero. Any help would be appreciated.


title: Logit with dummies when certain number of dummies must be used

text: In Stata 12, I'm using logit to model a process that is successful or fails and requires seven workers. I have a larger pool of potential workers and am using dummies to signify which seven individuals were working on the process, which resulted in a success or fail. So the coefficients will report the impact each individual worker will have on probability of success if they are included. However, (again) this process **requires** seven workers. All of the observations show exactly seven of the 30-something worker dummies activated. In reality, it doesn't make sense to imagine the addition or removal of a worker - the process will be fundamentally different if only six workers are present. My results can still show the relative impacts across workers, but I'm worried about the conceptual problem of having only six workers. I would be interested in comparing the impact of each worker relative to some "average" or "benchmark" worker, but I'm not sure how to do that.


title: How to create ROC curve to assess the performance of regression models?

text: I knew that, ROC curve are use to assess the performance of classifiers. But is it possible to generate ROC curve for the regression model? If yes, How?


title: What is iterative bootstrap? How is it used?

text: I have recently stumbled across a mention of "double/triple bootstrap" or "iterative bootstrap". As I understand, each bootstrap sample is bootstrapped again. What is the point? How is it used?


title: Multiple paired samples t-tests?

text: I'm trying to determine the effect of taking a training course on various aspects of being a good quality worker. I have pre-test and post-test data on 7 different ways job performance can be viewed. I'm currently performing 7 different paired samples t-tests on the data, to see the effect on each individual aspect. Is there a test to calculate the overall effect on job performance (by comparing all 7 aspects at once)? I'm trying to use SPSS but not sure how to proceed.


title: Categorical coefficients

text: Using the following table, is it possible to run a linear regression (or some other alternative in R) to calculate the value (intercept?) of "size" and "color" (coefficient?) This is a mock example and I may not in general know which is an intercept and which is a coefficient. Ideally I want some output that suggests               **Size**:     Small = 3     Med = 4          **Color**:     Red = 2     Yellow = 3          **Feature #2**      **Feature #3**... etc      I created this data by using price = (size == small ? 3 : 4) * (color == 'red' ? 2 * 3) In other words, I have a base price and then a single feature of "color" I've written some basic prolog which can solve for all the values in this limited example, but I'm looking for a more general solution               price,size,color     6,small,red     9,small,yellow     8,med,red     12,med,yellow      If I run a lm, I get unhelpful information               summary(lm(price ~ color + size))          Call:     lm(formula = price ~ color + size)          Residuals:         1     2     3     4       0.25 -0.25 -0.25  0.25           Coefficients:                 Estimate Std. Error t value Pr(>|t|)       (Intercept)    8.250      0.433   19.05   0.0334 *     coloryellow    3.500      0.500    7.00   0.0903 .     sizesmall     -2.500      0.500   -5.00   0.1257       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.5 on 1 degrees of freedom     Multiple R-squared: 0.9867,     Adjusted R-squared:  0.96      F-statistic:    37 on 2 and 1 DF,  p-value: 0.1155       running relimp confirms the importance:               calc.relimp(price ~ color + size, q)     Response variable: price      Total response variance: 4.779412           2 Regressors:      color size      Proportion of variance explained by model: 98.67%     Metrics are not normalized (rela=FALSE).           Relative importance metrics:                       lmg     color 0.6533333     size  0.3333333


title: CV of data after nonparametric regression

text: I want to do cross validation after running nonparametric regression on my data. Unlike parametric regression where I can first find my parameters and then easily handle the CV set with these parameters to find my MSE, in the non-parametric version this should be done differently. Given the training and CV set, I decided to do it as follows, take one data point at a time from CV set, look at its neighborhood in the training set using a window function, average the values in the window which gives the new value of the dependent variable. My question is whether the data point itself has to be included when averaging the neighborhood or it is just the neighbors which has to be considered. Your comments are appreciated since I am relatively new to this stuff.


title: Are there any formal definitions or statistical theories about this kind of events?

text: Say a basketball team wins a game. Obviously, at some point in the future, the basketball team will lose a game. So hence, after each game that the basketball team wins, the team is edging closer to their next lost. Are there any statistical tools that can help predict or explain this kind of situations? Another example would be the stock market, which can either close higher or close lower after each day. If a stock market closes higher today, eventually one day it will close lower. That day would be approaching nearer. Are there any statistical tools or tests that can help model or predict these kind of events? For example, the stock market closes up for the past 4 days in a row, so what is the probability that it will close up tomorrow? To get a better idea of what Im trying to explain. Lets say for example the past 5 game results (game 5 being the most recent) of a basketball team are: > Game 1: Win > > Game 2: Win > > Game 3: Win > > Game 4: Win > > Game 5: Lose So from Game 1 to Game 4, after each win, the team's next lost is approaching closer. Are there any statistical tools for these kind of situations?


title: How to: Prediction intervals for linear regression via bootstrapping

text: I am having trouble to understand how to use **bootstrapping** to calculate **prediction intervals** for a linear regression model. Can somebody outline a step-by-step procedure? I searched via google but nothing really makes sense to me. I do understand how to use bootstrapping for calculating confidence intervals for the model parameters.


title: Implementation of the cross validiation

text: I'm attending a course in computational statistics, which should be an applied course. We study different methods, which are important in "reality". One of these topics is Cross Validation. I'm faced with the following problem coming from a homework. We are given a dataset and suppose that the model is of the form $$ Y_i=m(X_i)+\epsilon_i $$ i.e. a nonparametric regression. We want to compute the generalized error using the leave one out Cross Validation score. This should be done by using kernel estimator `ksmooth`, local polynomials and smoothing splines. My first question is very general: Given the data how can I choose by eye a reasonable bandwith for the kernel estimator? See for example the following picture: ![enter image description here](http://i.stack.imgur.com/Fs3BD.png) Since this data looks rather "wild" it is for me not clear how to choose a bandwidth. My first attempt was just to run ksmooth playing around with different bandwidth. But as I said, here the data is wild, so it is (for me) hard to determine a reasonable bandwidth. The second problem is more concrete about the problem described above. So far I have the following code:                  cv <- function(data,used.function)     {       n <- nrow(data)       cv.value <- rep(0,length(n))       for (i in 1:n){           new.data <- data[-i,]           cv.value[i] <- used.function(new.data[,1],new.data[,2],data[i,1])       }       ## MSE       return(1/n*sum((new.data[,2]-cv.value)^2))     }               ### kernel estimator usind nadaraya-watson:     fcn1 <- function(reg.x, reg.y, x){       return(ksmooth(reg.x, reg.y, x.point = x, kernel = "normal", bandwidth = h)$y)     }     ### CV-score for kernel estimator:     (cv.nw <- cv(real.data, fcn1))      the function cv should be general such that I can apply local polynomials and smoothing spline too. The variable real.data contains the data. It is a $n\times 2$ matrix which store all the $x$ values and $y$ values. In the function body of cv I perform a leave one out cross validation. However using this code gives for cv.nw a NA. What is wrong with my code? I'm very thankful for your help.


title: R, getting data.frame factor name by index

text: I have a data frame, and I'm trying to get name of its factor (as string) by index. For example,               var1=factor(c(0,2,2,6,7),levels=0:7)     var2=factor(c(2,2,0,1,6),levels=0:7)     obs1=data.frame(var1,var2)     obs1[1]      returns                 var1     1    0     2    2     3    2     4    6     5    7      Is there a way to return just "var1"?   Thanks.


title: Multidimensional Differential Entropy

text: I am looking for a measure of entropy over multiple random variables, each with values between 0 and 1. Intuitively, it seems possible to talk about the expected value of information of several variables, which is entropy, but I am not sure how to go about. I know that in the continuous case, for a single value, we should look at the differential entropy, given by $$ - \int_0^1 f(x) log(f(x)) \partial x $$ , where $ f(x) $ is the probability distribution function, which in my case is a Beta Distribution, fitted to the data. I am fitting a Dirichlet distribution for the multiple variable case, since it is the multivariate generalization of the beta distribution. How do we measure the entropy of such a joint probability density function?


title: Spearman's Rho - from partial ranked variables

text: I have two variables which represent two performance measures.   I ranked a finite set of elements according these two variables.   Therefore, I have to ranks. Suppose the ranks are performed in descending order (the highest the measure the highest the value in the decision process of each element). For instance, an ordered set $\Omega$ of $10$ elements according the two measures A and B.               R#A = [5 3 1 9 2 10 6 7 4 8];       R#B = [9 7 4 8 5 6 1 3 2 10];      Suppose that now I truncate R#A and R#B in order to select the "top 5" elements:               R#A_5 = [5 3 1 9 2];     R#B_5 = [9 7 4 8 5];      In your opinion it is still possible to get the Spearman's rho correlation coefficient with these two partial orders? I know that   1) We are in the second step of the Spearman's because we are alrady dealing with ranks.   2) The Sample size is very low but it is just for explanation.


title: Fitting a parameterized 2D distribution

text: I frequently deal with datasets that rely on a binary classification based on a Pulse-Shape-Distribution (PSD) discrimination value. This is based on a 2d fit of the PSD value vs the amplitude of the data, a typical graph is show below ![psd](http://i.stack.imgur.com/sTRyw.png) The fits you see are the means and +/- 1 sigma for the two bands, where the mean and sigma is determined by chunking the 2d distribution along the x axis and performing a two gaussian fit along each slice. This works in most cases but sometimes the fits can bounce around a lot and on the lower end are less accurate. What I would like to do is to somehow parameterize the means of the gaussian to a functional form. For reasons based on the underlying physics of the data, it can be assumed that the mean of the lower distribution is constant and the mean of the upper is a monotonically decreasing curve (a cubic would probably be easiest to model). The integral of the two distributions relative to each other is not important, nor is the distribution of values along the x-axis. This changes for each dataset and is not important. What is needed is just the fits you see above. It seems that there are two options, both of which I am having a hard time thinking about how to implement: * Continue doing a 1d slicing fit, but constrain the values of the means in such a way that they follow a functional form * Fit a 2D distribution directly, but unsure of how to deal with the bin value and amplitude distribution not mattering. Is there a specific name, class of algorithms, or recasting of the problem that can help me on the path to implement a sort of 1.5-dimensional fit?


title: What is a conditional quantile function?

text: What is > conditional quantile of $Y(t)$ given $Y(t-1)$ where $Y(t)$ is a univariate > time series (they call it conditional value at risk in risk management). Can anybody explain this? Thanks


title: Estimating the response variable distribution in bivariate simple linear regression

text: So let's say I have a bunch of $IID$ samples $(X_n,Y_n)_n$ from some unknown distribution $p(x,y)$, where $Y$ is the response variable and $X$ is the explanatory variable. And suppose I plotted these points and saw a pretty obvious positive linear relationship, (for instance maybe $X$ is amount of sunlight and $Y$ is sunflower height, or something like that). Then suppose I wanted to try to get a good idea of what the distribution for $Y$ might be. Well I could just look at the $Y$ coordinates from my samples and plot a histogram and then fit a curve to it, or something like that, but I believe that would only make sense under the assumption that the $Y_n$ are $IID$, when all I know is that they are identically distributed, and not necessarily that they are independent. So instead what I do is plot the least-squares regression line $y=f(x)$ for my data, and then try to calculate the marginal distribution for $Y$: $$p(y)=\int_{x\in\Omega}p(x,y)dx=\int_{x\in\Omega}p(y|x)p(x)dx.$$ Now what I do is assume there's some standard normal variability about my regression line, and so I set $p(y|x) = N(f(x),1)$. However the problem in finding $p(x)$ is the same problem I have in finding $p(y)$, namely, my samples $X_n$ could very well not be independent. Assuming I'm understanding the theory correctly, how do I go about resolving this? **Edit** : Now I'm thinking maybe the $Y_n$ are independently sampled, but how would one go about proving this?


title: For what p-values can you say the data approached significance?

text: I am in the field of clinical psychology, where statistical significance is commonly said compared to p = .05. For what values of p can I say data _APPROACHED_ significance?


title: What test to use here?

text: I have a model created from a prior dataset and I want to know if it can be used on my current dataset. The original model gives the expected percentage of times an email is opened x hours after it is sent. I have data for another type of email showing how many emails were opened x hours after it was sent. What should I do to determine whether I can use the first model for the second dataset?


title: Expectation of density ratio of two iid variables

text: Let $X \sim N(0,1)$ and $Y \sim N(0,1)$ be independent RVs and let $f$ be their density function. I'd like to compute the expectation of the density ratio \begin{align} \mathbb{E}\left[\frac{f(X)}{f(Y)}\right] & = \mathbb{E} \left[e^{\frac{-X^2}{2}} e^{\frac{Y^2}{2}}\right] \\\ & = \mathbb{E} \left[e^{\frac{-X^2}{2}} \right]\mathbb{E} \left[e^{\frac{Y^2}{2}}\right]. \hspace{5mm} (1) \end{align} My initial thought was that we would simply have $\mathbb{E}\left[\frac{f(X)}{f(Y)}\right] = 1$ whenever $X,Y$ are i.i.d. But my Monte Carlo simulations are not converging. Seeing that, for example, $Y^2 \equiv Z \sim \chi^2_1$, I tried using the fact that $\mathbb{E}\left[g(Z)\right] = \int_{-\infty}^\infty g(z)f_Z(z) dz$ (Law of the Unconscious Statistician) with, $g(Z) = e^\frac{Z} {2}$ and $f_Z(z)$ being the Chi-squared density with 1 df, to evaluate the factors in $(1)$. However, I getting the integral: \begin{align} \mathbb{E}\left[e^{\frac{Z}{2}}\right] = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty \frac{1}{\sqrt{z}}dz \end{align} Am I going in the right direction? Thanks!


title: In RVM, is the kernel allowed to depend on the full dataset?

text: I want to use Relevance Vector Machines and need to define my custom made kernel. I was wondering if it is allowed for the kernel to depend on the full dataset. For exampe, I can calculate a certain weight vector from all the data and call this vector $w$ and then use the kernel K(x1,x2|w) (i.e. $w$ is a vector of parameters for the kernel). Are there any potential problems with this idea? Mbk


title: Obtaining Cohen's d from z-test statistic

text: After an ANOVA, I run post-hoc tests with the `glht` function in the `multcomp` package. `glht` returns a table similar to this one:               Linear Hypotheses:                                                                 Estimate Std. Error z value Pr(>|z|)       Control.Level1.Level1-Control.Level2.Level1==0              0.053125   0.019228   2.763   0.0477 *     Experimental.Level1.Level1-Experimental.Level2.Level1==0   -0.018750   0.027192  -0.690   1.0000       Control.Level1.Level2-Control.Level2.Level2==0              0.003125   0.019228   0.163   1.0000       Experimental.Level1.Level2-Experimental.Level2.Level2==0   -0.025000   0.027192  -0.919   1.0000      **If I need to calculate Cohen's d for the significant comparison, how can I do that, since I have z-value rather than t-test?**


title: Determine linear relationship between input and output with multiple input variables

text: Say that your input variables are in $\mathbb{R}^2$ with a univariate output variable. Say that you want to determine whether the output variable is a linear combination of the input variables. If you just had one input variable you could easily just graph input vs output and see how linear the relationship is. But if you have two input variables I don't think that you can simply say that that $y$ is a linear combination of $x_1$ and $x_2$ iff $x_1$ vs $y$ is linear and $x_2$ vs $y$ is linear. It seems like the reverse should be true, but I'm not sure if the forward argument is valid. e.g. you could have two non-linear polynomials whose linear combination is linear (couldn't you?) I suppose I could just fit a linear model and look at the residuals, but I was wondering if there was a more a priori approach to determining linearity (i.e. some sort of visualization).


title: Explaining the difference between error propagation and "grand variance"?

text: I have the following data concerned with simple radiation detection. I collect a background dose $B$ for 100 seconds and find the standard deviation and get the uncertainty $E_B$. I then take 10 measurements with a source for 20 seconds and find the standard deviation for each trial. $$ X_1 \pm E_1 $$ $$ : $$ $$ X_{10} \pm E_{10} $$ It's easy enough to find the total dose for each measurement using error propagation: $$ T_1 =X_1 - B \pm \sqrt{E_1^2 + E_B^2} $$ Finding the mean is easy enough now with all of my total values, but I'm getting different answers using different error propagation techniques, first I'm using the following formula: $$ s_f^2 = \frac{\partial f}{\partial x_1}^2E_1 + .. +\frac{\partial f}{\partial x_n}^2E_n $$ I get $$ E_{\bar T} = \frac{\sqrt{E_{T_1}^2+..+E_{T_{10}}^2}}{10} $$ I've also found this link which details how to find the "grand variance." I've used both methods and found different answers. Are there certain circumstances where the error propagation formula is not valid and the more rigorous method mentioned in the link has to be used?


title: How can the residuals and the fits of a regression model be correlated?

text: In the context of statistical regression and OLS, I understand that $Y$ lives in a n-dimensional space and has variance $\sigma^2$. The Hat Matrix, $H$, corresponds to an orthogonal projection of a given observation of $Y$ onto the column space of $X$, providing an estimate of $Y$, called $\hat{Y}$. Similarly, $(I-H)$ corresponds to an orthogonal projection of $Y$ onto a hyperplane perpendicular to the column space of $X$. Now, if $Y$ (and, by consequence, the error terms) are uncorrelated, how can its orthogonal projection onto two different planes become correlated? A geometric intuition would be helpful.


title: Decomposing total sum of squares

text: Consider the general linear regression model: $$y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \cdots + \beta_px_{ip} + \epsilon_i = \mathbf{x}_i^t \beta + \epsilon_i$$ where $\textbf{x}_i = (1,x_{i1},x_{i2},\cdots,x_{ip})^T$, $\beta=(\beta_0,\beta_1,\cdots,\beta_p)^T$ and $\epsilon_i$ are iid N(0,$\sigma^2$). I would like to see a complete proof of the following identity from first principles: $$\sum_{i=1}^n(y_i - \bar{y})^2 = \sum_{i=1}^n(\hat{y}_i - \bar{y})^2 + \sum_{i=1}^n(y_i - \hat{y}_i)^2$$ where $\hat{y}_i= \mathbf{x}_i^t \hat{\beta} $ ($\hat{\beta}$ is the least square estimator, $\bar{y}$ ia the sample mean of $y_i$). I know that the two terms on the right can be obtained by subtracting and adding $\hat{y}_i$ on the left side. But this introduces a "cross term": $$\sum_{i=1}^n2(\hat{y}_i - \bar{y})(y_i - \hat{y}_i)$$ Many texts claim that this is zero, but I have not seen a general proof of this statement. How can this be shown?


title: Probability of ($x\le 2Y$)

text: $X$ and $Y$ are independent and their probability density functions are $$f_X(t)=f_Y(t)=\left\\{\begin{array}{l} e^{-t},\:\text{if $t \geq 0$;} \\\ 0,\:\text{otherwise.}\end{array}\right.$$ $P(X \leq 2Y)$=? (The probability of $X \leq 2Y$)


title: Meaning of a perfect sequence of a graph

text: I'm reading Steffen L. Lauritzen's _"Graphical Models"_ and I'm having trouble fathoming the meaning behind the concept of a perfect sequence. Here's how it's defined in the book (page 14/15): > Given an undirected, marked graph $\mathcal{G} = (V,E)$. Let > $B_1,\ldots,B_k$ be a sequence of subsets of the vertex set V. Let $$H_j = > B_1 \cup \ldots \cup B_j$$ $$R_j= B_j \backslash H_{j-1}$$ $$S_j = H_{j-1} > \cap B_j$$ The sequence is said to be _perfect_ if the following conditions > are fulfilled: > > (i) for all $i > 1$ there is a $j < i$ such that $S_i \subseteq B_j$; > > (ii) the sets $S_i$ are complete for all $i$; > > (iii) for all $i > 1$ we have $R_i \subseteq \Gamma$ or $S_1 \subseteq > \Delta$ Here, $\Gamma$ and $\Delta$ are subsets such that $V = \Delta \cup \Gamma$ with $\Gamma \cap \Delta = \emptyset$, which represent quantitative and qualitative variables, respectively (I don't think this is really relevant to my question but I'm mentioning it just in case). The thing is, I can make up a simple graph and work out an example but the meaning behind this is not clear to me. What is a perfect sequence supposed to mean, intuitively? Also, should I interpret $R_1$ as $R_1 = B_1 \backslash H_0 = B_1 \backslash \emptyset = B_1$? This is not mentioned in the book. **NB** I've also asked this question on math.stackexchange.net but I thought that, since graphical models are statistics related, I would also ask here. If this does not sparkle with the moderators then feel free to close this question.


title: Understanding R output in Logistic Regression

text: I am following an example here on using Logistic Regression in R. However, I need some help interpreting the results. They do go over some of the interpretations in the above link, but I need more help with understanding a goodness of fit for Logistic Regression and the output that I am given. For convenience, here is the summary given in the example:               ## Call:     ## glm(formula = admit ~ gre + gpa + rank, family = "binomial",      ##     data = mydata)     ##      ## Deviance Residuals:      ##    Min      1Q  Median      3Q     Max       ## -1.627  -0.866  -0.639   1.149   2.079       ##      ## Coefficients:     ##             Estimate Std. Error z value Pr(>|z|)         ## (Intercept) -3.98998    1.13995   -3.50  0.00047 ***     ## gre          0.00226    0.00109    2.07  0.03847 *       ## gpa          0.80404    0.33182    2.42  0.01539 *       ## rank2       -0.67544    0.31649   -2.13  0.03283 *       ## rank3       -1.34020    0.34531   -3.88  0.00010 ***     ## rank4       -1.55146    0.41783   -3.71  0.00020 ***     ## ---     ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1     ##      ## (Dispersion parameter for binomial family taken to be 1)     ##      ##     Null deviance: 499.98  on 399  degrees of freedom     ## Residual deviance: 458.52  on 394  degrees of freedom     ## AIC: 470.5     ##      ## Number of Fisher Scoring iterations: 4        1. How well did Logistic Regression fit here?   2. What exactly are the Deviance Residuals? I believe they are the average residuals per quartile. How do I determine if they are bad/good/statistically significant?   3. What exactly is the `z-value` here? Is it the normalized standard deviation from the mean of the Estimate assuming a mean of 0?    4. What exactly are Signif. codes? Any help is greatly appreciated! You do not have to answer them all!


title: Missing data - best way to resolve?

text: As part of a school project I have to evaluate the impact of pre-school on the performance of children in academic tests (cognitive skills) as they enter kindergarten. I'm given a dataset with entries of 17'000 children, however, there are some problems with the data given: in 1696 cases, I do not have any data on their performance in these tests, in some more cases the gender, mother's education level, the father's education level etc. are not reported. All in all I do have almost 9000 children left with a complete set of entries. What is the best way to deal with this data? Just deleting these variables seems to violate the assumption of "random sampling" as there could be self- selection bias. Should I use multiple imputation?


title: Cox regression problem with time varying covariates

text: I was working on a cox regression problem and was refering to this link: http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_phreg_sect050.htm I have few questions regarding to this: As in the example mentioned, can we feed missing values to model? Also how do we handle multicollinearity in such a data setup.


title: Model with non-linear transformation

text: I don't understand this concept well and need help. I was choosing whether to use a linear model or apply a non-linear transformation in my model formula. To do a diagnostic, I quickly plotted my data: `plotalldaily <- ggplot(amsd, aes(ImpressionsA, Leads.T)) + geom_point(color="orange")+geom_smooth()` ![enter image description here](http://i.stack.imgur.com/xppJ2.jpg) From the plot, I guessed that a cubic polynomial transformation of my x variable should give me a better model fit. I referred this: http://www3.nd.edu/~rwilliam/stats2/l61.pdf . On Page 5, there is an explanation of polynomial models with cubic terms. So I checked the model fit using two formulae- one with the non-linear transformation & another simple linear: `test1 <- lm(Leads.T~ImpressionsA, amsd)` `test2 <- lm(Leads.T~I(ImpressionsA^3), amsd)` Strangely, the linear relationship is giving me a better model fit: Lower Standardized Error, Higher R-squared and Better Residuals Distribution.               TEST 1 Residuals      ![enter image description here](http://i.stack.imgur.com/VjYuO.jpg)               TEST 2 Residuals      ![enter image description here](http://i.stack.imgur.com/BM2nQ.jpg) I don't know what to make of it. Which model should I fit and what other kind of transformations should I try?


title: Using density estimate to look for problems in data collection - with R

text: I am analysing some data coming from different media (tv, press, ...). I would like to test if the data I have are good (if the collect was flawless), identify big issues and eventually understand what I should not rely on. One thing I observe, is that the density estimates (performed on the number of items produced by each media per day) for 4 out of 5 media closely resemble each other - even if the number of items collected for each media varies a lot (I mean a lot). Instead the fifth media (the radio) produces a weirder density estimate. The number of items collected for this media is not too small, if compared with the other media in my database. I use the geom_density() function in the R package ggplot2. Is this a good approach to observe if there were problems collecting the data? How can I investigate further? Any help appreciated!


title: Statistical issues with aggregating annual survey data from multiple years?

text: I am using a national telephone survey conducted every year by the CDC called the Behavioral Risk Factor Surveillance System (BRFSS) to answer a question about breast cancer screening rate in low-income counties vs. higher income counties in one state. The number of respondents in one year from the low- income counties is rather low, so I am interested in aggregating the survey data from 2011, 2012, and 2013. The survey methodology did not change substantially in that time. Here is a CDC document on the BRFSS methodology: http://www.cdc.gov/brfss/data_documentation/PDF/UserguideJune2013.pdf My question is: are there statistical issues if I aggregate the three years of data to have a larger sample of data? I will need to calculate crude screening proportions but also eventually use logistic regression as well to control for socioeconomic confounders. I appreciate your help!


title: Controlling for respondents' propensity to rate high/low or moderate/extreme

text: I am doing a simple linear regression to test the relationship between voters' self-reported policy stands and their approval ratings for two hypothetical presidential candidates in a simulated election campaign. Participants first report their own position on a range of issues (e.g. 'The government should increase taxes to fund better public education') on a 0-100 scale where 0='strongly disagree' and 100='strongly agree'. Later in the experiment, they also rate the two candidates (Candidate A & Candidate B) on similar issues, using the same scale where 0='strongly opposed' and 100='strongly in favour'. My question is: Is there a 'standard'/'common' way to control for the way respondents use the scales in general (i.e. across all questions in the experiment). Some respondents are much more cautious in their ratings, not rating anything outside of 40-60, while others are more extreme, entering values towards the ends of the 0-100 scale. Although in some cases this might reflect respondents' 'true position', in many cases it seems it is more a reflection of the way they understand the scale. Similarly, some respondents rate almost everything low, while others tend to rate most things high. One suggested possibility might be to use the rating of Candidate B on issue X as a control in the model for Candidate A's rating on issue X (so Candidate B's rating on issue X is one of the independent variables along with age, income, respondents' position on issue X, Y, Z etc with the dependent variable Candidate A's rating on issue X). Would this be a reasonable solution? Or perhaps a Multivariate linear model? With the two candidates' positions as the 2 dependant vars? Would this be a better way of doing this sort of analysis in any case? Many thanks in advance


title: Low accuracy for training in image classification

text: I'm a newbie using LinearSVM to train the classifier.   I labelled the images of 'buildings' as 1 and the others as -1. The training result is as follows : ![enter image description here](http://i.stack.imgur.com/oXmNA.png) and ![enter image description here](http://i.stack.imgur.com/5vduY.png) As you can see in the image some of the buildings have positive scores and others are having negative scores. Also the images labelled -1 are getting a positive score. Why is this happening?


title: Textbook with list of hypothesis tests and practical guidance on use

text: I found a list of statistical tests along with practical guidance on Wikipedia. Can anyone point me to something similar, but in a textbook form? I'm interested in particular in practical guidance (along the lines of "should have n>30 for Z-test")


title: Plotting a piecewise regression line

text: Is there a way of plotting the regression line of a piecewise model like this, other than using `lines` to plot each segment separately, or using `geom_smooth(aes(group=Ind), method="lm", fill=FALSE)` ?               m.sqft <- mean(sqft)     model <- lm(price~sqft+I((sqft-m.sqft)*Ind))     # sqft, price: continuous variables, Ind: if sqft>mean(sqft) then 1 else 0          plot(sqft,price)     abline(reg = model)     Warning message:     In abline(reg = model) :       only using the first two of 3regression coefficients      Thank you.


title: Understanding regressions - the role of the model

text: How can a regression model be any use if you don't know the function you are trying to get the parameters for? I saw a piece of research that said that mothers who breast fed their children were less likely to suffer diabetes in later life. The research was from a survey of some 1000 mothers and controlled for miscellaneous factors and a loglinear model was used. Now does this mean that they reckon all the factors that determine the likelihood of diabetes fit in a nice function (exponential presumably) that translates neatly into a linear model with logs and that whether the woman breast fed turned out to be statistically significant? I'm missing something I'm sure but, how the hell do they know the model? Thanks


title: Neural Networks - Calculating delta in Backpropogation

text: I'm developing software to create a neural network. I have the forward propogation code done, but when I started working on the algorithm for back- propogation I ran into a problem. I'm having problems calculating the delta in my backpropogation method. Let's say in my neural network I have 4 layers. I have 4 nodes in layer 1, 4 nodes in layer 2, 3 nodes in layer 3, and 1 node in layer 4. And let's say I am trying to find delta2 [delta2 = Theta2(T) * delta3 .* g'(z3)], I need to multiply Theta2 (4x4) by delta3 (3x1). This clearly shouldn't work because of the size of the matrix and vector, where am I making my mistake? Is my array for the weights, Theta, wrong? I have it of the dimension (# of nodes in current layer X # of activations from last layer).


title: How many failed both classes

text: I think this might be too easy compared to other questions on this site but I really can't find the solution.               60% of students passed class A     70% of students passed class B     80% of students passed at least one of both classes      a. What % of students passed both classes? b. What % of students passed class A but failed class B? I actually have the answers for both questions but no idea how to get there (a: 50% and b: 10%). Also, if anyone has a link to some theory about this kind of problems or more similar exercises I'd appreciate that a lot! Not sure how to look this up. And another question (but I am not sure if I have all the numbers necessary for this one): If 51.8% of population has obesity and 48% of the population are male,.... what % of the male population has obesity if there's 10% more females (than males) without obesity? Can this be solved with these numbers? I might have forgotten some part of the question so I'm not sure about this one.


title: Why use d-prime instead of percent correct?

text: In signal detection theory, people often use $d'$ to assess performance. Apart from the fact that $d'$ is in $z$ units (units of measurement transformed to standard deviation units, i.e., $z$ scores), making it comparable regardless of the original units of measurement. I can't see what the advantage in analysing $d'$ instead of proportion correct is. Don't both account for bias? Would they follow the same shaped ROC curves?


title: Importing time series from SQL base into R

text:
I have a number of regular daily measurements in a MySQL database that I'd like to manipulate using R. When it's returned from RMySQL, it looks like this:               > memdata             date    vsize     1  2011-04-22 3535.178     2  2011-04-23 5680.516     3  2011-04-24 5468.914     4  2011-04-25 4761.044     5  2011-04-26 4403.515     6  2011-04-27 4459.155      7  2011-04-28 4889.884     8  2011-04-29 5290.908     9  2011-04-30 5370.952     > str(memdata)     'data.frame':   9 obs. of  2 variables:      $ date : chr  "2011-04-22" "2011-04-23" "2011-04-24" "2011-04-25" ...
      $ vsize: num  3535 5681 5469 4761 4404 ...      Since many of the time series libraries expect `ts` objects, I'd like to convert the data frame into one, but there doesn't seem to be a straightforward way of doing that. This site has a lot of good examples, but none that deal with daily data. Any help would be greatly appreciated!


title: High frequency data series cleaning in R

text: I am looking at timeseries data in foreign exchange and bond markets (to test for reversion on extreme moves). Unfortunate "tick" data, namely high frequency data, is prone to many problems, and they obviously can significantly mess with the analysis. I'd like to know which R library can help with the following type of fairly frequent data cleaning problems: 1) one spike: ![enter image description here](http://i.stack.imgur.com/GCGja.png) This is typically created when one market maker prints a wrong quote in one tick, but there would have been no tradability at that price because it lasted for a split second. I'd like to eliminate the spike (but only if there is only one (or maybe 2) prints) 2) bid ask gapping: ![enter image description here](http://i.stack.imgur.com/bAN4d.png) In this case the market is fairly illiquid and the data algorithm is jumping between bids and asks (in this case 2bps wide) causing this weird cloud. Where should I start to clean this stuff, obviously trying to throw out the least amount of real data. I realise that the maxim of "look at the data" applies here, but when you're looking at 1000 series each with 100 days of data, you can see how this will become quickly impractical so I need some automated help. I'll also look at Python language methods if they're available or better.


title: Computing descriptives statistics for sites and locations based on literature search with sites having varying numbers of time points

text: I've collected a data set from the literature. These data are largely binomial: they are clutch sizes of a bird that lays only 1 or 2 eggs in most areas, and in some areas will lay 3. However the only information I have compiled are means, standard errors and sample sizes (# nests) from various sites. Some sites have multiple values reported from them (from different years), and I have split the sites into 3 categories based on location (A, B and C). I want to do 2 things that are essentially the same:   1. Produce one value (± SE or SD) for each site (so a site mean)   2. Produce a location mean ± SE or SD (so an average of all the sites within a given location category). Is it possible to do this relatively simply while keeping the information contained with the reported standard errors from each site? Currently, I've taken a mean of the means of each site within a location category, and just calculated the SE from the mean of means, but that's losing all the reported variation from each site.


title: Tests on binomial distribution

text: I have the observations $X(n)$, where $X(n)$ is the realization of a binomial random variable with probability of success $p(n)$, and with $Y(n)$ trials. The observations are independent across $n$. I would like to test the null hypothesis H0: $p(1)=(2)=\cdots=p(N)=0.5$. Is there a standard recommended test? An approach would be to perform a multiple comparison test with a correction of the significance level, but I wonder if othe methods would be possible. If $Y(n)=$const, I could have used a goodness-of-fit test, but this doesn't apply here. Suggestions welcome!


title: Explain the xkcd jelly bean comic: What makes it funny?

text: I see that one time out of the twenty total tests they run, $p < 0.05$, so they wrongly assume that during one of the twenty tests, the result is significant ($0.05 = 1/20$). xkcd jelly bean comic - "Significant"   * Title: Significant   * Hover text: "'So, uh, we did the green study again and got no link. It was probably a--' 'RESEARCH CONFLICTED ON GREEN JELLY BEAN/ACNE LINK; MORE STUDY RECOMMENDED!'" ![xkcd comic 882 - "Significant"](http://i.stack.imgur.com/QZsY3.png)


title: Bayesian estimation of multivariate Gaussian from noisy observations with known error variances

text: I have a dataset $\mathbf{D} = \\{ (\tau_i, \Gamma_i) : 1 \le i \le n \\}$ of observations $\tau_i = X_i + \epsilon_i$ from a $p$-dimensional Gaussian $X_i \sim \mathcal{N}(\mu, \Sigma)$ contaminated by additive noise $\epsilon_i \sim \mathcal{N}(0, \Gamma_i)$ where the $\Gamma_i$ are diagonal and observed. This gives a likelihood function of $$ \begin{aligned} L(\mu,\Sigma) &= \prod_i \int_{\mathbb{R}^p} \mathcal{N}(x_i | \tau_i, \Gamma_i)\, \mathcal{N}(x_i | \mu, \Sigma) \,dx_i \\\ &= \prod_i \mathcal{N}(\tau_i | \mu, \Sigma+\Gamma_i). \end{aligned} $$ Normally I would use a Normal-inverse-Wishart prior to get a Multivariate Student's t predictive distribution, but it's not clear how to do that with this likelihood. Is it at least possible to get a "sensible" approximation to the predictive distribution of $X$? Is there there some other prior I should consider?


title: Creating summary tables / graphics in R for categorical data

text: I am trying to summarise data from a household survey and as such most of my data is categorical (factor) data. I was looking to summarise it with plots of frequencies of responses to certain questions (e.g., a bar plot of percentages of households answering certain questions, with error bars showing confidence intervals). I found this excellent tutorial which I had thought was the answer to my prayers (http://www.cookbook-r.com/Manipulating_data/Summarizing_data/) but turns out this is only going to help with continuous data. What I need is something similar that will allow me to calculate proportions of counts and standard errors / confidence intervals of these proportions. Essentially I want to be able to produce summary tables that look like this for each of the questions asked in my survey data:               # sex condition  N(count) proportion SE of prop.  ci of prop     #   F   aspirin  4        0.9979145  0.4989572    1.5879046     #   F   placebo 12        0.5247655  0.1514867    0.3334201     #   M   aspirin  7        1.0674848  0.4034713    0.9872588     #   M   placebo  3        0.5291503  0.3055050    1.3144821      I would then want to plot barplots in ggplot (or similar) using these summary data. ![enter image description here](http://i.stack.imgur.com/X69Mj.png) I had thought to amend the code provided in the tutorial above to calculate the columns above, though as a relative newcomer to R, am struggling a little! I wondered if anyone could help or knows of packages / code that would do the job for me!


title: Small data set and validity of Multilevel Regression results

text: I want to estimate a hierarchical multilevel regression (two levels). My data set is fairly small: 25 groups and 255 individual observations distributed unequally among groups (group member size varies between 4 and 55). Several books point out that the small sample size could lead to biased estimates of the random effects. There was no clear conclusion which reasons are responsible for the bias. My question is as follows: which of the statements are true? \- Adding (more) regressors to Level 1 (individuals) will reduce the validity of the estimation results \- Adding (more) cross-level intaraction effects will reduce the validity of the estimation results \- Adding (more) random slopes will reduce the validity of the estimation results Could a MCMC approach help in my example? Thanks in advance, Linus holtermann@hwwi.org


title: Significance of correlation with only 11 datapoints

text: I am currently working on my thesis and calculating some simple correlation coefficients. I only have 11 datapoints available, but I would like to assess the statistial significance of my results. The asymptotic p-values are not significant at the 10% level. Because of the small sample, I additionally computed bootstrapped confidence intervals, which are, to my surprise, different from zero at the 10% level. My question: Can I claim statistical significance of my results according to my bootstrapped confidence intervals? Or is this not valid because of the small sample? Cheers and many thanks in advance! Chris


title: Three way anova

text: I am really hoping you can shed some light on a topic that I am struggling with. I have a 2x3x4 repeated measured anova. I have a significant 3-way interaction, and I want to make sure that I am using the correct post hoc comparisons and not violating any key statistics theory. If we consider my 2x3x4 RM ANOVA to be AxBxC. With the significant 3-way interaction I have performed the following post-hoc analysis: AxB at each level of C AxC at each level of B BXC at each level of A Is it typical for post-hoc tests on a three way interaction to perform all of the above 3 comparisons or would one typically just pick one or 2? They are all an important part of the story. Also regarding adjustment for multiple comparisons, I have used SIDAK to adjust for for this within each of the above 3 tests. I am assuming that I do not need to adjust each p value further to account for the total number of comparisons in all 3 of the above? Finally, most stats textbooks state that if you have an interaction term you can't reliably interpret the main effects. If I only have an interaction term for AxB, not AxC or BxC, is it still ok to interpret the main effect? Thanks for you help,


title: Hypothesis testing with latent variables

text: **Summary:** I created a model that tries to explain success factors in crowdfunding initiatives. I collected the data by scraping it directly from the crowdfunding platform. That of course limits my research in terms of data availability. In order to come up with conceptual model I adopted a theory, that assumes latent variables (such that I do not observe). What I observe is measures or proxies (directly collected form the website) that to some extent depict that latent variables which are embedded in my conceptual framework. Of course I set few hypothesis based on that latent variables. **Example:** I have adopt a theory about the online community reputation. For measuring reputation I set 2 measures - N of successfully funded projects, N unsuccessfully funded projects. That example is quite simple, however, but for some other concepts I have multiple measures on different scales. **Questions:** 1) How should I test my model reliability or/and validity? Is that necessary? I think I may not consider Cronbach's alpha in that case, since its not a survey approach, and often the measurement scales are different. Another option I was thinking of is CFA, but SPSS have only EFA. Would setting an EFA with defined number of factors equal to the number of my latent variables is an option? 2) How I should test (accept/reject) my hypothesis? Should I refer to the measures themselves, or should I refer to the latent variables? How I can create the latent variables if necessary?


title: Online model training in R, without a static data

text: Let's say I have a model in R, a regression tree created by "glm",using "data1" dataframe:               Model1 = glm(DepVar ~ . ,data=data1,family="binomial")      Is there a way to train the model with additional new data that will be coming, say data2 , without having data1 ?In other words, is there a way to create a new model Model2 that is based only on Model1 and data2, without keeping data1?