title: Contradiction between significant effect in multiple regression, but non-significant t-test on its own

text: I ran a multiple regression using 10 independent variables and the single dependent variable (consumer complaining behaviour). One of those independent variables was gender. The $R^2$ for the model itself was $.157 (F= 20.50, p = .000)$ which whilst not the highest $R^2$ score was at least significant. Down in the coefficients table Gender $(\beta = -.083, p = .006)$. As my supervisor explained it is a significant score that accepts the alternative, and has a negative relationship with CCB. Interpretively, it would mean men are more likely to complain than women (men = 1 women = 2). Now I got a bit curious and did a t-test to test the difference in means and as it turns out there is no significant difference between the gender groups. This is where I'm getting a bit confused... I'm not sure how I'm meant to interpret these results. It just seems like maybe they contradict each other?


title: Appropriate cluster method?

text: What cluster method would be appropriate for a study with 16 variables both categorial and quantitative? It is a repeated measures observational study. I read in a text that k-means would be best if the data was mostly quantitative although I may have interpreted that wrong. Also, what methods are best for throwing out certain variables prior to performing cluster analysis. I read similar articles to the study I am analyzing, and some opted to do a PCA first, or a discriminant function analysis afterwards to confirm the number of chosen clusters. Thoughts and comments?


title: Estimating 95% confidence interval for vertex of quadratic fit

text: Can anyone tell me how to obtain the 95% confidence intervals for the x and y coordinates of the vertex (peak in my case) of a quadratic fit? Is there a way to get these CIs in R? Thanks.


title: Determine the communication classes for this Markov Chain

text: Say we have a Markov Chain with probability matrix $$ P = \begin{pmatrix} 0.25 & 0.25 & 0.5 & 0 & 0 \\\ 0 & 0.66 & 0 & 0.33 & 0 \\\ 0 & 0.25 & 0.25 & 0.25 & 0.25 \\\ 0.5 & 0 & 0 & 0.25 & 0.25 \\\ 0 & 0 & 0 & 0.2 & 0.8 \end{pmatrix} $$ I'm confused, it may seem really basic but I don't want to leave it to chance in the exam. What would be the communication classes. I would have thought that there would be only one as all the states communicate, hence it is irreducible. But states 4 and 5 could be a class? I'm sorry for the basic question we had a poor lecturer and we've had to learn this mostly from books.


title: How can I get the percentile rank of a given score if I know the scores at standard percentiles?

text: For a distribution, let's say I know the scores at the 10th, 30th, 50th, 70th and 90th percentiles. I also know the interquartile range and range. Is it possible to get the exact percentile rank of a given score? I'll be using R if that makes things any easier. As a more concrete example, I have the following information: Lowest: 5.10, Highest: 6.00, iqr=0.12 10% <= 5.25 30% <= 5.71 50% <= 5.84 70% <= 5.95 90% <= 5.96 Let's say I want to know the percentile rank of the score 5.33. I can see that it's somewhere between 10 and 30. Is there a way to calculate this outright, or maybe somehow I can recreate (estimate) the distribution based on those percentiles and estimate it that way?


title: Modelling slopes over time

text: I have data of the price of a product before a newer version came out and after a newer version came out. I'd like to model the slope of the product pre the new product, and post the new product. ![enter image description here](http://i.stack.imgur.com/ffsHQ.png) Looking at the data, it is obvious when this point is and the negative slope over time for the product increases in magnitude. Linear models did not make much sense, since intercepts are different and not realistic? PRE:               Coefficients:                          Estimate Std. Error t value Pr(>|t|)         (Intercept)         668.12155    2.25824  295.86   <2e-16 ***     pre                  -0.23071    0.01968  -11.72   <2e-16 ***      POST (after day 150):               Coefficients:                            Estimate Std. Error t value Pr(>|t|)         (Intercept)           821.00351   10.96838   74.85   <2e-16 ***     post                   -1.13929    0.04899  -23.25   <2e-16 ***      Any advice on how to deal with this problem would be helpful.


title: Accounting for post-intervention bias following propensity matching

text: Using a large, cross-sectional survey of victims of violence, I am interested in testing the effect of alcohol intoxication (let's call it the 'treatment') on victims' subjective rating of the seriousness of the assault (outcome). I have run regression analyses to explore factors associated with seriousness ratings and, as you might expect, the amount of injury done is the strongest predictor of rated seriousness (more on this later). Victims who were drunk at the time of their assault had somewhat different characteristics from the other victims (e.g. more men, less educated...) and that that assaults on drunk victims resulted in more serious injury than assaults on sober victims, so I thought it might be a good idea to match victims with just the 'drunk at time of assault' as a distinguishing factor ('treatment') so that I can compare (drunk) apples with (sober) apples. I have run propensity score matching using 14 theoretically-relevant variables to identify a matched sample of sober and drunk victims. My problem is that, because important confounding factors like injury happened after the 'treatment allocation' (drunk at time of assault - intervention group; sober at time of assault - control group), I couldn't logically include those in my matching model. So, when I compare the groups, the biasing effect of injury on the outcome still remains. As far as I can see, I'm faced with three options but I'm probably missing something:   1. Include the statistically relevant, but logically problematic, 'post-treatment' factors in the matching (they are logically problematic because they happened after the 'treatment' was delivered but are likely to affect the outcome).   2. Control for the propensity score and the post-treatment factors in a regression of the outcome on victim intoxication using just the matched sample.   3. Stop messing around with treatment evaluation on a cross-sectional data set and just stick with regression. Thanks in advance


title: Are estimates of varying coefficients in ODE unique?

text: Suppose I have a function which is a solution of an ODE with varying- coefficients.Based on the simulation from the function plus Gaussian noise, I use 2nd order ODE to model the simulated data. But even I got good estimate of solution function, I have a very different estimates of varying coefficients. For example, from Jeff Cash Problem 9, let $\epsilon = 0.5$, the 2nd order ODE $$y''+\frac{4x}{0.5+x^2}y'+\frac{2}{0.5+x^2}y=0$$ $$y(-1)=y(1)=2/3$$ has the solution $$y(x)=\frac{1}{0.5+x^2}$$ I simulate data from $y(x)=\frac{1}{0.5+x^2}+ noise$ and use model $$y''+b_1(t)y'+b_0(t)y = 0$$ to estimate the $b_1(t)$ and $b_0(t)$. I am guessing the estimated $\hat{b_1}(t)$ and $\hat{b_0}(t)$ should approximate the true parameter numerically, that is $$\hat{b_0}(t)\approx \frac{2}{0.5+t^2}$$ $$\hat{b_1}(t)\approx \frac{4t}{0.5+t^2}$$ But from my simulation, the curves of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ are very different from the curves of $\frac{4x}{0.5+x^2}$ and $\frac{2}{0.5+x^2}$. Also, as I change the basis (bspline,fourier, monomial) to expand $b_1(t)$ and $b_0(t)$, they also have very different shape. But based on the numerical solution of estimated ODE $$y''+\hat{b_1}(t)y'+\hat{b_0}(t)y = 0$$ where $\hat{b_1}(t)$ and $\hat{b_0}(t)$ are estimates in term of any one basis system. The numerical solutions are pretty close to the true function $y(x)$. I am wondering if coefficient functions $b_1(t)$ and $b_0(t)$ are identifiable in the ODE model given a set of data. Since it gave me the expression that different estimates $\hat{b_1}(t)$ and $\hat{b_0}(t)$ can contribute the same numerical solution which is though close to the true solution. ![Estimated and true b0\(t\)](http://i.stack.imgur.com/TnssF.jpg) ![enter image description here](http://i.stack.imgur.com/UNUj6.jpg) ![enter image description here](http://i.stack.imgur.com/qwUyp.jpg) From figures above, we can see estimated coefficient functions $\hat{b_1}(t)$ and $\hat{b_0}(t)$ have different functional form from the 'true' functions $\frac{4x}{0.5+x^2}$ and $\frac{2}{0.5+x^2}$. But the estimated numerical solution of ODE(green line) is well matched with true solution $y(x)=\frac{1}{0.5+x^2}$ (blue line) even we have a very different coefficient functions. Any help is appreciated!               # simulate data     t = seq(-1,1,by=0.001)     nn = length(t)     yt = 1/(0.5+t^2)+rnorm(nn,0,0.2)          bbasis = create.bspline.basis(c(-1,1),norder=8,nbasis=30)     xfd = smooth.basis(t, yt, bbasis)$fd     xfdlist = list(xfd)     bfdPar = fdPar(bbasis)     bwtlist = list(bfdPar,bfdPar)          # perform the principal differential analysis     result = pda.fd(xfdlist, bwtlist)          # display coefficient functions     bwtlistout = result$bwtlist         beta0.fd = bwtlistout[[1]]$fd     beta1.fd = bwtlistout[[2]]$fd          plot(beta0.fd,main='beta0')     lines(t,2/(0.5+t^2),col='blue')     legend('bottomright',col=c('black','blue'),legend=c('estimated beta0','true beta0'),lty=1)          plot(beta1.fd,main='beta1')     lines(t,4*t/(0.5+t^2),col='blue')     legend('bottomright',col=c('black','blue'),legend=c('estimated beta1','true beta1'),lty=1)          # define the vary-coef function in terms of t     fbeta0<-function(t)eval.fd(t,beta0.fd)     fbeta1<-function(t)eval.fd(t,beta1.fd)          # BVP     yinit<-c(eval.fd(t[1],xfd),NA)     yend<-c(eval.fd(t[length(t)],xfd),NA)      tvec<-seq(t[1],t[length(t)],0.01)     col<-bvpcol(yini=yinit,yend=yend,x=tvec,func=fun2,parms=c(fbeta0,fbeta1),atol=1e-5,islin=T,nmax=3000)          #plot     plot(col[,1],col[,2],col='green',type='l')     plot(t,yt,col='darkgray')     lines(xfd,col='red',lwd=2)     lines(col[,1],col[,2],col='green',type='l')     lines(tvec,1/(0.5+tvec^2),col='blue')     legend('topleft',col=c('green','darkgray','red','blue'),legend=c('Numerical ODE solution','raw data','basis curve fitting','True ODE solution'),lty=1)


title: Preferred method for identifying curvilinear effect in multi-variable regression framework

text: Say some previous findings identified a curvilinear effect of X on Y, (specifically that X had a positive effect on Y, and that X^2 had a negative effect). You want to see if the same holds for your entirely different sample (although everything else between studies, constructs/measures, are exactly the same). Neither the previous study nor my study are experimental (so I am not manipulating X, merely observing it). There is no explicit theoretical reason why a curvilinear effect would occur. I will give some examples of how I would currently go about seeing if this is true, but I would like some input and peoples opinions on what they believe would be a preferrable method, if any of my suggestions are innapropriate, and of course if their are other alternatives. Examples: 1) Simply examine the bivariate scatterplot and fit some type of smoothing line to the mean of Y over bins of X (e.g. LOESS). Although confounding could be an issue, if a curvilinear effect exists their will likely be some evidence in its distribution. 2) Examine the partial regression plot or other visualization techniques to identify the effect of X on Y independent of other confounding variables. 3) Use some type of model selection criteria (e.g. BIC), and determine whether a model including X^2 is preferred over a model without X^2 4) Include a model with X^2 and see if X^2 has a statistically significant regression coefficient. Like I said any other suggestions are welcome as well. Edit: In this context, my main concern is to identify if the effect of X on Y is best represented in a similar manner to the previous study. While their could be other substantively interesting points to compare between studies (like the magnitude of the effect of X on Y), this is not my main concern.


title: Resources for learning to create data visualizations?

text: I'm interested in learning how to create the type of visualizations you see at http://flowingdata.com and informationisbeautiful. EDIT: Meaning, visualizations that are interesting in of themselves -- kinda like the NY Times graphics, as opposed to a quick something for a report. What kinds of tools are used to create these -- is it mostly a lot of Adobe Illustrator/Photoshop? What are good resources (books, websites, etc.) to learn how to use these tools for data visualization in particular? I know _what_ I want visualizations to look like (and I'm familiar with design principles, e.g., from Tufte's books), but I have no idea _how_ to create them.


title: Multiple Linear Regression

text: Some Regression doubts.. Say, I want to quantify the difference in the effect of BMI in onset age of diabetes for Hypertensive and non-hypertensive patients. I have done Multiple Linear Regression model with Age at onset of diabetes as outcome variable and BMI during diagnosis (continuous variable) and gender (binary) as predictor variables. ie. predicting onset age of diabetes from BMI adjusted for gender. (age_at_diag ~ bmi + gender) This is done separately for Hypertensive diabetes patients and non- hypertensive diabetic patients. Now my questions are,   1. Instead of doing it separately and comparing the coefficients can I add the hypertension variable (binary) in the above equation ie `age_at_diag ~ bmi + gender + hypertension`   2. What is the difference between `age_at_diag ~ bmi + gender + hypertension` and `age_at_diag ~ bmi + gender * hypertension` (ie hypertension as an interaction term)   3. How can you interpret the following MLR outputs              Call:     lm(formula = aad ~ diab_bmi * hyper + gender * hyper, data = t2d)          Residuals:         Min      1Q  Median      3Q     Max      -31.431  -6.605  -0.479   6.511  41.083           Coefficients:                       Estimate Std. Error t value Pr(>|t|)         (Intercept)       57.14978    0.99077  57.682  < 2e-16 ***     diab_bmi          -0.25695    0.02781  -9.239  < 2e-16 ***     hyper1            11.58879    1.40013   8.277  < 2e-16 ***     genderMale        -1.99808    0.36447  -5.482 4.35e-08 ***     diab_bmi:hyper1   -0.19536    0.03877  -5.039 4.79e-07 ***     hyper1:genderMale -1.52097    0.52443  -2.900  0.00374 **       and               Call:     lm(formula = aad ~ diab_bmi + gender + hyper, data = t2d)          Residuals:         Min      1Q  Median      3Q     Max      -30.629  -6.667  -0.480   6.475  39.449           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 60.84250    0.70204   86.67   <2e-16 ***     diab_bmi    -0.35801    0.01941  -18.44   <2e-16 ***     genderMale  -2.76840    0.26246  -10.55   <2e-16 ***     hyper        4.34875    0.25437   17.10   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: What is regularization in plain english?

text: Unlike other articles, I found the wikipedia entry for this subject unreadable for a non-math person (like me). I understood the basic idea, that you favor models with fewer rules. What I don't get is how do you get from a set of rules to a 'regularization score' which you can use to sort the models from least to most overfit. Can you describe a simple regularization method? I'm interested in the context of analyzing statistical trading systems. It would be great if you could describe if/how I can apply regularization to analyze the following two predictive models: Model 1 - price going up when:   * exp_moving_avg(price, period=50) > exp_moving_avg(price, period=200) Model 2 - price going up when:   * price[n] < price[n-1] 10 times in a row   * exp_moving_avg(price, period=200) going up But I'm more interested in getting a feeling for how you do regularization. So if you know better models for explaining it please do.


title: Degrees of freedom for post-hoc tests

text: I have conducted a one-way ANOVA on 8 groups. I understand that the degrees of freedom for that test is 7, but I don't understand how to ascertain the degrees of freedom for the post hoc tests: I conducted a series of t-tests as the multiple comparison post-hocs using the Holm-sidak method of familywise error rate reduction. How does one figure out what the degrees of freedom are for such tests? Is it simply `n1 + n1 - 2`, as it would be for regular t-tests? If so, since my samples sizes are 893 and 169 respectively for the post-hoc t-test in question, are the sample sizes getting a little high, forcing my Type I error rate up? Thanks!


title: Combining randomForests in R, why are the err.rate, mse and rsq components NULL

text: From `?combine` in randomForest: > The confusion, err.rate, mse and rsq components (as well as the > corresponding components in the test compnent, if exist) of the combined > object will be NULL. Why are the err.rate, mse and rsq components for the combined NULL? Is there an efficient way to re-calculate these metrics? I ask because I would like to figure out a way to use the `oob` re-sampling method with the `parRF` model in caret.


title: Drawing Chernoff graph in R with only four face features?

text: I have used **aplpack** package `faces()` method in R for drawing a Chernoff graph of the following data. I would like to draw a face with only four variables `a`, `b`, `c`, and `d` in order to explain why program `A` is better than program `B`. However, the faces() method is also using these variables for several other face features. Could you please help me to draw a Chernoff face having four features with these four variables. For example: length: a, width: d, eyes: b, and nose: c. Just I wanted to distinguish face graphs of programs A and B only with these four face features, and rest of the face features should be the same between programs A and B face graphs.               prog  a   b    c    d     A    0.8 7900  70  27     B    0.3 1920 393  43


title: Using Covariance Estimator to Perform Linear Regression?

text: Suppose you had a method for estimating the population covariance of a vector- valued random variable given observations of that random variable, say $f(Z) \rightarrow C$, where the rows of $Z$ are observations of the random variable. Can one abuse this process to perform a least squares regression $y = x^T\beta + \epsilon$, for $n$-dimensional vector $x$? The idea is that one would have a vector of observations of $y$, call it $Y$, and a matrix of paired observations of $x$, call it $X$, then one would compute the covariance of $Z = [X\; Y]$ (concatenate the matrix next to the vector) via $f$, call it $C$, then let $\hat{\beta} = C_{1:n,1:n}^{-1} C_{1:n,n+1}$. A few questions:   1. Will this work under optimistic conditions? (a simple simulation in Matlab shows that the precision is not great, but the results are within 4 sig figs, so I am guessing it will.)   2. Is this a known trick? if so, does it have a name I can google search for, or is it so trivial that it doesn't require a name?   3. Most importantly, if $f$ can deal with input where some values are missing (say, MCAR--missing completely at random), under what conditions will this technique behave reasonably for regression with missing values? **edit** I _am_ assuming that $x$ is drawn from a zero mean process and that the regression has no intercept term.


title: Combining many datasets to increase confidence

text: I have a few (5-6) data sets, each is a function of time, with the time span the same between datasets. These datasets are all statistics of various perspectives of something (the partitions of a graph), and I am trying to find points in time which show interesting changes across the different datasets. Normally these interesting features are local minima, but this may not always be the case. So for example, dataset1. is the variance of some measure and at time 0.4 it has a local minima. dataset2. is the number of maxima and it also has a local minima at 0.4. The problem is in reality sometimes a minima will show on all datasets and sometimes only on one or a few datasets. So my goal is basically to combine the information from all datasets to increase my confidence that a local minima at a particular time, is an interesting and significant point. A simple method would just be to count the number of datasets which have a local minima at that time and and use some threshold value, above which I define the minima to be significant. But I am wondering if there are more intelligent ways, i.e. ways which can see which datasets are more robust and reliable etc.


title: How to teach students who fear statistics?

text: I am about to help teach statistics to medical students this semester. I've heard many horror stories about the fear of these students from learning statistics. Can anyone suggest what to do with this fear? (Either links to people who are discussing this, or offer suggestions from your own experience)


title: Visualization of connections between groups

text: I have around ten groups (of companies). Each group is connected with each other group. The data I have is representing the strength of the connection. Imagine it's the number of times someone from group A sent an email to group B. The strength of a connection can be 0. There are two connections between two groups, A-B and B-A. a) What would be a good way to visualize this? I could imagine, for example, that each group is a circle. Lines connect the circles and the thickness of the line represents the strength of the connection. Being able to indicate the size of the groups would be a plus but not required. b) Do you know of any software tools to visualize this? The tool shouldn't be too expensive or should be available as a trial version as for now it's just a proof of concept. It doesn't have to be web-based. Thanks!


title: P-value for two variables

text: Let's say I have a point $P_1$ on a 2D scatter plot. The coordinates of $P_1$ are $(x_1, y_1)$. $x_1$ has a standard deviation of $\pm u_1$ and $y_1$ has a standard deviation of $\pm v_1$. I also have a point $P_2$. The coordinates of $P_2$ are $(x_2, y_2)$, with standard deviations of $\pm u_2$ and $\pm v_2$ respectively. How can I decide if $P_1$ and $P_2$ are different in a statistical significant manner?


title: How to subdivide a classification problem

text: I have a problem with a number of inputs and one binary output, which I have tried to train several classifiers to solve. Unfortunately, none of the classifiers (MLP, SVM, bagging) have achieved the required level of accuracy. I am thinking of subdividing the problem by splitting the dataset into two, and using a different classifier on each half, or the same classifier but with different parameters (or the same parameters, but it will be learning a different problem). One approach is to choose the attribute which has the highest correlation with the output, and split the data so that the lowest values are in one set, and the highest in another. I fond that one set gave very good accuracy, whereas the other had a lower accuracy than the combined data set. Is there an alternative approach to subdividing the problem that anyone can suggest?


title: Pearson Correlation with symmetric pairs?

text: I am correlating two features in the leaves of a dendrogram for a given tree distance (2,3,...). When i find the pair "A-B" at distance x I will also have the duplicate "B-A" in the same set for distance x along the tree branches. The result is that for every two leaves i will have for example a pair of (x=2,y=-3) and (x=-3,y=2).. If I do a scatter this symmetry along the xy will be very prominent. Can this be biasing my correlation values? Should i get rid of the redundant pairs, but which one to choose? Or it doesn't matter? Any ideas how to perform these tasks given a tree structure?


title: Statistical test on pairs of percentages

text: I am quantifying a measure of correlation in a graph/network under two different experimental conditions to determine whether there is a significant difference in the measure given by the change of condition. I am repeating this computation for a number (N=5) of heterogeneous graphs/networks in the sense that they may have different number of nodes/links. To avoid any dependency on these graphs/networks I am normalizing each measure by the summation of both measures under each condition. I am now simply using a t-test between both sets of proportions because of the limited number of samples, but I have concerns given the nature of my computations. Which statistical test would be more convenient in this case?


title: Is there an ideal number of measures in PCA ?

text: Working on a problem of extracting a minimal sub-set of criteria in a siting problem, I resorted to PCA, I have only 26 individuals (measures), I naturally thought it would be wise to ask whether that number is sufficient or not. I do obtain results, that's not the issue, I just wonder if I can trust them or not. So I guess, a complementary question would be : **_what is the ideal number of measures to have in order of insuring a solid PCA analysis ?_** Thanks for your replies, have a great day !


title: Plotting a "ragged array" - a list of vectors in R

text: I have a list of numeric vectors, not necessarily all the same length. I'd like to get a plot where each vector in the list is plotted vertically above a label for that vector on the x-axis. Is there an "automatic way" to do that in R? Here's an example set of data:                yl = structure(list(                A = c(4.232886, 3.338263, 1.765557, 1.438516),                 B = c(3.046518, 4.500074, 1.508675, 3.226517, 5.982371, 1.519653),                 C = c(5.788882, 7.644319)              ), .Names = c("A", "B", "C"))      In this case I'd like labels A, B and C on the x-axis, and the corresponding data values plotted vertically parallel to the y-axis.


title: Appropriate weighting for averages

text: What is the appropriate (or the advantages of disadvantages for each) measure for weighting mortality rates across countries: total population or total number of deaths? For example, if you have age-standardized cardiovascular disease mortality rates for 5 different-sized but otherwise comparable countries and want to find the mean, what would be the reasons for or appropriateness of weighting by population or number of deaths?


title: Calculating percentiles

text: Lets say I have multiple table representing 10 countries, and tables contains `user_id`, and `user_score`. Example:               Table CountryA                    Table CountryB       user_id | user_score             user_id | user_score     ---------------------            ---------------------       001   |    1245                  001   |    1023       002   |    1563                  002   |    950        :           :                     :           :       1000  |     850                  1000  |     1600      What I am trying to do is build a percentile table where each row is for the country and each column is 10th, 20th..                        10th, 20th   .... 100th     countryA score, score ....  score       I have started to do this on python, what my question is more related to calculating percentile. At the end I want to know what's the score you need to be getting to be on given percentile. Based on whats given in wikipedia (http://en.wikipedia.org/wiki/Percentile). Using `n = P/100 * N + 1/2` I would get the rank but not the score. Coming back to my problem, IF I were to use this I have to sort the table by score and then use `1000` as `N` and then look which row is tally with the answer that I get from the formula. Is this the correct way to do it? If the approach is correct, would be better that I write user_id, user_score to csv file and load it in to R, and may be R has a nice function to handle this?


title: Alternative to Wilcoxon signed-rank test, for data where many pairs are equal

text: I have several discrete datasets with a small range (-2 to +2) which I want to test against the null hypothesis "the distribution is centred around zero". The Wilcoxon signed-rank test was suggested to me as an appropriate choice, and some answers on this site back this up. The issue is that, for some datasets, a large portion of the pair deltas are 0, meaning most of the samples are discarded. The result is that some datasets pass the test at 95% confidence while others with nearly identical distributions fail it. What would be a good alternative test for such a scenario?


title: What is the objective Scikit-learn's Random Forest classifier is optimizing at each node?

text: I would like to ask what is the specific objective function that Scikit- learn's Random Forest classifier is optimizing at each node for the "Entropy" option. My understanding is that entropy is used as a measure here to split the training data reaching a particular node to the left and right subtree by ensuring that the split is done in such a way that the training data within each split is as "pure" as possible. However, what function is used to split the data? A linear separator? Axis-aligned linear separator? Conic-section?


title: What properties of a text makes it a spam/bad question?

text: I'm trying to identify numeric properties of a text message that make it a spam or, more specifically, a bad question on sites like this one. For example, would things like capital letter density matter? Can anyone make any suggestions based on experience?


title: Joint distribution of column sums when row sums are fixed

text: Suppose I have an $m$ by $n$ table $X_{ij} \in \\{0,1\\}$, where in each row, $r$ randomly chosen entries are set to 1 (the rest are 0), i.e. $\sum_j X_{ij}=r$. I know that e.g. the column sum $\sum_i X_{i1} $ has a binomial distribution $\mathrm{Bin}(m, r/n) $. But what is the joint distribution of all column sums, i.e. of $(\sum_i X_{i1},\; ...,\;\sum_i X_{in})$ ?


title: Modified Poisson binomial distribution

text: In the Poisson binomial distribution each trial is either yes (1) or no (0). Is there a distribution where each trial is either yes (N) or no (0)? I'd like to model a situation where I have, for example, 10 workers. Each worker can produce $N$ units of good in one day or fail with known probability and produce nothing. $N$ can be different for each worker. I need to know distribution of goods after a day of work. Poisson binomial distribution gives me an answer in case each worker produces 1 good.


title: Evaluating clusters of first-order Markov chains

text: I clustered my dataset of several thousand first-order Markov chains into about 10 clusters. Is there some recommended way how I can evaluate these clusters and find out what the items in the clusters share and how they differ from other clusters? So I can make statement like "Processes in cluster A tend to stay in state Y once they get there, which is not true for processes in other clusters." The transition matrices of those Markov chains are too large to just "look and see". They are relatively sparse, if that can help. My idea was to take all the transition matrices in a cluster, sum them and plot it as intensity in a picture (in a scale from 0 to 255). Is there something more "professional" I should try out?


title: Dealing with nuisance variables in a MANOVA situation

text: I have a study in which participants facebook profiles were rated by a group of students. The participants were then randomly allocated into one of three groups: the **uncertain condition** (where they could have been assigned either a liked-best or liked-least rating but they are not told which one), the **Liked-Best condition** (where they assigned the highest rating) and the **Liked-Least condition** (where they are assigned the lowest rating). Once they saw their ratings they completed a measure of mood and then in turn were able to rate how attractive they found the facebook profiles of those who rated them. After this they completed a measure of attachment style (ECR-S), self-esteem (RSES), the mood measure for a second time, and a measure of how much they had thought about the rating group in the last 15 minutes. The **dependent variables** in the study are **attraction levels** and **reported thoughts** of the participants to the rating group of students, as well as an **overall measure of mood** (mood across both time points summed and averaged). **Please bear with me as I am a complete novice at statistics, but the analysis of the data has me extremely confused.** My **IV** in the study is **Condition** , in that I want to see whether there was a significant difference between groups on the dependent variables. However I also have the measures of attachment style and self-esteem, the data for which my supervisor is planning on using in future studies but also wants me to include in mine. The problem is, I do not know where they fit in! Originally I thought that any effect of uncertainty on the DVs may be mediated by attachment style and self-esteem, BUT because the measurements were taken after the ratings had been given, they aren't exactly IVs as they may well have been affected by the ratings given. Please any help on how to analyse this would be much appreciated. I have been told I should use separate ANOVAs, but I do not feel this is justified due to the familywise inflation of the Type I error rate. I was thinking of using MANOVA to analyse Attraction and Reported thoughts (correlated at r = .34, p < .01) by Condition, with a one-way ANOVA for mood by condition. But then I am not sure how to deal with attachment style and self-esteem. Please Help!


title: Adding training examples to Bayesian classifier reduces accuracy

text: I'm working on a problem to predict/classify overall sentiment of a large amount of text, which I can verify on the next day. Each data point is a day and is composed of multiple articles. I bin the words by frequency in those articles to train a Bayesian classifier. I tried predicting the sentiment using all previous data. So on day N I would use all data from 1:N-1. However the accuracy of prediction over 110 test examples was only 51%. I noticed that the predictions were getting worse overtime, because the first 20 test examples were accurate at 70%. When i used only 1:20 examples in my training set for all subsequent predictions, the model was accurate at a rate of 57%. So my question is... How do determine the best features in that initial 20 days worth of data. Should I just use the top X features? Or is this 57% just randomness tempting me to follow an unfruitful path? Might there be a better model to use to classify binary sentiment of a large amount of text as positive or negative (multinomial bayes, sum)? Thanks, Ricky


title: What is the canonical example which show advantage of robust linear regression over LS linear regression?

text: What is the canonical example which show situation when robust linear regression has advantage over least square linear regression ? I was trying to simulate situation when some errors (20% of them) are generated from t-student distribution and 80% are from normal - both distribution with the same variance ! on datasets with 50 observation, and I cant see clearly that robust regression is better, here is my R code for this experiment :               library(MASS)     n=50 # size of datasets     N=1000 # number of regressions     wynik=matrix(0,N,2) # matrix with estimated coefficients     v=5  # parametr of t-student distribution     Sd=(v/(v-2))^.5 # standard deviation of gaussian distribution     a=1 # coefficient           for(i in 1:N){          x=rnorm(n,mean=1,sd=Sd)          e_norm<-rnorm(n,sd=Sd)     e_t<-rt(10, df=v )          y_norm=a*x+e_norm     y_t=a*x+c(e_t,rnorm(40,sd=Sd)) # wariant 2 część to outliery          Zm1=lm(y_t~x)$coef[2]         Zm2=rlm(y_t~x)$coef[2]          wynik[i,1]=Zm1     wynik[i,2]=Zm2     plot(1,1,main=paste(i))     }          plot(density(wynik[,1]),main="density of LS estimator(black) and robust estimator (green)")     lines(density(wynik[,2]),col="green")     # average values of LS and robust estimator     colMeans(wynik)


title: Why is variable selection necessary?

text: Common data-based variable selection procedures (for example, forward, backward, stepwise, all subsets) tend to yield models with undesirable properties, including:   1. Coefficients biased away from zero.   2. Standard errors that are too small and confidence intervals that are too narrow.   3. Test statistics and p-values that do not have the advertised meaning.   4. Estimates of model fit which are overly optimistic.   5. Included terms which can be meaningless (e.g., exclusion of lower-order terms). Yet, variable selection procedures persist. Given the problems with variable selection, why are these procedures necessary? What motivates their use? Some proposals to start the discussion....   * The desire for interpretable regression coefficients? (Misguided in a model with many IVs?)   * Eliminate variance introduced by irrelevant variables?   * Eliminate unnecessary covariance/redundancies among the independent variables?   * Reduce the number of parameter estimates (issues of power, sample size) Are there others? Are the problems addressed by variable selection techniques more or less important than the problems variable selection procedures introduce? When should they be used? When should they not be used?


title: Weighted harmonic mean vs weighted geometric mean

text: > **Possible Duplicate:**   >  Weighted geometric mean vs Weighted mean I searched for the differences between WHM and WGM. When to use each of them? when not?


title: Demand with different Measurement

text: Weekly's demand of a product has a probability function               Demand, x   Probability (fx)     ------------------------     0           0,1     1           025     2           0,4     3           0,15     4           0,10     5 or more   0      I have difficult to find out the   * expected value,    * variance and    * standard deviation  for the weekly demand. Where should I start?


title: Question in fitting regression

text: Consider the plot attached as the histogram of the outcome ($Y$) that is going to be the outcome in a linear regression. Clearly, the histogram shows the outcome is not normally distributed. How can I come up with a transformation that makes the data be normal so that I can fit a linear regression? My goal is to compare the effect of TRT (treatment) vs. CTRL (control). One obvious regression is: $Y = \mathrm{TRT} + \text{other covariates}$ Since $Y$ is not normal, do you think I can assess the TRT effect by considering $Y$ as a coefficient and TRT as outcome and fit a logistic regression? ![enter image description here](http://i.stack.imgur.com/WYw6D.png)


title: Origin of the Naïve Bayes classifier?

text: I've looking around Google Scholar for the earliest mention of this particular classifier and have not had much luck finding a definitive source. I've seen some sources cite as late as the 1980s and other as early as the 1930s. Does anyone know when the Naïve Bayes Classifier was developed and/or first used as a classification technique?


title: Magnitude of Pearson correlation after adjustment by simple linear regression

text: Let ${(x_1,x_2,z_1,z_2)}$ be real-valued vectors of equal length with $${\hat{r}_p(x_1,x_2) \approx \hat{r}_p(z_1,z_2) > 0}$$ $${\hat{r}_p(x_1,z_1) \approx \hat{r}_p(x_2,z_2) > 0}$$ where $\hat{r}_p$ denotes the estimated Pearson correlation coefficient. Now perform two simple linear regressions $$x_1 = \hat\beta_{0,1}+\hat\beta_{1,1}z_{1}+\hat{e}_1$$ $$x_2 = \hat\beta_{0,2}+\hat\beta_{1,2}z_{2}+\hat{e}_2$$ (The coefficients have two subscripts to make clear that these two regression equations do not share their coefficients.) The purpose of the regressions is to extract the residuals $\hat{e}_i$ as versions of the $x_i$ adjusted for/uncorrelated with the corresponding $z_i$. My question is whether general statements can be made about the magnitude of ${\hat{r}_p(\hat{e}_1,\hat{e}_2)}$ in relation to that of ${\hat{r}_p(x_1,x_2)}$, that is, is the estimated correlation between the $x_i$ expected to grow or shrink by "adjusting for" (removing linear correlation with) the $z_i$?


title: Box-Ljung test on white noise series

text: I generate this data in R:               set.seed(111)     ds=rnorm(1000)      When I perform Box-Ljung test to test the independency:               Box.test(ds,type='Ljung',lag=log(length(ds)))      it gave me p-value=0.5957, which is reasonable. However, when I perform this:               Box.test(ds[180:299],type='Ljung',lag=log(length(ds[180:299])))      it gave me p-value=0.00162045, which means the autocorrelation of this sub- sequence is significantly different from zero. This also happens for some other subsequences. However, I generated the time series using `$rnorm$` function in R with 1000 data points. Could anybody please explain this? Thanks very much in advance. For the convenience of you readers who are not using R, the sub-sequence ds[180:299] can be obtained here: http://ykang.hostoi.com/ds120.pdf


title: How to detect "bad inputs" when designing a classifier?

text: For example, warning casual drawings when recognizing hand-written digits. Some quick thoughts:   * A bad input may be an outlier for the training data   * A bad input may cause noticeable ties across different labels


title: Principal component regression analysis using SPSS

text: I have done multiple regression analysis (MLR) of my data and find out $R^2$ and $r$, and then to remove multicollinearity problem I used PCA. This analysis generated PC equal to my variables, I selected only those factors having eigen value greater than 1 after rotation (varimax rotation). Now I'm confused that if I put factor scores of only selected PC as independent variables in MLR I obtained $r$ and $R^2$ less than those obtained by simple MLR but it should be higher. Where I am wrong?


title: Reference on Path Analysis

text: This is not basically a research question rather an ask for help to provide me with some useful references on Path Analysis. This seems to be a trendy technique, but looks a little bit complicated. Is there any downloadable complete free book or PDF where I can find details on its theory and usage? I need to learn the technique for mixed type of variables. Which is the easiest software for performing path analysis? I have heard about AMOS. Can't we do it in SPSS? I would also like to have suggested some references where I can learn how to use the referred software for path analysis with mixed types of variables.


title: The best way to map correlation matrix from [-1, 1] space to [0, 1] space

text: This question was originally posted on SO and Qnan suggested me to post it here :) My question: let you have a correlation matrix; you would like correlations which are next to 1 and -1 go towards 1, while those next to 0 stay there. The simplest way is to use absolute values, e.g. if `Rho` is you correlation matrix then you will use `abs(Rho)`. Is there any way which is theoretically more correct than the one above? As an example: what if I use Normal p.d.f. instead of absolute value?               Adjusted Rho = N(Rho, mu = 0, sigma = stdev(Rho))      where N is the Normal p.d.f. function. Have you any better way? Update: Roman Luštrik suggested me to use the following R code:               x <- runif(min = -1, max = 1, n = 100)     tr <- (x - min(x))/diff(range(x))     plot(x)     points(tr, col = "red")      What are strengths and weaknesses of each method? Thanks,


title: The use of xtmixed with the quickicc command

text: I am attempting to do a icc using the `quickicc` command following the use of `xtmixed`. I'm attempting to check repeatability in a sample of 24 individuals with QT measures (a measure of ventricular depolarization). Each individual has various readings of QT; some have 31 whilst others have 5 measured at various times by one reader. My question: is it possible to do an analysis using the `xtmixed` command when obs vary in numbers when testing test-retest reliability or is it even appropriate? Furthermore, I'm unclear if I am using the commands properly as when I use the command `xtmixed qt || study_id:` I get a iteration string that goes on endlessly. Any assistance would be appreciated.


title: MCMC for structured matrix

text: I'm working with a Gaussian random field, which can be described by joint pdf as follows: $$X\mid \sim N\left ( \mu, \sigma^{2} \left ( I-C \right )^{-1}\right )$$ where $C$ is a structured matrix, $c_{ij}=0$ if $i$ and $j$ is not neightbours. I'm working with the Bayesian inference on such fields. So the main question is about how to sample such structured matrix? Every draw has to be of the same structure as previous, i.e. zeros for not neighbours and the $I-C$ has to be positive definite. Although there are papers on how to generate positive definite matrix by MCMC, but I have not found any about how to draw a structured matrix.


title: What statistical tools does a analyst (financial) require?

text: I just graduated in engineering where I have a very strong background in pure maths, but unfortunately I didnt study any statistics in my course. I know calculus, probability theory, measure theory, analysis etc. But I have not done any courses on regression etc. I am planning to apply for several "analyst" positions at various banks/investment firms in the coming months. Amongst some of them, they specifically require the applicants to have a good knowledge in statistics. What statistical techniques/concepts do I need and should I learn (self study), that an financial analyst would require? Basically I will be analysing and computing statistical and financial data, what techniques should I learn? Should I start with regression, time series etc.? As taken from wikipedia, this is something I would probably be doing: > Financial analysts use spreadsheet and statistical software packages to > analyze financial data, spot trends, and develop forecasts. On the basis of > their results, they write reports and make presentations, usually making > recommendations to buy or sell a particular investment or security. Senior > analysts may actually make the decision to buy or sell for the company or > client if they are the ones responsible for managing the assets. Other > analysts use the data to measure the financial risks associated with making > a particular investment decision. So what statistical concepts/tools/techniques/methods do I need to study, in order to do the above (in the quotes)? I just got the book George Snedecor: Statistical Methods, but what main concepts/techniques are most important/most used by analysts? I'd like to quickly start into those topics first.


title: Skewness in dependent variable (OLS, Gauss Markov, non-normality)

text: Consider the time series to be estimated with OLS: $Y_t = a + \bf{X}_t\bf{b} + e_t$ where $Y_t$ is skewed, $\bf{X}_t$ is a vector of regressor values at time $t$, $\bf{b}$ is a vector of coefficient estimates. Scenario (1): The design matrix is such that the skewness in the response is perfectly accounted for by the design matrix - i.e., there is no non-normality in the residuals. Are our coefficient estimates minimum variance? Scenario (2): The design matrix is such that the skewness in the response flows into the residuals (suppose for this scenario that the dependent has major skew and the residuals has smaller skew but still comes up on a hypothesis test of non-normality over the estimates of the residuals). Are our coefficient estimates minimum variance? What estimator should I use (I only know of 1: a variant of the Wild bootstrap that accounts for skew)?


title: How to assess accuracy of phone GPS in measuring distances?

text: I am a complete stranger to statistics (apart from mandatory courses in college), but lately I ran into an interesting real-world scenario. Recently I started jogging. I take my GPS phone with me to measure time and distance. I picked a route (about 4km, but I don't know exactly) and started running. I have observed that the distance measured is not constant, although the route is the same. Sometimes, the distance measured is a little more than 4km and sometimes a little less. I am interested in the precision of my phone GPS. But there is one small problem: when I arrive, I usually (but not always) check the phone, and if the distance is less than 4Km, I continue running to reach 4Km. So even though I have all the measurements saved in the history, I cannot use them to calculate the mean value, because sometimes I run a little bit more to reach 4km. What I do know, however, is that during my 20 runs, the minimum distance was about 3.9km and maximum 4.05km. So what I am interested in is how precise my GPS is and how long the path is. I expect the answer to be something along the lines "with confidence X%, the path is Y km long and the GPS has Z% distance variation". Or would I have to know the exact mean to calculate that?


title: Can complex sampling be ignored for exploratory patient-centered analyses?

text: Hello fellow StackExchange users, **Preamble:** I have been tasked with performing a cluster analysis (or possibly a latent class analysis, as I am pondering) to find non-overlapping groups of like subjects presenting with similar psycho-social traits, as measured in the dataset. There are _a priori_ hypotheses about the number of groups and their characteristics, but the analysis is to be exploratory rather than confirmatory nonetheless. The dataset contains several thousand subjects and about two-dozen variables of interest. So far, so good. However, there is a catch: the survey comes from a complex **multistage stratified design** , which obviously violates the SRS assumption and its resulting traditional standard errors. Such designs often require bootstrapping to provide valid inferences. Observation weights have been provided by the architects of the survey. **My question:** given that the statistics are exploratory and not confirmatory, and that no formal inferences of any kind are to be produced in this patient-centered analysis (not even precision estimates such as confidence intervals), can the complex survey design be ignored in good faith? The only caveat may be one of _external validity_ ; certain demographic groups have been oversampled. But using the observation weights when calculating the distance/covariance matrices should account for this, no? Many thanks,   Brenden


title: Bayesian model averaging in R

text: I have a logistic model that I've built with the `nls` function in R. I want to use Bayesian model averaging for variable selection, but I can't find a package for that in R. Are there any suitable packages? If not, is it possible to make a not too complicated script for it? Data example:               y<-sample(c(1,0),100,replace=T)          var1<-sample(c(1,0),100,replace=T)          var2<-sample(c(1,0),100,replace=T)          var3<-sample(c(1,0),100,replace=T)      The model:               Sw<- function(y1, N1,N2,N3) {             SA <- nls(y1~exp(c+(a1*N1)+(a2*N2)+(a3*N3))/(1+exp(c+(a1*N1)+(a2*N2)+(a3*N3)))         ,start=list(a1=-0.2,a2=-0.2,a3=-0.2,c=0.2))        SA        }                    model <- Sw(y, var1,var2,var3)      How would I do Bayesian model averaging on this? I have 190 observations, where about 70 are `1`s and 120 are `0`s. I have 13 variables in total.


title: How can r-squared be very high in all groups except one where it is very low?

text: I have run a regressions relating metal revenue to a metal index in each of 25 different factories. In 24 of the factorias $R^2$ is greater than .75. In one factory the $R^2$ is .061. Why is that $R^2$ is so low in one factory?


title: why Hierarchical Clustering pvclust vs. hclust got different result?

text: I am performing the hierarchical clustering analysis on a dataset of 25 viral populations using 3 viral components (variables) to construct a dendrogram with average method and correlation distance calculation. We firstly used the `hclust()` method generate a plot, but we still need to support the dendrogram by statistical analysis. So after transposing the data, we choose `pvclust()` to generate the dendrogramm. However, the plots constructed by `pvclust` and the one generated by `hclust` are totally different. We used the same data and same parameters (average method and correlation distance), but the results are so different. Why might this be? Here is the dataset https://github.com/yingfengisu/RSHOP/blob/master/TVA1.csv # hclust > ![enter image description here](http://i.stack.imgur.com/PaHfr.png)               ######### hclust method #############     sd.data=scale(tav.data)     dd=as.dist(1-cor(t(sd.data)))  # correlation-based distance     plot(hclust(dd, method="average"), main="Average Linkage with Correlation-Based Distance",xlab="", sub="", labels=tav.labs)      > ![enter image description here](http://i.stack.imgur.com/i02XB.png) * * * # pvclust > ![enter image description here](http://i.stack.imgur.com/EDRjX.png)               ######### pvclust method ################     tav.data0 = tav.data[,c(1,2,3)]     rownames(tav.data0)<-tav.labs     tav.data0 = as.data.frame(t(tav.data0))     sd.data0 = scale(tav.data0)     library(pvclust)     result=pvclust(sd.data0,method.hclust="average", method.dist="correlation",nboot=100,r=seq(0.7,1.4,by=.1))     plot(result)      > ![enter image description here](http://i.stack.imgur.com/CH3Fl.png)


title: Correlated varaibles and bootstrap forest possible problem?

text: I have been using this method to reduce the number of variables for a model but i have just had a bad thought. Due to the nature of bootstrap is it possible that by randomly using variables you are actually promoting correlated and collinear variables? I have 2 very dominant X variables that demote others to near noise. I found that bootstrap promoted other variables and basically gave me, what i thought, was a better model for pulling out these minor X variables No i am wondering if these other variables are simply being promoted due to correlation with the dominant variable? x1 is correlated to x2 but x1 is a better predictor of y run a bootstrap every time x1 is not present x2 is now promoted to be a predictor over many individual models x2 now makes a significant contribution ONLY because x1 is absent on enough occasions Is this a valid error i have made?


title: Probability density estimation with Dirichlet process

text: I'm trying to estimate a joint probability density function in a quite high dimensional space (around 15 or even more). I've heard about Dirichlet process methods to do this kind of things. Do you know a good paper/phd thesis/book to learn about this ? Thank You for your help


title: Has anyone publicly shared an implementation of RUSBoost in R?

text: There's no package available on CRAN, so I was hoping someone in the community had written their own function/package. I see it's been done in MATLAB, so I may just have to start with that and write the code myself... Algorithm description - "RUSBoost: Improving Classification Performance when Training Data is Skewed" Example of MATLAB code (which includes code for boosting in addition to the undersampling)


title: Create list from for loop if value in sequence is greater than .5

text: Trying to create a list of significant variables for a regression model. I have 82 variables, so in order to only include the significant ones, I created a list of the correlations and sorted them. I want to include the variables with correlations > 0.5 and <-0.5. I created a for loop, but am having difficulty getting a list. composite is a data frame               sortedcor <- sort(cor(composite)[,1])     regvar = 0     for (i in sortedcor){       if (sortedcor[i] < 0.5){         regvar[] <- i       }     }      getting this error: Error in if (sortedcor[i] < 0.5) { : argument is of length zero In addition: Warning messages: 1: In if (sortedcor[i] < 0.5) { : the condition has length > 1 and only the first element will be used


title: Multilevel analysis of dyadic data with all predictors categorical & on level 1

text: I’am trying to get a new perspective on some data on sleep deprivation and counter measures and thanks in advance everybody for investing some brain power in my problem. Here comes a brief description of my research I hope it helps to understand my question. I investigated how different counter tactics help to overcome the negative effects of sleep deprivation in a competitive situation. Therefore I paired participants and let them play a competitive reaction time task against each other. The dyads always consisted of one sleep deprived and one well-rested participant. That’s my within factor and I would say it’s level one in terms of multilevel analysis. My second iv are the three counter-tactics the sleep- deprived participant were provided with, resulting in three different types of dyads which in an ANOVA design (with dyads a focal variable) would be my between factor. But thinking about the multilevel analysis this second iv is a mixed variable hence all my predictors are categorical and on the lower level. Only the random effect of the dyad (which I was advised to address with a multilevel analysis) is on the upper level. My question(s) now: How can I do that in SPSS (I add the syntax I’ve produced so far below) and does my multilevel reasoning makes sense so far (doing a multilevel analysis wasn’t my idea and given the design and my limited knowledge of the topic I’m not 100% sure it’s even the right approach). Again: Thanks in advance for reading and thinking about my problem. * * *               MIXED time     /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1) SINGULAR(0.000000000001) HCONVERGE(0,     ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)     /FIXED=| SSTYPE(3)     /METHOD=REML     /PRINT=CPS DESCRIPTIVES G SOLUTION TESTCOV     /RANDOM=INTERCEPT | SUBJECT(Dyad) COVTYPE(VC)     /REPEATED=Sleep*Counter_Tactic | SUBJECT(Dyad) COVTYPE(CSH).


title: Is there a closed-form maximum-entropy analogue of the exponential distribution, but restricted to integers?

text: I've used exponential distributions to estimate maximum-entropy probability distributions before. Is there an analogue to the exponential that is restricted to integers? If so, does anyone have pointers to closed-form estimators for that distribution?


title: Detecting whether website visits are automated

text: I am trying to detect automated visits to a website. A typical data set for an automated client is of the form:               userid: visit_time1, visit_time2, ...     94562:  5, 10, 15, 25, 30                       ^                     missed       So the period of visit times is `5`. Occasionally (as I marked above), the automated client will miss a period but the next visit time will still be some integer multiple of the period. Please note that the data is somewhat noisy as well. So a more realistic data set would look like:               userid: visit_time1, visit_time2, ...     94562:  5, 11, 14, 25, 30                       ^                     missed       Whereas a non-automated client would have a data set that looks more like:               13345: 5, 21, 34, 89      with no clear pattern. I tried making a vector with `1's` at the visit time indexes and `0's` everywhere else. I then did autocorrelation and it worked quite well for subjectively determining whether a client is visiting at a regular interval. This was based on this question from the DSP stackexchange However, I still don't know how to programmatically detect this behavior. I also thought of testing the variance of the intervals between visits, but some intervals are multiples of the base period and this skews the mean.


title: Using BTL2 : Can we fix the scale of the ability-scores?

text: Using the Bradley Terry Luce 2 model: We want to rank 41 scripts by paired-comparison (package BTL2). **But our dataset is incomplete**. We can generate ability scores, but the scores at the end of the ranking have huge estimates (s.e). Also when we make dummy-items who are always winning, or always losing. Therefore we want to fix the scale: the best item gets an ability-score of 10, the worst item a ability-score of -10 (for the dummies). Is this possible with BTL2, if so, how can we do/ fix this?


title: regression for angular/circular data

text: I have supervised learning problem where targets are angles. If I would do simple regression then numbers 360 and 1 would be far away for my model, but actually they are close and predicting x and y coordinates doesn't feels right, since I am trying to predict just one number here. What is the proper way of doing such a problem? ![blue dots represent targets](http://i.stack.imgur.com/7t7aw.png)


title: Hyper-parameters for pretrain and fine-tuning

text: When doing deep learning, in particular dnn, it's shown that pre-train each layer in a unsupervised fashion, then fine-tune the weights using the labeled data in a supervised fashion. My problem is, since we're pre-trian each layer in a unsupervised fashion, so how could I pick those hyper-parameters, for example, each layer is a sparse autoencoder, so hyper-parameters includes regularization $\lambda$, sparsity $\rho$, sparsity pernalty control $\beta$ and the number of units in each hidden layer $n$, number of layers $l$. Am I supposed to do it this way: pick a set of hyper-parameters $(\lambda, \rho, \beta, n, l)$, pre-train each layer (until convergence or exceeding a number of iterations), then fine-tune (until convergence), and finally do the job (classification for example) and check out the accuracy. And then try out another set of hyper-parameters? This could be very time-consuming, any better idea?


title: Estimate population quantiles from subpopulations' quantiles

text: Suppose there is a population partitioned arbitrarily into a set of subpopulations that completely cover the original population. Assume that for some variable, we know each subpopulation's _quintiles only_ (or some other fixed set of quantiles), but not its distribution. However, the entire population's percentiles (and of course distribution) are unknown. Is there a way to estimate the quantiles of the whole population? This may sound like an odd situation, but it happens for many published economic statistics where access to the full data is not possible, only summaries of it. EDIT (CLARIFICATION): You may _not_ assume that each subpopulation is drawn randomly. In fact, each is likely to have its own bias. EDIT2: Additional details... (1) You may assume the variable is nonnegative, and further that we have a decent estimate for its minimum value across the population. (2) Rough or "unrefined" estimates are OK. FINAL UPDATE: Expanding on @whuber's suggestions and code snippets, I thought I'd add a nice example with some real-world data (hourly wages, if you're curious) that shows the final population upper/lower bounds, along with the _known_ population quantiles, which happened to be available for this case. I've also made some conservative assumptions about the min and max values of each subpopulation, which helps tighten up the intervals a bit. ![interval example](http://i.stack.imgur.com/x20AT.png)


title: Comparing three means

text: I'm performing a meta-analysis but I have a question: I'm evaluating three different treatment strategies (A,B,C) for each of those treatment strategies I have around 20 papers describing outcomes. Now here is my question: Each of those papers describe the mean age, however, some give a range, some give a standard deviation, some only a mean. Is there a statistical test how I can compare these means? If so, what test should I use. My second question is about comparing three groups with each other. For example, the percentage with hypertension is 65% in group A, 70% in group B and 75 % in group C. I do have the total number of patients for these groups (around 400 to 500 for all the groups). Data is extracted from retrospective studies, so the patients are not matched. What test should I use? Thank you very much for your help!


title: Population parameters of a regression

text: So this has really been bothering me and I was hoping for a (simple!) explanation if possible. Suppose I've specified a linear regression model: $$ Y = \beta_0 + \beta_1 X + \epsilon $$ And an alternative: $$ Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \epsilon $$ And I'm trying to estimate the $\beta$s, say perhaps through OLS (the exact method I don't think is relevant). My question is: what is the exact interpretation of the $\beta$s I am trying to estimate? The confusion arises from the fact that the population values of $\beta_1$ in either specification are presumably different, and this doesn't couch with my understanding of the population coefficients. I had always interpreted the $\beta$s as the partial derivative of $X$ on $Y$ 'in reality'. That is, if you were to change X holding other regressors constant the change in the expected value of Y. By providing a better and better specified model, one ensured that the _estimate_ of $\beta_1$ became more accurate (by separating out correlated variables in the error term). This was important to my understanding; $\beta_1$ was not contingent on the specification of my model -it remained an invariant feature of the population - but rather the estimator we had for $\beta_1$ (b1) changed and became more or less accurate depending on the model. All well and good, but this interpretation doesn't quite work in the example above. Suppose that the relationship between $X$ and $Y$ is curvilinear. If you were restricted to only include $X$ and not any higher order polynomials, then presumably the $\beta_1$ that would best describe the change in $E[Y]$ given a change in $X$ would be different than if you were to allow for higher order polynomials (in specification 2). So say, for arguments sake, the DGP was $$ E[Y] = 1 + 10 X - 2 X^2 $$ where $0<X<2$ to ensure the polynomial doesn't influence too heavily. In this case should the true value of $X$ in specification 1 be 10? Or, to fit it to that DGP when $X^2$ is not specified should it be ~6? It seems if it is the latter my understanding that the population coefficients do not depend on the specification go up in smoke! Please help!


title: Gather insights from quarterly financial forecast data

text: I am trying to analyze a quite large (~25,000 rows) dataset of financial forecasts. The forecasts are usually not derived from algorithms, but come from a large number of analysts who forecast the income and expenditures for their respective entities. The forecast starts 5 quarters in advance and is revised once per quarter, giving me a time series of 5 forecasts. The data also includes the actual value, which allows me to calculate the percentage error of the forecast. However, the data aggregates earnings and expenses and does not differentiate between currencies. Thus, one forecast may be revised from -10^6 in t-2 to 10^7 in t-1, and the actual may end up being 10^3, resulting in a huge PE. So far, we have been able to determine that a low volatility in forecasts ends up providing a better final forecast (measured as mean/median absolute percentage error). This came from 'manually' crunching the data, but I suspect that clustering or other ML techniques could reveal more useful patterns. So, if you have any quick ideas on how to   * find patterns in the forecasts   * find the commonalities of good predictions   * solve the issue of a skyrocketing PE for small actual values – I'm eager to listen and learn. I have been somewhat successful in 'manually' splitting the data based on quantiles of the APE, volatility, or on the order of magnitude of the actual, but I suspect that this could be "automated" by clustering the data, but I'm not sure how to go about doing this with the kind of data I have on my hands here. Thanks a lot, Henning P.S.: I have R, MatLab, and Weka at my disposal (though all previous calculations were done in R).


title: Centrality Measures for a directed multigraph

text: I have a directed multigraph that is used to represent an online discussion forum, where each actor is able to contribute comments to a discussion and have their comments commented on by other actors. Which centrality measures are most appropriate to use for a directed multigraph? Can I use the standard centrality measures used for non-directed binary graphs (e.g closeness, betweeness, eigenvector), or will this cause problems?


title: Determining sample size needed in a 4-way mixed ANOVA to reach significance

text: I have a data set for a 2 (between) x 3 (within) x 3 (within) x 3 (within) mixed ANOVA. I need to find the sample size necessary to find a significant result for one of the 3-way interactions including the between subjects factor and two of the within subjects factors. I'm trying to figure out the easiest/best way to simulate the addition of additional participants based on the current data. I currently have 24 subjects total (12 in each of the between subjects groups). I'm also familiar with R, MATLAB and Python. I've found some similar questions on the site but none of them seem to be trying to model the simulation based on existing data.


title: Maximum likelihood of function of the mean on a restricted parameter space

text: I've been trying to teach myself some of the fundamentals of statistics by trying to work through old qualifying exams. Here's a problem: Suppose $X_1, \ldots, X_n$ are a random sample from a normal distribution with mean $\theta$ and variance $\sigma^2$, where $\sigma^2$ is fixed and $\theta>0$ is a parameter. Find the maximum likelihood estimator of $\sqrt{\theta}$. **My work so far:** I have the likelihood function $$L(\theta|\mathbf{x})=(2\pi \sigma^2)^{-n/2} \text{exp}\Big(-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\theta) \Big).$$ Then, differentiating $L(\theta|\mathbf{x})$ with respect to theta, and equating the result to zero yields $$\sum_{i=1}^{n}(x_i-\theta)=0 \implies \theta=\frac{1}{n}\sum_{i=1}^{n}x_i \hspace{3mm} (=\bar{x}).$$ The second derivative is negative here, hence the MLE of $\theta$ is $\hat{\theta}=\bar{X}$. So, by the invariance property of MLEs, the MLE of $\sqrt{\theta}$ is $\sqrt{\hat{\theta}}=\sqrt{\bar{X}}.$ **Question:** My question is about the $\theta>0$ assumption (I'm pretty sure $\sqrt{\bar{X}}$ would be fine if there were no restrictions on $\theta$). So, for $\theta>0$, would the MLE for $\theta$ be something like $\max\\{\bar{X},0\\}$? I would greatly appreciate any feedback, corrections, etc.


title: How to find the best input value for this simple problem?

text: Suppose I have a neural network, with input variables $a, b, c, d, f, g$ and output variables$ m, n, o, p, q$. Given different input values, the neural network will output corresponding $m, n, o, p, q$. Now I want find out the best input values which can maxmize $m, n$, while minimize $o,p,q$ with different weights as well. So how can I find the best $a, b, c, d, f, g$? Currently I use a simple way, which calculate $x= w_1 m + w_2n+w_3 \frac{1}{o}+w_4 \frac{1}{p}+w_5 \frac{1}{q}$, then find the input to get maxmization of x. However this simple method assume $m, n, o, p, q$ are independent, which is not the case. So how should I solve this problem? Many thanks.


title: How to model the sampling distribution of the sample sum

text: I'm stuck on a stats problem and am wondering if someone can point out the error in my logic. Imagine you are planning for a camping trip for 50 people. Each person consumes 2.0 lb of food per day on average, with a standard deviation of 0.7 lb. What is the least amount of food you need to bring to be 97.5% sure that you will have enough food? **First line of reasoning:** The standard deviation for the average amount consumed by one person is **0.7/sqrt(50)** The standard deviation for the total amount consumed by all 50 people is **sqrt(50)*0.7** The desired answer is thus **2.0*50 + 2*sqrt(50)*0.7.** **Second line of reasoning:** SD(n*x) = n*SD(x). Thus, the standard deviation for the amount consumed by 50 people is **50*0.7** The desired solution is then **2.0*50 + 2*50*0.7** I know the second line of reasoning is wrong. Could someone explain why? Thanks! ** EDIT ** I think the error is that the standard deviation for the amount consumed by 50 people is not 50*0.7, because the people are independent. VAR(amount consumed by 50 people) != VAR(50 * amount consumed by one person)


title: How best to graphically represent an r-square effect-size measure?

text: It is easy to graph effect size (r) for a correlation in the form of a scatterplot. However, how do I graph the variance (r-squared)? Any thoughts on best practice?


title: Scaling count data by population size and area

text: Let's say I want to count the number of people in 10 different counties that have $x$ characteristic. As each of the 10 counties have different population sizes and areas, I want to scale the count data by population size and area. Is this the correct formula to use: $$\frac{\text{count of people in the county with x}}{\text{population size of county} \cdot \text{area of county}}$$


title: Displaying regression results in MATLAB

text: Along the same lines as this question, is there a nice way to display regression results in MATLAB from a single or many regressions in table or graph form?


title: Repeated measures ANOVA for circular / angular / directional data

text: I am looking for a test for circular data that is equivalent to linear repeated measures ANOVA (I have an experiment using human participants where the same sample of participants perform multiple experimental conditions and the dependent variable is an angle). I have noticed 3 tests for performing ANOVA on circular data:   1. Watson & Williams (1956),    2. Harrison & Kanji (1986),    3. Anderson & Wu (1995) http://fmatoolbox.sourceforge.net/API/FMAToolbox/General/CircularANOVA.html Do these tests require independent samples (as I suspect is the case for the Watson-Williams test) or can they be used on repeated measures data? If the former, is there another test that does what I desire?


title: Computing regression coefficient values in OLS

text: I have a question on testing the hypothesis that a particular regression coefficient in a simple OLS scheme with all the good assumptions is zero or not. In particular eq 3.12, in the book by Tibshirani and coll., _Elements of Statistical Learning_ defines the z-score as $$z = \beta/\hat\sigma\sqrt{v_j}$$ My question is that given the regression coefficients are all jointly normal, how can we separate one coefficient out like that? Is that a conditional z-score assuming all others betas have a particular value?


title: SVM parameter selection and model testing with cross-validation

text: I've read:   * Model selection and cross-validation: The right way   * Crossvalidation and/or testdata. Always use both or can one exclude the other? but I still don't get it. My problem is to construct a simple SVM, tune it's parameters and calculate the generalization error. Assume I have a dataset with e.g. 10 features and 200 samples. I don't want to waste data, since it's relatively small. My approach would be:   1. Split the dataset (e.g., 70/30) with holdout method into training / validation and test set.   2. Make a repeated n-fold cross-validation on the training set. I calculate the error rate (misclassification rate) after the folds are done (simply counting all misclassified samples). I try to minimize the error rate (or any other loss function?) each time the complete n fold is done. I store the error rate and choose the parameters with the lowest error rate. Then I retrain the model on the complete training set with the estimated parameters.   3. Calculate generalization error with the test set. The problems:   * The larger my test set is, the smaller gets the train set, so I discard potential information. Can this be solved via a "stacked" n-fold cv?   * Do I really have to make a REPEATED n-fold cv? Are there other possibilities?   * Is the error rate an appreciate loss function or should I choose another one (eg. the empirical error function or mse, but then I'd need a probability output, right?)?


title: Expectation of conditional normal distribution

text: I have two jointly normally distributed variables $s_1$ and $s_2$. I am now searching for the conditional expectation $$ E(s_1|s_1>r_1,\ s_2>r_2) $$ where $r_1$ and $r_2$ are constants. An idea of how to arrive at this and sources where it can be read up would also be greatly appreciated.


title: How to calculate a Gaussian kernel effectively in numpy

text: I have a numpy array with _m_ columns and _n_ rows, the columns being dimensions and the rows datapoints. I now need to calculate kernel values for each combination of data points. For a linear kernel $K(\mathbf{x}_i,\mathbf{x}_j) = \langle \mathbf{x}_i,\mathbf{x}_j \rangle$ I can simply do `dot(X,X.T)` How can I effectively calculate all values for the Gaussian Kernel $K(\mathbf{x}_i,\mathbf{x}_j) = \exp{-\frac{\|\mathbf{x}_i-\mathbf{x}_j\|_2^2}{s^2}}$ with a given _s_?


title: Bayesian vs frequentist Interpretations of Probability

text: Can someone give a good rundown of the differences between the Bayesian and the frequentist approach to probability? From what I understand: The frequentists view is that the data is a repeatable random sample (random variable) with a specific frequency/probability (which is defined as the relative frequency of an event as the number of trials approaches infinity). The underlying parameters and probabilities remain constant during this repeatable process and that the variation is due to variability in $X_n$ and _not_ the probability distribution (which is fixed for a certain event/process). The bayesian view is that the _data_ is fixed while the frequency/probability for a certain event can change meaning that the parameters of the distribution changes. In effect, the data that you get changes the prior distribution of a parameter which gets updated for each set of data. To me it seems that the frequentist approach is more practical/logical since it seems reasonable that events have a specific probability and that the variation is in our sampling. Furthermore, most data analysis from studies is usually done using the frequentist approach (i.e. confidence intervals, hypothesis testing with p-values etc) since it is easily understandable. I was just wondering whether anyone could give me a quick summary of their interpretation of bayesian vs frequentist approach including bayesian statistical equivalents of the frequentist p-value and confidence interval. In addition, specific examples of where 1 method would be preferable to the other is appreciated.


title: Statistical test to find variable causing failure

text: I have a database with tests that are either a pass or fail and around 25 description fields. Each description field is either a boolean or has between 3 and 10 options. I'd like to hook into my Postgres database using R and extract a subset of data that's hopefully representative of the entire database. Then I'd like to find out which fields correlate to test failures. Thanks for your help!


title: Does the sign of the adjusted residuals matter in a crosstable?

text: In a cross table, if the adjusted residual value for a cell is less than -1.96 or greater than 1.96 then we could say that this cell is determinant in the dependency (suppose pearson is 0.03). However, is this sign really important? I mean, if I get -2.5 can I say that the variables are negatively correlated? I'm creating a cross table with two variables: cluster membership (generated by kmeans) and a few company characteristics. So, If I get this -2.5 in a cell (cluster 3 x industrial companies), could I say that cluster 3 is "poor" of industrial companies? Thank you.


title: Can I validly expand the odds ratio analogously to the relationship between relative risk & the incidence rate ratio?

text: I am involved in a nested Case-Control study that involves cohorts of cases and controls entering a program, with the outcome of failure by six months. Cases and controls are individually matched on month / year of entry. I have a lot (>10,000) of both. A primary risk factor of interest is receiving a medical diagnosis (psychiatric, musculoskeletal, respiratory) after entry and before 6 months. I have the person-time at risk for cases and controls to develop this risk factor, and can directly measure the incidence density of onset of risk factor in each cohort. **Now, for background:** The relative risk (RR) is defined as: $$ \frac{\left(\frac{N_\text{Cases exposed}}{N_\text{population exposed}}\right)}{\left(\frac{N_\text{cases unexposed}}{N_\text{population unexposed}}\right)}. $$ The incidence rate ratio (IRR) is defined as: $$ \frac{{\left(\frac{N_\text{Cases exposed}}{\text{Person-time exposed}}\right)}}{\left(\frac{N_\text{cases unexposed}}{\text{Person-time unexposed}}\right)}. $$ The odds ratio (OR) is defined as: $$ \frac{N_\text{cases exposed}*N_\text{controls unexposed}}{N_\text{controls exposed}*N_\text{cases unexposed}}. $$ Under the rare disease assumption the OR approximates the RR. **My Questions:** Can I validly expand on the odds ratio using person-time: $$ \frac{N_\text{cases exposed}*\text{Person-time unexposed}}{N_\text{controls exposed}*\text{Person-time exposed}} $$   * to approximate the IRR?    * If not, why not, and what alternative would you suggest?    * If I can do this, would logistic regression would be appropriate?    * What other analytic approaches would you suggest?


title: Extrapolate forecasts to present data

text: The question title doesn't help, but i'll try to explain it: I have historical data (weekly points) from 2001/01/01 to yesterday and my variable has a strong seasonal pattern of 12 months (weather data). I'm using the stlf() function (from the R 'forecast' package) to make weekly forecasts - and getting very good results by the way. The problem is that for some models i don't have the 2 last years of data, so for these models my data goes from 2001/01/01 to 2010/12/31. What can i do to bring these forecasts 'to the future'? In other words, i need to extrapolate these forecasts to the current week. What i have: data from 2001/01/01 - 2010/12/31 What i need: make forecasts for the next weeks. To add confusion, i do have some points from the last month, but i don't know how to take advantage of this. Thanks for any help!


title: Intra- and inter- observer error, and the intra-class correlation coefficient (ICC)

text: I have 3 observers, who each perform multiple observations in which they count a phenomenon. I'm trying to quantify both inter- and intra- observer error for this study. I'm new to this and am finding the analysis options quite confusing. Someone recommended using the ICC. I've mocked up some data below with "number of words" on a page as the count variable - this mimics my data structure. I read Shrout and Fleiss (1979), and since I want to generalize to a population of observers it looks like I want the Two-way random effects model, ICC(2,1). My issues are:   1. for intra-observer error (or agreement), should I only include one observer per ICC? In other words, get the ICC for the repeated observations within observers?   2. for inter-observer error (or agreement), should I include all observers and the repeated observations, or should I aggregate the repeated observations so that there is only one for each observer?   3. I'm also not clear on whether I'd want the ICC for a single observation or averaged observations [i.e., ICC(2,1) vs ICC(2,k)]. Code:               # data with repeated observations     dat <- data.frame(         pages = c("page1", "page2", "page3", "page4", "page5"),          Peter1 = c(236, 179, 203, 218, 195),         Peter2 = c(232, 181, 198, 212, 197),         Peter3 = c(225, 176, 205, 216, 191),         Paul1 = c(227, 177, 204, 218, 196),         Paul2 = c(236, 179, 209, 222, 199),         Mary1 = c(234, 184, 201, 216, 191),         Mary2 = c(238, 182, 200, 218, 193),         Mary3 = c(240, 176, 206, 221, 197)         )          library(psych)          # intra-observer reliability?     ICC(x = dat[, paste0("Peter", 1:3)])     ICC(x = dat[, paste0("Paul", 1:2)])     ICC(x = dat[, paste0("Mary", 1:3)])          # inter-observer reliability?     ICC(x = dat[, -1])


title: Recommendations for a short, yet theoretically rigorous, seminar or course on principal components analysis

text: I am looking for a short course (1 week or less) that teaches the theoretical and applied aspects of principal components analysis in a mathematically rigorous way. I've seen university courses on multivariate statistics in general or the linear model, but not just PCA. Travel is not an issue, so good international courses are welcome, as long as they are in English. Also, I do not want to just learn how to run a PCA using some software package and interpret the results; I want to learn the detailed underlying theory as well. I cannot attend a semester-long course, nor do I want to do a web-only class.


title: Evaluate collisions in classification problems

text: **Short Introduction** In a classification problem, the objective is to identify to which of a set of categories a new observation belongs, on the basis of a set of examples whose category membership is known. The bottom-line of any solution to this problem is to find a set of features that describes the classification of the examples into the corresponding categories (hypothesis). When the size of the set of categories is large, we need more information to distinguish among the different categories (e.g., we could look for more features). Also, if the examples for a particular category have high variance (intra-class variance), it is more difficult to find the discriminant due to the existence of collisions between members of different classes. **Question** I'd like to know if there exists a quick exploratory method to evaluate how difficult a learning problem is. For example, I'd like to evaluate the number of collisions between different categories. The problem of training a classifier and performing some sort of cross-validation is that the solution is dependent on the choice of the learning model.


title: Back Transformation

text: If I had a response variable that was square-root transformed, and an explanatory variable that is log transformed, and I wished to back transform the model using the summary statistics below, such that Y ~ (X)^2, how would I interpret the meaning of the relationship between X and Y using the estimated Beta coefficient?   I thought it was interpreted as: "If there is a 1% increase in X, there is approximately a change of sqrt(2.1014))/100 units increase in Y."                   Residuals:             Min      1Q  Median      3Q     Max          -37.051 -12.096  -4.908   9.701  68.071               Coefficients:                     Estimate Std. Error t value Pr(>|t|)             (Intercept)  -3.0147     2.0827  -1.448    0.148             Dose.Back     2.1014     0.1679  12.514   <2e-16 ***         ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 16.28 on 1154 degrees of freedom     Multiple R-squared:  0.1195,    Adjusted R-squared:  0.1187      F-statistic: 156.6 on 1 and 1154 DF,  p-value: < 2.2e-16


title: Is there a way to calculate the probability of drawing X+ > y% of time given X-

text: Say I have a bag with 10 balls (either red or blue), and I know that there are 8 balls that are red. How do I calculate the probability of drawing a a red ball $>75\%$ of the time if the first ball I picked without replacements is red? Hope this is clear. Thanks!


title: Hypothesis testing for more than two samples?

text: We do a monthly fundraising letter at my nonprofit. I know how to set up tests in R to compare whether the proportion of response and amount received between two mailings is significantly different, but I would like to compare all twelve for the last year. I figure there is a test to do that, but I don't know what it would be. Can someone suggest a way to look at 12 samples at once?


title: Binomial regression with categorical explanatory variable

text: When doing a multiple binomial regression, is it possible to do the regression with a categorical variable as x at the same time as y is categorical? Most of my x-variables are continuous except for two, while my y-variables is all categorical.


title: Help with predicting binary sequence

text: I have 100,000 random binary sequences like this               01011101101001100111001100100001010...     11110111001110000110001100001010101...     00110010100011111000100010101010101...     10001001110101000111000101010011110...     11100100111100010100010010100001010...     ...      and so on. I also have a master binary sequence like this               ms = 00110101010100010101110010010010110.      My question is, how can I find the best random sequence to predict the next master sequence number?


title: Three open philosophical problems in statistics

text: I recently finished reading The Lady Tasting Tea, a fun book about the history of statistics. At the end of the book, the author, David Salsburg, proposes three open philosophical problems in statistics, the solutions to which he argues would have larger implications for the application of statistical theory to science. I had never heard of these problems before, so I am interested in other people's reactions to them. I am venturing into territory about which I have little knowledge, so I'm just going to describe Salsburg's portrayal of these problems and pose two general questions about these problems below. Salsburg's philosophical problems are: >   1. Can statistical models be used to make decisions? >   2. What is the meaning of probability when applied to real life? >   3. Do people really understand probability? > **Statistics and decision making** As an illustration of the problem presented in question 1, Salsburg presents the following paradox. Suppose we organize a lottery with 10000 unnumbered tickets. If we use probability to make a decision about whether any given ticket will win the lottery by rejecting this hypothesis for tickets with probabilities below, say, .001, we will reject the hypothesis of a winning ticket for all the tickets in the lottery! Salsburg uses this example to argue that logic is inconsistent with probability theory as probability theory is currently understood, and that, therefore, we currently do not have a good means of integrating statistics (which, in its modern form, is based in large part on probability theory) with a logical means of decision-making. **The meaning of probability** As a mathematical abstraction, Salsburg argues that probability works well, but when we attempt to apply the results to real life, we run into the problem that probability has no concrete meaning in real life. More specifically, when we say that there is a 95% chance of rain tomorrow, it is unclear to what entities that 95% applies. Does it apply to the set of possible experiments that we could conduct to obtain knowledge about rain? Does it apply to the set of people who might go outside and get wet? Salsburg argues that the lack of a means to interpret probabilities creates problems for any statistical model based on probability (i.e., most of them). **Do people understand probability?** Salsburg argues that one attempt to resolve the issues with the lack of a concrete means of interpreting probability is through the concept of "personal probability", proposed by Jimmie Savage and Bruno de Finetti, which understands probability as personal beliefs about the likelihood of future events. However, in order for personal probability to provide a coherent basis for probability, people need to have a common understanding of what probability is and a common means of using evidence to draw conclusions about probability. Unfortunately, evidence such as that produce by Kahneman and Tversky suggests that personal beliefs might be a difficult basis on which to create a coherent basis for probability. Salsburg suggests that statistical methods that model probabilities as beliefs (perhaps such as Bayesian methods? I'm stretching my knowledge here) will need to deal with this problem. **My questions**   1. To what extent are Salsburg's problems really problems for modern statistics?   2. Have we made any progress towards finding resolutions to these problems?


title: Statistical Significance of Population Change

text: How do you determine if the percent growth between two time periods is statistically significant? For example, the total population in 2000 was 180,000. It grew to 220,000 ten years later. How do you determine if this 22% growth was significant? Similarly, group A (6,000 people) earned £50,000 last year. Group B (9,000 people) earned £65,000. Are Group B's 50% earnings over Group A statistically significant?


title: Some questions about VAR-models, $\Phi$-matrix-coefficients and partial-(auto-)correlations

text: There is an abundance of literature about VAR-models, which teaches how to test preconditions, specify and estimate VAR-models for stationary and also cointegrated time-series. However, I'm still a bit puzzled by the meaning of the model's parameters. Namely the coefficients $\phi_{ijk}$ of the phi-matrices in the standard VAR(p)-model: $\boldsymbol{x}_t = \Phi_1 \boldsymbol{x}_{t-1} + \dots + \Phi_k \boldsymbol{x}_{t-k} + \dots + \Phi_p \boldsymbol{x}_{t-p}+ \boldsymbol{\epsilon}_t$ From what I understand, in the univariate case, the coefficient $\phi_{11,k}$ provides an estimate of the PAC (partial auto correlation) of the $k$-th Lag of $x_t$. I understand this to be the (auto-)correlation of $x_t$ with its $k$-th Lag $x_{t-k}$, when all other correlations $\phi_{11,t-l}$ with $l \ne k$ are controlled for.   1. What I am not clear about is how to interpret the coefficients $\phi_{ijk}$ in the more general multivariate case: Considering a bivariate VAR(p)-model, does $\phi_{11k}$ still have the same interpretation as in the univariate case or does it now, in fact denote the partial (auto-)correlation of $x_t$ with its $k$-th Lag $x_{t-k}$, when also all other correlations $\phi_{11,l}$, $\phi_{12,l}$ and $\phi_{21,l}$ ($l \ne k$) are controlled for?   If indeed, those other correlations for lags $l \ne k$ are controlled for, what about the correlation $\phi_{12,k}$ with the same lag as $\phi_{11k}$, is it controlled for, as well?   2. I am also interested in the partial correlations between all components of $\boldsymbol{x}_t$ without time-lag, so in the bivariate case between $x_{1,t}$ and $x_{2,t}$ after controlling for all the lagged correlations. This should in principle be given by $\Sigma_{\boldsymbol{\epsilon}}=\text{Cov}(\boldsymbol{\epsilon}_t)$. Did I miss something, here?   3. What is a good (standard) method of determining the significance of the partial-correlations within this setting (with and without time-lag)? Background to my questions is, that I have a dataset of 12 time-series, originating from couple-interactions (6 series per person). I would like to examine the series (total length = 180), which represent affectional variables which, hypothetically are closely related to one another. One simple way of obtaining an overview of dynamical structure of the dyadic system as a whole would be to estimate partial-correlations between the series (with and without lag). To avoid spurious correlations, I'm actually using cointegration techniques, thus gaining an ECM(p), which then can be reformulated as a VAR(p). I assume, that $p=2$ will suffice for the system in question.


title: Is the overlap between two gene expression samples significant?

text: I have performed an experiment to study the response of a yeast (that contains 5000 genes) to stress caused by heat shock. I have one list of 48 genes that are overexpressed at 37ºC and another list of 145 genes that are overexpressed at 42ºC. There are 38 genes that are overexpressed in both of them. By chance I expected only 1 gene overexpressed in both of them, how can I calculate if the overlap that I have obtained is significantly? How can I obtained the $p$ value? I know nothing about biostatistic or math software. Thanks you very much!!! Any help will be very welcome :)


title: Comparing the difference between two ranks

text: I have a list of teams in a local soccer league, ranked from 1 to 10, where the top (1st) ranked team has the most wins (and hence "points"). I'm trying to create a "power" list, so that if the 10th ranked team wins against the 1st ranked team, the game is more important that if the 1st ranked team wins against the 10th ranked team. Part of this formula involves calculating the opponent strength. For this I was using this formula:               opponent strength multiplier = ([Number of Teams] - [Ranking Position]) / 100       I was wondering if anyone had another (better) method to create a multiplier to compare teams? My final "power" calculation has a number of other factors, including points scored, importance of match, etc.


title: Looking for a proof that overfitting a model leads to greater variance estimates (under OLS)

text: So I've been trying to algebraically prove that overfitting a model leads to greater variance values for the parameter estimates. I've gotten close (reduced the problem to showing a certain matrix is positive definite) but I don't really like the approaches I've taken and I'd like to see if there's a simpler way. If anyone knows of a good strategy to use, or a pre-existing proof it would be appreciated, thanks!


title: Do all my dependent variables have to be normal for a non-parametric assumption?

text: Im currently attacking SPSS, and I'm just testing for normality... In the Shapiro-Wilk   Statistic test, my dependent variables have a wide range of significance levels. Obviously, I know that the sig. should be >0.05, which most the variables are...apart from 1 variable. **Does this mean all of my data is non-parametric and I have to use another test?**   * Here is my data below: **0.077** , 0.400, 0.169, 0.157.  _In bold is the variable which is not normal_