title: How to best represent median proportions using stacked bar charts?

text: I have a large dataset with several responses about an organization's labour costs. I would like to represent this data graphically. Typically what we do is provide a stacked bar chart (using R) of the percentages. However, there are some severe outliers in the data, and the mean becomes quite skewed by these outliers; obviously, the median is preferred. However, this becomes far more difficult to represent in a stacked bar chart. Here is an example of the data:                            Mean     Median      Cost A   51.39647 49.5389570       Cost B   15.07913  4.2676470       Cost C    5.98685  0.2238874       Cost D   27.53755 16.9540230      --------------------------------      SUM      100.00000 70.9845144      As you can see, the mean is not a good representation of the median. What I had originally done was make a stacked bar chart of the mean; however, I am not satisfied with this, as I don't believe it's a good representation of the data. I would like to represent the median. I have thought about normalizing the median data, so that, for example, Class A would become 49.54/70.98, but that would probably be even less representative than reporting the means. Using the median, it would be easy to say "The average company spends 4% of its labour costs on Cost B." This works well for reporting one statistic. However, I want to report the statistic as a _fraction_ of _all_ the costs for the _average_ company. Then I can no longer say what I said above, since the medians don't sum to 100. I need a good way to report, "The benchmark spending on Cost A is w%, on Cost B is x%, on Cost C is y%, and on Cost D is z%", so that sum(w,x,y,z)=100 and so that it is well representative of the population. My question is, does anybody know a good way to report this data, preferably in a stacked bar chart showing how the average company divides its spending among these four costs? The big idea is to get an understanding of how the "average" company divides its spending between these four costs.


title: Relationship between test stat and p-value in relation to t-test

text: I was wondering why it is necessarily true that if a test statistic exceeds the critical value of t, then it will also be true that the p-value will not exceed the level of significance.


title: Relationship between inverse gamma and gamma distribution

text: I have the following posterior distribution for $v$ $$f(v)\propto v^{-p/2}\exp\left(-\frac{1}{v}\frac{s}{2}\right)$$ and so clearly $$v\sim\text{Inverse-Gamma}\left(\frac{p}{2}-1,\frac{s}{2}\right)$$ Now can I say that $$v^{-1}\sim\text{Gamma}\left(\frac{p}{2}-1,\frac{s}{2}\right)$$


title: Correlate registered voters with local election results

text: I have a list of voters that participated in my local election with their age, gender, and precinct. I also have a list of the results of this election by precinct. Is there any way to find relationships between the two sets of data?


title: How to test the statistical significance for categorical variable in linear regression?

text: If in a linear regression I have categorical variable... how do I know the stastical signifance of the categorical variable? Let's say the factor $X_1$ has 10 levels... there will be 10 different resultant t-values, under the umbrella of one factor variable $X_1$... It seems to me that the statistical signifiance is tested for each level of the factor variable? No? @Macro: Following your suggestion, I've built the following example: It seems that x3 is useful and must be included in the model, from the below model comparison. But actually that's wrong... Any thoughts on this? Thank you!               n=100         x1=1:n     x2=(1:n)^2      x3=rnorm(n)     ee=rnorm(n)     y=3*x1-2*x2+x3+3+ee     lm1=lm(y~x1+x2+x3)     summary(lm1)          lm2=lm(y~x1+x2)      summary(lm2)          anova(lm1, lm2)          > anova(lm1, lm2)     Analysis of Variance Table          Model 1: y ~ x1 + x2 + x3     Model 2: y ~ x1 + x2       Res.Df     RSS Df Sum of Sq      F    Pr(>F)         1     96  82.782                                       2     97 146.773 -1    -63.99 74.207 1.401e-13 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


title: Regression model with categorical values

text: > **Possible Duplicate:**   >  How to test the statistical significance for categorical variable in linear > regression? As we know we can use linear models for numeric dataset(independent variables are numerical only), but what type model is applicable here when I have numeric + categorical dataset(independent variables are combination of numeric and categorical). for example I have two datasets 1.numeric dataset 2.numeric dataset + categorical dataset               1.numeric dataset (Prediction of price of home)          Independent variables     x1 =  numbers of bedrooms     x2 =  size of home in sq. feet          dependent variable     x3 =  price of home          here     dependent variable is numerical     independent variable is with numerical values               2.numeric dataset + categorical dataset(prediction of web visits)          Independent variables     x1 =  search time     x2 =  search query     x3 =  browser     x4 = country          dependent variable     x3 =  visits          here      dependent variable is numerical     independent variable is with combination of numerical and categorical values      I assume here that for dataset 1 linear model with lm() is applicable, but its not possible for second dataset. can any one suggest best technique for dataset 2 to be implemented with model for prediction.


title: Effective ways to display confusion matrices from different predictors in an academic publication?

text: I want to display the results of two different predictors' performance on a dataset. I have a confusion matrix for each of the predictors' results on the test cases. I want to present these confusion matrices in an academic publication (focusing on the false positives difference between the two). These matrices are 2x2. Which options are best suited for presenting the two confusion matrix information?   1. As a subfigure along with the dataset trajectory figure?   2. Within the caption of the dataset trajectory figure?   3. Two matrices in the equation environment in the discussion of the results?   4. As a list of numbers within a sentence in the discussion?


title: What is the difference between "Unweighted", "Weighted" and "Least Squares" Means?

text: These concepts have always confused me. In terms of DOE, I think the problem arises when the design is "unbalanced". So,   * what are their differences?    * Why do you choose one over the other/reasons?   * In one-way ANOVA, does using (can you?) LS means makes sense or LS means is only for two-way ANOVA and onwards? For ex: In a simple study, say, I am comparing if there is a diff. in height between M and F. I have an unbalanced design, say, I have more M than F (I selected equal sample, but some ppl died and this is what I have to use) Does having LS means means anything in this scenario? Or LS means= unweighted in this case? Can anyone please clarify in a simple and intuitive manner?


title: LME() error - iteration limit reached

text: In specifying a crossed mixed effects model, I am trying to include interactions. However, I get the following error message:               Error in lme.formula(rate ~ nozzle, random = ~nozzle | operator, data = Flow) :      nlminb problem, convergence error code = 1     message = iteration limit reached without convergence (10)      The model has the following: 1\. 3 nozzle types (fixed effect) 2\. 5 operators, each with 3 repeat measures on fuel flow from the 3 nozzle types. I was asked to include the interaction between nozzle type and operator in the model. This is my code for the model:               flow.lme <- lme(rate ~ nozzle, error= nozzle|operator, data=Flow)      Why would I get this error message??


title: How to perform randomization test for categorical data and/or chi-square?

text: I am very newbie in statistics and I my problem is as follows: I have a total set of 20,000 categorical observations (this is my whole dataset, say DATA1) which has a certain proportion of "yes" and "no" Ex: (id-1: yes, id-2: yes, id-3: no, id-4: no, ... id-20,000: yes). From the whole dataset, I have a particular subset of 3,000 observations of my interest (say DATA2) (Ex: id-1: yes, id-2: yes etc.). I would like to test if the proportion of "yes"/"no" in this particular sample is significant higher/lower than the whole dataset. The first approach I suppose to do is a chi-square test, is it right? If so, how the table should look like? Another approach that I am considering (and I think it is the most correct), is to random sample from the whole dataset (DATA1) 3,000 observations for N = 1000 rounds and count for each round the number of "yes" and "no". I suppose that I cannot plot the distribution and get z-scores because the variables are categorical. So, I'm stuck in this point because I don't know how to compare the simulated data with the observed data (DATA2) and finish with the statistical significance or computed p-values. Anyone could help me out on this problem? The easiest way to me is a solution in R. Thank you very much!


title: Best process for randomly drawing an odd and an even number

text: I want to draw a random number from 1 to 100. If this number is even, I want to then draw a random number from the odd numbers between 1 and 100. If the first number is odd, I want to draw a random number from the even numbers between 1 and 100. Are the following three strategies the same? (1) Draw a random number from the urn. Sort all odd (or even) numbers out and then draw a random number from the remaining numbers. (2) Draw a random number from the urn. Draw a second number from the urn and put it back until the second number is not odd (or even). (3) Draw two numbers from the urn and put them back until not both are odd (or even). As a real world example, imagine I want to draw a proband from a pool of probands. I then want to pair this person off with a second proband of the opposite sex. I want the sex of the first proband to be random (e.g. because of sequence effects), therefore I cannot draw first a man and then a woman, or vice versa. I could draw first a man before a woman, then a woman before a man, then again a man before a woman etc. But I want to avoid the predictable regularity of that. Contrary to the notification that pops up when I enter the title for this question, this is not a subjective question asking for opinions or guesses. I am sure the "relative randomness" of the events can be calculated: I would guess that you'd have to calculate the probabilities to draw a specific number under all three conditions, and the best "relative randomness" will be achieved, where this number has the lowest probability. Unfortunately, I don't know how to calculate this. Therefore I ask here. _I did not know the proper tags for this question. Please change them, if you feel others fit better._


title: How to set SMOTE parameters in R package DMwR?

text: For different imbalanced data-sets which rare class' proportion differ from 30% (rare) to 5% (rare), what is the best way to define the `Perc.Over` and `Perc.Under` in `SMOTE` method of `DMwR`? Currently I'm using this formula:               Over = ( (0.6 * COMMON_NO) - RARE_NO ) / RARE_NO          Under = (0.4 * COMMON_NO) / (RARE_NO * Over)      and finally transform it to percentage proportion:               Over_Perc = round(Over, 1) * 100          Under_Perc = round(Under, 1) * 100      and then feed it to SMOTE as following:               newData <- SMOTE(Label ~., df, perc.over=Over_Perc, perc.under=Under_Perc)      The thing is that, if the there is 30% rare data in the training set, performing classification without `SMOTE` gives a better performance than using `SMOTE`! While by having 5%-20% of rare data in the training set, using `SMOTE`, the classification algorithm gives a much more better performance, p.s. I'm using Matthew's correlation and AUC for performance measurement.


title: generate an exact pivotal quantity for theta

text: Given a random sample $x_1,...,x_n$ from $f(x|\theta)=\frac{\theta}{(x+1)^{\theta+1}}$, $x\ge0$ How to generate an exact pivotal quantity for $\theta$?


title: granger causality - using a different lag than the one selected by ACF/PACF

text: I'm testing for causality between two series. PACF and ACF select lags that, when used, do not yield significance in the causality tests. Playing around with the lag I find that I get causality significance for other lags...What are the problems associated with NOT using the lag selected by ACF/PACF? For the record, the range between ACF/PACF values was -3.07 to -3.05.


title: What is the level of measurement for years?

text: The years I won the gold medal in Olympics   1. 1990/1991   2. 1995/1996   3. 1998/1999 I understand that quantitative data can be used on numeric scale. Perhaps, this can. Also I believe this is a ratio data since it could have a natural zero. Am I correct?


title: Statistical analysis of sequential questions in SPSS

text: I am analysing a questionnaire using SPSS, but I do not know how to analyse ranking questions or questions where response from one question is used to answer the next. Can anyone help please?


title: binomial test for testing significance of classification very sensitive to hits count

text: I'm doing a multi-class classification task and I wonder if I'm doing the binomial test correctly since it is very sensitive to the count of correctly classified trials (hits count). Say, there are 20k trials and 3 classes. I calculate the p value as `p(#hits) = 1-binocdf(#hits,20k,3)` Then it seems to me, that the p value is very sensitive to comparibly small changes in #hits compared to total trials of 20k. `p(6700) = 0.31` while `6700/20k = 33.50%` `p(6750) = 0.10` while `6750/20k = 33.75%` `p(6800) = 0.02` while `6800/20k = 34.00%` Is it correct this way and I shouldn't be confused that this somehow is against my intuition or am I doing somethin wrong? So if I assume a 5% significance level, 6750 hits are not signifikant, but 6800 are? I guess this behaviour becomes even more prominent when there are more classes.


title: How do I cross error terms in Stata's xtmixed?

text: Let's say I have an experiment with three within-subject factors, A, B, & C. The data looks like this.                s  a  b  c      1  1  1  1      1  1  1  2      1  1  2  1      1  1  2  2      1  2  1  1      1  2  1  2      1  2  2  1      1  2  2  2      Simple enough. I have 49 subjects. Now, to do this ANOVA in R, I use               m1 <- aov(score ~ a*b*c + Error(subject/(a*b*c)), data)     summary(m1, type=3)     ... (clipped) ...     Error: s:b:c                       Df Sum Sq Mean Sq F value  Pr(>F)        b:c                1  4.608   4.608   8.121 0.00643 **     Residuals         48 27.236   0.567        That looks fine, and matches SPSS's repeated measures GLM. All is well. We can also do a mixed model in R using `lme4` and get the exact same results as this, as well as a mixed model done in JMP.               m2 <- lmer(score~a*b*c + (a*b*c|s), data)     library(car); Anova(m2, type=3, test.statistic"F")     ... (clipped) ...                          F Df Df.res    Pr(>F)      b:c             8.1206  1 48.000  0.006430 **       I can do the same thing in SAS and get similar results.               proc mixed data=mixedexample method=reml covtest;         class a b c s;         model score = a|b|c;          random intercept a|b|c/sub=s;     run;      I can do the ANOVA, by hand in Stata:               anova score a / s|a ///                 b / s#b ///                 c / s#c ///                 a#b / s#a#b ///                 a#c / s#a#c ///                  b#c / s#b#c ///                  a#b#c       I leave off the full interaction error term so it is the residual. Though, for some reason, the `df` of the main effect of `a` is twice the size of the others. But that's not the question I have. My question is: how do I do the `lmer` and `proc mixed` version of the full LMM in Stata? The simple version is               xtmixed score a##b##c || s:, reml      But how do I add the fully crossed error terms to `xtmixed` in the same way I do in `lmer` by adding `+ (a*b*c|s)`? The data is balanced with no missing values, so the LMM should be the same as the repeated-measures ANOVA, right? Why can't I do this in Stata? I barely know the basics of LMMs, but this is a Stata question. I'm just trying to figure out all the different ways of performing these two models. Also, if anyone knows a simpler way of doing the univariate ANOVA in Stata without specifying every single error term by hand? This may not even be the "right" way of doing this procedure, but because I get the same output everywhere else, how do I get Stata to do the same thing as R's `lmer`, SAS's `proc mixed`, and JMP?


title: Inconsistency in cross-validation results

text: I have a set of dataset recorded from subjects as they perform some particular cognitive task. The data consists of 16 channels and a number of sample points per channel and I want to classify this data according to the cognitive task being performed (everything is labelled). The issue is that I do not have a large amount of data (approximately 60 trials per session, 30 for each cognitive task) and I have 2 sessions. I am trying to train a Linear Discriminant Analysis (LDA) classifier to classify this data. The classifier is later to be used in **real-time** to give some form of output every number of samples. I made use of a 5-fold cross-validation to measure the generalization error of my classifier. The problem is that when I run this 5-fold cross validation a number of times, the results I get are not constant at all. Instead there's a significant variation in the overall accuracy (for example, first 5-fold cross validation may yield an average accuracy of 80%, 2nd yields an accuracy of 65%, 3rd yields an average of 72% etc...). Is this normal? If not, what could be the causes?


title: Do components of PCA really represent percentage of variance? Can they sum to more than 100%?

text: When reading O'Reilly's Machine Learning For Hackers, it says each component represents a percentage of the variance. I've quoted the relevant part of the page below. Speaking to another expert they agreed it is the percentage. However the 24 components sum to 133.2095%. How can that be? * * * (Excerpt from chapter 8, p.207) Having convinced ourselves that we can use PCA, how do we do that in R? Again, this is a place where R shines: the entirety of PCA can be done in one line of code. We use the princomp function to run PCA:               pca <- princomp(date.stock.matrix[,2:ncol(date.stock.matrix)])      If we just type pca into R, we’ll see a quick summary of the principal components:               Call:     princomp(x = date.stock.matrix[, 2:ncol(date.stock.matrix)])     Standard deviations:     Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7     29.1001249 20.4403404 12.6726924 11.4636450 8.4963820 8.1969345 5.5438308     Comp.8 Comp.9 Comp.10 Comp.11 Comp.12 Comp.13 Comp.14     5.1300931 4.7786752 4.2575099 3.3050931 2.6197715 2.4986181 2.1746125     Comp.15 Comp.16 Comp.17 Comp.18 Comp.19 Comp.20 Comp.21     1.9469475 1.8706240 1.6984043 1.6344116 1.2327471 1.1280913 0.9877634     Comp.22 Comp.23 Comp.24     0.8583681 0.7390626 0.4347983     24 variables and 2366 observations.      In this summary, the standard deviations tell us how much of the variance in the data set is accounted for by the different principal components. The first component, called Comp.1, accounts for 29% of the variance, while the next component accounts for 20%. By the end, the last component, Comp.24, accounts for less than 1% of the variance. This suggests that we can learn a lot about our data by just looking at the first principal component. * * * Code and data can be found here: https://github.com/johnmyleswhite/ML_for_Hackers/tree/master/08-PCA


title: What is the state-of-the-art regarding sampling from discrete distribution?

text: After struggling with auto-Poisson model (a.k.a. Random Markov Network with conditional Poisson distributions) trying to force Gibbs sampler to obtain discrete sample of the network (since I know conditional distributions) I realized that I have never seen discrete version of this sampler. I wonder what is the state-of-the-art regarding sampling from discrete distribution (not confined to Bayesian posterior sampling)? Is there something like Gibbs sampler to sample from conditional discrete distributions?


title: Measuring change of participation level in a population after a policy change

text: I would like to measure the change of a policy in a community, that is, check whether there was a change after the new policy was applied. The old policy required that people had to give away a right in order to participate in a community (similarly as academic journals work, where authors usually are asked to give away their copyright). The new policy dropped that requirement. The old policy could have refrained some people to contribute (or participate) to the community. The hypothesis would be:   * $H_0$: The change in the policy does not produce any effect in contribution (participation)   * $H_1$: The change in the policy produces an effect in contribution (participation) It seems to me I have to use Chi-square (with a contingency matrix) because I would like to measure changes in the population (contributions or participants). The following aspects makes me hesitate on how to approach the problem:   * There are participants that contribute before and after the policy change. Therefore, the groups before and after are not independents. (Does this matter for the research design?)   * The participation level is not homogeneous, neither before nor after. Because these contributions happen in a period of time (years), before the change there were participants who contributed for a short period of time, others for a longer period and other that still contribute. Questions:   * ¿Is Chi-Square the most appropriate approach for this problem? If not, what would be an adequate approach?   * Does it makes sense to be concerned about the previous aspects listed? If so, how should they be approached?   * Is there anything else I should consider? (something that could invalidate my analysis or something that could improve it) For what it is worth, after analyzing one community and if the null hypothesis ($H_0$) were rejected, I would try to analyze if there are changes in the level of contributions of the active participants in both before and after periods. After that, I would analyze other communities with similar policy changes and compare them (if applicable). Finally, it is an _ex post facto_ study. I have looked into another questions, but they seem different:   * How to detect a significant change in time series data due to a "policy" change? It is about a change in one person behaviour, not level of participation in a population.   * What test should I use to determine if a policy change had a statistically significant impact on website registrations? It might be, but the problem measure there is only registrations, not actual participation or contributions.


title: Multiple Hypothesis testing correction on non-independent data

text: Consider a vector $X_i$ with $i \in [1,N]$. Let us say that you divide the vector into two $X_{i}'$ with $i \in [1,j]$ and $X_i''$ with $i \in [j+1, N]$ for several $j$; then for each pair of sub-vectors associated with this $j$ you apply a test that delivers a p-value (like a Student's t-test or a Kolmogorov-Smirnov test). As you change your separating index $j$ systematically and repeat the test, you are comparing datasets that overlap strongly. Is there a way to easily calculate the number of _independent tests_ as required by the Sidak or the Bonferroni corrections, so that one can apply one of these? One can assume the $j$ index covers a consecutive series of indexes $j \in [lim_1, lim_2]$. Thank you very much.


title: How to control for a variable

text: The goal of the analysis is to find a correlation between birth context and stature, while controlling for age. Each person has an ID number, age, stature, and birth context. That last part is what is giving me trouble. There are three type of birth context: foreign, which means they were born in another country, immigrant, which means they were born within 9 years of their mother moving to the US, and established, which means that they were born after 10 years of the mother being in the US. I'm not sure how to control for age in this situation. Would I have to make separate dummy variables for each age, such as age 10, in the established group? Anything that points me in the right direction is appreciated.


title: How do you Interpret RMSLE (Root Mean Squared Logarithmic Error)?

text: I've been doing a machine learning competition where they use RMSLE (Root Mean Squared Logarithmic Error) to evaluate the performance predicting the sale price of a category of equipment. The problem is I'm not sure how to interpret the success of my final result. For example if I achieved a RMSLE of $1.052$ could I raise it the the exponential power $e$ and interpret it like rmse? (ie. $e^{1.052}=2.863=RMSE$)? Could I then say that my predictions were $\pm \$2.863$ on average from the the actual prices? Or is there a better way to interpret the metric? Or can the metric even be interpreted at all with the exception of comparing to the other RMSLEs of other models? Thanks for all your help!


title: what if I know my time series are cointegrated over a long period but not over a short period?

text: I am regressing time series on time series. I have tested for cointegration on the entire time sample (3 years) and the series are cointegrated. I need to make a rolling window of regressions (to get a time series of betas) over this entire time sample. My rolling window will probably be 50 or 100 days. What if my time series on this short period will not be cointegrated (even if I know it's cointegrated over the entire period)? Can I use a specific estimator (GMM or other) to estimate this regression? What is the bias on my estimated beta if I would regress using OLS?


title: 3-way measure of association, nominal variables

text: Cramer's V does a great job of capturing the degree of association between two nominal variables. What measure is best suited for capturing a similar association between three nominal variables?


title: Establishing that the population sampled of two separate surveys is the same

text: I have two surveys of two separate populations (I don't know that they are necessarily distinct, but they are from two different databases) that ask a similar set of questions. Some questions are basic demographics (e.g. age, income), while other questions are a bit more detailed or about their opinions (e.g. brand preferences, spending habits). How do I prove statistically that the two populations are "the same," or at least comparable? I know that I can do a t-test for individual questions, but is there a way to establish similarity on more than one dimension? The goal is to combine the surveys from these two populations into one series of survey data. For example, we may run survey A every six months, but we run survey B every month, except when we run survey A. I would then like to combine the results from survey A and survey B to have a monthly series of survey data.


title: Log normal distributions - particle sizes in an aerosol

text: Statistics were never my strong point and it's my first question, so please be gentle :) I'm doing some research using Computational Fluid Dynamics, CFD, to model the flow of an oil aerosol through a fibrous filter. The aerosol has a droplet size distribution that is log-normal. The existing code only allows the specification of a single size, I have to add the log-normal support. What I need to be able to do is if someone describes an aerosol flow as having, say, 10,000 droplets, an average droplet size of 425nm and a standard dev of 15, be able to, within my code, calculate that there are: aaa particles of size 100nm   bbb particles of size 200nm   ccc particles of size 300nm   ddd particles of size 400nm   ...   ...   ...   iii particles of size 900nm   jjj particles of size 1000nm when we are looking at particles in the range 100nm - 1000nm with 100nm step sizes. In turn I need to push this information, the particle counts, along with a lot of other stuff into the CFD solver and see what happens to the droplets. The big questions is, how do I calculate the number of particles for each step size. As I said, my stat's is quite limited, so I don't even know if saying _the average is 425nm..._ is the right way to frame the input data. My Masters supervisor has told me one thing, using the pdf, but my reading of the definition of log-normal, yes on Wikipedia, leads me to think he is wrong. Any thoughts on how I can work this out are greatly appreciated. Andrew Updated:   Thanks for the answers, but this is where my minimal stats shines through. OK so we rephrase that we have a GM of, say 328nm, and a GSD of 14.8nm. I then look at the definition of the cumulative dist func. and see: cdf = 1/2 + 1/2 * erf[(ln(x) - mu) / sqrt(2 * sigma^2)] - From wikipedia Am I correct in that mu = GM = 328nm and sigma = GSD = 14.8nm? Then in turn, what do I use for my x value? If I need to get values at 100nm, 200nm, 300nm,....., 900nm, 1000nm, do I use these values, lets just drop the "nm" part for now, as X, one at a time, or do I use a ln/log of them or some other magic number based upon the value in question? I don't need to actually code the erf() as c++ has a function call for it, it's just how I calculate the actual value that I pass into the c++ function that is causing me angst. Am I correct in that when I can calculate the input to the erf function that I calculate the cfd values which will be in the range 0 -> 1 which I then multiply by my sample size to get the number of particles below the point in question. I then just subtract successive values to get the particles in a range - correct/sort of correct/wrong??? Once again, thanks for the help, Andrew


title: Econometrics with proofs

text: Which econometrics book, or any other reference, do you believe provides the most rigorous view of the subject, preferably with proofs and full-blown discussion of mathematical topics? I mean a book that, for example, has the proofs of the Frisch Waugh, Gaus Markow, Cramer lower bound and the like. The thing is that econometrics books, including the handbooks, are mere statements of results with no demonstration whatsoever. One has to look at the original papers for proofs. This holds for recent developments (LATE, MTE etc) or for very established results such as the unbiasedness of the estimator of the error variance.


title: Reconcile difference in answers for conditional probability

text: Disclaimer: I am a mechanical engineer, not a statistics guy, so I sometimes miss things. **Background:** I have a copy of "Applied Probability Models with Optimization Applications" by Sheldon M. Ross (link) and I am working through the content, making sure that I understand it. On page 8 it gives the following as an example of conditional probability: > A prisoner is placed in a cell containing three doors. The first door leads > immediately to freedom. The second door leads to a tunnel which returns him > to the cell after one days travel. The third door leads to a similar tunnel > which returns him to his cell after three days. Assuming that the prisoner > is at all times likely to choose any one of the doors, what is the expected > length of time until the prisoner reaches freedom? It then goes on to say: $$ E\left[ Y | X =1 \right ] = 0 $$ $$ E\left[ Y | X =2 \right ] = 1 + EY $$ $$ E\left[ Y | X =3 \right ] = 3 + EY $$ From this it concludes: $$ EY = \frac {1}{3} \left[0+ \left(1+EY\right)+ \left(3+EY\right) \right ] $$ or $$ EY = 4 $$ I took this, and implemented it in MatLab and ran it a million times. Here is the code:               function conditionalProbEx1          clc;     close all;          y=myfun(1000000);          cdfplot(y)          [mean(y) median(y)]          %my game     function [t]=myfun(n)          if nargin==0         n=1;     end          t=zeros(n,1);          for i=1:n              %initialize t         t(i)=1;         loops=1; maxloops=200;     test=0;              while (test==0)&(loops<=maxloops)             %pick random number             ind=rand;                  %compute length of time for number             if ind<=(1/3)                 t(i)=t(i);                 test=1;                 %break             elseif (ind<=(2/3))&(ind>(1/3))                 t(i)=t(i)+1;             elseif (ind<=1)&(ind>(2/3))                 t(i)=t(i)+3;             end                  loops=loops+1;              end              %iterate          end          return      My values are all integers, but there are a million samples there, so it shouldn't be off much. When I plot it, I get the following: ![enter image description here](http://i.stack.imgur.com/wmRWs.png) The range of cumulative probabilities of remaining in the cell for only four days goes from 48.06% to 60.42%. It envelopes the median. The range of cumulative probabilities of remaining in the cell for exactly three days ranges from 48.06% to 44.38%, or 3.68%. **Question:** Why are you nearly 4x more likely to remain in the cell 4 days than you are to remain in the cell 3 days? Is it because of the "3-days" part of the loop? It feels fractal-ish, if that makes sense - the recurrence of the asymmetry in the step-size. How do I train my mechanical intuition to find this reasonable?


title: Predictive Modeler: How can learning Python and/or Java benefit me?

text: On a daily basis, I build predictive models (namely, logistic regression and credit scorecard models) using fairly large datasets (typically ~500k records and ~1k candidate variables) to predict outcomes such as payment vs. no payment, defaulting on a loan vs. not defaulting on a loan, fraud vs. no fraud, etc. Currently, I use both SAS and R to perform feature selection using methods such as Least Absolute Shrinkage and Selector Algorithm (LASSO) regression, Random Forest, and Gradient Boosting as well as a number of clustering algorithms and principal component analysis. As such, I am constantly searching for more efficient and effective ways of building better models, which has brought me to this post. I am very interested in learning about how Python and/or Java can make my life easier/more interesting in the realm of predictive modeling, and I am hoping that someone can provide me with the following: 1) What are some reasons to learn Python and/or Java to become a better predictive modeler? 2) Are there any references that you would recommend for me to learn more about this topic? 3) Which program would be more useful for me to learn, Python or Java (or some other programming language)? I really appreciate any help and/or insight that you can provide! Thanks!


title: Fitting non-normal data in lme4 with a family distribution

text: I'm currently working on fitting a model where we predict the level of some biomarker as a function of time (see image at bottom). I have two difficulties:   1. Each person contributes 2-3 datapoints because they're repeatedly measured before diagnosis (time=0).   2. The errors are not normally distributed.   3. The growth of the biomaker happens non-linearly  In an ideal world, I would love to do a mixed model with a gamma distribution and include a squared term to account for the non-linearity. For example:                   best.model = glmer(biomarker~time + I(time^2) + (time + I(time^2)|ID),                       data=data, family=Gamma)      This barks at me because it doesn't have enough df. (Because each person contributes three datapoints, a polynomial will fit each person perfectly. See Perfect fitting individual model in mixed model). So then I try to remove the polynomial term:               second.best.model = glmer(biomarker~time + (time|ID), data=data, family=Gamma)      This barks at me too, with the following error:               Error in initializePtr() :        function 'dataptr' not provided by package 'Rcpp'      After some experimentation, it seems I can either model the data with a non- mixed Gamma distribution (ignoring the dependence between datapoints), or fit with a mixed model and Normally distributed residuals. As a side note, I used to model this with a poisson and had no problems until I updated to the newest version of lme4. So, I have a list of questions:   1. Will modeling the data with a Gamma distribution "fix" the fact that the data has a non-constant rate of growth? (Forgive me if this is a dumb question. Generalized Mixed Models are not my forte!)   2. Would a different distribution work better (e.g., beta, log-normal, poisson)?   3. Am I way off in how I think I should analyze these data? Any other suggestions on how to better approach it? Thanks in advance for the help! ![enter image description here](http://i.stack.imgur.com/EFg7K.jpg)


title: What is French data analysis?

text: Some statistical methods - I do not remember if it is principal component analysis or something like that - are sometimes called "French data analysis". What is it exactly ? And some people say that this name is ironic, is it true, and why ?


title: Mean of a product of Gaussians

text: If I have two normally distributed random variables, X and Y, and I want to find the mean of the distribution that results from multiplying them together, which of the following two formulas should I use, and why? (and when would I use the other?):   1. E[XY] = (computational formula for the variance)                = Cov(X,Y) + E[X]E[Y]         2. E[XY] = (product of two Gaussian PDFs)                = (E[X]Var[Y] + E[Y]Var[X]) / (Var[Y] + Var[X])      If I take the example of X = Y, both Gaussian with mean=5 and variance=2, then the two formulas above lead to different answers:   1. E[XY] = 2 + 5*5 = 27   2. E[XY] = (5*2+5*2) / (2+2) = 5 I know I'm missing something - is it that the product of two Gaussian PDFs does not describe the distribution of the product of the corresponding variables? If it doesn't, then what does it (the second formula) describe and how should it used? Is it possible to explain the difference in meaning between the two formulas in plain English? :P Thanks in advance!


title: Hierarchical Bayesian Modeling

text: How does hierarchical bayesian modeling incorporate the correlation present between the observed variables (y1,....,yn) Thank you


title: All nearest neighbors in high dimensional space

text: Suppose I have a very large binary matrix representing $n$ customers and the $m$ products they bought, with $n$ and $m$ both rather large (in the order of millions). The matrix is also pretty sparse. As computing all (or even just the $k$) exact nearest neighbors for every customer is clearly an intractable problem, I'm looking for a way to approximate it. One idea I had is to use minhashing/LSH to first segment the customers into approximately $\frac{n}{k}$ buckets (this can be done by tweaking the similarity threshold I assume) and then simply consider that the $k$ neighbors of any given point are the other members of its corresponding bucket. Another slightly less efficient variant that might produce better results would be to segment into $\frac{n}{2k}$ buckets (i.e. bigger ones), and then search for the closest $k$ neighbors by brute force, inside a single bucket. These schemes would obviously not yield exact results, but I wonder if they might be acceptable nonetheless. I'd also be interested in any other method or idea to tackle this problem.


title: Coding of categorical variables in logistic regression

text: I have to do a binary logistic regression. I have a set of 7 independent variables. 4 of them are binary variables and the other 3 are categorical variables. The categorical variables are divided into 4 levels. For example, mode of payment: cash, internet banking, debit card, credit card. I am confused whether I should make 4 columns with 0-1 dummies or just to code them 1, 2, 3, 4. What would be the right method for these variables


title: Pre-stimulus baseline removal in R

text: I have the following scenario:   * trials were conducted where participants were exposed to multiple stimuli during the course of a trial   * a specific physiological response was continuously recorded during each trial For every set of stimulation phase response data, I would like to apply baseline removal based on the pre-stimulus data, i.e., data gathered about 500ms before the introduction of the stimulus. My questions is: how can I do that in R? Even a purely theoretical description would be very helpful.


title: Wilcoxon rank sum test methodology

text: I'm comparing the performance of two algorihtms $A$ and $B$ based on a metric. For each algorithm, I have $31$ independent samples that represent its performance. This samples are grouped in two sets: $X$ and $Y$ (algoritmh $A$ and algorithm $B$, respectively). For this task, I'm using a Wilcoxon rank-sum test with a level of significance $\alpha=0.05$ on R I have selected this nonparametric test because it makes no assumption about data distribution. This is my methodology:   1. **Hypothesis** Let $E(X)$ and $E(Y)$ be the means of $X$ and $Y$ respectively. Then, the one- tailed test is defined as follows: -$H_0$: $E(X) = E(Y)$   (the performance of both algorithms is similar) -$H_1$: $E(X) > E(Y)$   (the performance of $A$ is better than $B$)   2. **Means** $E(X) = 36.87548$, $E(Y)=37.72585$ `> summary(X) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 35.09 45.34 36.88 46.63 48.05` `> summary(Y) Min. 1st Qu. Median Mean 3rd Qu. Max. 9.332 34.860 41.120 37.730 42.510 48.110`   3. **Compare $A$ vs $B$** `p = wilcox.test(X, Y, alternative="greater", conf.level=0.95)$p.value` If $p=0.0170 < \alpha = 0.05$, then $E(X)$ should be greater than $E(Y)$ (It's what I think, but it's not greater!), We can reject $H_0$, so we can conclude that: _algorithm $A$ outperforms algorithm $B$ with a significance level of $\alpha=0.05$_.   4. **Compare $B$ vs $A$** `p = wilcox.test(Y, X, alternative="greater", conf.level=0.95)$p.value` If $p=0.9835 > \alpha$, then $E(X)$ and $E(Y)$ are similar (Again, It's my guess). We can not reject $H_0$, therefore we can conclude that: _$A$ and $B$ have similar performance_. **Data**               X = c(  45.51885768, 35.65081119, 44.60124311, 15.39979541, 48.05143243, 47.90604081,          7.58163868,  0.00000000, 40.94718019, 45.34194687, 28.55451125, 46.15113458,         48.03542321, 47.91413840, 45.38912357, 47.10730083, 47.42726563, 47.80316539,         34.51956662,  0.05853162, 45.29245167, 48.00199937, 45.28839538, 44.89017125,         45.47222435,  7.32111177, 43.35755055, 45.52413737, 45.53528261, 45.45233121,         3.04515436)          Y = c(  41.603296, 43.005620, 38.345339, 28.483733, 41.342548, 44.344933, 34.030309,         42.012604, 35.718175, 45.203532, 40.482022, 45.345594, 41.155187, 41.141522,         48.111677, 41.117034, 34.713158, 44.972073, 35.091889, 34.206018,  9.332199,         39.776291, 28.236449, 13.792789, 35.016681, 41.699073, 41.203090, 47.988765,         36.496740, 45.346652, 30.186485)      ![enter image description here](http://i.stack.imgur.com/Po8jS.jpg) **Questions**   1. is The Wilcoxon rank sum test based on means (as I think in steps 3 and 4) or is based on medians?    2. Is it possible that $A$ can be better than $B$ and at the same time $E(X) < E(Y)$?   3. Is the following methodology correct?              for x in (X, Y)             for y in (X, Y)             if x == y: continue                 p = wilcox.test(x, y, alternative="greater", conf.level=0.95)             if p < alpha:            x is better than y             else             if p > alpha:                p = wilcox.test(x, y, alternative="less", conf.level=0.95)                 if p < alpha:                x is worse than y             else            x and y have similar performance        4. If I wish to test $A$ against several algorithms $B$, $C$, $...$, what could be the best approach to take?


title: Obtaining a formula for prediction limits in a linear model

text: Let's take the following example:               set.seed(342)     x1 <- runif(100)     x2 <- runif(100)     y <- x1+x2 + 2*x1*x2 + rnorm(100)     fit <- lm(y~x1*x2)      This creates a model of y based on x1 and x2, using a OLS regression. If we wish to predict y for a given x_vec we could simply use the formula we get from the `summary(fit)`. However, what if we want to predict the lower and upper predictions of y? (for a given confidence level). How then would we build the formula? Thanks.


title: confidence intervals for predicted values in linear regression

text: I have made this linear regression model:               mtcars_lm <- lm(mpg ~ drat + hp, mtcars)      Using the effects package, I can predict values of `mpg` for every value of `hp` between 70 and 150, plus get a confidence interval for each value of `mpg`:               library(effects)     as.data.frame(effect("hp", mtcars_lm, xlevels = list(hp=seq(70, 150, 10))))             hp      fit        se    lower    upper     1  70 24.06201 0.9066002 22.20781 25.91622     2  80 23.54415 0.8355302 21.83530 25.25300     3  90 23.02628 0.7691351 21.45322 24.59934     4 100 22.50841 0.7087300 21.05890 23.95793     5 110 21.99055 0.6559716 20.64894 23.33216     6 120 21.47268 0.6128381 20.21929 22.72608     7 130 20.95481 0.5814753 19.76556 22.14407     8 140 20.43695 0.5638509 19.28374 21.59015     9 150 19.91908 0.5612604 18.77118 21.06699      And I can plot the result:               plot(effect("hp", mtcars_lm, xlevels = list(hp=seq(70, 150, 10))))      ![enter image description here](http://i.stack.imgur.com/KEmDi.jpg) My question is: what is the maths behind the calculation of the standard error and confidence interval for every value of `mpg`?


title: Calculating error bars for Excel Linear Regression

text: I've ben sent a forecast of sales from a consultancy. It uses Excel's LINEST function, taking 4 factors that seem to have affected sales in the past, and used them to make a prediction. How do I go about adding error bars to this prediction? I believe LINEST provides a standard error, but using that would give me constant error bars, and surely the error would increase the further inot the future we go? (Please note, I am not a statistician, so be gentle!) Thanks!


title: R multcomp: contrasts for Tukey

text: I want to use the `glht()` function from package `multcomp`. In the `linfct` argument, I provide my linear functions to be tested. I provide contrasts there, such as (for example) `c("Group1 - Group2 = 0", "Group1 - Group3 = 0")` to compare group no 1 against the other groups. I want to do Tukey's test in this fashion. I saw I could do (for example) `mcp(Group='Tukey')` for this, but this does all the pairwise comparisons and I want only specific contrasts. How can I do it?


title: Correlation between ordinal variable an nominal variable

text: I want to find correlation between job satisfaction (by likert-type questionare) and sex. Should I use Kruskal-Wallis test to find correlation between two?


title: Test to match different Likert scaled variables

text: I have a set of answers to a number of 7 point Likert scaled questions that ranged from very bad to very good. I have about 200 answers for each question. > Q1(type1): Very bad, ..., ..., ..., ..., Very good > > Q2(type2): ... > > Q3(type1): > > Q4(type1): > > Q5(type2): > > Q6(type1): > > ... What I'd like to do now is check which type1 question matches with which type2 question in such a way that I will know that a question of type 1 is as bad or as good as a question of type 2. I'm trying to figure out which test to use but I'm stuck. It can't be that I have to run lots of pretests but a variance analysis also doesn't seem right. I'm assuming my answers are interval scaled.


title: Non parametric cointegration test in R

text: I am not very comfortable with the (parametric) assumptions made in cointegration analysis. I have recently found that there is a way to carry out nonparametric conintegration analysis (a google search with the term 'nonparametric cointegration' yields several results). Is anyone aware of an R package that allows for nonparametric cointegration analysis?


title: Variance components estimates

text: In terms of estimating variance components: the random effect and the residual variance, both the Anova estimates and restricted likelihood estimators can be negative and result in zero for estimating the random effect. I am wondering whether there is an estimating approach which can guarantee positive estimators and whether the corresponding estimators can be computed using some existing R functions.


title: Estimate a model by minimising the sum of the one-, two, ... and h-step ahead forecasts?

text: When fitting (stationary) time series models, such as ARIMA models, the standard approach is to minimise the one-step ahead forecasting error, which is equivalent to performing maximum likelihood estimation of the Gaussian likelihood. Sometimes however, the main purpose of a model is for it to make good predictions on a considerable timescale (larger then one step, for example 4 or 12 steps). **I am thus wondering if approaches exist that estimate a model by minimising the sum of the one-, two, ... and h-step ahead forecasts.** I can imagine one would like for such a sum to be weighted, such that one could adhere more importance to 'early' errors. I have searched for quite some time now, but can only find some papers on "Multi-Step Ahead Estimation methods", which minimise (only) the h-step ahead forecasting error. Maybe I am using the wrong search keywords? Besides the fact that the approach I describe here seems intuitive to me, I also feel it could be used to put more emphasis on the data trend (and less on the noise) when estimating, potentially offering an alternative to an initial (and sometimes undesirable) denoising/smoothing of the data.


title: Can I have a bimodal likelihood function to represent a mix of two populations?

text: I came across the following toy example and am lacking a final answer/step to finish the analysis. Imagine a surgery or medical procedure where we don't know the success- probability. It can be thought of as a Bernoulli processes (hit or miss) with unknown theta. Now for the first patient the first try of the surgery is a success. I know how I can calculate the posterior when I start with a uniform prior. But know I want to include a 2nd patient. For him the surgery only works at the tenth trial. So he has 9-misses before the first hit. I also know how I would get a likelihood function from that and with a uniform prior that would be my posterior. But what if my prior is not uniform, but I have strong belief that for parts of the population the operation has a high success-probability but for another part of the population it only works in one out of ten or even twenty attempts. I know how to produce a prior distribution (or grid) that represents that success-probability beliefe (and bimodal distribution, as a mix of the population of easy and hard patients). What I don't know is how to calculate the likelihood function that strengthens that prior. I have two independent observations here (1,1) for instant hit and (9,1). Both have equal weight. If I update my bimodal prior with the first patient the bump of the pdf to the left (low-success-theta) would decrease and the posterior would increase in the middle and to the right. If I then use that posterior as a basis for the second patient the right-hand side of the distribution (high success) would decrease and the middle again would grow. So I was wondering if I actually need a bimodal likelihood function to update my belief. Am I making a mistake in my thinking here or do I miss some Bayesian fundamentals?


title: How and when to use the Bonferroni adjustment

text: I have two questions regarding when to use a Bonferroni adjustment:   * Is it appropriate to use a Bonferroni adjustment in all cases of multiple testing?   * If one performs a test on a data set, then one splits that data set into finer levels (e.g. split the data by gender) and performs the same tests, how might this affect the number of individual tests that are perceived? That is, if X hypotheses are tested on a dataset containing data from both males and females and then the dataset is split to give male and female data separately and the same hypotheses tested, would the number of individual hypotheses remain as X or increase due to the additional testing? Thank you for your comments.


title: Parametric and non-parametric tests?

text: I am trying to do tests to see if two numbers were derived from similar distributions. I am not a statistician so please excuse my use of words. I will try to be as accurate as possible. I have a two sets of numbers. If i do a Student's T-Test (unpaired, two way), i get a p-value that is much greater than 0.05. However, for the same sets of numbers if i do Mann-Whitney U test, or the Kolmogorov-Smirnoff test, the p-value is much smaller than 0.05. Is such a behavior expected. And should i consider the differences a significant or not. EDIT: After reading the answer posted, I guess I have to add one more subquestion, What do these tests test?


title: fastbw with rule="p" in R's rms package: why do results depend on number of covariates?

text: I've been trying to use the fastbw function from the rms package in R to perform logistic regression with backward selection, with p-values as exclusion criterion (I am well aware of the arguments against using p-values for this as opposed to e.g. AIC). However, the results are not in agreement with what I would get if I perform the backward selection manually, as fastbw often drops more factors in comparison. The results also seem to depend on the number of factors considered, even with the option               type="individual".      I created some simple example data in order to prove my point, which give the following result:               > fastbw(lrm(y~x1+x2+x3+x4),rule="p",type="individual")           Deleted Chi-Sq d.f. P      Residual d.f. P      AIC        x3      0.37   1    0.5412 0.37     1    0.5412 -1.63      x1      1.82   1    0.1771 2.20     2    0.3336 -1.80      x4      2.58   1    0.1082 4.78     3    0.1889 -1.22      x2      3.56   1    0.0591 8.34     4    0.0799  0.34          [...]          Factors in Final Model          None      I.e., x2 is dropped as the last of the factors considered, resulting in a model without factors. However, if I consider x2 only, I get the following result.               > fastbw(lrm(y~x2),rule="p",type="individual")          No Factors Deleted          Factors in Final Model          [1] x2      The same is true if I do the backward selection manually, as x2 considered separately has a p-value of 0.045. What might cause this behavior? Since x2 is the last remaining variable in the backward selection, the results shouldn't depend on associations with other model covariates.


title: How would you report (in publication) the results of a linear model fit using the poly function in R?

text: @John recently pointed out to me that R's `poly` function produces less correlated values (more orthogonal) to fit polynomial predictors, i.e. the transformed predictors have a lower correlation with each other than if I just mean centered the 1st order predictor prior to generating the nth order predictor(s). R's poly accomplishes this using an algorithm from Kennedy, W. J. Jr and Gentle, J. E. (1980) Statistical Computing Marcel Dekker, pg 343-344. I checked the book out from the library, but I have to say that I am not quite following what is going on. I've also made some plots comparing these values to the raw values and I sort of can see what is going on, but I still don't I fully understand it. It seems like the data is being processed and subjected to some transformation that involves the mean of the values and some prediction coefficients. I can always use `predict.poly` to turn the values from real values to `poly` values. My question is how should I report the results of a model fit with `poly`? The raw slope coefficients refer to a transformed value, if I report them, then I assume I should also report some of the coefficients from `poly` itself as well. Which ones? In any particular format? Is anybody going to really understand them? Is this a compelling reason to stick to raw mean-centered coefficients?


title: How to compare the accuracy of two different models using statistical significance

text: I am working on time series prediction. I have two data sets $D1=\\{x_1, x_2,....x_n\\}$ and $D2=\\{x_n+1, x_n+2, x_n+3,...., x_n+k\\}$. I have three prediction models: $M1, M2, M3$. All of those model are trained using samples in data set $D1$, and their performance is measured using the samples in data set $D2$. Let say the performance metrics is MSE (or anything else). The MSE of those models when measured for data set $D2$ are $MSE_1, MSE_2, $ and $MSE_3$. How can I test that improvement of one model over another is statistically significant. For example, let say $MSE_1=200$, $MSE_2=205$, $MSE_3=210$, and total number of sample in data set $D2$ based upon which those MSE are calculated is 2000. How can I test that $MSE_1$, $MSE_2$, and $MSE_3$ are significantly different. I would greatly appreciate if anyone can help me in this problem. Thanks in advance.


title: Testing hypotheses suggested by the data

text: It is widely known that one cannot validly test a hypothesis using the data in which the hypothesis-generating feature was observed. So one needs at least two samples, one for generating hypotheses and one to test them. But the two samples need to be similar in some way, eg. the finding that shoe size correlates with intelligence in a database of 5-15 year olds will likely reproduce in samples of this age group from other countries and times, but not samples of 45-55 year olds. I wonder: How different do these two datasets have to be? How similar? Are there general formal criteria for their difference/similarity? Is it valid to examine the same population at different times? Is it valid to partition a dataset into one for hypothesis generation and one for hypothesis testing?


title: PCA scores in multiple regression

text: A while ago a user on R-help mailing list asked about the soundness of using PCA scores in a regression. The user is trying to use some PC scores to explain variation in another PC (see full discussion here). The answer was that no, this is not sound because PCs are orthogonal to each other. Can someone explain in a bit more detail why this is so?


title: Response Range in Random Forest

text: I am using the random forest method (regression) to generate a predictor. Normally, I provide a response/answer along with each feature vector when I train my forest. However, for some of my responses, only a range is know and not the exact value. For example (hypothetically), say, the response is the selling price of a house and the feature vector for each house are things like size, shape, location, etc. What happens if, for some of the houses, only a price range is known (i.e., less than \$10,000) and not the actual value (sold for exactly \$20,000). How can/should I encode this information into the response? Maybe these houses should be omitted from training. But what if all houses that sell for less than \$10,000 are only provided as ranges (i.e., "less than \$10,000", "less than \$4500", or "less than \$90"). If I omit them then my predictor will be poor at assessing houses that sell for less than \$10,000.


title: How important are poststratification weights in regression analysis?

text: I am trying to decide whether I need to use weighted regression or not when estimating regression coefficients. Surely if I include age, gender, region, etc. as predictor variables in the model, then this should matter less but I read conflicting things in the literature. The reason I ask is that some of the models that I am interested in have no weighted estimators available. For example, if I am interested in the effect of education upon income, knowing that women are oversampled, could I just include sex as an extra predictor or do I really need to use weighted regression?


title: Probability that a component belongs to a giant component

text: In a configuration model where each node has a degree 2. What is probability that a node belongs to a component of size n.


title: Significance test for two mean differences? Like a t-test but for comparing differences

text: I'd like to find out if there is a statistical test that I could use to determine if a mean difference is significantly greater than another mean difference. To determine if one mean is different from another mean, I would use a t-test, but that would not apply here because there are not actual distributions; instead, there are differences. Below is a replicable example in R to help explain.               x1 <- c(1,1,1,1,0,0,0,0,1,1,1,1,1)     x2 <- c(1,1,1,1,0,0,0,0,0,0,0,0,0)          x1_mean_diff <- mean(x1) - mean(x2)     x1_mean_diff #sig different with a t test          y1 <- c(1,1,1,1,0,0,0,0,1,1,1,1,1)     y2 <- c(1,1,1,1,0,0,0,0,0,1,1,0,0)          y1_mean_diff <- mean(y1) - mean(y2)     y1_mean_diff #NOT sig different with a t test      Basically, I'd like to see if x1_mean_diff is significantly greater than y1_mean_diff. Is there a way to do this? Please let me know if I'm being unclear. Thanks!


title: How to perform correspondence analysis on ordinal data in SPSS?

text: I have 500 rows and three columns. Two of the variables (that i want to analyze) are in ordinal form and the other in scale. However, when i run the correspondence analysis wizard in SPSS i can only select the scale variable and not the ordinals. How is it possible to perform correspondence analysis on the ordinal variables? I have something like this: Variable1 | Variable2 ..Home . .| Strong ..Work . . | Strong


title: Correlated variables in a math model

text: Let's say you have 8 variables in a regression model. If some of them are correlated, what degree or percent of correlation should you considering removing some of the variables from the equation? How about covariance - what effect does that have? Lastly, how many variables are too many to inform the model? 5, 15, 50, 500? The model I have is a profit/financial type model and is fairly straightforward - profit is the output.


title: Two-state model with transient states for panel data

text: I need to analyse unbalanced panel data where the variable of interest is a dummy with two non-absorbing states (i.e. switch from 0 to 1 and 1 to 0 is possible at any time period). The goal is to predict the probability of being in state 1 for both out-of- time and out-of-sample data and comparing the results to the actual outcome in these data. Using this comparison, I want to identify cases which the actual outcome deviates from the prediction (i.e. outcome of 1 is predicted with a very high probability, but the outcome is 0; or outcome of 1 is predicted with a very low probability, but the outcome is 1). My approach would be to use survival analysis (e.g. Cox proportional hazards or random survival forests), but I am not sure this is the right way, since in survival analysis one state is usually an absorbing state, which probably induces the need for a different statistical approach. Questions:   1) Which approach would be correct to reach my goal (survival, some special case of survival, or something different)?   2) Can you recommend some theoretical background readings on the subject?   3) Can you recommend some practical examples/guides on using R to accomplish this goal?


title: Meaning of 'Data Adaptive' and 'Non Data Adaptive' time series representations

text: I was reading a paper and came across a tree-figure which split time series representations into "Data Adaptive" and "Non Data Adaptive" representations. In the Data Adaptive branch were representations like SVD, Sorted coefficients, piecewise polynomial etc. While in the "Non Data Adaptive" branch were representations wuch as wavelets, fourier/cosine transforms, random mappings etc. So I was wondering, what is meant by 'Data Adaptive' and 'Non Data Adaptive'?


title: Order Statistics

text: What is the motivation behind the use of order statistics in parameter estimation. In a very general sense, the first order statistic is considered to be an initial estimate to the location parameter. I wonder, how this would be possible. A layman based explanation would help.


title: 3 salesman probability

text: Three salesmen a, b, c make house-to-house sales. Their probabilities of making a sale if they all go out for a sales call are $\frac{5}{8}$, $\frac{3}{8}$, and $\frac{1}{3}$. What is the probability that: 1) a sale will be made 2) only one of the salesmen will make a sale


title: Problem with response optimization with three variables using Response Surface in Minitab

text: I'm intending to do a response optimization of one response, $y$, having three predictor variables, $x_1$, $x_2$, and $x_3$. These variables are coded in the following manner:               A    B   C   y     -1.00000    -1.00000    -1.00000    66      1.00000    -1.00000    -1.00000    80     -1.00000     1.00000    -1.00000    78      1.00000    1.00000     -1.00000    100     -1.00000    -1.00000     1.00000    70      1.00000    -1.00000     1.00000    100     -1.00000    1.00000      1.00000    60      1.00000    1.00000     1.00000     75     -1.68179    0.00000     0.00000     100     1.68179     0.00000      0.00000    80     0.00000    -1.68179     0.00000     68     0.00000     1.68179     0.00000     63     0.00000     0.00000     -1.68179    65     0.00000     0.00000      1.68179    82      0.00000    0.00000      0.00000    113     0.00000     0.00000     0.00000     100     0.00000     0.00000     0.00000     118     0.00000     0.00000     0.00000     88     0.00000     0.00000    0.00000       100     0.00000     0.00000     0.00000     85      What I have tried is (in Minitab) Stat -> DOE -> Define Custom Response Surface And choosing 3 responses, 6 center points, 1 replicate and alpha 1,682 and subsequently choosing Stat -> DOE -> Optimize Response variable and choosing Maximize and lowest 100 and target 118 I got that the maximum for $y$ is 101.711 and D=0.41 but A=0.6, B = -0.11 and C=0.15 which, if I've understood correctly, is invalid for this design. When I've googled all I see is people who merely make surface plots and conclude in what path the response is optimized - however using three variables I do not know how to perform this.


title: MANOVA multivariate normality

text: I'm running some analysis of student data and I'm having some questions regarding multivariate normality assumption of MANOVA. **_Should it be done using values of dependent variables or their residuals_**? Most of the textbooks are not enough clear about this, they just use the term "multivariate normality of DVs". I'm reading " _Discovering statistics using R_ " by A. Field et al. and they used raw scores of DVs. However, on StackExchange and the internet I've found different opinions on this and I'm not sure what to do. I calculated the multivariate normality for my 6 levels of independent variable, using scores of the dependent measures and they were all highly significant.Then, given that all my DVs are counts, I used log2(n+1) transformation to normalize the data and then I've run the same analysis again(I've added 1 to each count to be able to calculate the log of zero count). Finally, I've got p values for each of the groups of independent variable: 1) 0.0911(n=22), 2) 0.0236(n=15), 3) 0.00033(n=6), 4) 0.556(n=21) and 5) 0.0376(n=14), the number of participants in each group are in the parenthesis. The thing is that I'm not sure whether this is sufficiently enough for my MANOVA analysis. Box's M(test for homogeneity of covariance) was not significant, there is no multicolinearity (all my DVs are moderately correlated), only this multivariate normality is giving me problems. I've also read that MANOVA using Pillai trace statistic is robust for small data sets, when there is no violation of homogeneity of variance and when there is a problem with violations of other assumptions, but I'm not sure about this. Beside that small group 3 that has clearly significant p value, other two significant p values seem marginal to me, given the size of the data(p=0.037 and p=0.0236). Is this good enough or there is something else I can try do to fix this normality problems? I've also used Mahalanobis distance to check for the multivariate outliers. When I run it on the all the data, I don't get any outliers, and when I run it per group of independent variable, there are a lot of them, but not a single on in that group with low p value. Also, when I remove one of those outliers that are found, the p values are actually getting much worse. Should I run the test for univariate outliers on all the data and then try to remove those ones or maybe even per group of independent variable? I'm not sure about removing the outliers, especially since the groups are already small and how should I justify that decision.


title: Statistical distribution of a scaled complex Wishart matrix

text: Given that $\mathbf{x}\in\mathbb{C}^{N\times 1}$, where $\mathbf{x}\sim\mathcal{CN}(0,\mathbf{R})$ then $\mathbf{x}\mathbf{x}^\mathrm{H}$ is complex Wishart, i.e, $\mathbf{x}\mathbf{x}^\mathrm{H}\sim\mathcal{CW}(1,\mathbf{R})$.   * If the PDF of $\mathbf{x}\mathbf{x}^\mathrm{H}$ is known, then what is the distribution of $\mathbf{x}\mathbf{x}^\mathrm{H}\mathbf{A}$, given that $\mathbf{A}$ consists of constants, Hermitian, singular and rank-1?   * What happens if $\mathbf{x}\in\mathbb{C}^{N\times m}$, i.e., $\mathbf{x}\mathbf{x}^\mathrm{H}\sim\mathcal{CW}(m,\mathbf{A})$?


title: Estimating absolute difference between two conditions

text: A clinical trial examines the effect of two treatments upon lung function. Patients in each group contribute data at baseline and at the end of the trial. Baseline data are provided for each group and then end of treatment data. ![Data](http://i.stack.imgur.com/QLyge.jpg) I need to calculate the mean change (SD) over the time of the trial for each drug so that I can use the data in a meta-analysis - how should I approach this?


title: Understanding the output of a bootstrap performed in R (tsboot, MannKendall)

text: I have a question regarding the interpretation of the tsboot call in R. I checked the documentation of both the Kendall and the boot package, but am no smarter than before. When I run a bootstrap using e.g. the example in the Kendall package, where the test statistic is Kendall's tau:               library(Kendall)     # Annual precipitation entire Great Lakes     # The Mann-Kendall trend test confirms the upward trend.     data(PrecipGL)     MannKendall(PrecipGL)      which confirms the upward trend:               tau = 0.265, 2-sided pvalue =0.00029206      The example then continues to use a block bootstrap:               #     #Use block bootstrap      library(boot)     data(PrecipGL)     MKtau<-function(z) MannKendall(z)$tau     tsboot(PrecipGL, MKtau, R=500, l=5, sim="fixed")      I receive the following result:               BLOCK BOOTSTRAP FOR TIME SERIES     Fixed Block Length of 5      Call:     tsboot(tseries = PrecipGL, statistic = MKtau, R = 500, l = 5,      sim = "fixed")               Bootstrap Statistics :      original     bias    std. error     t1* 0.2645801 -0.2670514  0.09270585      If I understand correctly, the "t1* original" is the original MKtau, the "bias" is the mean of the MKtau from the R=500 bootstrapped time series, and the "std. error" is the standard deviation of the MKtaus from the 500 samples. I have trouble understanding what this means - this basically tells me that all 500 MKTaus are lower than the original, and that the original t1* is in the range of 3 sd of the bootstrapped MKtaus. Ergo it is significantly different? Or would I say the MKtau for the data set is 0.26 plus/minus standard error? I'm sorry for the lengthy question, but I am a statistics novice and am learning via self-study, lacking someone to bounce this probably really simple problem back and forth with.


title: Adjusted R squared on a holdout set

text: The formula for adjusted $R^2$ is: $$ 1 - \frac{(n-1)}{(n-p-1)}(1-R^2) $$ where $r^2$ is the coefficient of determination, $n$ is the number of points, and $p$ is the number of parameters the model has. If I want the adjusted $R^2$ from a holdout set, is $n$ the number of points in the data the model was trained on, or the number of points in the holdout set?


title: Skewed variable - better log10 or ln?

text: I am building a logistic regression model and one of my independent variable sis heavily skewed. Is it better to use the ln or the log10? Why? And how to correct the skewness of a variable that contains a lot of zero values?


title: Applying Frisch-Waugh-Lovell theorem to IV regression in R

text: I am estimating an instrumental variables linear regression that has a large number of indicator (factor) variables. I don't particularly care about the coefficient estimates on those indicator variables. In Stata's ivreg2 package there is a "partial" option that applies the Frisch-Waugh-Lovell theorem to orthogonalize the dependent and exogenous variables to the indicator variables. After this transformation the indicator variables are not estimated because they do not affect the coefficients on the variables I am interested in. My question is, is there something like this in R? It doesn't have to be part of an IV regression package but I am looking to orthogonalize one set of variables to another set of variables. This seems like something that would have already been implemented. Thanks.


title: Type I error in Negative Binomial GLM

text: In the following R code, the expected type I error should be 0.05 (I think), but it is consistently greater than 0.05. I want to know why.               library(MASS)         nsim <- 1000     alpha <- 0.05     out <- rep(NA,nsim)     for(i in 1:nsim){       x <- runif(20,0,10)       y <- rnbinom(20,mu=3,size=1.285714)       out[i] <- summary(glm.nb(y~x))$coefficients[2,4]     }     mean(out<alpha)


title: What's the algorithm for finding sequences used by TraMineR?

text: I'm working an analysis about finding frequent sequences in a event-state dataset using the R package `TraMineR` (and `arulesSequences` too). In `arulesSequences` the algorithm used to find frequent sequences is the cSPADE algorithm. But what is the algorithm `TraMineR` uses for the function `seqefsub`? Does it use the GSP algorithm or the cSPADE algorithm as well? I couldn't find any information about that in the user guide and other information material.


title: Properly interpret the alpha / beta parameters in the Beta Distribution

text: For quite a while I believed that the proper interpretation of a Beta distribution with $\alpha$ and $\beta$ is: " _what is the most likely $P$ given $\alpha -1 $ success (heads), and $\beta -1 $ of failures (tails)_ ", which also made sense when $\alpha -1 $ and $\beta$ are equal to 1, which means that you don't have any heads / tails, i.e. a uniform distribution. Yesterday, after asking a different question Is the beta distribution really better than the normal distribution for testing the difference of two proportions? , this understanding was undermined. Can anyone help me to better understand the interpretation of the $\alpha$ and $\beta$ parameters? I tried to get help from this post: here but that didn't help either.


title: Theory on discriminant analysis in small sample size conditions

text: I see a similarity between a problem I'm working on and Linear (or Quadratic) Discriminant Analysis when the sample size is smaller than $p+1$. I'm interested in theory bounding the generalization error of either classifier when the data can be assumed to be dense (i.e. all of the data parameters are potentially informative) and the estimated class covariance(s) are singular (and so one has to regularize to get a working classifier). For example it seems to me that, in this setting, pseudo-inversion should (in the generalization error sense) be the worst possible thing to do, as one discards potentially useful discriminating information in the null space of the covariance matrix. Making this intuition more concrete doesn't seem straightforward though. Searches on Scholar using likely-looking strings e.g. "discriminant analysis" AND "small sample size" return thousands of papers, largely from the face recognition literature and, as far as I can see, propose different regularization schemes or LDA/QDA variants. My interest is in any _guarantees_ on the performance of such classifiers in the small sample size case, given that there are plenty of alternative ways in which one can solve the numerical issues. Can anyone point me, please, to some relevant theory?


title: Combine variables that are extremely lightly populated?

text: A similar question to my other question about mixed distributions. Here i have quite a few variables that are populated to less than 5%, many are even populated to less than 1% this 1% would represent 10 of 1000 entries If i was designing a survey for example these groups are too 'fine' in their detail These again are a collection of almost OR events (If variable 1 and 2 are populated 3 - 10 will be zeros) In the real world these are loosely 'types' of food eaten. So they all have a loose link they are important as 'weight'is an output My issue is that combining them all together would loose a lot of information but leaving them separate would seem to include both too many variables with any one not having any descriptive power. Is some form of combination my only real option? Ie, external information saying input variable 1,2,3 can be combined and 4, 5, 6 are similar enough to be combined. This may or may not be possible as there may not be any links like this


title: Hazard model with time varying covariates

text: I am trying to build a discrete time event model. The event occurs once in 20 years period of time. Most of the example i saw use only covariates which are same over period of time. e.g. sex,race etc. But i have observed values of covariates for different time periods.And there are around 10-15 covariates which are observed for 20 years. I want to model all the time varying covariates for each year. e.g. in SAS the option is like: model event=var1 var2 var3 number_of_years Here they have var1 var2 as race, sex etc, which are constant over period of time(time independent). But i have like different var1 values observed for number of years. So for second entity, its again different for number of years. The event occurs once for given number of years in 1/0 form I am using SAS, but can use R as well if that helps. Please let me know, if you have any questions. To clarify further: Please refer to: http://ccsr.ac.uk/qmss/summer/Paris09/documents/SingerWillet.pdf in the last you can see logistic regression code. I want to know how to include predictors like cash value along with dummy variables D01-D12 which are varying for 12 time periods.


title: forecast improvement using Kalman FIlter clearing

text: I have been facing a wall after doing a forecast of wind speed time series data using ARIMA with python. I have result with a nrmse growth going from 2% to 15% and now what I want is to use kalman filter to reduce my forecast errors and have better result. 1/Actually I don't know even after lot of readings what is doing what. My main goal is to improve my forecast results (that I obtained using ARIMA). Can you enlighten me on Kalman Filter for forecasting univariate time series? 2/ I also want to know how I can know if a time series is non linear or not ? The definition on various papers does not allow me to write a code based on them. Thanks for your considerations?


title: Meta-analysis with CMA

text: I wonder if I could ask you all a quick question regarding the input of effect size data into Comprehensive Meta-analysis(CMA)? I have tried contacting via the software portal with no response. I am currently performing a meta-anaylsis evaluating the role of a serum biomarker in the development of nausea. Most of the studies divide subjects into those that are nausea resistant to a particular stimulus and those who are nausea sensitive and report levels of the biomarker pre and post stimulus. For instance in the nausea sensitive group (n=10) it may be 1 pg/mL +/- SD 0.5 at baseline and then 2.5pg/mL +/- SD 0.5 after the stimulus whereas in the nausea resistant group (n=15) it may be 1 pg/mL at baseline +/- SD 0.6 and then 1.2 pg/mL +/- SD 0.8 after the stimulus. In CMA I am trying the enter this effect size data under two groups/correlation > unmatched group pre and post > means, SD pre and post, N in each group but how do I calculate the pre/post correlations and what do I standardize the results by (i.e. change score in SD or post score SD). Any help or advice would be gratefully received,


title: General linear model with interaction term in SPSS

text: Following my question here, I will be most grateful for further assistance to clear my confusion. The confusion is because I have read in many places that when you use ANCOVA in GLM (in SPSS), the effect of the covariate is controlled. I am unsure why it should be controlled in my case below (unless I am wrongly interpreting what a 'covariate' means). Or using the wrong test! **Question 1** **What is the best method to analyse the following problem in SPSS?** I am looking at how the relationship between A and B behaves or differs based on gender. In other words, how does gender impact on the relationship between A and B. A and B are interval variables. A is the independent and B is the dependent variable. Gender, of course, is a categorical variable. (The nature of the variables, as I understand it, makes factorial ANOVA inappropriate, leading to a choice of ANCOVA for this example.) My understanding is that B goes in the dependent variable box, gender goes in the factor box and A goes in the covariate box in SPSS GLM (univariate). **Question 2** _**_*What does the interaction mean?**** I understand I can get an interaction between A and gender, which then means that A and gender have a joint effect on B. If this is significant, does it mean that gender affects the relationship between A and B. **Question 3** **Should I look at the effect of A and gender separately or only focus on the interaction between A and gender?** I am very new to GLM and it is fascinating (if I can get it right!).


title: Johansen's cointegration test and multicollinearity

text: 1. Which variables to use when running Johansen's cointegration test? levels or the differenced?   2. How can I better detect and avoid multicollinearity before running Johansen's cointegration test?


title: Dynamic mean of randomly distributed events

text: The aim is to estimate an error on a stochastic event rate dynamically. This Question points in the same direction, I am interested in a theoretical extension. I read out the event counter second-wise, every black $1$ is a counted event (new events over time, see the plot below). During the measurement I am estimating the event rate, so as more statistics is accumulated, mean event rate (red) should asymptotically become more accurate. ![dyn_mean1](http://i.stack.imgur.com/lPFLq.png) As one can see, the mean value oscillates around true value of 0.5, ![dyn_mean2](http://i.stack.imgur.com/0VJKY.png) even after one order of magnitude more events collected. **Practical question:** How can one calculate the number of events needed to estimate the mean value to a maximum error ($0.5\pm \sigma$)? - **answered (?)here** **Theoretical question:** Can this oscillation be described analytically? Can you suggest further reading? The events are radiation counts, so they are uncorrelated, may by Poisson- distribution applied? **Addendum:** Idealized first approximation - every 10th event is non-zero: ![dyn_mean_reg1](http://i.stack.imgur.com/iRD6Y.png) May be this curve is superimposed with the realer-life example above, is any techniques of partitioning in arbitrary functions applicable here?


title: A question on no. of training examples and decision trees

text: I have a set of around 200,000 training instances. Each training instance consists of an attribute called $duration$, a discrete integer type and a time series of floating-point values, in form of a vector of length=$duration$. So if each training instance is to be considered a vector, it is of a length $duration+1$. The problem is is to build a classifier and while I was working on the problem, I tried the following approach: Because each training vector had a different length equal to the duration attribute, I decided to build different classifiers for different durations. In this way, I could do away with the problem of comparing 2 time series of different lengths. Now I had to divide the data into groups where each member would have the same duration. This led to segregating the data into 20 groups (for different durations) with 10000 instances each. Several have questioned me saying that 10,000 might not be enough to build a model at all. I have tried to convince myself that such a 10,000 may be enough by the following reasoning: Take for example a decision tree. It works by deciding on a particular attribute and building recursively, sub-trees as its children. When the subtrees are created, they are done so by using lesser number of training instances than the number used for training its parent node as the the training instances used to train the subtree is a subset of instances in the training set which have a specific value of the attribute on the parent node as per the branch. In this manner, as a decision tree grows down, the no. of training instances it uses to create a sub-tree reduces, but it doesnt mean that the subtrees down the tree by any means are incorrect or unjustified. Questions: When I divided my data into groups and built for models for each group, didnt I do a similar thing (and implicitly gave a decision tree like structure to the classifier)? Are my colleagues right in questioning my approach? Also note that duration is an attribute which contributes heavily towards learning the target function. By dividing into groups and building a classifier for each group, am I not able to capture the effect of duration on the target? Are there better ways to achieve results in this case?


title: Ridge, lasso and elastic net

text: How do ridge, LASSO and elasticnet regularization methods compare? What are their respective advantages and disadvantages? Any good technical paper, or lecture notes would be appreciated as well.


title: Are real and imaginary components of frequency element of fft correlated?

text: I want to use model-based clustering to classify 1,225 time series (24 periods each). I have decomposed these time series using the fast Fourier transform and selected the harmonics that explain at least a threshold percentage of time series variance for all time series in the sample. I want to do model- based clustering on the real and imaginary parts for each transform element of a give time series because it would potentially save me from having to account for temporal autocorrelation in model based clustering across periods of a time series. I know that each complex element of the fast Fourier transform is independent from other elements, but I do not know if the imaginary and real parts of the output for a given output element are independent. I would like to know because if they were, it would allow me to maintain the default assumption of the Mclust package in R for model-based clustering that the variables analyzed have a multivariate Gaussian distribution. NOTE: The input is real-valued, and I have converted from a two-sided to a one-sided spectrum by removing redundant frequency elements and multiplying the positive frequencies (other than the mean component) by two per the advice I got from another StackOverflow answer here: http://stackoverflow.com/questions/8264530/how-do-i-calculate-amplitude-and- phase-angle-of-fft-output-from-real-valued-in


title: How to simulate from a Gaussian copula?

text: Suppose that I have two univariate marginal distributions, say $F$ and $G$, which I can simulate from. Now, construct their joint distribution using a Gaussian copula, denoted $C(F,G;\Sigma)$. All the parameters are known. Is there a non-MCMC method for simulating from this copula?


title: Mean of one group vs. the rest of the sample

text: I have data from a survey with Likert-type questions (Rate how strongly you feel about this issue on a scale from 1 to 10). We know the respondents home town (as well as some other information, such as gender, age, educational level etc.). We found that the average (mean) score on one of the survey questions was a lot bigger for survey takers from one town than for those from other towns. Here is my questions: How do we test whether it is significantly so?   * Can we test town A vs. all the other towns together?   * Or do we have to test it vs. each other town? (do a Anova with all towns followed by post-hoc tests?)  Any pointers would be much appreciated!


title: Interpretation of results in non-inferiority study

text: I need some help interpreting some results from a non-inferiority study. The primary endpoint for assessment of non-inferiority is the change from baseline in Power of Attention Score. The study is looking to show that the new treatment is no worse than placebo in its impact on a patient's cognitive function. The non-inferiority limit is 121ms. Results from the analyses show a difference between treatment groups of 23.23ms in terms of increases from baseline in power of attention score (Placebo-Treatment), so the new treatment is showing a smaller increase from baseline than placebo. The 95% CI for this difference is (-145.1,191.59). Given that the mean increase from baseline was smaller in the new treatment group and the upper limit of the 95% CI is < 121ms (non-inferiority limit), then I assume we can conclude non-inferiority (i.e. any advantage in placebo is less than non-inferiority limit) as the CI lies to the left of the limit? I am asking this question as I am a little confused on whether to be using the upper bound of lower bound of the non-inferiority interval. The 95% CI does in fact include the lower limit. Also, the CI is very wide so I am wondering how valid this conclusion is anyway based on the number of patients included in the analysis. Thanks for your help with this.


title: Formal definition of random assignment

text: I am looking for a formal definition of random assignment. Let $\mathbf{Z}$ be a vector of treatment assignments in which each element is 0 (unit not assigned to treatment) or 1 (unit assigned to treatment). In a JASA article, Angrist, Imbens, and Rubin (1996, 446-47) say that treatment assignment $Z_i$ is random if $\Pr(\mathbf{Z} = \mathbf{c}) = \Pr(\mathbf{Z} = \mathbf{c'})$ for all $\mathbf{c}$ and $\mathbf{c'}$ such that $\iota^T\mathbf{c} = \iota^T\mathbf{c'}$, where $\iota$ is a column vector with all elements equal to 1. In words, the definition seems to be this: assignment $Z_i$ is random if any vector of assignments that includes $m$ assignments to treatment is as likely as any other vector that includes $m$ assignments to treatment. This definition seems unsatisfactory. What if I decide a priori that I want to rule out a particular vector of assignments, and choose one of the remaining vectors at random? This practice would not satisfy the AIR definition, but it would still be random assignment. Here is an example. Imagine a binary assignment to treatment for each of two subjects. Let $\mathbf{Z}$ be the vector of treatment assignments. Then $\mathbf{Z}$ has four possible values: {0, 0}, {0, 1}, {1, 0}, and {1, 1}. By the AIR definition, assignment is random only if $\Pr(\mathbf{Z} = \\{1, 0\\}) = \Pr(\mathbf{Z} = \\{0, 1\\})$. But why should this be the definition of random assignment, or even a necessary condition for it? What if I simply decide that I want to rule out {0, 1} and choose at random from the three remaining vectors? It seems that this practice is consistent with conventional understanding of random assignment but inconsistent with the AIR definition. So: is there a formal definition of random assignment that encompasses the idea that the experimenter may rule out some assignment vectors a priori?


title: Probability puzzle about zombies

text: I am thinking about writing a simple game about zombies. I got stuck trying to calculate how many people should become zombies. **Here are my conditions:** We have a small rural town of 700 people. One night 200 zombies come to town. Each zombie has a single 30% contact chance (per day) with 50% chance of infection of one human and turning him into a new zombie. How many people shall become zombies on the **first** , **second** and **third** day? I decided to leverage the chance of infection according to amount of people in town. The more people, the less chance for contact (which may look controversial, yet let it be so). So `(200/700)*0.3=0.086` or `8.6%` contact chance. Then I started off thinking that one zombie has `8.6%` chance of infecting one human and during one day there shall be 200 tries per 200 humans only. I tried to find a number of tries at which the probability of infection is 99.9% and it appears to be 60. Then every ~60th human shall be infected which means `700/60=11.6` people per day. I don't have any background in statistics so I am probably wrong, it would be nice if anyone would help me.


title: Simple similarity comparison over time

text: I've seen that Pearson's r is typically used for evaluation between non- experimental conditions. However, I've also seen Pearson's r used to indicate effect size - is this acceptable? Is it a useful descriptive statistic at least? For example, there are two groups of students. In one group students are put into pairs to do a task, in the other they aren't. (This is a task whereby improvement is to be expected over time.) There is a moderate effect size between time and similarity between partners after being grouped. Is this enough: r = 0.25 for paired students over time, and 0.67, 0.4, 0.3 after grouping by pairs. (It seems almost beside the point to indicate that nonpaired individuals showed no significant similarity to everyone else over time.)


title: Confusion related to semisupervised learning in random walk

text: I am trying to understand the semi supervised learning in random walk. Lets say I have 10 classes and I have some labelled and unlabelled points. Now, I need to find the labels for the unlabelled points using semi supervised learning in random walk. I can define the transition matrix P for the nodes/elements such that every entry $P_{ij}$ gives the probability of moving from node i to j. Now its given that I can propagate the labels. If P is transition matrix, I can have P resetted to P = $$P_{ll} P_{lu}$$ $$P_{ul} P_{uu}$$ and if Y represents a matrix of probability distributions over the label set, then I can use the following iterative algorithm to get the labels for the unlabelled points. Let say $Y_l$ be the set of labelled points given and $Y_u$ be the set of unlabelled points for which we have to find the labels. Lets says there are ten labelled points given for the 10 labels and I have to find the labels for the remaining 100 points lets say, then there is this iterative algorithm $Y^{0} \leftarrow Y$ $t \leftarrow 1$ repeat $Y^{t} \leftarrow PY^{t-1}$ $Y_{l}^{t} \leftarrow Y_{l}$ until convergence to $Y^\inf$ $\tilde{Y} \leftarrow Y^{inf}$ I didn't get how to initialize this Y vector at the beginning. Lets say I have 110 points given. I have label 1.2.3...10 for the ten points, then how am I going to initialize this Y matrix and in the end when I get $\tilde{Y}$ how will I know which class it belongs to. I mean I will just have some values. How am I going to know which class the unlabelled points from $\tilde{Y}$ belong to. If it had been binary I would have known, because if the value was greater than 0.5, I would have said it belongs to class 1 otherwise 0. But what in the case when I have ten labels.


title: Why are probability distributions denoted with a tilde?

text: What is the meaning of the tilde when specifying probability distributions? For example: $$Z \sim \mbox{Normal}(0,1).$$


title: Documented/reproducible examples of successful real-world applications of econometric methods?

text: This question might sound very broad, but here is what I am looking for. I know there are many excellent books about econometric methods, and many excellent expository articles about econometric techniques. There are even excellent _reproducible_ examples of econometrics, as described in this CrossValidated question. In fact the examples in this question come very close to what I am looking for; the only thing missing in those examples is that they are only _research_ reports, without any mention of how the results of the study fared in a _real-world application_. What I am looking for are documented/reproducible examples of real-world applications of econometric theory that ideally have the following characteristics:   1. They should be **reproducible** , i.e. contain detailed description of (and pointers to) data, econometric techniques, and code. Ideally the code would be in the R language.   2. There should be detailed documentation showing that the technique **succeeded in the real world** , according to a well-quantified measure of success (e.g. "the technique helped increase revenues because it enabled improved forecasting of demand, and here are the numbers involved") I am using the term _econometric_ quite broadly here -- I mean any sort of _data-mining_ , _statistical data-analysis_ , _predictiion_ , _forecasting_ or _machine-learning_ technique. One immediate problem in finding such examples: many successful applications of econometrics are done in a for-profit setting and are therefore proprietary, so if a technique worked well, it probably will _not_ be published (this is especially true in the case of proprietary trading strategies), but I am nevertheless hoping there are published examples that have at least property (2) above if not both (1) and (2).