title: Machine learning curse of dimensionality explained?

text: I'm having trouble understanding the curse of dimensionality. Specifically, I came across it while doing the `scikit-learn` tutorial in python. Can someone please explain the below in a simpler manner? Sorry I have been trying to understand for the longest time and can't understand how they came up with the calculation for number of training examples to achieve an efficient KNN estimator? Here is the explanation: > For an estimator to be effective, you need the distance between neighboring > points to be less than some value d, which depends on the problem. In one > dimension, this requires on average n ~ 1/d points. In the context of the > above KNN example, if the data is described by just one feature with values > ranging from 0 to 1 and with n training observations, then new data will be > no further away than 1/n. Therefore, the nearest neighbor decision rule will > be efficient as soon as 1/n is small compared to the scale of between-class > feature variations. > > If the number of features is p, you now require n ~ 1/d^p points. Let’s say > that we require 10 points in one dimension: Now 10^p points are required in > p dimensions to pave the [0, 1] space. As p becomes large, the number of > training points required for a good estimator grows exponentially. link here EDIT: also is the tilde (`~`) supposed to represent approximate in that example? or the python tilde operator?


title: What is wrong with teleoanalysis (if anything)?

text: I came across the term 'teleoanalysis' in relation to medical trials. There is no Wikipedia article on it yet, but Google search has a few links to scientific papers and their critiques. As far as I understand, one of the main concepts of teleoanalysis is: if A causes B, and B causes C, then we can conclude that A causes C. Here is a quote from the teleoanalysis founding paper: > _‘It may also be necessary to quantify the individual effects that relate to > separate steps in a causal pathway–that is, the effect of factor A on > disease C is determined from the estimate of the effect of A on an > intermediate factor B and the estimate of the effect of B on C, rather than > by directly measuring the effect of A on C. The exercise is like putting > together the pieces in a jigsaw puzzle.’_ One of the criticism mentions that this A->B->C logic is deployed when randomised trials either failed to prove relationship between A and C or are impossible to conduct, but there is some evidence of relationship between A and B as well as B and C. The objection is that this logic is unscientific. I was wondering if there is anything wrong with A->B->C reasoning from experiment design and general statistical science point of view?


title: Mixed model with 1 observation per level

text: I'm fitting a random effects model with `glmer` to some business data. The aim is to analyse sales performance by distributor, taking into account regional variation. I have the following variables:   * `distcode`: distributor ID, with about 800 levels   * `region`: top-level geographical ID (north, south, east, west)   * `zone`: mid-level geography nested within `region`, about 30 levels in all   * `territory`: low-level geography nested within `zone`, about 150 levels Each distributor operates in only one territory. The tricky part is that this is summarised data, with one data point per distributor. So I have 800 data points and I'm trying to fit (at least) 800 parameters albeit in a regularised fashion. I've fitted a model as follows:               glmer(ninv ~ 1 + (1|region/zone/territory) + (1|distcode), family=poisson)      This runs without a problem, although it does print a note: > Number of levels of a grouping factor for the random effects is _equal_ to > n, the number of observations Is this a sensible thing to do? I get finite estimates of all the coefficients, and the AIC also isn't unreasonable. If I try a poisson GLMM with the identity link, the AIC is much worse so the log link is at least a good starting point. If I plot the fitted values vs the response, I get what is essentially a perfect fit, which I guess is because I have one data point per distributor. Is that reasonable, or am I doing something completely silly? This is using data for one month. I can get data for multiple months and get some replication that way, but I'd have to add new terms for month-to-month variation and possible interactions, correct? * * * ETA: I ran the above model again, but without a `family` argument (so just a gaussian LMM rather than a GLMM). Now `lmer` gave me the following error: > Error in (function (fr, FL, start, REML, verbose) : Number of levels of a > grouping factor for the random effects must be less than the number of > observations So I'd guess that I'm not doing something sensible, as changing the family shouldn't have an effect. But the question now is, why did it work in the first place?


title: Citation for the weakness of sequential univariate regression?

text: I'm rebutting a paper in which the author has done linear regression in two steps:   1. Regression against a single predictor   2. Regression of the residuals from the first model against a second predictor The problem with this approach is that if the two predictors are correlated (they are) then it tends to force more of the variance into the first predictor. The results thus depend arbitrarily on the order of the two steps. In this case if the steps were reversed, or if a proper multivariate regression were performed, then the author would have reached the opposite conclusion. I've demonstrated all of this, but it would be nice to have a good citation for the problem too - it must be well know. Thanks in advance for your suggestions.


title: Chi-squared distribution as sum expressed

text: Let $X_1, \ldots, X_n$ be i.i.d. exponentially distributed random variables with density $$\eqalign{\theta^{-1} e^{-x/\theta}, &x \ge 0 \\ 0, &x \lt 0} $$ and let $Y_i = X_{(i)}$ denote the order statistics such that $Y_1 \leq \cdots \leq Y_n$. How to show that $$ 2\frac{\left(\sum_{i=1}^{r}Y_i\right) + (n-r)Y_r}{\theta} $$ has a chi-square distribution with $2r$ degrees of freedom? I wrote the the joint density of $(Y_1,Y_2,...,Y_r)$ but nothing became apparent.


title: How to predict rank in ordinal regression using R?

text: I have done an ordinal regression using `lrm` from the `rms` package.               load(url("http://www.produnis.de/R/NEXT-Q0-OR.RData"))     require(rms)     mydata <- NEXT.or     fit <- lrm(QuitJob~QD,data=mydata);fit      I would like to predict the rank in "QuitJob" (1=won't quit ... 5=will quit very soon) by giving a value for "QD" (how much work to do). For example, if "QD" is 78, what rank of "QuitJob" is likely? If I type in               predict.lrm(fit, data.frame(QD=78))      I get                       1      0.3191933       Do I understand it right, that the probability of being in rank 1 is 31,92% ?! If so, how can I get the probabilities for ranks 2-5 ?


title: difference between hierarchical dirichlet process and nested dirichlet process

text: There have some extensions to Dirichlet process. One is Hierarchical Dirichlet process, and another is Nested Dirichlet Process. What are the differences between these two? I once read the paper of Bayesian Nonparametric Inference -- Why and How. The author discusses their differenes as follows, but I am not clear how to understand this. Any more explanations is highly appreciated. ![enter image description here](http://i.stack.imgur.com/U9tzc.jpg)


title: Alpha correction for two-way ANOVA (and post-hoc) with unequal variance?

text: I have a multifactorial design (3 treatments x 6 time points, n = 6). One sample per subject. Because my sample size is small, I've used the KS test for normality - all groups pass. Still, I'm not confident I can assume normality due to small n. That said, my main issue (I believe) is lack of equal variance. The differences of SD between groups is substantial (up to 13-14x greater), all tests (Levene's, Brown-Forsythe, Bartlett's, etc) show very small p values, indicating unequal variance. Log transformation doesn't help. Variance does not follow any obvious trend with respect to the size of the mean. Is a normal two-way ANOVA followed by a multiple comparison test using an adjusted significance level (to be more conservative) in each analysis an appropriate way around this? I'm using JMP and GraphPad Prism - I'm not looking for sophisticated alternatives to the parametric two-way ANOVA or post-hoc to deal with unequal variance - as far as I know, no such alternative tests even exist for normal multi-factorial designs (unmatched, non-pair wise). Shall I run the ANOVA and post-hoc anyway, and report the unequal variance violation? If so, with or without an adjustment to alpha? Thanks, Statsnewb


title: Similarity Coefficient and Geochemical Correlation

text: I have a set of around 160 samples from different sediment bodies, or stratigraphic units, which I am attempting to correlate using geochemical analysis of 25 elements for each sample, and a similarity coefficient often applied in correlating volcanic ashes. The literature gives no mention of assumptions of normality for the coefficient. It is assumed that each sample is representative of the variability in that stratigraphic unit, as representative sections were sampled while in the field. To do these analyses, the sample is split to a smaller fraction that is loaded in to a plastic capsule, which contains the sample for the XRF irradiation. Some workers in my department who have done XRF on solid objects such as larger rock fragments, insist that 25-30 sub-samples must be analyzed per sample, So in other words, for each of my 160 samples, I would prepare 25-30 splits of the sample for XRF analysis. These ideas have propagated through the department. They are however using MDFA with its assumption of multivariate normality. Also, I am measuring sands, silts and clays, plus smaller granule gravels, of which hundreds if not thousands are visible to the XRF beam while in the capsules. Being that there does not seem to be any assumption of normality with the Correlation Coefficient being used, is there any basis for this idea that I would need to prepare 25-30 splits of each sample in my case? I apologize if this has been explained poorly.


title: How to modify post-hoc power analysis simulation to perform a-priori power analysis?

text: I am interested in doing a sample size calculation for exacerbation rates. I am drawing my data from a negative binomial distribution and my test statistic is based on a negative binomial glm. Currently, I have a code for a power calculation for a given sample size. How do I modify this code to calculate sample size? The following are the input parameters for the sample size calculation:   * Alpha level = 0.04,    * target power = 90%,    * event rate for reference group (r0) = `mu0 <\- 0.8/52*48`,    * Rate ratio (rr) under `ha = mu1/mu0` where `mu1 <\- mu0*0.6`,    * average exposure time = 48 weeks,    * sample size allocation ratio = 1 (so each group has equal sample size; i.e., `n1 = n0`), and    * dispersion parameter (`k = 0.4`).               library(MASS)     n     <- 400        # sample size     alpha <- 0.04     int   <- 1000     k     <- 1/0.4      # shape     mu0   <- 0.8/52*50  # control exacerbation rate     mu1   <- mu0*0.4    # exacerbation rate for trt1     mu2   <- mu1        # exacerbation rate for trt2     t0    <- rep(0, n)     t1    <- rep(1, n)     ta1   <- c(t0, t1)     ta2   <- c(t0, t1)     c1    <- NULL     c2    <- c1     test1 <- c1     test2 <- c1          for (i in 1:int){       pbo <- rnegbin(n, theta=k, mu=mu0)       d1  <- rnegbin(n, theta=k, mu=mu1)       d2  <- rnegbin(n, theta=k, mu=mu2)       r1  <- c(pbo, d1)       r2  <- c(pbo, d2)            x1  <- cbind(ta1, r1)       x2  <- cbind(ta2, r2)            dt1 <- as.data.frame(x1)       dt2 <- as.data.frame(x2)       fm1 <- glm.nb(r1 ~ ta1, data = dt1)       fm2 <- glm.nb(r2 ~ ta2, data = dt2)            test1 <- c(test1, summary(fm1)$coef[2,4]) #$       test2 <- c(test2, summary(fm2)$coef[2,4]) #$       c1    <- c(c1,    summary(fm1)$coef[2,3]) #$       c2    <- c(c2,    summary(fm2)$coef[2,3]) #$     }          #power     result <- (sum( (test1<alpha)&(test2<alpha) )/int)      result


title: Homogeneity of Variance, 2 Way Completely Crossed Design

text: I am looking for some advice as to how I might handle having an unfavorable Levene's test outcome, that is to say a highly significant value. **DOE** Various containers with water were microwaved for 30 seconds under various conditions, the response variable the max temperature. The experiment run order was randomized. There were 4 levels for container types; ceramic, plastic, paper, and foam. 3 Levels of microwave turntable conditions, dirty, wet, and clean. 12 treatment combinations, 4 replicates per treatment, i.e. a completely crossed design. The goal is to examine how the water temp is impacted by these 2 factors. ![enter image description here](http://i.stack.imgur.com/vc4jI.png) As one can see from this interaction plot, there is a pattern when it comes to microwave cleanliness, clean microwaves preform the best while wet ones preform the worst. Unfortunately, I ran a Levene's Test on the treatment combinations, and it comes out to be highly significant, which means the assumption of homogeneity is broken, and ANOVA is inappropriate. I believe this is caused by the greater deal of variation in the ceramic temps. I know Welch t-tests on all pairwise comparisons is an option, but then I loose the ability to say the different levels of microwave cleanliness are significantly different (which seems like a crime given the interaction plot). I also thought about doing a transformation of the response variable, but a box-cox shows that none of the possible transformations are particularly intuitive. ![enter image description here](http://i.stack.imgur.com/uGOod.png) Any ideas about how to handle this homogeneity issue I would be much appreciated.


title: Question on MLE

text: I apologize for the basic question. If $\{p_\theta(x): \theta\in K\subseteq\mathbb{R}\}$ is a smooth family of distributions, then the MLE $\hat{\theta}_n,$ under suitable regularity conditions satisfies $\sqrt{n}(\hat{\theta}_n-\theta)\to\mathcal{N}(0,I(\theta)^{-1})$ where $I(\theta)$ is the Fisher information. We also know from the Cramer-Rao bound that no estimator can outperform this estimator in the sense of having a smaller mean squared error (except on a set of zero Lebesgue measure, the superefficiency phenomenon). I am interested in knowing similar results for other families of distributions which do not satisfy the above conditions. The canonical examples are: $\{\mathrm{Uniform}[\theta,\theta+1]: \theta\in\mathbb{R}\}$ and $\{\mathrm{Uniform}[0,\theta]: \theta>0\}.$ In these cases, the MLE estimates the parameter to within $O(1/n),$ since in the former case, $\hat{\theta}_n = \min\{X_1,X_2,\ldots, X_n\}$ and $n(\theta-\hat{\theta}_n)\to -Z,$ where $Z\sim\mathrm{exp}(1)$ and in the latter case, $\hat{\theta}_n = \max\{X_1,X_2,\ldots, X_n\}$ and $n(\theta-\hat{\theta}_n)\to W,$ where $W\sim\mathrm{exp}(\frac{1}{\theta}).$ Do we know (1) If MLE is always an order-optimal estimator in the sense that it estimates $\theta$ to $O(\frac{1}{n^\alpha})$ where $\alpha$ is the best possible over all estimators? (2) If MLE is not just order-optimal but optimal in a suitable sense similar to that of the Cramer-Rao lower bound, in that, any other estimator can outperform MLE only on a set of Lebesgue measure 0? (3) Any general results on the distribution of convergence of $n^\alpha(\theta-\hat{\theta}_n)$ where $O\left(\frac{1}{n^\alpha}\right)$ is the rate of the convergence of the estimator?


title: What's the correct way to simulate a division winner?

text: So I'm trying to simulate one winner for each of the 4 NFL divisions, **assuming that all teams are of equal ability**. It's clear that the probability of each team being the division winner is $1/4$. I wonder if the following simulation strategy works. Instead of simulating each division, I randomly rank all 16 teams from 1 to 16, then compare them in group of 4. I did run the code, and the probability does converge to .25, but it takes 100,000 simulations to even get kinda close. Thus, I'm unsure whether the proposed simulation strategy is valid or not?               one_simu <- function() {       rank <- sample(1:16, replace=F)       groups <- split(rank, ceiling(seq_along(rank)/4)) # Split into 4 groups of 4       result <- rep(0, 16)       for (i in seq_along(rank)) {         group_of_i <- groups[[ceiling(i/4)] ] # The group that team 1 belongs to         if (rank[i] == max(group_of_i)) {           result[i] <- result[i] + 1           }       }        return(result)     }          apply(replicate(1000, one_simu()), 1, mean)     apply(replicate(100000, one_simu()), 1, mean)


title: How can I code in SAS a tree to bin continuos variables?

text: I don't want to use SAS-Miner, because I want to create the code to do the optimal binning, even if it's only using entrophy criteria. Is there a code to do so? Thanks.


title: oversampling glm in R: how to define `weights`

text: Let's say I have 10 positives out of 1000 observations. I'd like to run `glm` on the 10 positives and a sample of 10 non-positives (so a total of 20 records in the dataset going into the analysis). To maintain ratios, the weighing scheme is set to 1 for positives and 99 for non-positives, and `weights = c(rep(1,10),rep(99,10))` if the 10 positives come first in the data, followed by the 10 sampled non-positives:                glm( is_positive ~ . , data = D_sampled , family = binomial , weights = c(rep(1,10),rep(99,10)) )      Is this the right way to run an oversampled logistic with 'glm' (keeping in mind that the intercept will have to be corrected)?


title: updating clusters in NCut spectral clustering

text: Are there any methods to update the clusters formed in k-way NCut when new data points arrive or any change in the similarity matrix?? thanks


title: Conditional variance model with external regressors

text: I need to estimate the following conditional variance model $$ y_t = \sigma_t\,\varepsilon_t $$ where $y_t$ are the observed data, $\varepsilon_t$ are iid Normal shocks (zero-mean and unit-variance) and where the variance $\sigma_t$ is given by $$ \sigma_t^2 = \omega+\sum_{j=1}^m\zeta_j\,v_{j,t}+\sum_{j=1}^{q}\,\alpha_j\,\varepsilon_{t-j}^2+\sum_{j=1}^{p}\beta_{j}\,\sigma_{t-j}^2, $$ where $v_{j,t}$ with $j=1,...,m$ are $m$ external observed variables. I am asking whether I can find some efficient numerical routines, especially with Matlab/Octave.


I understand the difference between partial correlation and semi-partial and zero-order correlation. This terminology makes sense to me in that the partial correlation of X and Y partials out the effect of one or more variables from both X and Y, and semi-partial correlation only partials out the effect of one or more variables on either X or Y (i.e., it only partially partials; i.e., it semi-partials). For example, SPSS labels the semi-partial correlation as part correlation. This is the source of some confusion to students encountering this output for the first time.   * **Why is the semi-partial correlation sometimes called the "part correlation"?**   * **What is the etymology of the term "part correlation"?**


title: Not very known but effective data mining algorithms?

text: I recently came a cross Matrix Factorization for classification or recommendations. I was very shocked that it was one of the most effective techniques in the Netflix competition. But actually I think very few people have heard about it. Other techniques like SVM, ANN ... etc are very popular. I was wondering if there are other machine learning algorithms that are not popular or very known but showed far better performance than other popular techniques in some domains or problems.


title: Convergence of batch gradient descent in logistic regression

text: I am not really sure about how it behaves when using batch gradient descent in logistic regression. As we do each iteration, $L(W)$ is getting bigger and bigger, it will jump across the largest point and $L(W)$ is going down. How do I know it without computing $L(W)$ but only knowing old $w$ vector and updated $w$ vector? If I use regularized logistic regression, will the weights become smaller and smaller or any other patterns?


title: CART - Classification And Regression Trees

text: I try to prune one regression tree build with the rpart function in R. To decide where to prune the tree I used the plotcp function. But I noticed that if I use the same predictor variables and in the same order the plotcp graph always change? How is this possible? Thank you for your explanation in advance.


title: Is it possible to find the combined standard deviation?

text: Suppose I have 2 sets: **Set A** : number of items $n= 10$, $\mu = 2.4$ , $\sigma = 0.8$ **Set B** : number of items $n= 5$, $\mu = 2$, $\sigma = 1.2$ I can find the combined mean ($\mu$) easily, but how am I supposed to find the combined standard deviation?


title: Zero-inflation on steroids: choose among Poisson, negative binomial and zero-inflated regressions

text: I am struggling to fit alternative count models into my data. I guess my problem is just too many zeros. This is my data               > summary(smpl)         response        predict1          predict2              Min.   :0.000   Min.   :   1.00   Min.   :    22005        1st Qu.:0.000   1st Qu.:   3.00   1st Qu.:  4669705        Median :0.000   Median :   8.00   Median : 12540318        Mean   :0.017   Mean   :  23.27   Mean   : 20382574        3rd Qu.:0.000   3rd Qu.:  20.00   3rd Qu.: 25468156        Max.   :3.000   Max.   :1584.00   Max.   :145348049          > table(smpl$response)       0   1   2   3      987  10   2   1       I tried three regressions: basic Poisson, negative binomial and zero-inflated but the only formula returning coefficients without warnings is the Poisson:               > summary(glm(response ~ ., data = smpl, family = poisson))          Call:     glm(formula = response ~ ., family = poisson, data = smpl)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.3871  -0.2214  -0.1722  -0.1148   4.7861            Coefficients:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -3.472e+00  3.521e-01  -9.862  < 2e-16 ***     predict1     3.229e-03  7.271e-04   4.442 8.93e-06 ***     predict2    -6.258e-08  3.060e-08  -2.045   0.0409 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for poisson family taken to be 1)              Null deviance: 150.67  on 999  degrees of freedom     Residual deviance: 135.84  on 997  degrees of freedom     AIC: 170.06          Number of Fisher Scoring iterations: 8      The negative binomial returns a warnings on both the convergence and the alternation limit               summary(glm.nb(response ~ ., data = smpl))          Call:     glm.nb(formula = response ~ ., data = smpl, init.theta = 0.04901296596,          link = log)          Deviance Residuals:           Min        1Q    Median        3Q       Max       -1.28844  -0.17677  -0.14542  -0.09808   2.38314            Coefficients:                   Estimate Std. Error z value Pr(>|z|)         (Intercept) -3.899e+00  4.587e-01  -8.499  < 2e-16 ***     predict1     1.226e-02  2.144e-03   5.720 1.06e-08 ***     predict2    -5.982e-08  3.407e-08  -1.756   0.0791 .       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Negative Binomial(0.049) family taken to be 1)              Null deviance: 69.927  on 999  degrees of freedom     Residual deviance: 55.940  on 997  degrees of freedom     AIC: 152.37          Number of Fisher Scoring iterations: 1                             Theta:  0.0490                Std. Err.:  0.0251      Warning while fitting theta: alternation limit reached            2 x log-likelihood:  -144.3700      Warning messages:     1: glm.fit: algorithm did not converge      2: In glm.nb(response ~ ., data = smpl) : alternation limit reached      and the zero-inflated (from the `pscl` package) doesn't return anything at all               > summary(zeroinfl(response ~ ., data = smpl, dist = "negbin"))          Call:     zeroinfl(formula = response ~ ., data = smpl, dist = "negbin")          Pearson residuals:          Min       1Q   Median       3Q      Max      -0.45252 -0.08817 -0.05515 -0.04210 19.56118           Count model coefficients (negbin with log link):                   Estimate Std. Error z value Pr(>|z|)     (Intercept) -1.477e+00         NA      NA       NA     predict1     2.678e-03         NA      NA       NA     predict2    -1.160e-07         NA      NA       NA     Log(theta)  -1.241e+00         NA      NA       NA          Zero-inflation model coefficients (binomial with logit link):                   Estimate Std. Error z value Pr(>|z|)     (Intercept)  4.869e+00         NA      NA       NA     predict1    -1.329e-01         NA      NA       NA     predict2    -1.346e-07         NA      NA       NA     Error in if (getOption("show.signif.stars") & any(rbind(x$coefficients$count,  :        missing value where TRUE/FALSE needed      Then my questions are:   1. Is there anything I can do in terms of "formula tweaking" with the negative binomial (to avoid the warnings) and with the zero-inflated (to get the coefficients)?   2. Looking only at the results above (thus including problems with convergence and alternation limit) should I select the negative binomial model since it seems, looking at the AIC, to fit better than the Poisson in my data?


title: Is there an analytic formula for the kurtosis of a (student t) mixture distribution?

text: The mixture distribution should be composed of noncentral scaled student t components. The mean of the mixture distribution can be calculated easily just by weighting the component means according to the mixing law. I also could derieve a tractable analytic formula for the mixture variance by means of the "law of total variance" in a straightforward way. But for the kurtosis I found no low-level "shortcuts" to get a general analytic formula in this case. Does anyone have seen such a formula in a book? Any suggestions for a straightforward derivation of such a formula would also be very helping! Thanks a lot in advance!! Jo


title: Why features compression is good?

text: I'm reading about deep learning and that in principles it's a features compression technique and that is why it works. Now my question is why compressing features from 200 or so into 4 is better? How and why does it influence classifiers performance or the construction of the decision boundary? Look at the following picture for example: ![enter image description here](http://i.stack.imgur.com/BMpiP.png) Assume that blue points are class-0 and red points are class-1. A classifier might see that any instance that has an $x_2$ lower than $0.5$ is of class-0 and any instance that has $x_2$ above $0.5$ is class-1. At this time $x_1$ is redundant. But now look at this example: ![enter image description here](http://i.stack.imgur.com/P7zdm.png) In this example $x_1$ becomes vital for a good performance. So why features compression is better? Or could (if not always better) lead to a better performance?


title: Basic question about the filtration in the martingale formed by sum of iid random variables

text: It is known that if $X_i$ are iid, $E(|X_i|) < \infty$ and $E(X_i) = 0$ then $S_n = \sum_1^n X_i$ is a martingale. Suppose all $X_i$ are defined wrt sample space $\Omega$. I don't understand why $S_n$ is $\sigma(X_1, .. X_n)$-measurable? Specifically, I am confused about the sample space for $S_n$. To be $\sigma(X_1, .. X_n)$-measurable the sample space of $S_n$ should be $\Omega$. I do not see why. Any help would be appreciated.


title: Silhouette scores for different distance metrics

text: I clustered a data set using PAM with a euclidean distance metric and a pearson correlation distance metric. The average silhouette value of the correlation clusters is higher at most points than the euclidean one, but their maximas are at the same number of clusters, and peaks occur at similar cluster numbers. Is there any way to compare the two on the same scale?


title: Is it true that values of parameters do not vary from sample to sample?

text: This is something I was told. However, it seems to me that parameters of a population can be functions of time, in certain situations, in which case the values of parameters _could_ vary from sample to sample, if the parameter varies with time, unless the samples are all taken at the same time. Or, am I missing something? Are populations with time-varying parameters ill- defined, in some way? If so, what are the criteria that makes a population 'well-defined'? Thanks in advance.


title: Training error of SVM in small n large p case

text: There are n data points in R^p space, where p>n. We perform binary classification using linear SVM. As the VC dimension of linear SVM is p+1, it is capable to shatter n points with any label arrangement. Does that mean we should get a 0 training error? However, I rarely get such training error from my experience working with high dimensional data. Does that suggest an under- fitting?


title: Is a linear mixed model a special case of Gaussian process regression?

text: Are there any introductory references or resources that discuss this connection if there is one? I have not been able to find any ones suitable for beginners.


title: Large discrepancy between confidence intervals for correlation coefficient using bootstrap and Fisher's Z transform

text: I found that the default confidence intervals provided by cor.test (Fisher's Z transform according to the help) are quite different from the nonparametric boostrap confidence intervals for the below data. Is it normal this large discrepancy? I am doing something wrong?               dat<-data.frame(x=1:5,y=c(2.4,2.6,3.5,7,8.1))     cor.test(dat$x,dat$y)$conf.int     bootCor<-boot(dat,function(d,i) return(cor(d[i,1],d[i,2])),R=1000)     quantile(bootCor$t,c(.025,.975),na.rm=T)


title: Finding outliers on simple dataset?

text: I'm pretty inexperienced in statistics and am wondering the best way to tackle a couple of simple questions. I am given the following data on 1,000 hypothetical slot machines: Attempts made, wins, payout, payout per attempt (win rate * payout) Some slot machines have over 15,000 attempts, while others have as few as 50. Here is what I'm wondering:   1. What is the best way assess which machines have payout per attempt values that can be classified as outliers? Seeing as some machines have so few attempts.   2. How can I determine the statistical significance of a given payout per attempt value when comparing it to a specific number? For example if as a gambler I want to target machines with higher than X payout per attempts? I am working with these values in Excel also. Thank you very much!


title: Where can I use kernels other than Gaussian (like Cauchy, laplacian) in kernel methods in machine learning? Or maybe in kernel density estimation?

text: In few papers I read that - kernel used doesn't really matter for kernel density estimation but bandwidth of the kernel is the most important factor. But I did not see any mathematical explanation to this and can't completely believe in this.


title: Is minimum-k a distributive measure or an algebraic measure?

text: Minimum-k (ie, finding the K smallest elements of a dataset) seems like it should be a distributive measure. If I know the smallest k items for two separate datasets, I can compute the smallest k items for the merged dataset by just merging the two sub-results. However, some lecture slides categorized minimum-k as algebraic, the reasoning being that k numbers need to be sent from each sub-division to compute the global min-k, as opposed to something like max, where each sub-division just reports a single number (it's own max). This distinction doesn't seem particularly relevant though -- the data being sent is the smallest k items, not some collection of "sub-values" used to compute the smallest k items (as in the case of average being computed via sum and count).


title: Non-negative Matrix Factorization - Basic Question: Cluster Assignment

text: Given the result $V_{m\times n} \approx W_{m\times k} \cdot H_{k\times n}$, where columns of $V$ are data points and $m$ is data dimension, what is the function by which you assign the data points to the $k$ classes? Is the usual method to use the normalized columns of $W$ as centroids? It seems like that would be problematic since even a normalized version of $W$ is not unique. For example if $R$ is a rotation matrix $WR\cdot R^{-1}H = W \cdot H$ where $WR$ is still normalized. Thanks so much...


title: Question about EFFECT SIZE

text: I have a question. If I found the effects of a variable is insignificant (p>.05), but the effect size d is 0.2 (small effect). So I used GPower to calculate the power and I found it is 0.14, which means I need a larger sample size to get a significant results. So is the hypothesis rejected or supported? because it will be significant if i have a larger sample size. Thanks Let me make it clear. My hypothesis is that **priming improves creative performance.** So if p>.05, d=0, the hypothesis is rejected, **priming has no effect on creative performance**. My question is if p>.05, d=0.2 (small effect), is the hypothesis rejected or supported?


title: How did Gallup form this index?

text: I am reading this Gallup poll, and in it there is a "Positive Experience Index" that is used without really explaining how it was created, and therefore what, precisely, it means. From reading the first paragraph, it appears that five yes/no questions were asked, and the index was created from this data. Among other things, it isn't clear whether Gallup weighted any of these responses in its index. My questions are the following: (1) How could I find out how this index was created? (2) Is there a a general method that is regularly used by statisticians in such a situation that I just don't know about, and this is why I'm ignorant? (My background is not in statistics.) (3) If the answers to the previous questions are null, did Gallup just do a poor job presenting its findings? I.e. _What's going on here?_


title: How to capture a dataset characteristic?

text: I am working on microaggregation problem, where all records are clustered in groups of a minimum size. We must optimize the sum of the squared error (SSE). As I run my algorithm on forest and vlcensus (very large census) datasets, the results for vlcensus are promising, while for the forest they are not satisfying. If there is any way I can determine what is the main differences or characteristics for the datasets (for example if vlcensus is more uniform than the forest, or if it is more sparse or clustered than the forest, ...). Forest is 581012x10 and vlCensus is 149642x13 numerical datasets. Regards


title: If $X_{n+1}$ is a martingale subject to $Y_0,\ldots,Y_n$, then is it a martingale with respect to $Y_0^2,\ldots,Y_n^2$?

text: I don't have a very solid foundation in measure theory, and this always seems a bit confusing to me so I would appreciate any help. We are given $ E \left( X_{n+1} | Y_0,\ldots,Y_n \right) = X_n. $ Prove or disprove $ E \left( X_{n+1} | Y_0^2,\ldots,Y_n^2 \right) = X_n $ I am thinking that if $F=\sigma \left(Y_0,\ldots,Y_n \right)$ and $G=\sigma \left(Y_0^2,\ldots,Y_n^2 \right)$, I need to prove that F = G? Is this correct? Then I can do something like this: $ E \left( X_{n+1} | G \right) = E \left( E \left( X_{n+1} | F \right) | G \right) = E \left( X_{n+1} | F \right) = X_{n+1} $. Also is $E \left( X_{n+1}^2 | G \right) = X_n^2$ (a martingale) given $ E \left( X_{n+1} | Y_0,\ldots,Y_n \right) = X_n. $


title: Follow up a significant interaction in ANOVA with dichotomous variables

text: Hopefully this is an easily answerable question, though I've been searching for an answer without any luck. I've run a 2x2x2x2 mixed-design ANOVA with two within-subjects and two between-subjects factors. I'm trying to follow up the only significant interaction, which is between the two between-subjects factors. These are scary image vs. non-scary image, and high-fear vs. low-fear. The interaction is significant because the low-fear group barely differ between the scary and non-scary images, but the high-fear group have both a lower score than the low-fear group for the non-scary image, and a higher score than the low-fear group for the scary image. I would like to know whether each of these differences is significant, that is, whether the difference between fear groups is significant at each level of the image variable, but I'm not sure how to test it. Many thanks.


title: (Spline) curve fitting $(x,y)$ points when $dy/dx$ is also known for each point

text: I have an $(x,y)$ dataset consisting of points on an unknown curve and I've been using spline fits to generate a curve joining these points as a guide-to- the-eye (and nothing more). However, I also know $dy/dx$ at each of my points, and so I know that the spline fit is not ideal: e.g. I know that there's a maxima at one point and a minima at another, and the spline fit is misleading at these points. What kind of curve fitting should I use if I want to also incorporate this extra information? I have access to Matlab, Mathematica, R and anything open source. Suggestions for even just a place to start reading would be much appreciated. Thank you. Edit: the only thing that I have thought of is to manually make additional points either side of my "real" points at $x+dx$ and $x-dx$, and then use this new set of points to generate the spline, but this seems a little dubious to me.


title: How to calculate most congested roads from vehicle counts

text: I'm studying traffic flows, and I'd like to know what are the busiest/most congested roads. From this website http://www.dft.gov.uk/traffic- counts/about.php, it looks like vehicle counts are typically measured in AADF which is the Annual Average Daily Flow - the average number of vehicles that will travel between two checkpoints on an average day. However, these AADFs are measured on road sections of varying lengths. The website seems to think that "Traffic" ie AADF x road length (measured in vehicle miles) is a useful measurement, but what is this number actually useful for? If this is the case, a very sparse but long road would have the same value as a busy short road. It doesn't seem to make any sense using this definition of "traffic". Surely, to figure out how congested a road is, you need to _divide_ the AADF by the road length. Can anybody confirm this? Or make it clear how the AADF _is_ useful? Thank you.


title: Conditional cross-validation sampling

text: When applying cross-validation to evaluate the predictive performance of a binary classification model, **is it acceptable to separately sample from cases and non-cases to achieve class proportions in the training and testing partitions that closely match those in the original data set**? My application is a Cox model (yes this is not a binary classifier). My data consist of ~500 non-cases and ~100 cases. When I partition the full data randomly for ten-fold cross-validation, there is high variance in the number of cases that make it to the test partition. In the populations for which I'd like to make good predictions, however, the incidence of the disease I am studying is pretty similar to that in my data. I therefore feel as if my conventional cross-validation scheme provides estimates of the mean predictive performance over many parallel universes in which incidences vary hugely. I have a feel that the answer should depend on the performance metrics used. What I use is net reclassification improvement, integrated discrimination improvement, and difference in C Index, to compare two models. All of these are sensitive to class proportions. Should I sample separately? Should I not? And why?


title: Learning parameters of non-parametric Bayesian models

text: I have a sample of Chinese restaurant process which I want to model as Pitman–Yor process. How do I determine parameters of Pitman-Yor model from given sample? For Dirichlet process I would just use the fact that on average K=αlogN where is K is observed number of tables and N is total number of visitors but I see no similar way for the Pitman-Yor process.


title: Understanding Sequential Probability Ratio Test (SPRT) Likelihood Ratio

text: I am a software developer looking to develop an alternative for the simple hypothesis testing scheme described here. In short, the test works as follows:   * Two URLs are compared for their ability to _convert_ visitors.   * Discrete samples are captured. Each of them either _converts_ or fails to do so. A conversion is a desirable event.   * A higher _conversion rate_ is seen as desirable. I am looking into SPRT because while the users of my testing system are comfortable setting significance level and power ahead of time, they are not comfortable setting sample size ahead of time. They want a test with a stopping rule rather than one that has to capture a specific sample size. I understand most of the math in the linked Wikipedia article. However, I don't understand how to compute the likelihood ratio. Could someone provide a concrete example of how to compute the likelihood ratio and $S_{i}$ for the following events during the example test:   * 0.95 significance level   * 0.80 power   * Two URLs: URL1 and URL2 Events:   1. Visitor arrives at URL1. Visitor converts.   2. Visitor arrives at URL1. Visitor fails to convert.   3. Visitor arrives at URL2. Visitor fails to convert.   4. Visitor arrives at URL1. Visitor fails to covnert.   5. Visitor arrives at URL2. Visitor converts.   6. Visitor arrives at URL2. Visitor converts. If you could write out the likelihood ratio and $S_{i}$ at each step I would appreciate it!


title: Multiple imputation of time variables -- which step to impute?

text: Lets assume I have a survival analysis study with an exposure, two covariates, and two time related variables. Say date of diagnosis and date of death. Combined, the two time related variables will be used to make a survival time variable: `t = date of death - date of diagnosis `. The problem is that some of those date variables are missing. That should be addressable by multiple imputation, but I'm not sure _which_ part I should be imputing - the dates used to make t, or t itself, based on the values of the other covariates. My concern with using the date is the influence that will have on the imputed values - that if there's a small number of people diagnosed on say May 8th that they'll have a hugely disproportionate influence on the imputed values vs. other demographic characteristics for someone diagnosed May 8th but missing their outcome date. If it helps, the outcome in question isn't seasonal or time-related, so time _shouldn't_ - at least to my mind - have a heavy influence.


title: Negative Binomial Regression and Heteroskedasticity

text: What relationship does a negative binomial regression have to heteroskedasticity and if one still needs to check and/or correct for it how would this be done?


title: Good algorithm for processing positional estimates

text: We (my team) are building a robot which will navigate around an arena. The robot uses a camera to determine its position based on markers on the wall. We have tested this and found it can determine position well when close to the wall. However, at further distances, the markers appear to move randomly by single pixels on the sensor (due to vibration, sensor noise and other factors), causing the positional estimates to vary considerably. Initially we took a mean of the captured samples, but if we have a list of samples like this:               (4.45, 5.06)     (4.43, 5.07)     (4.42, 6.48)     (3.95, 5.67)      then the outliers considerably change the averaged result, making the calculated position inaccurate. So, we need to remove these data points - either replace the values with more reasonable estimates, or mark them as outliers and ignore them while taking the mean. Ideally, any such algorithm should be capable of handing only a few data points, as we will only be able to see a maximum of 7 markers at any one time, and typically only 4-5. I took GCSE Statistics about 3 years ago (and got an A) and I vaguely remember an algorithm which did just this - but my knowledge has since faded.


title: Least squares logistic regression

text: I have seen it claimed in Hosmer & Lemeshow (and elsewhere) that least squares parameter estimation in logistic regression is suboptimal (does not lead to a minimum variance unbiased estimator). Does anyone know other books/articles that show/discuss this explicitly? Google and my personal library have not helped me here...


title: How to obtained composite index using factor analysis on likert scale data

text: I am analyzing the influence of social capital on household income. I have collected data through a questionnaire survey for my study.Different dimensions (variables) of social capital were measured using 5-point multiple likert items.The scores on each item was summed up to obtain the single value for each variable.Now i have 6 such variables and I want to perform factor analysis on these 6 variables to obtain a single social capital index for each household.My question is that can I summed up the factor scores (obtained through factor analysis) to obtain social capital index? if my approach is not appropriate, please suggest me some alternatives. I have to use this social capital index in multiple regression for further analysis.


title: How to interpret type I (sequential) ANOVA and MANOVA?

text: My primary question is how to interpret the output (coefficients, F, P) when conducting a Type I (sequential) ANOVA? My specific research problem is a bit more complex, so I will break my example into parts. First, if I am interested in the effect of spider density (X1) on say plant growth (Y1) and I planted seedlings in enclosures and manipulated spider density, then I can analyze the data with a simple ANOVA or linear regression. Then it wouldn't matter if I used Type I, II, or III Sum of Squares (SS) for my ANOVA. In my case, I have 4 replicates of 5 density levels, so I can use density as a factor or as a continuous variable. In this case, I prefer to interpret it as a continuous independent (predictor) variable. In R I might run the following:               lm1 <- lm(y1 ~ density, data = Ena)     summary(lm1)     anova(lm1)      Running the anova function will make sense for comparison later hopefully, so please ignore the oddness of it here. The output is:               Response: y1               Df  Sum Sq Mean Sq F value  Pr(>F)       density    1 0.48357 0.48357  3.4279 0.08058 .     Residuals 18 2.53920 0.14107       Now, let's say I suspect that the starting level of inorganic nitrogen in the soil, which I couldn't control, may have also significantly affected the plant growth. I'm not particularly interested in this effect but would like to potentially account for the variation it causes. Really, my primary interest is in the effects of spider density (hypothesis: increased spider density causes increased plant growth - presumably through reduction of herbivorous insects but I'm only testing the effect not the mechanism). I could add the effect of inorganic N to my analysis. For the sake of my question, let's pretend that I test the interaction density*inorganicN and it's non-significant so I remove it from the analysis and run the following main effects:               > lm2 <- lm(y1 ~ density + inorganicN, data = Ena)     > anova(lm2)     Analysis of Variance Table          Response: y1                Df  Sum Sq Mean Sq F value  Pr(>F)       density     1 0.48357 0.48357  3.4113 0.08223 .     inorganicN  1 0.12936 0.12936  0.9126 0.35282       Residuals  17 2.40983 0.14175       Now, it makes a difference whether I use Type I or Type II SS (I know some people object to the terms Type I & II etc. but given the popularity of SAS it's easy short-hand). R anova{stats} uses Type I by default. I can calculate the type II SS, F, and P for density by reversing the order of my main effects or I can use Dr. John Fox's "car" package (companion to applied regression). I prefer the latter method since it is easier for more complex problems.               library(car)     Anova(lm2)                 Sum Sq Df F value  Pr(>F)       density    0.58425  1  4.1216 0.05829 .     inorganicN 0.12936  1  0.9126 0.35282       Residuals  2.40983 17        My understanding is that type II hypotheses would be, "There is no linear effect of x1 on y1 given the effect of (holding constant?) x2" and the same for x2 given x1. I guess this is where I get confused. **What is the hypothesis being tested by the ANOVA using the type I (sequential) method above compared to the hypothesis using the type II method?** In reality, my data is a bit more complex because I measured numerous metrics of plant growth as well as nutrient dynamics and litter decomposition. My actual analysis is something like:               Y <- cbind(y1 + y2 + y3 + y4 + y5)     # Type II     mlm1 <- lm(Y ~ density + nitrate + Npred, data = Ena)     Manova(mlm1)          Type II MANOVA Tests: Pillai test statistic             Df test stat approx F num Df den Df  Pr(>F)         density  1   0.34397        1      5     12 0.34269         nitrate  1   0.99994    40337      5     12 < 2e-16 ***     Npred    1   0.65582        5      5     12 0.01445 *                # Type I     maov1 <- manova(Y ~ density + nitrate + Npred, data = Ena)     summary(maov1)                    Df  Pillai approx F num Df den Df  Pr(>F)         density    1 0.99950     4762      5     12 < 2e-16 ***     nitrate    1 0.99995    46248      5     12 < 2e-16 ***     Npred      1 0.65582        5      5     12 0.01445 *       Residuals 16                                                 Sorry this has become ridiculously long for a fairly basic question. Any insight would be greatly appreciate. Thank you and happy New Year! Dan Hocking


title: How do I use the standard regression assumptions to prove that $\hat{\sigma}^2$ is an unbiased estimator of $\sigma^2$?

text: I'm working through an econometrics textbook and it's proving that $$ \sigma^2 = E(\hat{\sigma}^2) = \frac{SSR}{n-2} $$ I followed the proof (an example of which is shown on talkstats) until it reaches this step: \begin{align} E\left[\displaystyle\sum\limits_{i=1}^n \hat{u}_{i}^{2}\right] &= E\left[\displaystyle\sum\limits_{i=1}^n (u_{i}-\bar{u})^{2}\right] + (\hat{\beta}_{1}-\beta_{1})^{2} E\left[\displaystyle\sum\limits_{i=1}^n (x_{i} - \bar{x})^2\right] \\ &+ 2 (\hat{\beta}_{1}-\beta _{1}) E\left[\displaystyle\sum\limits_{i=1}^n (u_{i}-\bar{u})(x_{i}-\bar{x})\right] \\ &= \dots \end{align} The book says that:   1. The expected value of the first term should be $(n-1)\sigma^2$   2. The expected value of the second term should be $\sigma^2$.   3. The expected value of the third term should be $2 \sigma^2$ because the term itself can be written as $2(\hat{\beta}_1-\beta_1)^2 \hat{\sigma}_x^2$. For the first term, I can work backwards and get that: \begin{align} \hat{\sigma}^2 &= \frac{1}{n-1} \displaystyle\sum\limits_{i=1}^n (u_{i}-\bar{u})^{2} \\ &= \sigma^2 \end{align} but it seems like this only works if we _assume_ that $\hat{\sigma}^2$ is an unbiased estimator. For the second term: \begin{align} (\hat{\beta}_{1}-\beta_{1})^{2} E\left[\displaystyle\sum\limits_{i=1}^n (x_{i} - \bar{x})^2\right] &= (\hat{\beta}_{1}-\beta_{1})^{2} E\left[(n-1) \hat{\sigma}_x^2\right] \\ &= (\hat{\beta}_{1}-\beta_{1})^{2} (n - 1) E\left[\hat{\sigma}_x^2\right] \\ &= (\hat{\beta}_{1}-\beta_{1})^{2} (n - 1) \hat{\sigma}_x^2 \end{align} but this doesn't match what the text shows, so I'm a bit confused as to the final step for this term. For the third term: \begin{align} 2 (\hat{\beta}_{1}-\beta _{1}) E\left[\displaystyle\sum\limits_{i=1}^n (u_{i}-\bar{u})(x_{i}-\bar{x})\right] &= 2 (\hat{\beta}_{1}-\beta _{1}) E\left[\displaystyle\sum\limits_{i=1}^n u_{i}(x_{i}-\bar{x})\right] \\ &= 2 (\hat{\beta}_{1}-\beta_{1}) \displaystyle\sum\limits_{i=1}^n E\left[u_{i}(x_{i}-\bar{x})\right] \\ &= 2 (\hat{\beta}_{1}-\beta_{1}) \displaystyle\sum\limits_{i=1}^n E\left[u_{i}\right]E\left[(x_{i}-\bar{x})\right] \\ \end{align} I would have thought that because $E[u_i] = 0$, this term would cancel out, so I'm not sure how to derive $2\sigma^2$ from it. Can someone give me a hint in the right direction for these steps? I'm firmly trying to work through all of the proofs myself with only a little guidance from the book, but I'm not quite there yet. I noticed this question, which links to this question, so maybe with more reading I'll see where to apply that. (If that isn't it, hints are welcome, but I'll still work through those questions to see if it helps)


title: Outlier detection in time series data

text: I checked different questions on similar topics, but none were exactly the answer I wanted and I am confused. I am working with big data, the data has a bursty nature with high frequency. I considered features one by one with respect to time (equivalent to one time- series) and want to remove outliers in selected time series. I am implementing this in java using weka. I read a lot about this problem but I did not find which exact method would work best for any time series. It would be also great, if any optimal outlier method could also smooth time series because at the end of day, I need to give multiple time series data to PCA for finding correlation. As you know, PCA is sensitive to outlier, noisy nature and missing values. I know, I cannot find all three things in one algorithm but for me, outlier detection is tough part.   Please give your views on it.   Thanks


title: How to combine variances from sensors where each observation has its own variance?

text: I have a set of measurements $x_1$ ... $x_n$. These measurements are normally distributed, measuring the same value. However due to the way the data is measured, each $x$ has its own standard deviation: $s_1$ ... $s_n$. In other words I have a sensor which returns pairs (x,s). Now, I wan to estimate parameters of distribution of $x$ using $\text{N}$ samples. The common mean would be just a sample mean. What about variance? I could not apply "pooled variance" formula because I do not know how many samples were used to estimate each $s$.   Can I just use a mean of $s$ as sample variance? **Update:** I can not make any assumptions on how sensor produced (x,s) values. They might be based on some hidden iterations or perhaps something else. However it is safe to assume that each (x,s) pair returned to me is independent from others and measuring the same true value of 'x'.


title: How can I determine the best fit normal distribution from this information?

text: I'm trying to take a normal distribution of points, and force them to become a uniform distribution. I've had little success on S.O., so I thought I'd ask here. Basically, I have a hash function which takes an `X`, `Y`, and seed value and generates numbers in a normal distribution with a mean of 0. These numbers vary only slightly from the numbers on all sides of them. When mapped out with the `X` and `Y` and the hashed value being `Z`, it creates a terrain map. The map has some large peaks far past 1 and -1, but the middle 50% of the values lie between (-0.4,0.4). I'm trying to smooth this map out so that it retains its shape, but has a more-or-less uniform (rather than normal) distribution. As I said, the middle 50% of the values lie between (-0.4,0.4). The theoretical bounds for the hash function is (-2.25,2.25), though after generating a billion samples, the range of numbers found were about (-1.75,1.75). I think I need to take the above information to determine the best fit normal distribution, then use it to transform each value. As I found out on S.O., it's a really tough problem for me to explain. Hopefully someone here can at least point me in the right direction explaining it, or understands what I'm trying to do. * * * The probability densities of my distribution, in blue, and `Normal(0,.72)`, in red: [-3.8,3.8] ![alt text](http://i.stack.imgur.com/AAGtA.png) * * * Their cumulative probabilities: [-3.8,3.8] ![alt text](http://i.stack.imgur.com/7euTm.png) * * * The probability density after inverse probability transform: [0,1] ![alt text](http://i.stack.imgur.com/cJyta.png) * * * And the cumulative probability compared to `Uniform`: [0,1] ![alt text](http://i.stack.imgur.com/CHiQo.png)


title: Analytical solutions to limits of correlation stress testing

text: Suppose I have some initial correlation matrix. I want to stress each correlation in the matrix by the same constant simultaneously (except the diagonal; lets call this global parallel stress since it affects the entire matrix by the same constant at the same time). I start with adding 0.01% to each correlation and check if the matrix is still PSD. I continue increasing the stress levels by small increments. Eventually I encounter a stress level above which the matrix would no longer be PSD. Let's call this the upper stress boundary. I repeat the same procedure for negative stresses and find the lower stress boundary. My general observation was that the minimum eigenvalue initial correlation matrix is roughly equal to the upper stress boundary. **Is there an analytical solution for this problem?** Generally, I am more interested in the upper stress boundary. I stress the entire matrix at the same time because it is a simple way to reduce diversification benefit across all market variables. However, from a theoretical point of view it would be interesting to find out how to analytically calculate both lower and upper stress boundaries. EDIT: please see a previous thread, which might be relevant here (some brilliant idea's from kwak): http://stats.stackexchange.com/questions/2572/correlation-stress-testing


title: How to explain such a big difference between parametric and non parametric test (and other questions)?

text: I am advising a small medical study (two groups, treatment is a dummy variable), i.e. a 2x2 contingency table. I'm comparing the value of the Pearson's $\chi^2$ test and a non parametric competitor McNemar's $\chi^2$ test. **Edit** The answer brings an additional questions: they have matched each case in group1 with 4 controls in group2 (matched according to all the variables they deem important, such as sex, age, hemisf) except treatment. Setting aside the question of whether the matching is well done (i.e. whether all important variables have indeed been identified), doesn't the fact that to each case corresponds **four** controls (i.e. not one) artificially inflate the significance of the results? (they still have to send me part of the data, this is why the table below does not reflect this 4 to 1 ratio). I obtain the following (very different) results (n=116) (R CrossTable() function).:                    [,1] [,2]     [1,]   39    9     [2,]   49   19      Statistics for All Table Factors ## Pearson's Chi-squared test with Yates' continuity correction Chi^2 = 0.844691 d.f. = 1 p = 0.3580586 ## McNemar's Chi-squared test with continuity correction Chi^2 = 26.22414 d.f. = 1 p = 3.039988e-07 ## Fisher's Exact Test for Count Data Sample estimate odds ratio: 1.672924 Alternative hypothesis: true odds ratio is not equal to 1 p = 0.279528 95% confidence interval: 0.6356085 4.692326 The McNemar is the approximate version, but the exact version gives the same conclusions (strong rejection of the null). My question is: how can i understand such a large difference between $\chi^2$ and McNemar ?


title: How to combine confidence intervals for a variance component of a mixed-effects model when using multiple imputation

text: The logic of multiple imputation (MI) is to impute the missing values not once but several (typically M=5) times, resulting in M completed datasets. The M completed datasets are then analyzed with complete-data methods upon which the M estimates and their standard errors are combined using Rubin's formulas to obtain the "overall" estimate and its standard error. Great so far, but i'm not sure how to apply this recipe when variance components of a mixed-effects model are concerned. The sampling distribution of a variance component is asymmetrical - therefore the corresponding confidence interval can't be given in the typical "estimate ± 1.96*se(estimate)" form. For this reason the R packages lme4 and nlme don't even provide the standard errors of the variance components, but only provide confidence intervals. We can therefore perform MI on a dataset and then get M confidence intervals per variance component after fitting the same mixed-effect model on the M completed datasets. The question is how to combine these M intervals into one "overall" confidence interval. I guess this should be possible - the authors of an article (yucel & demirtas (2010) Impact of non-normal random effects on inference by MI) seem to have done it, but they don't explain exactly how. Any tips would be much obliged! Cheers, Rok


title: Choosing N when doing t-tests on bootstrapped samples

text: Specifically the Welch's T-Test, but probably any T-Test, requires a value for Mean, Variance and N to be used to calculate the T-statistic and the degrees of freedom. I am concerned about these values if the sample has been bootstrapped (re-sampled with replacement) before using the T-Test. Specifically, I'm concerned that the arbitrary choice of how many times to bootstrap (10K, 20K, 100K) will change the T-statistic even though the underlying distribution is the same. For instance, if I have 100 sample values with an unknown distribution and I re-sample with replacement 10,000 times, then I use those 10,000 samples to calculate a mean, variance, and N for a t-test. **N:** Should I use the sample size _before_ the bootstrap or _after_ the bootstrap? I'm concerned that my choice of how many times to run the bootstrap (10,000 is an arbitrary number) will make the distribution look more accurate than it is. Should I use n=10,000 or N=100? _Full disclosure: I have a fairly good knowledge of statistics, enough to use the terms, but perhaps not enough to ensure I have used them 100% correctly. Please excuse confusions of terminology._


title: What is a good academic citation for cross-validation?

text: I performed a holdout cross-validation analysis on a multilevel model fit. The purpose of this was to show that we didn't have a problem with over-fitting, for which it worked just fine. Now we are writing it up for publication and I need a citation to support my methodology. I am looking for a good canonical statistical reference, ideally a book, that does a nice job explaining why holdout cross-validation is a real thing that people do and makes sense in this application. The paper will be published in a biological journal, so I am looking to point non-statistical types to a general reference. Somehow, none of my books seem to quite do it. The wikipedia entry would be perfectly adequate for my purpose, but I'd rather not cite wikipedia. Any suggestions?


title: Continuous and Categorical variable data analysis

text: I have three variables:   * distance (continuous, variable range negative infinity to positive infinity)   * isLand (discrete categorical/ Boolean, variable range 1 or 0)   * occupants (discrete categorical, variable range 0-7)  I want to answer the following statistical questions:   * How to I compare distributions that have both categorical and continuous variable. For example, I like to determine if the data distribution of distance vs occupants varies depending on the value of isLand.   * Given two of the three variables, can I predict the third using some equation?   * How can I determine independence with more than two variables?


title: Identifying outlier data in high-dimensional settings

text: I have a data set with high-dimensional feature space. Are there any pre- processing methodologies that can detect outliers from this data set? The outlier, I mean, are the ones that tend to be very different with other data points and can be the ones that generated by improper measurements.


title: Crossvalidation and/or testdata. Always use both or can one exclude the other?

text: I'm trying to build a two class classifier on a dataset of around 570 samples. Im evaluating several classificiation stratigies (LDA, QDA,RDA, logistic, logistic with some additional ellements like splines ...) I have difficulties if I should take out a 25% of my 570 samples as a test set or if should just rely on 10 fold crossvalidation to estimate my test error. I know this subject has been discussed before here on stackexhange but Im still confused when I read te answers. What I'm thinking right know is that I have to take out 25 %, keep it aside and build my models using the other 75% of the data (so trainingdata). Also perform 10 fold crossvalidation on this trainingdata and use this cross validation error to decide between several classification methods/models. After selecting a final model, do a last check with the test data (that has not been used otherwise) to do a last performance check of my data. I can not use this testdata to compare classification methods. Am I correct here? best regards


title: Computing the likelihood for multiplicative error in the independent variable

text: A problem has recently arisen for me which involves estimating relative strength of various mechanisms contributing to an overall quantity. These strength parameters $q_j$, $j\in[M]$ stays fixed between observations of the quantity $\kappa_i$, $i\in[N]$, but the weighting factors $p_{ij}$ change between observations. These weights are not observed directly, but are simulated. The dominant source of error in the whole process is in the simulations of the weight, and unfortunately the error here is not uniform AWGN, but "in the exponent", i.e. $p_{ij}=e^{\mu_{ij} + \delta_{ij}}$. Accordingly I think I have the following model \[\kappa_i = \sum_{j=1}^M q_j p_{ij} \epsilon_{ij}\] \[\log(\epsilon_{ij})\sim \mathcal{N}(0, \sigma^2) \] Given that I want to find an 'optimal' assignment for $\vec{q}$, I write the likelihood function \[\log\mathcal{L}(\vec{q}) = \sum_{i=1}^N \log\left(\int\bigg(\prod_{j=1}^M d\epsilon_{ij} P(\epsilon_{ij})\bigg)\delta\left(\kappa_i-\Sigma_jq_jp_{ij}\epsilon_{ij}\right)\right) \] Unfortunately, my attempts at putting this in a useful form haven't yielded results, and I can't seem to find any (readable to a non-professional statistician) references for this type of problem. I think it might be a good start to replace the log-normal with a more tractable gamma distribution. Any suggestions? I would like to know if there is perhaps some way of transforming this problem so it is amenable to usual techniques. I don't need an analytic solution, but without a likelihood function I'm not sure how I'd even do it numerically.


title: How to approach missing data...ignore or not to ignore?

text: What would be the best option when you come across missing data. Do you exclude that person completely from the data analysis or not. If not then how do you go about it? i have data on two groups at three different time points. At each time point they had to undertake a 6 minute walking test. Now some participants they stopped at stage 4 (because physically they could not go any longer) and others managed to complete all of the stages (1-6th stage)..some managed to complete only 2 stages and so on. Thank you in advance


title: Hierarchiqual prediction using R

text: I'm pretty new in R, and I couldn't find any information about a package who can do the following: supposing that I have a set of data (for instance, different text documents), which can have several classes. For example, a datum could be a Sport, a Sport with Ball, a Sport without Ball and a Car. I'd like to be able to predict to which category the data belongs, considering that I may not hit that the datum is a Sport with Ball, but I'd be happy if I correctly predict that it's a Sport. Which package can provide this kind of stuff? Some examples would be useful, if possible. Thanks in advance


title: Sample mean of independent, normally-distributed variables

text: My book gives the following lemma: if $X_1, X_2, \dots, X_n$ are independent normally distributed random variables, such that each $X_i$ has mean $\mu$ and variance $\sigma^2$, then the mean of $\bar{x}$ is normally distributed with mean $\mu$ and variance $\sigma^2$ over $n$. Following this, I attempted a problem where each $X_i$ had mean 25 and standard deviation 1. I had to find the probability that the mean of $X_1, X_2, X_3, X_4$ was greater than 26. I ended up with the probability being zero - I do $P(X < 26) = P(Z < (26-25)/0.25) = P(Z < 4)$. I think I'm doing something wrong somewhere. Can someone point out where I'm going wrong please?


title: How to prepare interactions of categorical variables in scikit-learn?

text: **What is the best way to prepare interactions of categorical features before fitting with scikit-learn?** With `statsmodels` I could conveniently say in R-style `smf.ols(formula = 'depvar ~ C(var1)*C(var2)', data=df).fit()` (same in Stata with `regress depvar i.var1##i.var2`). Can `sklearn.preprocessing.PolynomialFeatures` (in v0.15, currently dev) be used with categorical variables?


title: Proof of Point-Biserial Correlation being a special case of Pearson Correlation

text: I have been examining the use of the Point Biserial correlation as a statistic to measure the relationship between a dichotomous variable and a continuous one. Wikipedia et. al. seem to concur that the Point Biserial Correlation is a special case of the Pearson Correlation, but I cannot find a proof for this, algebraic or otherwise, and it is making me wary of using this in the context of the research I am doing (I need to do some statistical confidence testing afterwards). I have tried deriving the truth myself, but have chased everything round in a circle. Any advice greatly appreciated.


title: Extracting variables chosen by glmnet or lars packages in r

text: I'm dealing with a large p small n problem (p>1000, n=150) and have decided to use glmnet and lars to select possible explanatory variables (as I understand _is_ their functionality) rather than do this manually. I have gotten out both a lars, and a glmnet object from the functions, but now can't discern how to extract from those objects, the names of the coefficients / the optimal model, that the functions have chosen. Any help would be greatly appreciated!


title: Comparing Multivariate Means

text: Apologies in advance if this is very basic; I am quite new to statistics. I have a large set of 7 dimensional data, for two groups. In my example, they are express relative preferences over a certain bundle of goods (for men, and for women). (Weightings of the parameters in a representative utility function, so for example X1=<2,3,1,4,5,2,7>). As such, one could interpret them as either two sets of multivariate data, or equivalently, two sets of 7-dimensional directional unit vectors, without any loss of relevant information. I want to see if the populations have different means. First question, should I try to use a multivariate t test, or delve into the literature on directional statistics? (the latter seems to represent the intuition behind the data structure better) Secondly, can you recommend R packages that might be able to do either of the above? The 'circular' statistics package, as far as I can tell, does not extend to N dimensional vectors. Thirdly, and I'm throwing this out there more if anyone could indicate which area of literature to examine, suppose the sets of data represent ordered variables? And supposing they were paired instead of independent (say representing twin 1 and twin 2)? What would be the best topic to investigate to help answer questions such as, for example, does a skew in the 'twin 1' data towards the first few ('lower order' variables) correlate with a similar skew towards the first few in the corresponding 'twin 2' data? Or if twin 1 highly prefers 'higher order' goods twin 2 prefers 'lower order' goods? And so on. Thanks in advance!               EDIT   (adding information given in comment)      So a little bit of context - the data arises from a simulation exercise where users are given a budget of say £10 and a menu of prices for different items. From each 'market' we assume their choice is their optimal bundle and we use a generic utility function = $U = \alpha(\ln (x_1))+\beta(\ln (x_2) +\dots$. For $k=7$ items to compute their relative weights for each item. That is, the solution can be represented by the 7D unit vector; magnitude is irrelevant. I'm trying to see if men and women differ systematically in their preferences.


title: How to get/calculate confidence intervals for slopes of the continuous*categorical interactions?

text: Not sure if this should be in Stack Exchange or here, but here goes. I have a linear mixed model in SPSS, and I want the confidence intervals for slopes of the continuous*categorical interactions. Below is the output from the analysis, with estimates of fixed effects and the estimated marginal means. I want to report the effects of my continuous variable (here ZlnSCLd) in each category (period). SPSS prints the slopes in the estimates table as relative to the reference category, so here, the slope is 10.49 for period 205, and the slopes for the rest of the periods are in relation to that - i.e., for period 201, the slope is ZlnSCLd+interaction = 10.49-20.18 = -9.69. This can be seen in the two estimated marginal means, the first with ZlnSCLd set to 1, the second with the default zero, where the difference of the 201-periods in the two tables equals -9.688. ![enter image description here](http://i.stack.imgur.com/ekwOT.png) So I know how to get the slopes, but what about the CIs? The slope for the ref category has the CI ready, but if I report the slope for period 201 as -9.69, I cannot report the CI for it as [-41.84, 1.48] that I have for the relative estimate -20.18. Can I simply add the estimate of the ref category also in the CIs so that [-41.8, 1.48] becomes [-31.3, 11.97], or do I have to take into account the CI of the reference category somehow (because the model - and the relative estimate - depends on the ref cat)? (On a side note, why is the CI shown by SPSS not +-1.96*SE? It is close, but not quite: [-41.62131513, 1.258071195]. This makes me hesitant to calculate it myself.)


title: How would you calculate the probability of matching the value of a sequenced list?

text: Suppose we have a sequence of items, from an alphabet of 2: {0,1}. For example, this 8-bit sequence: 10001010 You generate a random 8-bit string. The probability of you matching the first bit is $1/2$ or 50%, the last bit also $1/2$ or 50%, and of matching both: $1/2 * 1/2$ or 25%. Part 1 concerns how many attempts someone would have to make to have the first and last bit match.   1. Intuitively, it seems like 4 attempts would yield a high probability of achieving this match. But statistics is not always intuitive... Part 2 concerns any particular bit in the sequence, but trying to match the most number possible.   1. What is the probability of matching _any_ bit correctly? (Of matching any $n$ bits correctly?)    2. How many random generations would someone need to make to reach $n$ bit matches with probability=$1/2$? (With probability=$3/4$?)


title: Kalman filter observation model

text: I am using the following equations for Kalman filter. The state vector is $x$ and observations are $z$. Both $x$ and $z$ are multivariate vectors of same length. Forecast model: $x_{t+1} = 0.963 x_t + \mathcal{N}(2.19,3.34)$ Observation model: $z_t = -0.00542849 x_t + \mathcal{N}(0.4524,0.038952)$ These relations hold between state and observation vectors element by element (since the lengths are same). The observations are confined within an interval [0, 1] whereas $x$ values may range from [0,80]. I am not sure if the noise values I have determined are correct. I just fit a linear regression between $z$ and $x$ and used the MSE to find the standard deviation of noise. Is that correct? Somehow, the update step introduces a lot of noise. So much that I am in fact better off not updating at all. How does one estimate the noise covariance matrices given a set of $(x_{t+1},x_t)$ and $(z_t,x_t)$ values? Is my methods of fitting a linear regression and using $\sqrt{MSE}$ as standard deviation of noise correct? EDIT: I should mention here that the relationship between $z$ and $x$ seems linear element-by-element and that's why I am learning a regression model.


title: Prediction with not atomic features

text: I would like to use another type of data, not atomic data, as a feature for a prediction. Suppose I have a Table with those Features: - Column 1: Categorical - House - Column 2: Numerical - 23.22 - Column 3: A Vector - [ 12, 22, 32 ] - Column 4: A Tree - [ [ 2323, 2323 ],[2323, 2323] , [ Boolean, Categorical ] ] - Column 5: A List [ 122, Boolean ] I would like to predict/classify ... Columns 2 ... for example.... I am making a Software to automatically respond questions... Any type...like "Where Foo was Born ?" ... I first make a query to a search engine --->>> then I get some Text data as a Result. So I do all the Parsing Staff... Tagging, Stemming, Parsing, Splitting... My first approach was to make a table, each row with a line of text.. and a lot of Features...like ... First Word ... Tag of First Word.. Chunks, etc.. But with this approach I am missing the relationships between the Sentences. I would like to know if there is an algorithm that look inside the Tree Structures... Vectors... and make the relations and extract whatever is relevant for predicting/classifying. I rather know a library that does that then an algorithm that I have to implement... Thank you very much !


title: Is there a version of multivariate multinomial logit?

text: I'm working with a data set with 2-3 response variables and 7 predictor variables. All the variables are categorical. If there were just one response variable, I think a multinomial logit would be the right model, but there are 2 or 3. So my question is - is there a multivariate version of the multinomial logit? I've looked at several books on categorical data, but haven't seen anything like this (mainly using Agresti 2002). I have about 2000 observations, though I'll probably need to split it up into 2 or 3 data subsets to really see what's going on. One thing I was thinking about is converting it to counts and use a model for count data. I could also combine the 2-3 response vars into 1 categorical with a lot of categories, but I think that will lower the chances of anything showing up for any of the categories. I could also do 2-3 separate models, one for each variable, which is obviously not as good. I might also be able to get rid of some of the predictors (I think 3 of the 7 have the most explanatory power). I'm not opposed to using machine learning methods, I've found some interesting stuff already with decision trees. thanks, -paul


title: Experimental design considerations for Randomized Block design

text: Is there a use of using a Randomized Block design when you only have one treatment and you only want to test the effects of this single treatment? The experiment in question is to see whether a training exercise (the treatment) effects the grades of a number of students. Students are grouped in a control and experimental group and are subject to a test. The experimental group are then given training and then asked to re take the exam.


title: Is there an R equivalent for the SPSS custom Tables feature?

text: I have used the custom tables feature in SPSS for my job almost every day. I am trying to learn R . How can I recreate the custom tables functionality in R. Is there some specific package that I need to use?


title: Is there an easy way to combine two glm models in R?

text: I have two logistic regression models in R made with `glm()`. They both use the same variables, but were made using different subsets of a matrix. Is there an easy way to get an average model which gives the means of the coefficients and then use this with the predict() function? [ sorry if this type of question should be posted on a programming site let me know and I'll post it there ] Thanks


title: How to increase the proportion of explained variance in PCA?

text: When I was performing a PCA on a large data set, I only get the PC1 explaining about 25% of variance and PC2 about 20%, PC3 about 15%, PC4 about 10%... So I wonder if there is a way to increase the Proportion of Variance Explained in PCA.


title: Big Data Use Cases

text: I am going to study about big data use cases in the real world. I want the stuff to be preferable in book format(not article or blog posts or website articles) and high level, industry based and vendor neutral. Up to now I came to this cases which are concise books up to 100 pages that describe how Industry can solve its problems with big data analytics. Can anyone recommend some similar items?   * How Data Science Is Transforming Health Care   * The Information Technology Fix for Health the first two available here and Book: Big Data - The Future of Education


title: Best option to cluster variables (not cases) in R

text: In SPSS, if I use the hierarchical clustering procedure, I have the ability to cluster both variables and cases using a variety of methods and distance measures. For this task, I would like to use R to cluster my variables. For context, my data come from a survey and the respondents were able to select multiple items from a block of options. In my datafile the data are coded as 0,1. I am trying to learn R so any help you can provide will be greatly appreciated. Best regards, Brock


title: How do I deduce the SD from regression and ANOVA tables?

text: ![enter image description here](http://i.stack.imgur.com/GTB1c.png) This is a Minitab printout. I want to find the value of A5, or S. I think S is supposed to be the sample standard deviation, but I don't know how to calculate it. Any tips on how I should go about calculating it?


title: How to check for bivariate Gaussianity without the use of regression?

text: What steps could be taken to check for bivariate Gaussianity without using regression based check? Can we somehow employ the use of definition of variogram measure for assessing spatial variability?


title: How do I calculate single variable confidence intervals assuming different distributions with R?

text: I'm trying to convert a bunch of Stata commands to their R equivalents, and I'm struggling with how R does handle confidence intervals for inferential statistics of single variables. In Stata, I can use `ci variable` to calculate normal confidence intervals, `ci variable, b` to calculate binomial intervals, and `ci variable, p` to calculate the intervals assuming the variable is distributed as Poisson. R seems slightly more complicated, though (`confint()` only works with models…). It appears that I have to run the respective tests to get the confidence intervals (`t.test`, `binom.test`, and `poisson.test`). I'm fine doing that, but these functions seem to be less generalized than Stata's `ci` command. For example, to calculate a normal 95% confidence interval I use `t.test`:               set.seed(1234)     x <- rnorm(100)     t.test(x)          # 95 percent confidence interval:     #  -0.35605755  0.04253406      However, it seems to be a lot more complicated to calcuate binomial or Poisson intervals. For example, if I have a column of binary data (yes/no; 1/0) like `x` below, Stata appears to convert it into count data automatically when running `ci x, b`. In R, I have to convert it to count data on my own with `sum()`:               set.seed(1234)     x <- sample(0:1, 100, replace=TRUE)  # lots of 0s and 1s     binom.test(sum(x), length(x))          # 95 percent confidence interval:     # 0.3503202 0.5527198       Is manually feeding in the sum and the length of the variable the official R-esque way of calculating the confidence interval assuming a binomial distribution, or is there a better way? Likewise, what's the most R-esque way to calculate confidence intervals for a variable assumed to follow a Poisson distribution--the equivalent of `ci x, p`:               set.seed(1234)     x <- rpois(100, 5)  # random Poisson distribution     # ... magic R voodoo ...     # Confidence interval!      So, I guess in summary, what's the best R equivalent for Stata's `ci x`, `ci x, b`, and `ci x, p` commands for single variables?


title: What's a good normality test for time-series / sleep data?

text: I have data that is, at its most basic level, activity counts per minute of fruit flies. Flies are considered to be asleep when they have no activity for at least 5 minutes. The data essentially looks somewhat like a sinusoidal curve when plotted as activity vs. time or sleep vs. time. I have researched various normality tests and have so far been using Shapiro-Wilks (SW), Kolmogorov-Smirnov (KS), and the Anderson Darling (AD) test. I've ended up using the AD test mostly because its supposed to be a better form of the KS test (which I have seen used in the literature). If it matters, the biological reason for non-normality is that these flies are mutants and some of them tend to be grossly hyperactive. I hesitate to use a normality test just because it's commonly seen in the literature. So, are there better tests than the ones I've used so far? And are there specific reasons to use them with time-series / sleep data?


title: Graphical nominal model

text: Suppose I have a set of $k$ matrices. $$ \mathbb D = A_1,A_2,...,A_k $$ Each column of $A$ is categorical vector. $$ A = v_1,v_2,...,v_n $$ I want to find the mapping $$ f: A \rightarrow Q $$ Where $Q$ is a $(n,n)$ matrix such that $$ \frac{1}{k} \sum_k \text{MSE}\left(Q_{k,real}(ij) - Q_{k,predict}(ij)\right) $$ is minimized. MSE is chosen because $Q(ij) \in 0,1$ representing some continuous event occuring. Thus far several distance metrics $d(A(i),A(j))$ have been directly used to form the elements of $Q(i,j)$. They are not parameterized in anyway by the bins of the frequency vector or bins of the cross-tabs of the $v$s. Chi- squared gave bad results. I settled on symmetric uncertainty(information theory). Forming $Q_{predict}$ it easy to show that the error can be greatly minimized by Z-scoring it assuming a bivariate Gaussian distribution, however the results are less than spectacular. (This feels wrong and there is still a lot of error.) I'm looking for a well studied supervised machine-learning/fitted-model for the entire process, accounting for the graphical structure and cross- tabulation data. So far I've heard suggestions from computer science people on Belief Nets, recurrent neural networks, structured SVM and multiple correspondence analysis. Multiple correspondence analysis is the only method though expect nominal data, but it doesn't fit anything nor model the joint PMF. Ideally I would like the model to be interpretable and not a black- box(like SVM). I'm much better at linear algebra and programming than I am at statistics. I've thought about building a nested model and solving it with something like stochastic gradient descent. I'm willing to go experimental on any suggestions. I have an abundance of CPU cycles and training data but no good model.


title: Should estimates over a subsample be similar to the full sample when using panel data?

text: I am estimating the effect of $x$ on $y$; $y$ and $x$ are available both by city and by year. During the first period, $T_1$, both $x$ and $y$ generally increase. During the second period, $T_2$, both generally decrease. I have about 30 cities over 40 years.   1. How can I use regression to attribute the average effect of $x$ on $y$ to the overall average increase in $y$ during $T_1$ and the same for the decrease in $T_2$? For example, I want to say something like during period $T_1$, $x$ can explain 5% of the increase in $y$ over $T_1$. How would I compute $x$? Normally, I divide the coefficient by the average difference between the last year and the first year, but that approach does not sound sensible if $x$ and $y$ both increase and then decrease.   2. I find statistically significant results when I do the regression over the entire time period. However, if I break it separately to just subperiods, $T_1$ and $T_2$, neither are statistically significant and the coefficients change sign. Is that simply because the smaller sample makes the standard errors blow up? Should I find effects of $x$ and $y$ over subsamples of the panel?


title: Orthogonalized regression reference?

text: I'm not sure what this is called but I remember seeing a colleague of mine doing a multivariate regression much like $$Y \sim X_1 + X_2 + X_3$$ and then he said he would "orthogonalize" so he regressed each of the IVs against the other ones, and took the residuals..as in: $$Z_1 = residuals(X_1 \sim X_2 + X_3)$$ $$Z_2 = residuals(X_2 \sim X_1 + X_3)$$ $$Z_3 = residuals(X_3 \sim X_1 + X_2)$$ after which he redid the regression: $$Y \sim Z_1 + Z_2 + Z_3$$ I have never seen this before, is there a reference or name under which I can search read more about this and why/when would one want to do something like this?


title: Modification of error function and weight update rule of backpropagation algorithm

text: ![enter image description here](http://i.stack.imgur.com/kKvo3.jpg) I'm trying to understand the neural networks with a given example. On the graph, observed instances are shown. It states that the function is symmetric around x=3.5. So the question is : How would you modify the error function and the weight update rule of backprogation algorithm? The confusing point is does it whether asks me start computing step by step by forward propagation then calculate errors and update the weights or there is a tricky part that I'm missing. Because there is no information about the number of hidden layers etc. is given. When using a neural network, I have to use the information below: ![enter image description here](http://i.stack.imgur.com/uTaTD.jpg)


title: Most efficient order to learn LaTeX, Sweave, Beamer?

text: I am very interested in learning how to create reoccurring reports from my R code and ggplot2 visualization. I understand that LaTeX seems to be a possible answer, and to use it with R most use Sweave. And for presentations from R to LaTeX people use Beamer. My question is, what I should learn first, or what should I learn concurrently? I don't know Sweave or LaTeX. Should I learn LaTeX (at least a bit), then learn Sweave? Or would you suggest learning them at the same time? Links to tutorials supporting your answers greatly appreciated.


title: How does linear discriminant analysis reduce the dimensions?

text: There are words from "The Elements of Statistical Learning" on page 91: > The K centroids in p-dimensional input space span at most K-1 dimensional > subspace, and if p is much larger than K, this will be a considerable drop > in dimension. I have two questions:   1. Why the K centroids in p-dimensional input space span at most K-1 dimensional subspace?   2. How are the K centroids located? There is no explain in the book and I didn't find the answer from related papers.


title: CART analysis with multiple dependent variables

text: I want to do a CART analysis with **multiple dependent variables.** Which program is able to do that?


title: Can I run a regression where dependents are coefficients from other regression and independents are R squared from these regressions?

text: I calculate Beta risk for multiple assets by running multiple regressions: Return = intercept + Beta*Market_Return Then I want to compare the Betas. I run other regression where Beta is dependent variable against several independent variables: Beta = intercept + X + Y... May I include R squared from first regression as an independent variable in the second regression? Sort of variable that describes how well assets returns could be explained by market returns? Beta = intercept + X + Y + Rsqrt...


title: Hypothesis tests and confidence intervals

text: Suppose we have data set $X$ with $20$ observations. A confidence interval for the mean $\mu$ is always the same no matter what hypothesis we test? So if we test $H_a: \mu < \mu_0$, $H_a: \mu \neq \mu_0$ or $H_a: \mu > \mu_0$, the confidence interval will be the same?


title: How do I check for normality prior to using a Mixed Anova?

text: I have just conducted an experiment in which I have measured glucose concentration in patients (4 different treatment groups) at 3 different weeks ( pre, 6 weeks and 12 weeks). I want to use Mixed Anova to test if there is any significant difference between treatments in glucose concentration (also possible time and time*treatment interaction). However, Mixed Anova requires the data to be normally distributed. I am using SPSS. Must I perform for example a Kolmogorov‐Smirnov/Shapiro-Wilks test (glucose data (group1, time1,...,etc.) of every group separately (at each time point)? Or is there another way to look at normality before running Mixed Anova?


title: Using a general linear model to perform z-test

text: I am trying to perform a z-test with a brain-imaging software. I was told that the best way was to use a general linear model with one subject in one group and the rest in the other group. Using this method gives results that look identical to calculating the average and standard deviation of the control group and performing a specific z-test. Is it accurate to use a general linear model in this way? What statistical test might it be using? Thanks for your help!


title: Is "A person from Group X was 3 times more likely than a person from Group Y to have the disease" an accurate gloss of an odds ratio of 3?

text: All the interpretations I've seen have referred instead to the odds being 3 times higher in Group X. However, It is not clear to me that there is a meaningful difference between these formulations. Edit: In thinking that it was OK to express things in this way, I was influenced by http://sphweb.bumc.bu.edu/otlt/sparta/html/32a.html, which states: "An odds ratio is interpreted as if it were a relative risk. In this case an odds ratio of 3.46 indicates that people who ate food from the bakery had 3.46 times the risk of developing hepatitis compared to people who did not eat food from the bakery." Is that website perhaps wrong? Or am I perhaps interpreting it incorrectly?


title: ANOVA in R effect of variables on one another

text: I am seeing some thing very strange when I do an ANOVA in R. The actual data doesn't matter here, check this out - > summary(aov(t(freq)~t(mob)+t(tab)+t(comp))) and the F value is 4.4% for the comp variable, and >50% for the other two. It seems from the F value that there is a significant relationship between the t(comp) variable and the freq variable. But now - > summary(aov(t(freq)~t(comp))) and the f-value for the comp variable is now 77%. So now, there is no significant relationship between the freq and comp variables. What is going on here? The variables were independent. Why do the two ANOVA's give me seemingly contradictory results?


title: Violin plots interpretation

text: I am comparing the distrubution of different groups using the violin plots, however most of the online resources I found are just related to how to make the plots and very basic interpretation of the results (the median variation, the data is clustered or not). I am looking for detailed examples that I can follow as my guideline to correctly interpret the violin plots.


title: Testing parameter estimates

text: If I have derived parameter estimates for a distribution for some data, do I need to conduct a significance test on these to determine how valid they are? Or is the same effect done by just performing a goodness-of-fit test for the distribution?


title: Regression with correlated explanatory variables

text: I have variables of the following kind (coded in R):               set.seed(2)     dependent.variable = rnorm(12)     exp1 = c(1,3,4,8,3,4,1,5,6,6,7,9)     exp2 = c(1,3,6,2,1,1,3,4,6,4,1,1)     exp3 = exp1*exp2      my three explanatory variables are related by an equation of the kind: $exp1 = \frac{exp3}{exp2}$. For my real variables the relation might be slightly more complex and I haven't figured out the equation yet. I want to know if any of these three variables influence `dependent.variable`. Does it make sense? Can I simply run:               lm(dependent.variable~exp1*exp2*exp3)      or               lm(dependent.variable~poly(exp1,2,raw=T)*poly(exp2,2,raw=T)*poly(exp3,2,raw=T))      without having to care about the fact that the explanatory variables are not independent? If not, what should I do? Thanks a lot!


title: Confirmatory Factor Analysis with R lavaan

text: I tried to do a CFA with the `lavaan` package in R. Here is my model: ![enter image description here](http://i.stack.imgur.com/f1w70.jpg) I know how to define my latent variables like this:               MASTER.model <- ' Usability  =~ PUS1 + PUS2 + PUS3 + PUS4 + PUS5                         PerceivedUsefulness =~ PU1 + PU2 + PU3 + PU4 + PU5                        PerceivedEaseOfUse   =~ PEOU1 + PEOU2 + PEOU3 + PEOU4 + PEOU5                        BehavioralIntent =~ BI1 + BI2      But I don't know how to define the links between `Usability`, `EOU`, etc. Is this a covariance, written with e.g. `Usability~~PerceivedEaseOfUse` or do I write it like this: `PerceivedEaseOfUse~Usability` which means `PerceivedEaseOfUse` is regressed from `Usability`? My try was:               MASTER.model <- ' Usability  =~ PUS1 + PUS2 + PUS3 + PUS4 + PUS5                      PerceivedUsefulness =~ PU1 + PU2 + PU3 + PU4 + PU5                        PerceivedEaseOfUse   =~ PEOU1 + PEOU2 + PEOU3 + PEOU4 + PEOU5                        BehavioralIntent =~ BI1 + BI2                        PerceivedUsefulness ~ PerceivedEaseOfUse                        BehavioralIntent ~ PerceivedEaseOfUse                        BehavioralIntent ~ PerceivedUsefulness                         Usability ~~ Usability                         PerceivedEaseOfUse ~~ PerceivedEaseOfUse                         PerceivedUsefulness ~~ PerceivedUsefulness                                          PerceivedEaseOfUse ~ Usability'      This gave me proper results, but my RMSEA was really poor, which would mean my model doesn't fit my data. I just want to be sure it is not wrong because of me using R in a wrong way. Thank you.