title: after using AIC, how to determine the contribution or effect size of a individual covariate?

text: I am confused and looking for advise. I have found myself in this same situation repeatedly in the last few months. I want to know if covariate _X_ is influential or important. However, I also hypothesize that covariate _K_ , _L_ , and _W_ could be influential or are important. So, I generate lets say ~10 a priori hypothesized models. 1 full model and then variations of the covariates - thus creating subsets or nested models. (are these nested? ive no idea) I have run the AICc and find the top 4 models are "similar" and all within 2 delta. There is model#1 _K_ + _L_. model #2 _K_ + _X_ , model #3 of _K_ + _X_ + _L_ model #4 Intercept ONLY. Again, all I really care about understanding is if variable _X_. Because I did NOT set out to run all permutations of the variables (limiting myself to 10 model combinations) I cannot use the relative variable of importance calculation. So, how then can I evaluate the contribution of the variable X? I notice that variable K is in all of the top ranked models - but do not know how to evaluate the contribution of that variable either. How can one evaluate the effect of that 1 variable? I can run a univariate analysis to see if it improves AICc value from the null model - but then should I look at the odds ratio? Is there some sort of magnitude of contribution test? You just have to infer the strength by looking at the change in log likelihood values? Basically, I want to quantify how/why that variable _X_ is an uninformative parameter. Can I do this? if so - how? If you have one of the top models (or even the top model) being the Intercept Only model - how do you explain the performance of the rest of the models in relation to that? In my example it was ranked as 4th best model - so there seems to be "something" from the other variables, but... Oh - and the best model AIC weights never get above .28 - so they are really crap.


title: How to interpret these custom contrasts?

text: I am doing a one way ANOVA (per species) with custom contrasts.                    [,1] [,2] [,3] [,4]     0.5    -1    0    0    0     5       1   -1    0    0     12.5    0    1   -1    0     25      0    0    1   -1     50      0    0    0    1      where I compare intensity 0.5 against 5, 5 against 12.5 and so on. These is the data I'm working on ![enter image description here](http://i.stack.imgur.com/L7uVk.png) with the following results               Generalized least squares fit by REML       Model: dark ~ intensity        Data: skofijski.diurnal[skofijski.diurnal$species == "niphargus", ]             AIC      BIC    logLik       63.41333 67.66163 -25.70667          Coefficients:                 Value Std.Error  t-value p-value     (Intercept) 16.95 0.2140872 79.17334  0.0000     intensity1   2.20 0.4281744  5.13809  0.0001     intensity2   1.40 0.5244044  2.66970  0.0175     intensity3   2.10 0.5244044  4.00454  0.0011     intensity4   1.80 0.4281744  4.20389  0.0008           Correlation:                 (Intr) intns1 intns2 intns3     intensity1 0.000                           intensity2 0.000  0.612                    intensity3 0.000  0.408  0.667             intensity4 0.000  0.250  0.408  0.612           Standardized residuals:            Min         Q1        Med         Q3        Max      -2.3500484 -0.7833495  0.2611165  0.7833495  1.3055824           Residual standard error: 0.9574271      Degrees of freedom: 20 total; 15 residual      16.95 is the global mean for "niphargus". In intensity1, I'm comparing means for intensity 0.5 against 5. If I understood this right, the coefficient for intensity1 of 2.2 should be half the difference between means of intensity levels 0.5 and 5. However, my hand calculations don't match those of the summary. Can anyone chip in what am I doing wrong?               ce1 <- skofijski.diurnal$intensity         levels(ce1) <- c("0.5", "5", "0", "0", "0")         ce1 <- as.factor(as.character(ce1))         tapply(skofijski.diurnal$dark, ce1, mean)            0    0.5      5        14.500 11.875 13.000      diff(tapply(skofijski.diurnal$dark, ce1, mean))/2           0.5       5        -1.3125  0.5625


title: cross validation in ridge regression for classification. regularization issue

text: I perform ridge regression for classification. To find regularization parameter I do K-fold cross-validation with classification accuracy as a measure. This gives me some $\lambda$, which I then use in training of a final model on the whole available training data. The problem is that when I take 10*$\lambda$ my test accuracy on separate dataset is much better than with $\lambda$. I cannot see a reason for that. Tell me please, why this might happen? The lambda I get is 10^4 and beta coefficients are about 10^(-3). I have about 15000 features, which I standardize before doing regression.


title: How to analyze these data for HS science project

text: 10 kids had access to 4 computers, each with a different game (a-d). The test was 5 min.               kid 1 played game a for 238 sec     kid 2 played game a for 263 sec     kid 3 played game b for 280 sec     kid 4 played game a for 245 sec     kid 5 played game a for 220 sec, game b for 25 sec, and game d for 78 sec     kid 6 played game a for 280 sec     kid 7 played game a for 213 sec     kid 8 played game a for 47 sec, and game b for 211 sec     kid 9 played game b for 254, and game c for 21 sec     kid 10 played game a for 249      7 kids had a preference for game a, 3 kids had a preference for game b, no kids had preferences for games c or d. What is the best way to analyze if there is a statistical difference in these data?


title: Within-class covariance matrix

text: Let $\textbf{X}$ be an $n \times p$ matrix with the rows containing observations and the columns containing features. Also assume that the features are centered at $0$. Let $C_k\subset \{1, \dots n \}$ contain the indices of the observations that belong to class $k$. Why is an estimate of the within-class covariance matrix $$ \widehat{\Sigma_w} = \frac{1}{n}\sum_{k=1}^{K} \sum_{i \in C_k} (\textbf{x}_i- \hat{\mu}_k)(\textbf{x}_i-\hat{\mu}_{k})^{T}$$ For a particular class $k$, are interested in the covariance between observations. So we would have $\binom{|C_{k}|}{2}$ covariances for $k \in \mathbb{N}$. The matrix should be of size $\binom{|C_{k}|}{2} \times k$. **Added.** $w$ is just an abbreviation for within-class variance.


title: Finding a distribution by simulation?

text: Suppose you have one variable, $x$, with 8 data points, with a sample mean of 60%, and a sample standard deviation of 7%, and let’s also assume you know the sample comes from a lognormal distribution (or from a distribution with a heavier tail than lognormal). A sample with mean = 60%, SD = 7%, could be produced by an underlying distribution with, say, mean = 55%, SD = 10%; or it could be produced by a true mean of 53% and SD = 9%; or… Any number of other combination of mean and SD could produce the data we see. Now, I looked at true means running from 30% to 100% and SD’s running from 1% to 50% (all in 1% increments) – in other words, 3,550 combinations. For each combination, I simulated 10,000 groups of 8 data points from a lognormal distribution, calculated the mean and SD of those 8 points, and figured out how many of the 10,000 were “close” to my sample mean of 60% and SD of 7%. I defined “close” as being within 0.5% either way – in other words, if the simulated mean was between 59.5% and 60.5%, and the SD was between 6.5% and 7.5%, I assumed it “close.” Let’s suppose with a true mean of 53% and true SD of 9%, 10 out of the 10,000 simulations were “close” to my sample. After doing that for each of the 3,550 combinations, I had a distribution of the likelihood that each combination is really the true underlying combination. Next, given the likelihood of each combination, we can ask what is the probability that $x$ is, say, 70% or less. Again, we run through each of the 3,550 combinations, and get the answer. For example, if the true mean is 40%, with a 15% SD, that probability is, say, 0.903; if the true mean is 65% with a 12% SD, the probability is, say, 0.539. If those were the only 2 combinations, and they were equally likely, then we’d say the probability of $x$ being 60% or less is, say, (0.903+0.539)÷2 = 0.721. So my question is do you think that these calculations make sense? If no, please correct me, if yes, can you please provide me with any other way that is more neat than this one. Thanks.


title: Datasets for testing source code authorship attribution

text: I have been working on a research paper related to source code authorship attribution. In order to measure the performance of the technique, described in my paper, I need several data sets. Following papers described similar systems and tested those systems using several datasets. Unfortunately these datasets are not publicly available. Can somebody help me to find out, several publicly available datasets, suitable for testing source code attribution? (I’m particularly interested in datasets related to Java programming language) Thanks Shevertalov M, Kothari J, Stehle E, Mancoridis S. On the use of discretised source code metrics for author identification. Proceedings of the First International Symposium on Search Based Software Engineering, HarmanM, Di Penta M, Poulding S (eds.), IEEE Computer Society Press: Windsor, United Kingdom, 2009; 69-78. Burrows S, Tahaghoghi SMM. Source code authorship attribution using n-grams. Proceedings of the Twelfth Australasian Document Computing Symposium, Spink A, Turpin A, Wu M (eds.), RMIT University: Melbourne,Australia, 2007; 32-39. Frantzeskou G, Stamatatos E, Gritzalis S, Katsikas S. Source code author identification based on n-gram authorprofiles. Artificial Intelligence Applications and Innovations, Maglogiannis IG, Karpouzis K, Bramer M (eds.).Springer: New York City, New York, 2006; 508-515. Frantzeskou G, Stamatatos E, Gritzalis S, Katsikas S. Effective identification of source code authors using bytelevelinformation. Proceedings of the Twenty- Eighth International Conference on Software Engineering, OsterweilLJ, Rombach D, Soffa ML (eds.), ACM Special Interest Group on Software Engineering, ACM Press: Shanghai,China, 2006; 893-896. Frantzeskou G, Stamatatos E, Gritzalis S, Chaski CE, Howald BS. Identifying authorship by byte-level n-grams:The source code author profile (SCAP) method. International Journal of Digital Evidence January 2007; 6(1):1-18. Frantzeskou G, MacDonell SG, Stamatatos E, Gritzalis S. Examining the significance of high-level programmingfeatures in source code author classification. Journal of Systems and Software March 2008; 81(3):447-460. Frantzeskou G, Stamatatos E, Gritzalis S. Supporting the cybercrime investigation process: Effective discriminationof source code authors based on byte-level information. Proceedings of the Second International Conference on Ebusiness and Telecommunication Networks, Filipe J, Vasiu L (eds.), Institute for Systems and Technologies ofInformation, Control and Communication, INSTICC Press: Reading, United Kingdom, 2005; 283-290. Burrows S, Uitdenbogerd AL, Turpin A. Application of information retrieval techniques for source codeauthorship attribution. Proceedings of the Fourteenth International Conference on Database Systems for Advanced Applications, Zhou X, Yokota H, Kotagiri R, Lin X (eds.), Springer: Brisbane, Australia, 2009; 699-713.


title: If $X$ has a log-normal distribution, does $X-c$ also have a log-normal distribution?

text: Let's say $X$ has a log-normal distribution and there is one real positive number $c$. then is it right to say that $(X -c)$ also has some log-normal distribution? My feeling is that, it can't be, because $(X - c)$ may take negative value whereas a log-normal distribution is only defined on the positive domain. Can somebody disprove that?


title: Is Kolmogorov-Smirnov test valid with discrete distributions?

text: I'm comparing a sample and checking whether it distributes as some, discrete, distribution. However, I'm not enterily sure that Kolmogorov-Smirnov applies. Wikipedia seems to imply it does not. If it does not, how can I test the sample's distribution?


title: The difference between SVD and SVD++

text: What if any is the connection between SVD (the one you learn about in your linear algebra course) and SVD++ (the one from the Netflix prize)? I know they both want to find latent factor spaces. But that seems to be it. Am I missing anything?


title: Is there a radical difference in how bayesian and frequentist approaches treat nuisance parameters?

text: The wiki article on credible intervals has the following statement: > credible intervals and confidence intervals treat nuisance parameters in > radically different ways. What is the **radical** difference that the wiki talks about? Credible intervals are based on the posterior distribution of the parameter and confidence interval is based on the maximum likelihood associated with the data generating process. It seems to me that how credible and confidence intervals are computed is **not dependent** on whether the parameters are nuisance or not. So, I am a bit puzzled by this statement. PS: I am aware of alternative approaches to dealing with nuisance parameters under frequentist inference but I think they are less common than standard maximum likelihood. (See this question on the difference between partial, profile and marginal likelihoods.)


title: Non-parametric test for unequal samples with subsequent post-hoc analysis?

text: Is it okay to perform a Kruskal-Wallis on four unequal samples? Further to this, is there a subsequent pair-wise post-hoc suitable for two unequal sample sizes?


title: Seeking to understand asymmetry in hypothesis testing

text: I need to have one understanding on statistical hypothesis testing. In a typical hypothesis test, we have 2 opposite hypotheses; namely Null and Alternative. Here my textbook says that "those 2 hypotheses are not symmetrical in the sense that if we swap the hypotheses then the result will alter". Here I am unable to grasp the point which that textbook wants to say. Can somebody explain to me in detail? It would be helpful if someone can give some example of that asymmetry as well. Appreciate your help. Thanks,


title: Different prediction plot from survival coxph and rms cph

text: I've created my own slightly enhanced version of the termplot that I use in this example, you can find it here. I've previously posted on SO but the more I think about it I believe that this probably is more related to the interpretation of the Cox Proportional hazards model than with the actual coding. ## The problem When I look at a Hazard Ratio plot I expect to have a reference point where the confidence interval naturally is 0 and this is the case when I use the cph() from the `rms package` but not when I use the coxph() from the `survival package`. Is correct behavior by coxph() and if so what is the reference point? Also, the dummy variable in the coxph() has an interval and the value is other than $e^0$? ## Example Here's my test code:               # Load libs     library(survival)     library(rms)          # Regular survival     survobj <- with(lung, Surv(time,status))          # Prepare the variables     lung$sex <- factor(lung$sex, levels=1:2, labels=c("Male", "Female"))     labels(lung$sex) <- "Sex"     labels(lung$age) <- "Age"          # The rms survival     ddist <- datadist(lung)     options(datadist="ddist")     rms_surv_fit <- cph(survobj~rcs(age, 4)+sex, data=lung, x=T, y=T)      ### The cph plots This code:               termplot2(rms_surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05,               se.type="polygon", yscale="exponential", log="y",                xlab=c("Age", "Sex"),                ylab=rep("Hazard Ratio", times=2),               main=rep("cph() plot", times=2),               col.se=rgb(.2,.2,1,.4), col.term="black")      gives this plot: ![cph\(\) termplot2](http://i.stack.imgur.com/qnWsL.png) ### The coxph plots This code:               termplot2(surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05,                se.type="polygon", yscale="exponential", log="y",                xlab=c("Age", "Sex"),                ylab=rep("Hazard Ratio", times=2),               main=rep("coxph() plot", times=2),               col.se=rgb(.2,.2,1,.4), col.term="black")      gives this plot: ![coxph\(\) termplot2](http://i.stack.imgur.com/VWfNo.png) ## Update As @Frank Harrell suggested and after adjusting along suggestion in his recent comment I got:               p <- Predict(rms_surv_fit, age=seq(50, 70, times=20),                   sex=c("Male", "Female"), fun=exp)     plot.Predict(p, ~ age | sex,                  col="black",                  col.fill=gray(seq(.8, .75, length=5)))      This gave this very nice plot: ![Lattice plot](http://i.stack.imgur.com/RVXy3.png) I've looked at the contrast.rms again after the comment and tried this code that gave a plot... although there is probably much more that can be done :-)               w <- contrast.rms(rms_surv_fit,                        list(sex=c("Male", "Female"),                             age=seq(50, 70, times=20)))          xYplot(Cbind(Contrast, Lower, Upper) ~ age | sex,             data=w, method="bands")      Gave this plot: ![The contrast plot](http://i.stack.imgur.com/Q7AW2.png) ## UPDATE 2 Prof. Thernau was kind enough to comment on the plots lack of a confidence waist: > The smoothing splines in coxph, like the ones in gam, are normalized so that > sum(prediction) =0. So I don't have a fixed single point for which the > variance is extra small. Although I'm not yet familiar with GAM this does seem to answer my question: **this seems to be an issue of interpretation.**


title: Latent Semantic Analysis - Co-occurrence of words

text: Let $A[n\times m]$ represents the term-document matrix, where, $n$ is the number of terms and $m$ is the number of documents. This matrix can be composed into 3 matrices (SVD decomposition) such as, $A = U\times W \times V^T$ Its truncated version can be written as: $A_k = U_k\times W_k \times V_k^T$ Now the term-term matrix can be written as $T_k = U_k\times W_k \times W_k^T \times U_k^T$ According to the papers I read this $T_k$ represents the high order co- occurrence information. I understand that this matrix represents the first order co-occurrence but I did not understand how this matrix represents the high-order co-occurrence information. For me, I think that $T^N_k = U_k\times W_k \times (W_k^T \times U_k^T)^N$ will represents that $N^{th}$ (upto) order co-occurrence. Could anybody please explain it to me.


title: When should I update a recommendation engine?

text: [I asked this on StackOverflow and was told it would be a better fit here] I am including a basic recommendation engine in a very small project for my final exam. I understand the code and the math but I am not too certain as to when I should update the values. In this project it would be fairly simple to update in real time but in the real world would that not take an extreme amount of cpu time? Would it make more sense to run the algorithm on a schedule? In the project I am using linq2entities, storing the purchases in a table that I use to calculate the correlation coefficient. I only have a couple of thousand users in the db but I am trying to get an idea of what would be more efficient in practice.


title: Difference between a 2 factor ANOVA and mixed effects model

text: The `lme4` package in R includes the `cake` dataset.               library(lme4)     head(cake[,2:4], 20)        recipe temperature angle     1       A         175    42     2       A         185    46     3       A         195    47     4       A         205    39     5       A         215    53     6       A         225    42     7       B         175    39     8       B         185    46     9       B         195    51     10      B         205    49     11      B         215    55     12      B         225    42     13      C         175    46     14      C         185    44     15      C         195    45     16      C         205    46     17      C         215    48     18      C         225    63     19      A         175    47     20      A         185    29      I've analysed the `cake` dataset using two different models below. The first model is a 2 factor ANOVA:               summary(aov(angle ~ temperature + recipe, cake))                  Df Sum Sq Mean Sq F value   Pr(>F)         temperature   5   2100   420.1   6.918 4.37e-06 ***     recipe        2    135    67.5   1.112     0.33         Residuals   262  15908    60.7                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      ...and the second is a mixed effects model, with `temperature` as a random effect:               lmer(angle ~ recipe + (1| temperature), data=cake, REML=F)     Linear mixed model fit by maximum likelihood      Formula: angle ~ recipe + (1 | temperature)         Data: cake        AIC  BIC logLik deviance REMLdev      1893 1911 -941.7     1883    1877     Random effects:      Groups      Name        Variance Std.Dev.      temperature (Intercept)  6.4399  2.5377        Residual                60.2560  7.7625       Number of obs: 270, groups: temperature, 6          Fixed effects:                 Estimate Std. Error t value     (Intercept)   33.122      1.320  25.093     recipeB       -1.478      1.157  -1.277     recipeC       -1.522      1.157  -1.315          Correlation of Fixed Effects:             (Intr) recipB     recipeB -0.438            recipeC -0.438  0.500      Is someone able to provide a summary of what the mixed effect model has done differently to the ANOVA?


title: Distribution of random variable

text: The chloride concentration at depth $x$ mm at time $t$ is: $$ \text{C}(x,t)=\text{Cs}(1-\text{erf}(x/(2(D*t)^{0.5})) $$ where $\text{Cs}$ is chloride surface concentration with Lognormal distribution, $D$ is the diffusion coefficient with normal distribution, and $\text{erf}$ is error function. On the other hand, we have:   $$ \ln(1.08*\text{icorr})=7.98+0.7711*\ln(1.69*\text{C}(x,t))-\frac{3006}{T}-\frac{\exp(8.03)}{(1+1.69*\text{C}(x,t))^{0.549}}+2.24*t^{-0.215} $$   where $T$ is temperature (Kelvin) with normal distribution. As you can see, $\text{icorr}$ depends on $T$ and $C(x,t)$ and $C(x,t)$ depends on $Cs$ and $D$. We know the distributions of $Cs$, $D$, and $T$. **How can I determine the distribution of $\text{icorr}$?**


title: Pearson's or Spearman's correlation with non-normal data

text: I get this question frequently enough in my statistics consulting work, that I thought I'd post it here. I have an answer, which is posted below, but I was keen to hear what others have to say. **Question:** If you have two variables that are not normally distributed, should you use Spearman's rho for the correlation?


title: R package for fixed-effect logistic regression

text: I'm looking for an `R` package for estimating the coefficients of logit models with individual fixed-effect (individual intercept) using Chamberlain's 1980 estimator. It is often known as Chamberlain's fixed-effect logit estimator. It's a classic estimator when dealing with binary outcome panel data (at least in econometrics), but I just don't find anything related to it in the CRAN. Any clue?


title: Can I simply remove one of two predictor variables that are highly linearly correlated?

text: Using Pearson's Correlation Coefficient, I have several variables that are highly correlated ($\rho = 0.978$ and $\rho = 0.989$ for 2 pairs of variables that are in my model). The _reason_ some of the variables are highly correlated is because one variable is used in the _calculation_ for another variable. Example: $B = V / 3000$ and $E = V * D$ $B$ and $E$ have $\rho = 0.989$ Is it possible for me to just "throw away" one of the variables?


title: How to prove the exponential expansion recursive form?

text: Consider the following exponential expansion form: $$ \exp\left[\sum_{k=1}^\infty \gamma_k x^k\right] = \sum_{j=0}^\infty\delta_j x^j $$ where $\gamma_k$'s are known, and $\delta_0=1$, $$\delta_{j+1} = \frac{1}{j+1}\sum_{i=1}^{j+1}i\gamma_i\delta_{j+1-i}$$ Does anyone know how to prove it?


title: For a classification problem if class variable has unequal distribution which technique we should use?

text: for eg. if I have a class variable credit scoring with two classes good and bad, where #(good) = 700 and #(bad)= 300. I do not want to shorten my data. which technique I should use? I was using SVM but it is giving all bad as good in predictions.


title: How should I select winners in this Alternative Vote result set, where there are no majorities?

text: In the simulated Alternative Vote below, I had 10,000 people vote for 100 candidates using AV style voting (also known as instant runoff voting). Most people want person #57 or person #78 to win, but neither was able to reach a majority.   * What math/stats function will allow me to eliminate those with the least support? I could just take the mean, and eliminate those below that line, but there may be a smarter way of doing this. **More detail on my use-case:** I plan to use AV-style voting to allow people to choose which charities to automatically donate to. I want to allow multiple winners, and to payout with a ratio of some kind that evenly distributes the donations.               Voting start.     Count Ballots start.     Output: (PersonID)  (Votes)     1 1     70 3     3 1     73 3     78 4932     11 1     81 1     82 15     16 21     85 2     87 1     88 3     21 1     89 13     90 1     23 18     91 1     93 1     96 1     100 1     33 1     34 2     35 1     39 1     46 3     47 1     54 1     56 1     57 4963     63 4     65 1      PS - I'm not a stats guy, so retagging would be appreciated


title: "finely" labeled classification problem

text: I came across a kind of classification problem. Suppose I have a dataset with labels, where each label is one of a1, a2, a3, b1, and b2. I want to make a classifier that estimates "coarse" label, a or b. Here is a concrete example. Imagine a set of animal images labeled by the its species. We want make a classifier that predicts big classes, like mammal, reptile or fish. **Q1.** Does anybody know the name of this kind of classification problem setting? * * * To solve the problem, we can think of several approaches.   1. First make the label coarse (a1->a, a2->a, b1->b,... ) and solve it as a standard classification problem.   2. First solve the standard classification problem, predicting a1, a2,... Then make the label coarse.   3. First solve the standard classification problem, predicting a1, a2,... (level-0 classifier) Then use the output scores of level-0 classifier as input, and make a level-1 classifier which outputs a, b,.. (stacked generalization approach) **Q2.** What kind of approach is good for this kind of problem? I guess the approach 1 is not optimal since part or label information is thrown away. As Marc Claesen has pointeded out, this approach 1 is good if there is an example that is obviously a member of a but not obviously a member of one of a1, a2,...


title: Defining a threshold for percentage data with high amounts of exact zeros

text: I am trying to analyse my biological data which is derived from flow cytometry and describes the percentage of cells with a certain property. Data is available for two conditions: once the cells areuntreated and once the cells were stimulated. For the untreated condition I analysed 57 samples, for the stimulated condition 39, which can be considered as more or less independent random variables. In both conditions, I have a lot of samples with exactly 0.0000% cells with the property of interest, however in the treated condition there are some samples that show "high" percentage (up to 0.35%). (exemplary data posted at the end of the post) I now want to use the samples of the "untreated" condition to define a threshold for random background activation in order to say in the "treated" condition from which % on, I can treat the sample as "really" positive for the property. Initially, I tried to define the mean of the untreated condition, add 2 times the standard deviation (calculated on Poisson model) and set this as the threshold. However, it turned out that the standard deviation as defined by Poisson distribution is getting quite big when calculating with these small numbers. This led to the threshold cutting off a number of events I think should be real. The problem is now, that I don't know which model I should apply to solve this problem. I tried transformations however, since the 0 values are exactly 0, I never got rid of them. Further, I read something about zero-inflated Poisson distributions but I don't know whether that helps here and I also don't really know how to apply that. **So to summarize, the problem is: How can I say about a single event in the treated condition whether it is significantly higher than the mean of the untreated condition?** If someone has some advice on how to deal with this data, I would be really grateful. Thanks a lot. Cheers, Florian **Data** Treated: 0.0000, 0.0000, 0.0000, 0.0000, 0.0130, 0.0000, 0.0000, 0.0000, 0.0021, 0.0010, 0.0000, 0.0000, 0.0016, 0.0021, 0.0054, 0.0056, 0.0067, 0.0000, 0.0022, 0.0052, 0.0000, 0.0004, 0.0000, 0.0000, 0.0021, 0.0000, 0.0000, 0.0283, 0.0142, 0.0215, 0.3570, 0.0000, 0.0000, 0.0093, 0.0109, 0.0000, 0.0000, 0.0011, 0.0000 Untreated: 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0021, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0043, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0039, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0020, 0.0010, 0.0015, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000


title: Most suitable algorithm for optimizing Maximum likelihood function

text: What is the most suitable optimization algorithm for optimizing maximum likelihood estimator? In excel I used GRG non linear optimization algorithm, is that good enough? I want to write my own code to understand optimization better so I want to choose most suitable non linear constrained optimization algo for MLE. Cheers


title: Alias Structure of one-fourth replicate of a $4^2$-Factorial Design with interaction $\text{A}{B}^3$ confounded

text: For finding the Alias of main effect _A_ , i started as the following : $\text{A}$$\times$ $\text{A}{B}^3$$=\text{A}^2\text{B}^3={(\text{A}^2\text{B}^3})^2=\text{A}^4\text{B}^6=\text{B}^2$(mod $4$)$={(\text{B}^2)}^2=\text{B}^4=0$(mod $4$) Where am i doing mistake?


title: Calculate prediction error for a multinomial output like in the "effects" package

text: After fitting a multinomial model to my data with the "multinom" function (package nnet), I want to show the effect of selected variables controlling for other variable values. I know that the "effects" package do mainly what I want, but I want to be able to calculate the prediction error (confidence interval) by myself. Does someone could tell me the methodology and if possible the R code? I think we should use the delta method, but I'm not sure how to apply it in this case. Here is a small example code (based on data available in the effects package)               library(nnet)     library(effects)     mod <- multinom(vote ~ age + gender, data=BEPS)     summary(mod)          # Call:     # multinom(formula = vote ~ age + gender, data = BEPS)          # Coefficients:     #                  (Intercept)         age gendermale     # Labour             1.2241862 -0.01562320  0.1682676     # Liberal Democrat   0.4979706 -0.01551381  0.1240998          # Std. Errors:     #                  (Intercept)         age gendermale     # Labour             0.2277826 0.003830006  0.1204621     # Liberal Democrat   0.2694373 0.004578836  0.1436882          # Residual Deviance: 3186.266      # AIC: 3198.266           plot(allEffects(mod))      ![output of the effect function](http://i.stack.imgur.com/TbYtq.png) The only thing I need is to be able to calculate the values of errors shown in this graph ! Thank you in advance,


title: Is a data point significantly larger than a certain distribution average?

text: I have a simulated distribution with mean 12.53% and standard deviation 11.83%. The sample size is big enough (10,000) to assume it is a Normal distribution. How do I properly test if the value "26.05%" is significantly larger than the mean 12.53%? Can anyone please help me to write the null hypothesis, as well as the test, or just give me any reference that I'm not being able to find (or most probably to recognize) on the web?


title: Minimum Number of Variables for Effective use of the Genetic Algorithm

text: From my understanding of the genetic algorithm the population consists of individuals, where each individual is a potential solution made up of "genes", and each gene is a variable. So for a cost function that takes in 3 variables the individuals will have 3 genes. The optimization happens by selection (picking the fittests from the population), crossover (swapping genes between pairs in the population), and mutation (randomly changing certain genes). Because of this it seems to me that the genetic algorithm is best suited to high dimensionality optimisation problems, i.e ones with lots of variables, so that there are more genes that can be swapped and mutated. Is this true? and what is a good number of genes for use in the genetic algorithm? Thanks!


title: PCA realization in python

text: Input matrix M x N, M rows (number of samples), N cols (dimesionality of data). But I don't understand why we must transpose matrix when put it in np.cov?               def pca():         #M x N         xs= np.loadtxt("data_3d.txt",delimiter=" ", skiprows=1, usecols=(0,1,2))         print xs.shape         # print xs              #get mean         mean= np.mean(xs,axis=0)         # print mean.shape         # print mean              #N x M         data= (xs-mean).T # why need transpose?         print data.shape         # print data              #N x N         covData=np.cov(data)#calculate covariance matrix         print covData.shape              eigenvalues, eigenvectors = np.linalg.eig(covData)         print eigenvalues.shape # N long         print eigenvectors.shape # N x N         print eigenvalues         print eigenvectors      And how then project and reconstruct data? By projection we need to obtain matrix M x k where k< N. By reconstruction we must obtain data before we put it in np.cov , so it must be N x M.               #sort and get k largest eigenvalues     k=2     idx = eigenvalues.argsort()[-k:][::-1]     print idx          eigenvalues = eigenvalues[idx] # k long      eigenvectors = eigenvectors[:,idx] # N x k     print eigenvalues.shape     print eigenvectors.shape     print eigenvalues     print eigenvectors          #projection and reconstruction     pr= np.dot(data.T,eigenvectors) # (M N) * (N k) => (M k)     rec= np.dot(eigenvectors, pr.T) #(N k) * (k M) => (N M)     print (data-rec) # test reconstruction error


title: How to calculate the next year's January hourly forward prices for an asset given the historical hourly prices?

text: I have an exercise given to me during the recruitment process which asks to calculate the next year's January hourly forward prices for an asset given the historical hourly prices of that asset and next year's January average price (average of all hours in that period). The exercise asks to take account of the seasonality (modelling seasonality), hence there is a 9 year hourly price dataset given. I have studied statistics and econometrics at school but have forgotten most of the techniques that could be used for price estimation. Any idea how this exercise can be solved? I can provide a sample dataset on request.


title: Mean Absolute Error and Data Distribution

text: I use a memory-based learning model to predict human scores in a [0, 10] range (quiz results). As a forecast error metric I use Mean Absolute Error. I was wondering what is the relation between MAE and data distribution. How can i use the data (or the error) distribution to explain the error and possibly argue whether it is acceptable. For example, can a MAE of 60% be characterized as "good" for a non-normal data distribution rather than a normal one? Since i am new to statistics, I would dearly appreciate any suggested reference to related bibliography. Thank you in advance.


title: Trend analysis with three time points in a repeated measures design

text: Assume that there are 15 subjects and that a quantitative feature is measured for all subjects at three equally-spaced times points. The research of interest is to find out wether there is a linear trend across time. Admittedly the number of time points is quite limited, but what are the options for trend analysis? I can think of two approaches: 1) The traditional trend analysis is usually performed with weights assigned based on orthogonal polynomials. With three time points, this would lead to the weights of -1, 0, and 1. The problem with this approach is that the data at the midpoint is essentially ignored (under the equal space assumption), leading to a contrast between the two end time-points. 2) Alternatively I can try linear mixed-effects modeling with lme from R package nlme using the three time points as an explanatory variable (1, 2 and 3?). Would this approach be reasonable? Any suggestions? Thanks!


title: Measuring goodness of fit for mixed logistic regression model - inconsistent results from R squared and AUC

text: I am trying to assess the goodness-of-fit or accuracy of 6 generalised linear models. I first assessed this using AUC (calculated from function auc1 described here), and got results ranging from 0.65 to 0.82, and therefore concluded that my models explain a reasonable amount of the variation in my datasets. However, I then tested them again using Nakagawa and Schielzeth's (2012) method for obtaining R^2 from GLMM, and Nagelkerke's modified R^2 based on likelihood ratio (both in R's MuMIn package). I would expect slight differences in the results, but the results I have are totally inconsistent. Anyone have a suggestion for why this might be? Nakagawa and Schieleth's R^2 function in MuMIn comes with a warning that it is in the experimental stage and should be used with caution. Should I therefore just rely on the results from the AUC calculations? Details about the models: I am using my models to try to identify predictive factors in the presence or absence of particular viruses in several wildlife species. The fixed effects are various individual traits (eg. sex, age, month captured, presence of parasites, scaled mass index (SMI), and random variables are site (nested in region), year, and observer(because of potential bias in measurements). Models differ by virus or host species, and explanatory variables not relevant to certain species are not included in those models, but the models are on the whole very similar. Variation in goodness-of-fit measures:               Model        AUC        Nakagawa & S's           Nagelkerke's                              r.squaredGLMM            r.squaredLR                              (R^2m/R^2c)*             (adj.R^2)      Model 1      0.78       0.96/0.96                0.32     Model 2      0.72       0.22/0.22                0.09     Model 3      0.81       0.96/0.96                0.18     Model 4      0.78       0.09/0.37                0.07     Model 5      0.84       0.28/0.54                0.16     Model 6      0.65       0.08/0.08                0.04      note: * R^2m refers to variation from fixed effects only, R^2c includes random effects (so the full model)


title: Relation between variational Bayes and EM

text: I read somewhere that Variational Bayes method is a generalization of the EM algorithm. Indeed, the iterative parts of the algorithms are very similar. In order to test whether the EM algorithm is a special version of the Variational Bayes, I tried the following:   1. $Y$ is data, $X$ is the collection of latent variables and $\Theta$ is the parameters. In Variational Bayes we make can make a approximation such that $P(X,\Theta|Y) \approx Q_X(X)Q_\Theta(\Theta)$. Where $Q$s are simpler, tractable distributions.   2. Since the EM algorithm finds a MAP point estimate, I thought that Variational Bayes can converge to EM if I use a Delta Function such that: $Q^1_\Theta(\Theta)=\delta_{\Theta^1}(\Theta)$. $\Theta_1$ is the first estimate for the parameters as usually done in EM.   3. When $Q^1_\Theta(\Theta)=\delta_{\Theta^1}(\Theta)$ is given, $Q^1_X(X)$ which minimizes the KL Divergence is found by the formula $$Q^1_X(X)=\frac{\exp(E_{\delta_{\Theta^1}}[\ln P(X,Y,\Theta)])}{\int\exp(E_{\delta_{\Theta^1}}[\ln P(X,Y,\Theta)])dX}$$ The formula above simplifies to $Q^1_X(X)=P(X|\Theta^1,Y)$ , this step turns out to be the equivalent of the Expectation step of the EM algorithm! But I can't derive the Maximization step as the continuation of this. In the next step we need to calculate $Q^2_\Theta(\Theta)$ and according to Variational Bayes iteration rule this is: $$Q^2_\Theta(\Theta)=\frac{\exp(E_{P(X|\Theta^1,Y)}[\ln P(X,Y,\Theta)])}{\int\exp(E_{P(X|\Theta^1,Y)}[\ln P(X,Y,\Theta)])d\Theta}$$ My question is, is really VB and EM algorithms are connected in this way? How can we derive EM as a special case of the Variational Bayes, is my approach true? Thanks in advance.


title: Intra- and inter-rater reliability on the same data

text: I've read some of the answers here (specifically the 'Inter-rater reliability for ordinal or interval data' one). I'm still somewhat perplexed though! I have data for 4 raters, who each rated every subject's CT scan, on two occasions, for the presence or absence of 9 signs. The two readings were separated by a washout period. As such, my data looks like this (for each sign/variable):               A1  A2  B1  B2  C1  C2  D1  D2     0   1    1   1   0   1   0   0     1   1    0   0   1   1   1   1     etc      With A1 representing the first reading by rater A, and A2 the second, and so on. Calculating the intrarater reliability is easy enough, but for inter-, I got the Fleiss Kappa and used bootstrapping to estimate the CIs, which I think is fine. Except, obviously this views each 'rating' by a given rater as being different 'raters'. I thought about using the ICC with the WinPepi function which can provide interrater agreement for fixed raters, but I suspect it might not be appropriate for binary data? Would you just get the Fleiss and report it as such? PS: The answer provided below by StasK is brilliant, and did clarify some things, but because of an ambiguity in the original question it isn't quite what I was looking for...


title: Linear regression with factors in R

text: I'm trying to understand how exactly factors work in R. Let's say I want to run a regression using some sample data in R:               > data(CO2)     > colnames(CO2)     [1] "Plant"     "Type"      "Treatment" "conc"      "uptake"        > levels(CO2$Type)     [1] "Quebec"      "Mississippi"     > levels(CO2$Treatment)     [1] "nonchilled" "chilled"        > lm(uptake ~ Type + Treatment, data = CO2)          Call:     lm(formula = uptake ~ Type + Treatment, data = CO2)          Coefficients:      (Intercept)   TypeMississippi  Treatmentchilled              36.97            -12.66             -6.86        I understand that `TypeMississippi` and `Treatmentchilled` are treated as booleans: For each row, the initial uptake is `36.97`, and we subtract `12.66` if it's of type Mississippi and `6.86` if it was chilled. I'm having trouble understanding something like this:                > lm(uptake ~ Type * Treatment, data = CO2)           Call:      lm(formula = uptake ~ Type * Treatment, data = CO2)           Coefficients:                      (Intercept)                   TypeMississippi                             35.333                            -9.381                   Treatmentchilled  TypeMississippi:Treatmentchilled                             -3.581                            -6.557        What does it mean to multiply two factors together in an `lm`?


title: What is the best way to learn the fundamentals of probability required for machine learning algorithms?

text: I took a probability course in university a few years ago, but I'm going through some machine learning algorithms now and some of the math is just befuddling. Specifically right now, I'm learning the EM algorithm (expectation maximization) and it seems that there is a large disconnect between what is required and what I have. I'm not asking for a book or a website, but what is the way to learn enough of these topics to be able to get a thorough understanding of the algorithms that use them? Is it required to go through a book and do hundreds of exercises? Or is that overkill in this sense? edit: If this is the wrong location for this question, please vote to migrate :)


title: Approximation of logarithm of standard normal CDF for x<0

text: Does anyone know of an approximation for the logarithm of the standard normal CDF for x<0? I need to implement an algorithm that very quickly calculates it. The straightforward way, of course, is to first calculate the CDF (for which I can find suitably simple approximations on Wikipedia), and then to take the logarithm of that. Obviously I'd like to avoid the time-cost of having to calculate two special functions, not to mention that the tiny intermediate value (in the tail) rules out using fixed-point arithmetic which is much faster than floating-point arithmetic. I know there are hundreds of approximations for all kinds of statistical functions, but the fact that this is the logarithm of one makes it harder to find one. I'd be very grateful if anyone could point me to one, or to a source where I might find one.


title: half-cauchy prior for scale parameter

text: I am looking for a prior for a scale parameter for which I have prior knowledge such that: "$\sigma$ typically does not exceed 100." ("typically" meaning that occasionnally this can happen). In such a context, I notice in the paper "Prior distributions for variance parameters in hierarchical models" of Andrew Gelman the following recommandation: > [...] When more prior information is desired, for instance to restrict σ > away from very large values, we recommend working within the half-t family > of prior distributions, which are more flexible and have better behavior > near 0, compared to the inverse-gamma family. A reasonable starting point is > the half-Cauchy family, with scale set to a value that is high but not off > the scale." As I understand it, a Cauchy (thus half-Cauchy) distribution has an infinite variance and I am not confortable with the idea of building an informative prior with an infinite variance density. Have you some insight on why my interpretation is bad/unsuited ? Moreover, have you some alternative proposals for my prior ?


title: Basic 2SLS IV Questions in Stata

text: **(1)** If I believe my instrument is exogenous conditional upon a few exogenous variables, do I include them only in the first stage? I.e. would the command be:               ivregress 2sls Y (X= inst1 inst 2 exog1 exog2) exog3 exog4      Where `Y` is the dependent variable, `X` is an endogenous variable, `inst1` and `2` are the instruments for `X`, `exog1` and `exog2` are what make the instrument exogenous to the error in the second stage, and `exog3` and `exog4` are just second stage variables that enhance precision? Would I need to include `exog1` and `exog2` outside the parenthesis also? **(2)** If I want to interact endogenous variable `X` with `exog3`, how should I instrument `X` and interact it? As follows? `ivregress 2sls Y (X*exog3= inst1 inst 2 exog1 exog2) exog4`


title: Correlation tests – multivariate correlation matrix?

text: I just got comments from a reviewer to a submitted article and didn't understand what I should do very well. Here are the tests I performed: We first used principal components analysis (PCA) to ordinate the sites according to topological and physical variables, and to bring out potential shifts in physical conditions following anthropogenic disturbance (`ref`: reference, `log`: logging, `fog`: formerly gold mining and `cug`: current gold mining). Prior to analysis, continuous variables were log-transformed, and discrete variables expressed in percentages were arcsine-transformed. Plots of the first two ordination axes usually capture most of the variance and consequently contain most of the information that is likely to be interpretable (Waite et al., 2000). Neighbouring sites in the scatterplots were expected to define areas with similar physical environments. Conversely, sites having a large distance to each other were expected to be distant in the feature space, according to environmental characterization. To bring out relationships between water chemistry, disturbance, and local environments, significant differences among a priori groups were tested using Kruskal–Wallis tests using site coordinates on the two first axis of the PCA. Then, significant differences in physical characteristic between a priori groups of sites (`Ref`, `Log`, `Fog`, `Cug`) were further assessed using Wilcoxon tests. **Comments:** > At an unbalanced design (Ref=57, Log=15, Fog=9, Cug=14), the PCA can be a > good ordination method but not an adequate comparison method! Maybe > correlation tests (or multivariate correlation matrix) could be much more > informative due that several physical variables are nested into each stream.   * Does he meant to use raw data instead of PCA coordinates?   * What should I do?   * How can I perform a correlation test using environmental data (quantitative) and categorical labels (qualitative data)? My purpose was to check if environmental variables differ among my 4 categories of environmental quality.


title: coupon collector's problem - labelled or not?

text: I have come across the famous coupon collector's problem. Linking the problem to the Twelvefold way, my first reaction has been to use the case with `unlabelled balls with labelled boxes`. However, it seems that the question is answered typically (correctly?) with `labelled balls with labelled boxes` case. Please see this question question for the approaches. So I wonder why this is the case.


title: Alternatives to ratios of counts if denominators can be zero

text: This is probably a really dumb question but I can't find the answer anywhere. I have been given a bunch of data and asked to calculate proportions between two variables. The variables are both counts (this is social science stuff). The problem is that there are a number of zeroes dotted althrough my dataset which of course present a problem when calculating proportions. Is there a way that I can transform the data to get rid of the zeros? I'm not going to analyze it, so I don't need it to satisfy any specific assumptions.


title: Trying to determine the similarity of the distributions of labels

text: I asked this question on stackoverflow first since it is more of an implementation question, but I got no answers yet, therefore I'm asking it here again. I have a dataset which is annotated by 2 different sources. This means that for every sample, 2 different sources give a label. I want to see how similar the two sets of labels are. Therefore I thought of using the Spearman correlation and the test for homogeneity with chi squared. I made this question because it is the first time I'm doing this and also in Python. So I'm not sure I'm doing them right, and would like some feedback. For example I have my 2 sets:               set1 = [1,5,7,3,2,4,...]     set1 = [1,3,7,2,2,1,...]     rho, pval = scipy.stats.spearmanr(set1,set2)      The results are:               rho = 0.456498145813     p value = 0.0      What do these results mean? (I think there is very little correlation, though why do I have a $p$ value of 0?) And for the test for homogeneity I calculate the frequencies and use the chi squared test:               chisq , p = scipy.stats.chisquare(np.array(distributionSet1),np.array(distributionSet2),6)      The results here are:               chisq = 47611.7764023     p = nan      How can the $p$ value be `nan`? And what does such a high chi squared value mean? For the degrees of freedom, I calculated 6 (I have 2 sets, so 2 populations, and 7 categories, so $(2-1)\times(7-1)=6$).


title: Can Someone Explain How Factor Multiplication Works with Factor Graphs?

text: I'm taking the Probablistic Graphical Model course here: https://class.coursera.org/pgm-003/ This class uses the concept of Factors extensively with regards to graphical models: http://en.wikipedia.org/wiki/Factor_graph I'm learning about making Clique Trees, http://en.wikipedia.org/wiki/Tree_decomposition, and part of that includes multiplying factors together to create new clique potentials, which will be used for inference later when you want to find marginal probabilities. So, if you have: F(A,B) and F(C | D) , you can make a clique with scope C(A, B, C, D) by using factor multiplication with F(A,B) and F(C | D) And that's where I get lost. Can you just multiply conditional and joint probabilities like that? And then end up with a new joint probability? I would really appreciate if someone could explain how that's done.


title: Cluster analysis and chi square test

text: I performed a cluster analysis and now want to consider if the variable “I have a smart phone” (yes/no) significantly differ between the cluster solutions by using chi square test. The p-value was greater than 0.05, so I have to accept the null hypothesis (there is no significant difference between the clusters and having a smart phone). Within all cluster solutions, the majority of the people has a smart phone. Can I therefore assume that all interviewees, no matter to what cluster they are assigned, use smart phones in their lives? Does this makes sense? Thanks!


I have an elementary question that I'm sure must be answered in textbooks somewhere, but I haven't found it yet. Why is the average the right way to deal with Gaussian noise? Let's flesh this out a bit. Here's my model. There is an unknown parameter $\theta$. We have iid random variables $X_1,\dots,X_n$, defined as $X_i = \theta + Y_i$ where $Y_i \sim \mathcal{N}(0, \sigma^2)$. In other words, the $Y_i$'s are Gaussian noise, and each $X_i$ is a noisy observation of the underlying parameter. Now, given observations of $X_1,\dots,X_n$, we want to estimate the underlying parameter $\theta$. I think I remember reading that the optimal way to estimate $\theta$ is using the average: i.e., $\hat{\Theta} = \frac{1}{n} (X_1 + \dots + X_n)$. Is this true? In what sense is this optimal, and why is it optimal?


title: Can AIC be used to compare an ARMA model to an ARMA-GARCH model?

text: Suppose I have one time series and two competing models that describe it. Model 1 is ARMA$(p_1,q_1)$, model 2 is ARMA$(p_2,q_2)$-GARCH$(r,s)$. I obtain AIC values of model 1 and model 2. I would like to conclude that the model with the lower AIC value is the prefered model. However, I am not sure whether I can compare AIC values coming from two different classes of models such as ARMA and ARMA-GARCH.      My intuition is the following: if we find a model that nests both model 1 and model 2, we should be able to use AIC to compare the former two. Indeed, ARMA$(p_1,q_1)$ and ARMA$(p_2,q_2)$-GARCH$(r,s)$ are both nested within ARMA$(\max\{p_1,p_2\},\max\{q_1,q_2\})$-GARCH$(r,s)$ and may be obtained from the latter model by applying zero restrictions on coefficients.      However, I am not sure whether my approach makes sense. Could someone please elaborate on this issue and answer the question: can we, or can we not, use AIC to compare ARMA models with ARMA-GARCH models?


title: Time series Data Analysis and Forecasting by country and time factor

text: cty year qtr tl Argentina 2009 Q4 3 Argentina 2010 Q1 2 Argentina 2010 Q2 7 Argentina 2010 Q3 7 Argentina 2010 Q4 10 Argentina 2011 Q1 7 Argentina 2011 Q2 7 Argentina 2011 Q3 1 Argentina 2011 Q4 7 Argentina 2012 Q1 5 The data set has around 40 countries with each country having data for 5 years and four quarters. cty=country, year=Year, qtr Quarter, var=count of people. As I am new to R, can someone assist me with how to prepare the dataset for time series and fit a ARIMA/ARMA models etc. followed by the forecast in R? Some basic codes to prepare the data, to perform analysis and forecast for this series as an example would be helpful.


title: What do people call the groups into which quantiles divide the population?

text: Is there a correct technical name for the group of observations between two quantiles? For example, if you have the values of the four cut-points that divide a population into five groups of equal size, my understanding is that these levels are the quintiles. However, there is often a reference to “quantile means” as the means of the subgroups of the total population when divided by the quantiles. My impression is that is not considered correct usage among statisticians, as the quantile means should refer to the mean value of the quantiles, across samples, say, not to the groups of observations bracketed by them. Is there a proper name for these groups? I am interested in learning more about the statistical properties of the groups defined by the quantiles, rather than the quantiles themselves, and am having trouble finding even basic information – whether means for large populations will be approximately normal, how to estimate mean and variance when all you have is parameter estimates for the underlying distribution from which the sample is drawn, whether the sample means are unbiased estimators of the population means for these groups, construction of confidence intervals for means, etc. If I had a correct name for the groups, I think I should be able to find all these things rather quickly. In addition, I would be grateful if anyone could point me toward their favorite entry to literature that addresses these questions. Thanks!


title: Normalized likelihoods

text: AIC (BIC) model selection methods are widely used. These methods can select non-nested models unlike likelihood ratio type selection that requires model to be nested. The AIC has definition $\text{AIC}=-2\text{log}L+2V$ where $L$ is the likelihood of the model and $V$ is the effective number of paramters in that model. Take for instance making a decision whether to use a binary or an ordinal logistic regression model. You receive the response variable in an ordinal form with labels say $1,2,3 $. But it is required from you to predict observations with label $\{1,2\}$ as one group and observations with label $3$ as another. Therefore you have two ways to go 1)Group the data into two groups with labels $0$ for observation with labels $\{1,2\}$ and $1$ for for observations with label $3$ and use a binary logistic model. 2)Use the original $3$ labels $1,2,3$ in an ordinal logistic regression model and end up using the boundary from your ordinal model that separates observations with labels $2$ from observations with labels $3$. Since it is an ordinal model it will automatically be the boundary that separates group $\{1,2\}$ from group $3$ I will represent the two scenarios![Blue and red lines are boundries from an ordinal model, black line is the boundry from a binary model](http://i.stack.imgur.com/a7CWT.png) with plots. Blue and red lines are boundaries from an ordinal model, black line is the boundary from a binary model. Black points are labelled with 1, red points 2 and blue points 3. Therefore as described above I am either using the blue line when I choose to use an ordinal model or the black line when I choose to use the binary model. Using the AIC to choose either the binary or the ordinal model is possible. But the two models differ in two things. 1) Likelihood scale 2) Parameter scale. Is there a normalization that can be done so that the two models can be compared in terms of their AICs?


title: How to judge if a supervised machine learning model is overfitting or not?

text: Can anyone tell me how to judge if a supervised machine learning model is overfitting or not? If I don't have an external validation dataset, I want to know if I can use ROC of 10 fold cross validation to explain overfitting. If I have an external validation dataset, what should I do next?


title: Ripley's K Function and L Function for Point Patterns

text: The following is a spatial point pattern: ![enter image description here](http://i.stack.imgur.com/3NXRG.jpg) and these are the corresponding Ripley's K function and L function for this data: ![enter image description here](http://i.stack.imgur.com/SsZAs.jpg) ![enter image description here](http://i.stack.imgur.com/tOFv2.jpg) How are these functions interpreted?


title: Model validation and verification for Markov Chain switching model

text: Assuming I have a discrete-time Markov chain with only five states. The chain will be used for the prediction of the macroscopic states which are observable and coming from a timeseries. I use maximum likelihood estimation to calculate the transfer probability matrix. Now I hvae the TPMs and I can plot various realizations in time. How do I validate, that the chain is "right"? If I had many states predicting a time series I could use the MSE or something like that. But I think this doesn't make sense here. Do you have any tips for me? Thank you!


title: 1sigma error on the mean

text: My x and y data (n) have a non normal distribution. I'm measuring the mean of the distribution (I'm using Python):               x_m=np.mean(x)     y_m=np.mean(y)      I would like to measure error that represent the standard deviation of the mean (from the 16th and 84th percentile). Do you think this is the right way to do it?               y1=(numpy.percentile(y,16))/n     y2=(numpy.percentile(y,84))/n      Are there any other possibility to measure such an error? Thank you for your help.


title: Conditions in law of large numbers

text: The (Strong) Law of large numbers states that $ \frac{1}{N}\sum_{k=1}^N h(X_k) \rightarrow \mathbb{E}\left[h(X)\right]$ a.s in $\mu$ as $N\rightarrow \infty$. but I can't find any conditions on $h(\cdot)$ other than it needs to be $\mu$-integrable. I am not overly comfortable with measure and integration theory so I might miss the point here but what conditions are set on $h$? The proofs I have found on LLN that doesn't involve measure theory, doesn't have any function $h$ and I don't really know what is included in being $\mu$-integrable. I need this in my thesis for explaining why my Sequential Monte Carlo estimators are consistent. Any reference to a book that has this proof might be really helpful but understanding is my biggest priority here.


title: Binary or Multinomial Logistic Regression?

text: I wish to analyze the following : Predictor Variable (IV): Satisfaction of sexual needs as important (4 items scale and respond based on 4 point likert scale. Sum up to get the item score.) Response Variable (DV): Condom usagae (2 options : never or sometime) **Questions:**   1. Should I use binary logistic or multinomial logistic? (some people tell me to use multinomial logistic but a book said to only use it when the DV has more than two levels, and my DV only has 2 levels - never or sometime).   2. How can I use SPSS to analyse this? I need step by step help.  The SPSS dialog box for logistic regression has three boxes:   1. Dependent : I put in Condom Usage    2. Factor(s): (I am not sure should whether I put Satisfaction of sexual needs as important here?)   3. Covariate : (I am not sure whether I should put Satisfaction of sexual needs as important here?) I am very sorry. Maybe my question sounds silly but i really need help, as i am beginner. I couldn't get any tutor within my town.


title: How to handle unequal sample size with multiple factors?

text: The following dataset shows observations (count data) of cedar trees growing in small, medium, or large forest gaps (i.e. openings). You can see that that there is an unequal sample size and normality is questionable. I would like to test the following hypotheses:   * Ho: Cedar trees grow on all substrates equally among all gap sizes   * Ho: Growth on individual substrates (e.g. rock) does not differ between gap size classes My initial thoughts are that ANOVA would be inappropriate due to unequal sample sizes; contingency tables may be able to handle the structure of the data, although probably cannot address the hypotheses I am interested in testing. Goodness of fit tests seem promising. What statistical methods would be appropriate to address these types of data and hypotheses? Data is in the following form: ![enter image description here](http://i.stack.imgur.com/CpBIF.jpg)


title: How can I reach the plot of outlier effetc of a TC outlier?

text: In this article there is a plot of the outlier effect. Detecting Outliers in Time Series (LS/AO/TC) using tsoutliers package in R. How to represent outliers in equation format? How can I have this plot for my tc outlier?


title: Hidden Markov Model segmentation of different proportions of binary data

text: I need to segment a sequence of 0s and 1s by their proportion at relatively large scales. As an example, let's define 5 different states that represent 5 different ratios of 1s & 0s.               Alphabet: 1 and 0          State      Definition                  emission prob.     state 0: 100% zeroes and 0% ones      0:0.999   1: 0.001     state 1: 75% zeroes and 25% ones      0:0.75     1: 0.25     state 2: 50% zeroes and 50% ones      0: 0.5     1: 0.5     state 3: 25% zeroes and 75% ones      0: 0.25    1: 0.25     state 4: 0% zeroes and 100% ones      0: 0.001   1: 0.999      Attempts: With all the transition probabilities that I've tried so far and the emissions of each state, the output of my model sequences of states is only either state `0` or state `4`. Example: data (binary): > > 00000000000001111111111110000000000101010101010101010000001000100010011001000010 output I get no matter how I change the transition probs. (in states): > > 00000000000004444444444440000000000404040404040404040000004000400040044004000040 output I need (in states): > > 00000000000004444444444440000000000222222222222222220000001111111111111111111111 I have the impression that I am missing some basic theory rather than an implementation problem. For instance, I smoothed the data by aggregating obtaining the ratio of `1` vs `0` in a arbitrarily defined window, and in this way I can see the intermediate states between state `0` and state `4`. Nevertheless, I don't want to smooth the real data as I need to justify then the smoothing window size. Is using HMMs a good solution for this problem?


title: Can a model of P(Y|X) be trained via stochastic gradient descent from non-i.i.d. samples of P(X) and i.i.d. samples of P(Y|X)?

text: When training a parameterized model (e.g. to maximize likelihood) via stochastic gradient descent on some data set, it is commonly assumed that the training samples are drawn i.i.d. from the training data distribution. So if the goal is to model a joint distribution $P(X,Y)$, then each training sample $(x_i,y_i)$ should be drawn i.i.d. from that distribution. If the goal is instead to model a conditional distribution $P(Y|X)$, then how does the i.i.d. requirement change, if at all?   1. Must we still draw each sample $(x_i,y_i)$ i.i.d. from the joint distribution?   2. Should we draw $x_i$ i.i.d. from $P(X)$, then draw $y_i$ i.i.d. from $P(Y|X)$?   3. Can we draw $x_i$ not i.i.d. from $P(X)$ (e.g. correlated over time), then draw $y_i$ i.i.d. from $P(Y|X)$? Can you comment on the validity of these three approaches for stochastic gradient descent? (Or help me rephrase the question if necessary.) I would like to do #3 if possible. My application is in reinforcement learning, where I'm using a parameterized conditional model as a control policy. The sequence of states $x_i$ is highly correlated, but the actions $y_i$ are sampled i.i.d. from a stochastic policy conditioned on the state. The resulting samples $(x_i,y_i)$ (or a subset of them) are used to train the policy. (In other words, imagine running a control policy for a long time in some environment, gathering a data set of state/action samples. Then even though the states are correlated over time, the actions are generated independently, conditioned on the state.) This is somewhat similar to the situation in this paper. I found a paper, Ryabko, 2006, "Pattern Recognition for Conditionally Independent Data," which at first seemed relevant; however, there the situation is reversed from what I need, where $y_i$ (the label/category/action) can be draw not i.i.d from $P(Y)$, and $x_i$ (the object/pattern/state) is drawn i.i.d. from $P(X|Y)$. **Update:** Two papers (here and here) mentioned in the Ryabko paper seem relevant here. They assume the $x_i$ come from an arbitrary process (e.g. not iid, possibly nonstationary). They show that nearest-neighbor and kernel estimators are consistent in this case. But I'm more interested in whether estimation based on stochastic gradient descent is valid in this situation.


title: Confidence intervals when using Bayes' theorem

text: I'm computing some conditional probabilities, and associated 95% confidence intervals. For many of my cases, I have straightforward counts of `x` successes out of `n` trials (from a contingency table), so I can use a Binomial confidence interval, such as is provided by `binom.confint(x, n, method='exact')` in `R`. In other cases though, I don't have such data, so I use Bayes' theorem to compute from information I do have. For example, given events $a$ and $b$: $$ P(a|b) = \frac{P(b|a) \cdot P(a)}{P(b)} $$ I can compute a 95% confidence interval around $P(b|a)$ using $\textrm{binom.confint}(\#\left(b\cap{}a),\#(a)\right)$, and I compute the ratio $P(a)/P(b)$ as their frequency ratio $\#(a)/\#(b)$. Is it possible to derive a confidence interval around $P(a|b)$ using this information? Thanks.


title: What hypothesis test to use for categorical variables?  Possibly in R?

text: Edit: I think this is a better question, Say, I have categorical characteristics such as gender, race. How should I use Fisher's test and chi- square test? I was looking at this: http://math.hws.edu/javamath/ryan/ChiSquare.html end edit. I want to do a hypothesis testing for cancer patients. Everyone is admitted to the hospital, a small portion went developed cancer, while the rest didn't. The dependent variable is categorical, not numerical, meaning something like severe, mild, and unknown. How should I give a hypothesis test? Is it possible? Thanks! Edit: The data analyst used Fisher's Exact or Chi-squared test. I don't know how I can do that, R asks the 2 test groups of equal length


title: How to arrive at a specific formulation of the relative median deviation?

text: I am an economist currently working with this book: Frank Cowell - Measuring Inequality On page 25 a formulation of the relative mean deviation is given as follows: $$ M = 2 \left[ F\left(\bar{y}\right) - \Phi(\bar{y}) \right] $$ $F$ is the CDF, $\Phi$ is the proportion of total income received by persons who have an income less than or equal to $y$ ( per the book's definition: $\Phi=\frac{1}{\bar{y}} \int_0^y zdF(z)$), and $\bar{y}$ is the mean. All this is also defined on page 152 in the appendix. The appendix also gives a definition of $M$: $$ M = \int \left| \frac{y}{\bar{y}} -1\right|dF $$ The book says that the former formulation can be derived from the latter, but I have no idea how to begin with this. How do I perform the integration here and get to the first formulation?


title: Removing parameters in AMOS when $\chi^2$ value and DF are 0

text: I am VERY new to SEM and AMOS. I ran a model where I obtained 0 DF, so it could not compute the fit. So, I am left to cut out parameters. HOWEVER, I'm not sure where I can cut them in my particular case. I have attached here a screenshot of my model...you can see there are not too many parameters. I am trying to test to see if Emotional Models of Attachment (measured by early relationships) affect level addiction (measured in the survey).   1. I'm not sure how I can cut down parameters without losing the whole meaning. For example, I can cut out the latent variable, and just look at the data as a path analysis from the early relationship measures to Addiction... but that would eliminate the theory that Attachment Models are involved, right?    2. Also, in researching more about constraining, I read that I can set the regression (number by the arrow) at 1. But in my case, this would sort of ruin the whole point of running my model, correct? When is good practice to constrain from a latent variable? Sorry if these are stupid questions...I'm just new to all this. THANK YOU SO MUCH for taking the time to answer! I wish I could buy you all a beer! ![SEM screenshot](http://i.stack.imgur.com/xFFJO.png)


title: Survfit function for gbm cox model

text: I just fitted a boosted regression coxph model:               cox=gbm(Surv(periods, event) ~ grade + fico_range_low + revol_util + dti, data=notes)      However, I want to obtain the survival curve from the model similar to the `survfit()` function in the `survival` package. Does anyone know how to obtain using the model from the `gbm` package?


title: Paired test: many repetitions, few subjects

text: I have a situation similar to this: 9 normal-sighted subjects are each tested under two conditions: say, with and without spectacles (of a fixed power) and asked to reach towards an object many times in each condition (approximately 5000 repetitions, but varies with condition and subject). The goal is to make statements about the effects of spectacles on the reach error (distance from the object). Let us say that each repetition in a given subject is independent from the other repetitions. a) I could do a paired t-test with 8 degrees of freedom (9 subjects). b) Since the sample-sizes are different in each condition, I looked for the possibility of a weighted paired t-tests, and found an implementation in R for a one-sample t-test in the weights package and also a related weighted mean in the HMisc package. However, I could not find any references discussing the concepts behind weighting by the sample-size. c) I think I could do separate 2-sample tests in each individual subject with approximately 5000+5000-2 degrees of freedom, and report those. But I am wondering what assumption would be needed in order to pool the repetitions across subjects for each condition and then do a simple two-sample t-test between the two conditions: I guess here one assumes that the within-subject variance is much larger than the between-subject variance and then one would have to pretend that each repetition from a single subject is like sampling a different subject. d) As a related question: Does the answer to c) change if we had an unpaired situation similar to the one above, where instead of 9 subjects tested on 2 conditions, we had 18 subjects, with 9 randomly assigned to condition 1 and the other 9 to condition 2. What assumptions would allow one to simply pool all the repetitions from the 9 subjects in each condition into 2 large sets of data that can then be compared with a 2-sample test with a huge number of degrees of freedom ? Thanks !!


title: Determining which theoretical prediction is more in line with actual data

text: I have a question concerning the analysis of experimental data. Imagine the following situation (which I simplify to make the point clearer):   * I ask my participants to give a numerical answer to 3 questions, say, A, B, and C. So in total there are 3 data points per participant.   * Now, I can compute several "expected" answers to question C according to competing theories. To simplify, imagine I consider theories X and Y, which yields two distinctive predictions for C: C-X and C-Y.   * These "expected" answers are computed by combining the actual answers to questions A and B. For instance, let's say that theory X predicts that C equals A+B, while theory Y states that C should amount to the average of A and B (this is just an example, not the real case of my dataset).   * Thus, note that my two "expected" predictions for C are based on a theory, but also on the actual answers given to A and B. This means that I have different predictions for each participant (since I also have different answers to A and B for each participant). And it also means that the predictions are not independent from the actual data of A and B (not the ones I want to model, i.e., answers to C).   * Finally, note also that the theoretical predictions may correlate between them, since they are relatively similar computations on the same values (A and B). My goal: As you may suppose by now, I would like to check whether the actual answers to C support theory X, or theory Y, or both, or none. A relative, rather than absolute, measure of support (e.g., "the evidence is 3 times more likely under theory X than uncer Y") would also be fine, but this sounds more like model selection and I'm not sure that this is what I want. Would you help me decide the type of data analysis that suits better this situation? I can easily come up with simple solutions (e.g., a correlation model), but maybe they have drawbacks or there are more informative and sofistified alternatives... Any suggestion is welcome and I'm willing to learn! EDIT: To further clarify, A, B, and C can be assumed continuous.


title: Comparing statistics of networks of different sizes

text: Hi This maybe a basic Stats question. Let say I have 3 networks of different sizes. Size in terms of number of nodes and links. Network n1, n2 and n3 have v1, v2 and v3 nodes and l1, l2 and l3 links. That is all we know about the networks. For each of these networks we compute some parameters based on some criteria. For example, param1 = no. of nodes with criteria 1 param2 = no. of nodes with criteria 2 param3 = no. of nodes with criteria 3 Result                  nodes links param1 param2 param3      n1   v1    l1     a      b      c        n2   v2    l2     d      e      f        n3   v3    l3     g      h      i      Question If I average (a+d+g/3) to get the expected param1 across networks, would that make sense knowing that the size of the networks is different and the parameters are a function of the size of the network? How to reconcile with the fact that different networks have different sizes and therefore cannot be compared directly.


title: Classification Tree Interpretation

text: I know how to do analysis in Decision tree but I am not sure how the data is classified in the output. Is there any criteria there in decision tree for the output? Can I assign the output at what basis we need the data. Please help me to know this.


title: Computing event rates given RR + CI and total sample size in each treatment group

text: I am looking at some data for the risk of mortality in patients undergoing treatment A vs treatment B and I am given the total number of patients in each treatment arm and the relative risk + confidence intervals of mortality. How would I go about finding the actual number of deaths in each treatment arm? Here are the hard numbers:               Study 1     Adjusted total # of patients for treatment A: 8366     Adjusted total # of patients for treatment B: 10251     Adjusted RR for mortality: 0.72 (0.66-0.78)          Study 1     Adjusted total # of patients for treatment A: 23113     Adjusted total # of patients for treatment B: 26819     Adjusted RR for mortality: 0.78 (0.72-0.83)      So I would like to compute the actual number of patients who died in each treatment arm for both studies 1 and 2. Not sure if this helps but the paper states that: > adjusted survival curves were estimated with the use of the inverse- > probability-weighting approach of Cole and Hernan. For each treatment group, > the survival curves adjusted with the use of inverse probability weighting > represent the expected rate of survival if the treatment of interest were > applied to all study patients. Using estimated rates of survival among > patients undergoing A and among those undergoing B, we calculated risk > ratios at specific time points and used bootstrap methods to obtain 95% > confidence intervals. I had a look at the Wikipedia article for RR (http://en.wikipedia.org/wiki/Relative_risk) but its not the most helpful. Thanks!


title: Replacing categorical variables with historic response rate

text: In Linoff and Berry's "Data Mining Techniques" they mention reducing the number of categorical variables in a classification model by replacing the variable with the historic response rate. "When building model sets for directed data mining, a powerful transformation is to replace categorical variables with the historical measure of what you are trying to predict. So, historical response rate, historical attrition rate, and historical average customer spend by ZIP code, county, occupation code, or whatever are often more powerful predictors than the original categories themselves." Is this valid from a statistical perspective? Are there any papers that discuss this technique?


title: Finding a function for data points with disturbances

text: I have a question regarding how to find an equation for data points that have disturbances. For example, I have the following data:   * P = 366 measured values    * T = 366 measured values   * t = [ 1 : 366 ], representing the days of the year (index) So at each t (day), we have value of P and corresponding value of T. When drawing P and T vs. t, it shows a continued curve for P with disturbances. These disturbances are caused by the change of T and it is clear that P is mainly affected by T, meaning that: P = f(T) This figure is shown here: ![http://store2.up-00.com/Sep11/VKq30182.jpg](http://i.stack.imgur.com/w8VGd.jpg) This figure shows P vs. T: ![http://store2.up-00.com/Sep11/esZ30295.jpg](http://i.stack.imgur.com/NPuhd.jpg) My aim is to find this relation between P and T or in another words: P = f(T). What I have tried and think so far is first to smooth the curve P and from the smoothed points and smoothed curve, some function may be established. Then, including the effect of disturbances (up and down) by some function, may be exponential. Could you please guide me how can I manipulate this problem to get a final model P = f(T).


title: Mediation analysis with curvilinear mediator and categorical predictors

text: I am currently conducting analyses where I wish to investigate how one variable (M) mediates the association between SES (3-level categorical variable) and a measure of health (continuous variable). The association between M and health is curvilinear, meaning that both low and high levels of M is associated with increased health problems. Any suggestions on how to go about achieving this? I have studied the literature on mediation, but haven't been able to find an approach which seems suitable.


title: Estimate the Kullback-Leibler divergence

text: I would like to be sure I am able to compute the KL divergence based on a sample. Assume the data come from a Gamma distribution with shape=1/.85 and scale=.85. > set.seed(937) > > theta <- .85 > > x <- rgamma(1000, shape=1/theta, scale=theta) Based on that sample, I would like to compute the KL-divergence from the real underlying distribution to an inverse gaussian distribution with mean 1 (mu=1) and precision 0.832 (lambda=0.832). I obtain KL = 1.286916. Can you confirm I compute it well ? Thank you !


title: Simulate Arellano-Bond

text: I have fitted a dynamic panel data model with Arellano-Bond estimator in gretl, here is the output:               Model 5: 2-step dynamic panel, using 2332 observations     Included 106 cross-sectional units     H-matrix as per Ox/DPD     Dependent variable: trvr                        coefficient   std. error     z       p-value        ---------------------------------------------------------       Dtrvr(-1)    0.895381     0.0248490    36.03    2.55e-284 ***       const        0.0230952    0.00226823   10.18    2.39e-024 ***       x1          -0.0263556    0.00836633   -3.150   0.0016    ***       x2           0.127888     0.0171532     7.456   8.94e-014 ***          Sum squared resid    605.9396   S.E. of regression   0.510180          Number of instruments = 256     Test for AR(1) errors: z = -4.29161 [0.0000]     Test for AR(2) errors: z = 1.62503 [0.1042]     Sargan over-identification test: Chi-square(252) = 105.139 [1.0000]     Wald (joint) test: Chi-square(3) = 2333.35 [0.0000]      I have 2 questions about the results:   1. How do I assess the fit?   2. How can I simulate from the model? Regards


title: How to calculate mean absolute error using cross-validation?

text: When I use 5-k cross validation, is the mean absolute error (MAE) equal to $$\text{mean}(\text{MAEs of each 5 steps})=\frac{1}{5}(MAE_1 + \cdots + MAE_5)$$ or equal to $$\text{mean}(\text{absolute errors of all predictions})$$?


title: ANOVA shows no significant interaction effect, but Tukey's Test shows significant differences between groups?

text: I am a PhD student but I am having a hard time understanding certain aspects of statistics. I just ran an ANOVA and I got no significant interaction effect between Density (0 and 20) and Species (species A and species B). However, Tukey's test shows that there is a significant difference between these variables. Species A at a density of 0 is significantly larger than the other combinations. Is this relevant given that ANOVA has already established that there is no significant interaction effect between Species and Density? Why would there be significant differences in Tukeys if ANOVA says there is no interaction effect in the first place?


title: Combining independent Gaussian probabilities

text: I am using three Gaussian distributions with which I generate random numbers to represent many candidate xyz points. I use some selection criteria (details not particularly relevant) to decide on the final position, but I would like to know what the probability of my final xyz point is, and indeed what the probability of getting something "better" (with respect to the mean values of the Gaussian) would be. The three Gaussians (corresponding to my x,y and z coordinate respectively) look like:   $ 1: \mu = -3.89764, \sigma = 0.0838747$   $ 2: \mu = 1.29137, \sigma = 0.0176682$   $ 3: \mu = -5.67947, \sigma = 0.0462101$ and the final position my algorithm settles on is: -3.87008, 1.27922, -5.71548 Currently my approach is to use the cumulative distribution function to calculate the fraction of the $x$ Gaussian between -3.9252 and -3.87008, which corresponds to the distance from the mean to the final position in the $x$ Gaussian on both sides of the mean, and I calculate this to be 0.257531. Similarly for the $y$ and $z$ Gaussians I calculate: Fraction of y between 1.27922 and 1.30352 is 0.508343   Fraction of z between -5.71548 and -5.64346 is 0.564178 If I multiply these together I end up with a value of 0.07385, which I take as the probability to settle on a position as good (or better) with respect to the three mean values of the input Gaussians. Is this procedure correct? Would that suggest that there is ~0.93 probability to settle on a worse result? There is tension between this and the fact that if I multiply the inverse of my above fractions together (eg (1-0.257531) * (1-0.508343) * (1-0.564178)) I get a value of 0.159. I expect this could be as I'm only considering the case where all xyz positions are "at least as good" or all are "worse" with respect to the means, and neglecting the cases where some are better and some are worse. Any clarity that you might be able to provide would be much appreciated! Ultimately I would like to calculate a p-value for the final position.


title: IV and Exogeneity Tests

text: Read something about the Hausman test that didn't sound right in some grad course handout online. It stated that the Hausman null of OLS and IV not being statistically different, if not rejected, means there is no evidence of endogeneity. It also states that under the null it is assumed that OLS is consistent+efficient and that IV is consistent+inefficient. That doesn't sound right to me. If OLS is biased and IV is correctly identified, wouldn't rejecting the null be evidence of no endogeneity? And why would we assume OLS is consistent for the null?


title: Can we use a coefficient of variation as a statistic for testing homogeneity of variance?

text: There are well-known HOV tests, for example Levene's, Brown-Forsythe's, Barlett's available in SAS and R. But $C_v$ also deals with dispersion, can we use it to understand the degree of HOV? July 15, 2014 UPDATE: For example, if $C_{v_1}$ > $C_{v_2}$, where $C_{v_1}$ is the coefficient of variation for sample 1 and $C_{v_2}$ is the coefficient of variation for sample 2, can we say that the assume that samples 1 and 2 are differently dispersed? P.S. by "dispersion" I mean variance (standard deviation).


title: Marketing Mix Modeling in R

text: I have 24 months of sales data by customer and when these customers were touched by different marketing channels (email, coupon, etc). I am trying to build a Marketing Mix Model to show me the baseline sales and the lift by different channels. I know SAS can do this but I am wondering if this can be done in R. Anybody knows what package in R can be used?


title: Why use vector error correction model?

text: I am confused about the **Vector Error Correction Model** ( **VECM** ). Technical background:   **VECM** offers a possibility to apply **Vector Autoregressive Model** ( **VAR** ) to integrated multivariate time series. In the textbooks they name some problems in applying a **VAR** to integrated time series, the most important of which is the so called spurious regression (t-statistics are highly significant and R^2 is high although there is no relation between the variables). The process of estimating the **VECM** consists roughly of the three following steps, the confusing one of which is for me the first one:   1. Specification and estimation of a **VAR** model for the _integrated_ multivariate time series   2. Calculate likelihood ratio tests to determine the number of cointegration relations   3. After determining the number of cointegrations, estimate the **VECM** In the first step one estimates a **VAR** model with appropriate number of lags (using the usual goodness of fit criteria) and then checks if the residuals correspond to the model assumptions, namely the absence of serial correlation and heteroscedasticity and that the residuals are normally distributed. So, one checks if the **VAR** model appropriately describes the multivariate time series, and one proceeds to further steps only if it does. And now to my question: If the **VAR** model describes the data well, why do I need the **VECM** at all? If my goal is to generate **forecasts** , isn't it enough to estimate a **VAR** and check the assumptions, and if they are fulfilled, then just use this model?


title: Hypothesis testing involving pathological distributions

text: Given a very large number of samples $n$ from the population, suppose that I am trying a standard hypothesis test on the means: $$\begin{align}H_0:&\overline{x}\leq\mu_0\\ H_1:&\overline{x}>\mu_0 \end{align}$$ I am interested in the power of my test as I increase $n$. If I am given nice probability measures associated with null and alternate hypotheses (s.t. Gaussians), my task is pretty easy: I can readily characterize the power as a function of $n$ and my significance level (and show how it increases as you add more observations). Now, suppose that the probability measure associated with $H_1$ is pathological in that it has no mean (e.g. Cauchy distribution) or the mean is infinite (e.g. appropriately-parametrized Levy distribution). The probability measure associated with $H_0$ still has a finite mean. I believe this hypothesis test can be written down as follows: $$\begin{align}H_0:&\overline{x}=\mu_0\\ H_1:&\overline{x}~\text{infinite or undefined} \end{align}$$ Is there a way to perform such a test, supposing that my number of observations from the population is very large ($n\rightarrow\infty$)? I think there should be some way, considering that as one increases $n$, if the underlying distribution has a mean, then the average should converge to it by Law of Large Numbers. And, I am assuming that if if the underlying distribution has no mean, the average doesn't converge to anything. However, can one make a precise statement about the power of a statistical hypothesis test in terms of $n$ in that case? What would a test that separates these two hypotheses look like? I wonder if one can prove that in such a case there exists a test that has power of 1 for a finite sample size... **EDIT** : Attempted to clarify what the second hypothesis test might look like.


title: What are Hommel Hochberg corrections?

text: I have recently been introduced to to Hommel Hochberg corrections. I am trying to find a simple explanation about what this actually is/does, but am having no luck. Can anyone please give a brief and simple description about Hommel Hochberg corrections?


title: Must there be "an effect to be mediated" in mediational analysis (i.e., must IVs & DVs be correlated)?

text: Baron and Kenny outlined several steps to aid in determining if a mediational analysis is appropriate to test a particular hypothesis. The very first step was "Show that the initial [independent] variable is correlated with the outcome [dependent]". This is referred here as an "effect to be mediated" ($X\rightarrow Y$). Since then, however, several authors have highlighted that an effect to be mediated is not necessary (e.g., "Reconsidering Baron and Kenny: Myths and Truths about Mediation Analysis"). At first I found this difficult to conceptualize, but now I understand the rationale. Perhaps it is best to consider it in the context of a mediational model in which the independent variable is positively associated with one mediator and negatively associated with another other, leading to a total effect of near 0 (e.g., $X$ causes $Z$, $Z$ causes $Y$, $X$ has a negative effect on $P$, $P$ has a negative effect on $Y$). Despite the lack of a total effect, mediation could be statistically significant for both mediators. From a statistical viewpoint this makes sense to me, and indeed I have even explained it to a few others; however, I have trouble understanding this from a logical or even philosophical standpoint. Start by considering that mediational models are inherently causal. How can a mediational model be theoretically possible without "an effect to be mediated"? Cause is defined by one thing leading to change in another. If changes in $X$ are not associated with changes in $Y$, how could mediation conceivably be present? In other words, if changes in $X$ lead to concurrent changes in mediators that lead to no effect on the dependent variables, how could this possibly be considered causation? What might be needed is a special instance in which changes in $X$ lead only to changes in certain mediators and leads to a total effect, but this seems like a different topic. Consider this example in a universe in which only four variables exist (i.e., all possible mediators are present):   * Intelligence ($X$) causes perceived need to have a healthy lifestyle ($Z$, mediator 1), which leads to weight loss ($-1Y$)   * Intelligence causes increased attraction to video games ($A$, mediator 2), which leads to weight gain ($+1Y$).   * Increasing or decreasing intelligence does not change weight Thoughts?