title: Why does MICE fail for one dataset and not the other?

text: I am getting this error in MICE               Error in seq.default(1, ncol(pred)) : 'to' must be of length 1      My dataset is very large but I have been able to create a reproducible example with smaller subsets, one of which works, and the other of which fails with the above error:               dt.fail <- structure(list(hsp = c(49, 48, 42, 49, 49, 49, 50, 50, 43, 41, 43, 50, 50, 43, 45, 41, 45, 43, 50, 49), age = c(76, 83, 76, 88, 73, 41, 83, 80, 66, 56, 50, NA, 84, 46, 64, 84, 77, 69, 80, 77), sbp = c(137L, 60L, 155L, 170L, 95L, NA, 140L, 95L, 122L, 166L, 123L, 127L, 59L, 55L, 120L, 207L, 118L, 112L, 150L, 185L), hr = c(70, 90, 84, 64, 80, NA, 110, 77, 105, NA, 50, 80, 37, 38, 87, 76, 88, 93, 90, 80), CTP = c(NA, 57, NA, NA, NA, NA, 100, 142, NA, 83, 56, NA, 66, 48, 46, 88, 43, NA, 79, NA), const = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Names = c("hsp", "age", "sbp", "hr", "CTP", "const"), row.names = c(NA, 20L), class = "data.frame")          dt.OK <- structure(list(hsp = c(2, 1, 5, 1, 5, 4, 1, 2, 4, 1, 1, 5, 5, 1, 2, 5, 1, 2, 5, 2), age = c(65, 100, 84, 52, 70, 83, 68, 81, 91, 80, 69, 54, 92, 66, 90, 42, 78, 55, 83, 78), sbp = c(120L, 116L, 166L, 145L, 58L, NA, 146L, 156L, 130L, 110L, 83L, 162L, 101L, 135L, 118L, 130L, 122L, 136L, 115L, 115L), hr = c(68, 92, 94, 75, 143, NA, 79, 128, 90, 55, 58, 63, 95, 70, 100, 80, 59, 74, 103, 65), CTP = c(76, NA, 196, 59, 78, NA, 111, 50, NA, 60, 54, 92, 113, 65, 54, NA, 102, 53, 80, 69), const = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Names = c("hsp", "age", "sbp", "hr", "CTP", "const"), row.names = c(NA, 20L), class = "data.frame")          require(mice)     require(pan)          PredMatrix <- quickpred(dt.OK)     PredMatrix ["CTP",] <- c(-2, 1, 1, 1, 0, 2)          imp.OK <- mice(dt.OK, meth=c("","","","","2l.pan",""), pred=PredMatrix, maxit=1, ,m=1)          imp.fail <- mice(dt.fail, meth=c("","","","","2l.pan",""), pred=PredMatrix, maxit=1, ,m=1)


title: Confused by Benjamini-Hochberg-like false detection rate calculation in paper

text: I'm trying to figure out where the false detection rate correction in a genetics paper is coming from and failing. The paper in question is looking at a sequencing method where every sample is read multiple times giving, each time, either a positive or negative result. They're trying to threshold the number of positive required before a sample is accepted to give an overall false positive rate of 1%. It's using something very close to the Benjamini- Hochberg method but, for a reason I'm struggling to figure out they've used the total negatives rather than the total trials for the denominator. Benjamini-Hochberg would be $p(k) \le \frac{k}{m}\alpha$ where $p(k)$ is the probability of getting $k$ successes in $m$ trials and $\alpha$ is the desired $q$-value. This paper is using $p(k) \le \frac{k}{u}\alpha$ where $u$ is the number of negative trials. Why would they be doing this? I can see that it is a more conservative approach but I can't find any theoretical justification for using this particular formulation. The papers in question are these: http://www.ncbi.nlm.nih.gov/pubmed/20436463 and http://www.nature.com/nature/journal/v462/n7271/full/nature08514.html


title: Why is naive Bayes a linear model?

text: People always said that naive Bayes is a linear model. I am not able to understand why, so can anybody explain?


title: Line chart vs Area chart - NOIR

text: Do you know NOIR classification of the data? NOIR - nominal ordinal interval ratio. I want to argue that a line chart can be used with both interval and ratio data whereas area chart should be used only with ratio data. This is because ratio data has a meaningful 0 and in case of an area chart the baseline of the graph is 0 (meaning that the area spans from the curve to the horizontal axis). Therefore it would be unreasonable to use an area chart with interval data. However this is my feeling only and I am not sure. What do you think about it? Do you know any publications that discuss this matter?


title: How to improve forecasting accuracy?

text: I got some users' history data and generated some sequences of real numbers. The length of each sequence is between 15 and 25. What's more, I do not know whether these sequences have patterns and the frequency is not known as well. My goal is using each sequence to predict its next value, and then I use auto.arima in R to do this. However, the accuracy of the prediction is low. Anyone have any good ideas to improve the accuracy? One of these sequences is:                   1.5959709882736206           0.7300914525985718           2.0011744499206543           3.6755871772766113           0.8066112399101257           1.3413848876953125           3.371157646179199           0.4400146007537842           2.637667655944824           2.1453769207000732           2.341433048248291           2.3429665565490723           1.1187453269958496           1.4169363975524902           3.328829050064087           4.157748699188232           3.9255290031433105           2.7843635082244873


title: How to model a five point rating item predicted by multiple variables some of which are factors?

text: I have carried a questionnaire, which has a variety of questions in which the responses are given in a rating form (e.g. 1-5, from bad to good). I now would like to test whether the response to question A is affected by several variables, which include some factors (e.g. gender), some continuous (e.g. age) and also some discrete (other 1-5 responses). My first thought was that I should use a GLM, but I am unsure if this is the best for this type of data and which would be the appropriate distribution. If anyone would have an idea of what is best to do, I would greatly appreciate to know it.


title: Creating contrast matrix for linear regression in R

text: Sorry for asking such stupid question, but I cannot understand the criteria behind contrast matrix, should be created for designing linear models in R. I have read limma user guide (P: 101) and there was an example for that. The example is about estrogen(present/absent) and its late (48h) and early (10h) effect on a cancer. I cannot understand in what sense they could take away the effect of time by assigning 1 for both e48 and E48. Here are the R commands:               treatments <- factor(c(1,1,2,2,3,3,4,4),labels=c("e10","E10","e48","E48"))     contrasts(treatments) <- cbind(Time=c(0,0,1,1),E10=c(0,1,0,0),E48=c(0,0,0,1))     contrasts          > treatments     [1] e10 e10 E10 E10 e48 e48 E48 E48     attr(,"contrasts")         Time E10 E48     e10    0   0   0     E10    0   1   0     e48    1   0   0     E48    1   0   1     Levels: e10 E10 e48 E48


title: How to analyse a ranking and rating scale together?

text: In some market research consumers are asked to rank the features of a product based on priority. For example, _Rank the following features for a device based on your priority (1 being the topmost priority)_               Storage capacity         6     Portability              5     Touch interface          1     Keyboard                 4     Long battery life        2     Entertainment on the go  3      _On a scale of 1 to 5, rate the features based on importance (1 being very important)_                                   1   2   3   4   5     Storage capacity    1                    Portability                 3            Touch interface     1                    Keyboard            1                    Long battery life       2                Entertainment on the go         4         Now based on the ranking and rating, I want to assign weights and in the end find out whether the consumer prefers a laptop or a tablet PC/ Say in this example the consumer has ranked touch as 1 and battery life as 2; these are the attributes of a tablet. But he has rated keyboard and storage capacity as the most important, which are the features of a laptop with a keyboard. How do I combine both these and probably assign some weights and come up with a score? If the score is above a certain level then the consumer prefers a tablet and below a certain level then the consumer prefers a laptop.


title: Understanding ordered logistic regression

text: I have the following toy data:               x <- structure(c(2L, 2L, 3L, 1L, 2L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 3L, 2L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 2L, 2L,      2L, 3L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 3L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 2L, 3L, 2L, 3L,      3L, 2L, 2L, 2L, 3L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 3L,      3L, 2L, 2L, 2L, 3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor")          y <- structure(c(2L, 2L, 3L, 1L, 2L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,      2L, 3L, 2L, 1L, 2L, 2L, 3L, 2L, 2L, 3L, 2L, 3L, 2L, 3L, 1L, 2L,      2L, 3L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 3L, 3L, 2L, 2L, 2L,      2L, 3L, 2L, 2L, 2L, 1L, 2L, 3L, 2L, 2L, 3L, 3L, 1L, 3L, 2L, 3L,      3L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 2L, 1L,      3L, 2L, 2L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 1L, 3L,      3L, 3L, 2L, 2L, 3L, 3L, 2L), .Label = c("1", "2", "3"), class = "factor")          z <- structure(c(1L, 1L, 3L, 2L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 3L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L,      1L, 3L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 3L, 1L, 3L,      3L, 1L, 1L, 1L, 3L, 1L, 3L, 3L, 1L, 1L, 1L, 1L, 3L, 3L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 3L,      3L, 1L, 1L, 1L, 3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor")      I have replaced the 1's in `x` with 2's in `z` and vice versa. Now when I do an ordered logistic regression in R with `polr` (library MASS), I get the following coefficients:               f1 <- polr(y~x,Hess=TRUE)     f2 <- polr(y~z,Hess=TRUE)          coef(summary(f1))            Value Std. Error  t value     x2  25.95727  0.3028808 85.70127     x3  30.21524  0.5463144 55.30742     1|2 24.02167  0.3480269 69.02246     2|3 27.77068  0.3432316 80.90944          coef(summary(f2))              Value   Std. Error       t value     z2  -21.495979 6.530398e-10 -3.291680e+10     z3    4.257964 8.119540e-01  5.244095e+00     1|2  -1.935599 3.567345e-01 -5.425880e+00     2|3   1.813399 3.411874e-01  5.314964e+00      It seems that something is not correct. Why relabeling the levels is changing dramatically the estimates for the SEs?


title: Omitted variable bias formula with more than one variable

text: For a general model $$y_{i} = \alpha + \beta_{1}X_{1} + \beta_{2}X_{2} + \epsilon_{i}$$ regressing $y_{i}$ on $X_{1}$ alone will result in $\beta_{1}$ being biased given by: $$plim \: \widehat{\beta}_{1} = \beta_{1} + \beta_{2}\frac{Cov(X_{1},X_{2})}{Var(X_{1})}$$ Now I stumbled across a paper which has a model of the form $$S_{i} = \alpha D_{i} + \beta P_{i} + u_{i}$$ without any other variables and no intercept. The true value of $\beta$ is supposed to be zero because the event $S$ preceds event $P$. The authors then state the omitted variable bias of the two coefficients as: $$\widehat{\alpha}_{OLS} = \alpha + \frac{Cov(D,u)Var(P) - Cov(D,P)Cov(P,u)}{Var(D)Var(P) - Cov(D,P)^{2}}$$ and $$\widehat{\beta}_{OLS} = \frac{Cov(P,u)Var(D) - Cov(D,P)Cov(D,u)}{Var(D)Var(P) - Cov(D,P)^{2}}$$ I have tried for a long time now to reproduce this result but I cannot figure out how they arrived at these two bias expressions. The problematic point is that I apparently cannot understand from which OVB formula they start because in their setting it cannot be the usual one for a single variable written above. The paper is by Adda et al. (2011) for which the posed problem can be found on PDF page 52.   If someone could point towards the right direction I would be most grateful.


title: Sum standard deviation vs standard error

text: I'm having difficulty in determining what exactly the difference is between the 2, especially when given an exercise and I have to choose which of the 2 to use. These is how my text book describes them: **Sum standard deviation** _Given is a population with a normally distributed random variable $X$. When you have a sample $n$ from this population the population is:_ $X_{sum} = X_1 + X_2 ... + X_n$ with $\mu_{Xsum} = n \times \mu_x$ and $\sigma_{Xsum} = \sqrt{n} \times \sigma_x$. **Standard error** _When you have a normally distributed random variable $X$ with mean $\mu_X$ and standard deviation $\sigma_X$ and sample length $n$, the sample mean $\bar{X}$ is normally distributed with $\mu_{\bar{x}} = \mu_X$ and $\sigma_{\bar{x}} = \dfrac{\sigma_X}{\sqrt{n}}$_ These 2 are awefully similair to me to the point I can't at all decide which to use where. Here are the problems where I discovered I couldn't: **Problem 1** _A filling machine fills bottles of lemonade. The amount is normally distributed with $\mu = 102 \space cl$._ $\sigma$ = $1.93\space cl$.   * Calculate the chance that out of 12 bottles the average volume is $100 \space cl$. The problem itself is easy, however the troublesome part is what to choose for the standard deviation of the sample. Here they use $\dfrac{1.93}{\sqrt{12}} $ which I can live with, until I encountered the second problem. **Problem 2** _A tea company puts 20 teabags in one package. The weight of a teabag is normally distributed with $\mu = 5.3 \space g$ and $\sigma = 0.5 \space g.$_   * Calculate the chance that a package weighs less than 100 grams. Here I thought they'd also use $\dfrac{0.5}{\sqrt{20}}$, but instead they use $\sqrt{20} \times 0.5$. Can someone clear up the confusion?


title: Mathematical background for neural networks

text: Not sure if this is appropriate for this site, but I'm beginning my MSE in computer science (BS in applied mathematics) and want to get a strong background in machine learning (I'm most likely going to pursue a PhD). One of my sub-interests is neural networks. What is a good mathematical background for ANNs? Like in other areas of machine learning, I assume linear algebra is important, but what other areas of mathematics are important? I plan to read Neural Networks: A Systematic Introduction or Neural Networks for Pattern Recognition. Does anyone have any input or alternative recommendations?


title: Weighted sum of chi distribution

text: How can I approximate a weighted sum of _chi_ distribution? I found several documents for the weighted sum of **chi-squared** distribution and approximation with gamma distribution but I didn't find anything about a weighted sum of _chi_ distributions.


title: Predicting the edges of a graph

text: I have a dataset of paired relations, indicating whether $a$ is in relation with $b$. It is better to consider this dataset as a graph where each node has a numerical value as its feature. Let's say this feature could possibly be varying between $-10$ and $10$. Now the question is: Are two nodes related one to each other, or not? If they are then the value of this expression is $1$: $(\text{node}_1,\text{node}_2) = 1$, and there will be an edge between them in the graph, otherwise $(\text{node}_1,\text{node}_2) = 0$, which means that there will be no edge between $\text{node}_1$ and $\text{node}_2$. Let's put the problem in context this way. The data is about medicine. Each medicine has one feature (solubility), and some medicines are more effective when used together. So the dataset is about testing the two medicine altogether. Now if the medicine (node) gets more effective for, say, an illness then the value of $(\text{node}_1,\text{node}_2)$ would be 1, or in other words there will be an edge in the graph. We've done this experimentally, but now what I want to do is to leave-one-node-out learning and trying to guess it's edges, for I could add more support toward my experiment. Now that I have the data, I want to predict the edges of a new node in a way that each time I am leaving one node out or leave-one-node-out. I am leaving the node and all its connections(edges) out and from the information I have about the other nodes in dataset I want to predict the edges of the node which i left out. My question is what is the best method/algorithm and solution for my problem?


title: What are some well known improvements over textbook MCMC algorithms that people use for bayesian inference?

text: When I'm coding a Monte Carlo simulation for some problem, and the model is simple enough, I use a very basic textbook Gibbs sampling. When it's not possible to use Gibbs sampling, I code the textbook Metropolis-Hastings I've learned years ago. The only thought I give to it is choosing the jumping distribution or its parameters. I know there are hundreds and hundreds of specialized methods that improve over those textbook options, but I usually never think about using/learning them. It usually feels like it's too much effort to improve a little bit what is already working very well. But recently I've been thinking if maybe there aren't new general methods that can improve over what I've been doing. It's been many decades since those methods were discovered. Maybe I'm _really_ outdated! **Are there any well known alternatives to Metropolis-Hastings that are:**   * **reasonably easy to implement,**   * **as universally appliable as MH,**   * **and always improves over MH's results in some sense (computational performance, accuracy, etc...)?** I know about some very specialized improvements for very specialized models, but are there some general stuff everybody uses that I don't know?


title: Measure of goodness-of-fit in Ordinal Logistic Regression with continuous independent variable

text: In case of the ordinal logistic regression, both of the goodness-of-fit statistics, Pearson and Deviance goodness-of-fit measures, should be used only for models that have reasonably large expected values in each cell. If you have a continuous independent variable or many categorical predictors or some predictors with many values, you may have many cells with small expected values. SPSS warns you about the number of empty cells in your design. In this situation, neither statistic provides a dependable goodness-of-fit test. So if I have some continuous independent variables then is there any measure of goodness-of-fit? How do I know if my model fitted well or not?


title: True meaning of $\pi_0 = 1$ in false discovery rate

text: I am not a statistician so apologies if my terminology is wrong I have a dataset of >14 million p values derived from Fisher's exact testing on genome scale sequencing data. Benjamini-Hochberg correction of these p-values turns out only ~1500 significant p-values so I have been attempting to use qvalue R package to relax the FDR and therefore gain a greater number of significant p-values. Applying the qvalue package with default options returns a $\pi_0$ of 1               qobj<-qvalue(p.fish.unadj)     qsummary(qobj)     pi0:    1      I presume this is a "good thing" and means that:   1. I have good power in my dataset to detect _true null_ results and therefore _true alternative_ results also   2. $\pi_0=1$ is an approximation because my true alternatives i.e. significant tests is such a small number compared to the overall size of the dataset Am I right?


title: confusion related to maximum a posteriori estimation

text: I was reading this article in wikipedia related to MAP http://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation. However, I had this confusion when it says MAP estimation is a limit of Bayes estimators (under the 0-1 loss function), it is not very representative of Bayesian methods in general. What does it mean MAP estimation is a limit and it has presented a criticism of MAP stating that the mean/median estimates will be better. I didn't get why it is so? Any guidance will be appreciated.


title: Representing experimental data

text: I have an argument with my advisor over data visualization. He claims that when representing experimental results, the values should be plotted with " _markers_ " only, as presented in the image bellow. While curves should only represent a " **model** " ![Markers.png](https://dl.dropbox.com/u/6484897/markers.png) I on the other hand believe that a curve is unnecessary in many cases in order to facilitated readability, as shown in the second image bellow: ![Lines.png](https://dl.dropbox.com/u/6484897/lines.png) Am I wrong or my professor? If the later one is the case, how do I go around to explain this to him.


title: arima regression Improving the model

text: Given the following model which relates the full year home sales to the unemployment rate (observed or estimated) I get a projected increase of 14% for 2013 over 2012... last year the same approach over projected by 6% (the 2012 projection was for 41,992 and the actual is coming in about 39,535) So I think the model is over projecting but I'm at a dead end thinking of a good (valid) way to modify it? BTW the projection)s) is identical to what a simple straight line regression in a spreadsheet yields. I want to also thank the fine people at Stack overflow who got me this far in my 1st foray into R http://stackoverflow.com/questions/14032768/csv-input-to- r-forecast-with-dates-via-r-studio#14032768 Pointers appreciated.                # load the base data as presented in the question      Workbook1 <- structure(list(Year = structure(1:10, .Label = c("31-Dec-04",      "31-Dec-05", "31-Dec-06", "31-Dec-07", "31-Dec-08", "31-Dec-09",      "31-Dec-10", "31-Dec-11", "31-Dec-12", "31-Dec-13"), class = "factor"),      total = c(51439L, 59674L, 58664L, 55698L, 42235L, 37918L,      36234L, 36965L, 39535L, NA), UnemplRt = c(5.7, 4.7, 3.8,      3.7, 4.3, 8.5, 10.9, 10, 8.3, 7.1)), .Names = c("Year", "total",      "UnemplRt"), class = "data.frame", row.names = c(NA, -10L))          # Make a time series out of the value     dependent <- ts(Workbook1[1:9,]$total, start=c(2004), frequency=1)          # load forecast package     require(forecast)     # load independent variables in variables.     unemployment <- ts(Workbook1[1:9,]$UnemplRt, start=c(2004), frequency=1)         unemployment_future <- ts(Workbook1[10:10,]$UnemplRt, start=c(2004), frequency=1)          # make a model that fits the history     fit2 <- auto.arima(dependent, xreg=unemployment)          # generate a forecast with the already known unemployment rate for 2013.     fcast2 <- forecast(fit2,xreg=unemployment_future)     fcast2          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95     2013       45168.02 38848.92 51487.12 35503.79 54832.25      This yield exactly the same result a a simple spreadsheet liniear regression. And I belive the suggested increase in home sales is too high, last year the projection turned out to be about 6% too high. So I'm trying to torture the numbers in some statically valid method to come up with a somewhat lower number for 2013


title: Factor analysis on multiply imputed data

text: I have a data set with approximately 500 observations on eight key variables. There are a lot of missing data; only about 1/12 of the observations are complete. I am using `PROC MI` and `MIANALYZE` in SAS to run various regressions on multiply imputed data, and this is working well. (There are about 200 variables in total, and there are high correlations among them which helps multiple imputation.) However, I would also like to do factor analysis on the imputed data. This does not appear to be easily done in SAS, and it poses some interesting challenges:   1. The signs of the factors are arbitrary, so different imputations could yield opposite signs;    2. What was factor 1 in one imputation could be factor 2 in another imputation (although in this case, with so few key variables, it is likely that one factor is enough), and probably other issues as well. I could do exploratory factor analysis on each imputed data set, of course, and I could then average them on a sort of ad-hoc basis, but this seems very sloppy. Some Googling did not reveal any solutions, but ... well, Google doesn't always find everything. Any help, ideas, references appreciated. I have access to SAS and to R.


title: ARCH and GARCH Models

text: I found excellent notes on ARCH and GARCH models here. On page 3 it is given that: Standard time series models: \begin{eqnarray*} Y_{t} & = & E\left(Y_{t}|\Omega_{t-1}\right)+\epsilon_{t}\\ E\left(Y_{t}|\Omega_{t-1}\right) & = & \mu_{t}\left(\theta\right)\\ E\left(Y_{t}|\Omega_{t-1}\right) & = & E\left(\epsilon_{t}^{2}|\Omega_{t-1}\right)=\sigma^{2} \end{eqnarray*} What does this $\Omega$ means here? Any help will be highly appreciated. Thanks


title: Where and why does deep learning shine?

text: With all the media talk and hype about deep learning these days, I read some elementary stuff about it. I just found that it is just another machine learning method to learn patterns from data. But my question is: where does and why this method shine? Why all the talk about it right now? I.e. what is the fuss all about?


title: Do M-estimators and L-estimators overlap?

text: In my self-study, I have read the wikipedia entries and some books in regards to M-estimators and L-estimators. I understand that M-estimators are so called "M" because they "Maximize" the likelihood function, and/or functions of the data itself. I also understand that L-estimators are so called because they use "Linear" combinations of order-statistics of the data. It seems to me though that there are overlaps between M-estimators and L-estimators, is this not correct? (For example, the mean seems to me that it can be derived as an M-estimator, or as an L-estimator). Put another way, if we were to draw a Venn Diagram of all M and L estimators, how much would they overlap with each other, if at all? **Edit:** Thanks to the feedback, I understand that such a Venn Diagram might have infinite elements. What I would like however, is very rough percentages of overlap, _based on ones' experience_ , since I know that this experience is far greater than mine! I am just after a rough estimate based on what experts and practitioners in the field have seen and come across. Thank you.


title: Stata - Help for switchr command

text: I'm trying to perform a markov switching regression model in stata using the command _switchr_. The command syntax is the following one: _switchr eq1 eq2 [weight] [if exp] [, cluster(string) strata(string) sigequal tol(real) tout(integer) noisyn(integer) ]_ where _eq1_ and _eq2_ are the regime equation (the equation in which the variables coeffients change with respect to the regime) and the main equation (or classification equation). The stata help explains that equation as saying " It must be defined by the user so that its dependent variable is a variable containing an initial guess of the partition of the observations into the two components. This initial guess variable will be overwritten by switchr to contain the probabilities that a given observation adheres to the first component regression." I'm not able to understand what the "main equation" means. Which are the regression variables of the main equation? I hope I was pretty clear asking this question. Thanks for helping.


title: Measuring correlation or dependence between two data sets

text: Is there any statistical test or measure to evaluation the degree of correlation or dependence between two sets of data-points ?


title: Hierarchical recommender

text: I have a hierarchical feature matrix, By that I mean that each item may belongs to one or more category, so my data will look something like that |  User |  Categ  | Item  |   |  1       |   12      | 120   |   |  2       |   15      | 411   |   |  9       |   35      | 411   |   |  1       |   15      | 321   |   where the numbers represent ids for users, categories and items respectively. As you can see there is a hierarchy here, a set of categories for each user and a set of items within that category, where one item can belong to more than one category. My goal is to recommend new items to the new incoming users. I already have a user similarity algorithm in place and I can find the top n similar users. Now I need to use the above table to recommend one or more items to the new user.   * Easy way is to take the categories from the most similar user, pick some items from these categories and suggest them to the new user  But I am looking for a better approach that will consider the combination of **_categ-item_** simultaneously. This is a very large data set, so the algorithm should scale. Any suggestion?


title: Is it appropriate to use SE in place of SD in rnorm()?

text: I have wood density data for a number of tree species given a tree's state of decay. I am presented with mean density and standard errors (SE). I am NOT provided with sample sizes. For example, a decay class 1 fir tree has a mean density of 0.305 with a standard error of 0.024. I am using R. I would like to generate a density value for a given decay class 1 fir tree similarly to using `rnorm(n, mean, standard deviation)` when pulling from a normal distribution. In this case, is it appropriate to substitute the SE for SD, such as `density<-rnorm(1, mean=0.305, sd=0.024)`?


title: what does this +/- of "average" mean?

text: I have a set of data points and I was told that for it to be valid, every data point should be within +/- say 5% of the average of the entire dataset. I have no problem of checking each data point and reach a conclusion if it is valid. My question is, how can I relate above statement to a statistic measure so that I can say, how valid it is. Put in another way, is there any relationship between the +/-5% deviation from average and the dispersion measure such as sd etc? Thanks Oliver


title: Dispersion measure - probability density function

text: I am wondering whether someone has a tip on this potentially very basic question. i have done some grid-based bayesian analysis and ended up with a non-standard discrete posterior density function over an interval of 100 observational values (per graph) below and as shown in the vector format per below. I can easily find the corresponding value of the distribution in R (which.max) but am struggling to find a way to calculate any measure of dispersion of the probability distribution. Ideally - i would like to have a quick way of calculating in R that gives me the answer that x% of the density lies between e.g. the index values of 20 and 80; any tips very welcome... ![enter image description here](http://i.stack.imgur.com/TfGml.png)                 [1] 4.615412e-02 7.685555e-02 7.936106e-02 9.355243e-02 9.390333e-02 9.415255e-02 8.812764e-02 8.094180e-02       [9] 7.111677e-02 6.096649e-02 5.067963e-02 4.113239e-02 3.252751e-02 2.514941e-02 1.900486e-02 1.405895e-02      [17] 1.018350e-02 7.230007e-03 5.033092e-03 3.437899e-03 2.305284e-03 1.518504e-03 9.831581e-04 6.261243e-04      [25] 3.925111e-04 2.424243e-04 1.476553e-04 8.878700e-05 5.277273e-05 3.104750e-05 1.810726e-05 1.048550e-05


title: priors for strictly positive index or score types of variables

text: Is there a prior that's commonly used for "index" or "score" type variables that are user-defined as a weighted sum of a small number of variables (sometimes with pre-defined interaction contributions)? I know the gamma, inverse gamma, binomial, negative binomial, poisson, etc distributions are strictly positive, but they correspond to processes that don't describe something like an index/score type variables well. Distributions of index/score type of values are usually overdispersed relative to these distributions. I could construct my own mixture distribution, but I was wondering if there is are standard practices for this kind of data.


title: Chi-square independence in political sciences

text: I've been assigned to do a cross-class project between political science and math about the behavior of voters in a democratic state. More specifically, I'm focusing on a party from a national election in 2007, where I have made a table consisting of two rows – mandates and news about that party – and one column for each election day:                   Day        1   2   3   4   5   6...         Mandates  11   8  11  13  12  12         News      30  46  12  33  48  42      I would like to check whether the appearence of news in the media about this party have affected the number of mandates given to the party in the daily opinion surveys. Would it be _legitimate_ to perform a chi-square-independence test on this data or would I not be able to conclude anything? Why/why not? -- To get a little more advanced, I was thinking if it is possible to evaluate if the news have affected some days particularly than others with the chi-square- independence test ( **yes it has to be this method** ). I would then calculate the change in the mandates for each day which would make my table look like this:                   Day        1   2   3   4   5   6...         Mandates  11  -3   3   2  -1   0         News      30  46  12  33  48  42      Then for each column, I would manually calculate the expected result. Example for column "Day 2": The observed results are that there has been 46 news and the party has lost three mandates. We assume that the news haven't affected the number of mandates and thus they would stay equal and the change will be zero: Observed data:                   Day        2         Mandates  -3         News      46      Expected data:                   Day        2         Mandates   0         News      46      If I do so for every column, would I then be able to evaluate that the news has affected the number of mandates one day, but not another, mathematically? I would really like to test for this, because then I would be able to pin point some days to focus on from which I can extract articles and look upon the positiv/negative rhetorics made by the party. Sorry for the bad tables, but I was not able to post images. Thanks in advance. /Brinck10


title: Appropriate statistical test to test if probabilities are accurate

text: I have some data that looks like this:               Prob    Outcome     0.09    0     0.10    0     0.10    0     0.11    1     0.84    1     0.99    1     0.86    1     0.78    1     0.86    1     0.00    0     etc.      i.e. a bunch a probabilities each with a single test. What statisitcal test should I use to test the hypothesis that the probabilities are correct? **Further details** : The data points are _combat probabilities_ from the game Civilization IV, and I have over 3000 of them in my set. Thus, each probaility is generated using some unknown formula from different input data, depending on the relative strengths of the units in that battle. It has been suggested that the outcomes do not accurately reflect the probabilities given: for instance, the computer player wins more often that it should, based on the probabilites displayed, which is what we want to test. So there is a link insofar as we assume the probabilities displayed are generated using the same formula for each line. It's this unknown formula that we want to test for consistency with the actual results.


title: How to handle samples that belong to multiple classes in supervised learning?

text: I'm currently facing a problem in which some of my training examples belong to more than one class at the same time, say, sample $y_i$ pertains to class $A$ and $B$. I was thinking that a solution to it would be to consider that sample as two-fold, i.e., consider it as two samples, one for class $A$ and one for class $B$. However, my problem is that I'm performing a one-vs-all strategy, in which I think this solution may cause numerical errors (the feature matrix would have identical rows)! Do you know of any reference to this kind of problem (or the technical name for it)? Thanks in advance!


title: Integrating unknown parameter out of the likelihood function

text: Suppose we wish to model two variables $x$ and $y$, as having an underlying linear relation with added errors. That is, with data $(x,y)_i: i = 1,...,n$, we model: $$ \begin{pmatrix} x_i \\ y_i \end{pmatrix} \sim \mathcal N\left(\begin{pmatrix}u_i\\v_i\end{pmatrix},\Sigma\right) $$ and $v_i = a + bu_i$ with some parameters ($a,b)$ and $\Sigma$ is a diagonal. Assume that the value $u_i$ follow a normal distribution with mean $\mu$ and variance $\tau^2$. Write the likelihood of data given the parameters; you can do this by integrating over $u_1,...,u_n$ or by working with the multivariate normal distribution.   Given density function: \begin{align} f(x,y|u,a,b,\Sigma) &= \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left[\frac{(x-u_i)^2}{\sigma_x^2}+\frac{(y-a- bu_i)^2}{\sigma_y^2}\right]\right) \\ p(u_i)&= \frac{1}{\sqrt{2\pi}\tau} \exp\left(-\frac{1}{2}\frac{(u_i-\mu)^2}{\tau^2}\right) \\ f(x,y|a,b,\Sigma) &= \int ...\int f(x,y|u,a,b,\Sigma)p(u_i)du_1...du_n \end{align} My question is how to get the integration above? What distribution does it follow?


title: Determining which AdWords have the highest amount of user click throughs

text: Possible Duplicate First, I am not a statistician (though I'd like be one) but I am trying to understand how different tests can be utilized to examine samples. Let's say that I have five of Google AdWords ads and I'm trying to develop a test to show which has been more effective at getting users to click on that ad.                    January   February      Ad A - 850       5000            Ad B - 900       5300           Ad C - 880       5100           Ad D - 880       5100             Ad E - 800       5000        In terms of frequency, Ad B has the highest amount of user click throughs in both January and February. (Is it fine just to say that one has a higher freqency and is therefore the more effective ad without relying on any significance test) However, for each month, I want to use a statistical test to determine which ad was more effective. Initially, I assumed a chi-square (pearsons) test would be appropriate, but with five seperate ads, I'm not sure if that's the appropriate plan of action. Can anyone help me understand what action to take and why, or is it enough just to rely on the frequency distribution of the five ads.


title: Simple question about notation

text: In the context of likelihood-based inference, I've seen some notation concerning the parameter(s) of interest which I've found a little confusing. For example, notation such as $p_{\theta}(x)$ and ${\mathbb E}_{\theta}\left[S(\theta)\right]$. What is the significance of the parameter ($\theta$) in subscript notation above? In other words, how should it be read? My first assumption was that it simply meant "with parameter $\theta$"; for example, for $p_{\theta}(x)$, it would read: "The probability density of $x$ with parameter $\theta$." However, this probably isn't correct because $p_{\theta}(x) = L(\theta)$ and, in general, $L(\theta)$ is _not_ a distribution (i.e. it does not integrate to unity); hence it can't be a density, can it? In addition, in the case of ${\mathbb E}_{\theta}\left[S(\theta)\right]$, I'm not sure what it changes relative to ${\mathbb E}\left[(S(\theta)\right]$ (i.e. with the subscript $\theta$ omitted). _In the above $S(\theta)$ and $L(\theta)$ signify the score function and likelihood function respectively._


title: Weighted linear regression

text: I have a set of $n$ events. Each event has $m$ variables. At least 1 event produces an observation. It is possible for several events to occur simultaneously. E.g. 4 Events. 3 variables each: $E_{11}=0.4, E_{12}=-0.3, E_{13}=-0.4, E_{21}=0.3, ... , E_{42}=-3.0, E_{43}=2.1$   produces $y_1 = -5$ Here all events occurred simultaneously. However, it is also possible for this to occur: $E_{11}=0.2, E_{12}=0.2, E_{13}=-2$ and all other events are 0 (they did not occur). So I want to scale up all the coefficients to match the observation. Effectively I want to be able to weigh the events with 4 weights; $W_1$ to $W_4$ so that the total contribution for the first vector (all events) would be: $W_1\cdot$(regression coeffs of $E_1) + W_2\cdot$(regression coeffs of $E_2) + ... + W_4\cdot$(regression coeffs of $E_4)$. Where each $W$ here has been normalized (and all $W$ are $>0)$. In solving this regression, you have the usual matrix of observations $A$, a sparse weight matrix $w$, the vector of unknown coefficients $x$ and the resulting observations $y$. I.e. \begin{align} \displaystyle \left( \begin{array}{lllll} E^1_{11} & E^1_{12} & ... & E^1_{42} & E^1_{43}\\ E^2_{11} & E^2_{12} & ... & E^2_{42} & E^2_{43}\\ ... \\ E^{k-1}_{11} & E^{k-1}_{12} & ... & E^{k-1}_{42} & E^{k-1}_{43}\\ E^{k}_{11} & E^{k}_{12} & ... & E^{k}_{42} & E^{k}_{43} \end{array}\right)\\\cdot \left(\begin{array}{lllll} \frac {W_1}{\sum_4 W_i} & \frac {W_1}{\sum_4 W_i} & ... & \frac {W_4}{\sum_4 W_i} & \frac {W_4}{\sum_4 W_i}\\ \frac {W_1}{W_1+W_2} & \frac {W_1}{W_1+W_2} & ... & 0 & 0\\ ... \\ 0 & 0 & ... & 0 & 0\\ 0 & 0 & ... & 1 & 1\end{array}\right)^T \left(\begin{array}{l} a_{11}\\a_{12}\\a_{13}\\...\\a_{41}\\a_{42}\\a_{43}\end{array}\right)\\= \left(\begin{array}{l} y_1\\ y_2\\ ...\\ y_{k-1}\\ y_{k} \end{array} \right) \end{align} In the above equation, the first observation had all 4 events happen, the second one only events 1 and 2 occurred etc. The last observation only had event 4 occur. The $E$s and $y$s are known. I need the vector of coefficients (the $a$s) and the weights $(w_1$ – $w_4)$. How do I solve this efficiently? Thanks.


title: Each letter has a $40\%$ chance of being replied to; how many letters to send to be $99\%$ sure that you get 200 replies?

text: A researcher knows that the probability that a questionnaire will be reponded by mail is 40%. He wants to be 99% sure he will get at least 200 responded questionnaires back. How many questionnaires must he send by mail to be 99% sure he will get at least 200 answers?


title: Percentiles of a distribution

text: I'm measuring error on the median of my data using bootstrapping. I would like to obtain I sigma error bars on my data, so I'm measuring the 16th and 84th percentiles of my data. Should I divide these percentile values by the square root of the number of data or not? Thanks in advance.


title: Dixon Q Test - can I (all by myself) calculate the value of the Q critical?

text: Im just wondering .. Everywhere I see only tables with the Qcritic as the results. Lets say I have a big set of data samples (like 1000 or more). Is then a way of mathematically calculate a Qcritic value? Is there any mathematical formula for Qcritic?


title: MIMIC fail to converge

text: I am trying to fit a MIMIC model measuring poverty in a household sample. The model looks like this: (HUMCAP -> education job) (HOUSINGQUAL -> floor wall dwelling tv refrigerator radio electricity lighting heating cooking washingmachine) (POV->HUMCAP HOUSINGQUAL) However, when trying to fit this model it iterates forever. I don't understand why bc according to the k(k+1)/2 - rule my model should be identified. I would be very thankful if someone had any ideas. Thank you !! cross-posted: http://www.statalist.org/forums/forum/general-stata- discussion/general/3006-mimic-fail-to-converge


title: Testing if a vector is in a population of vectors

text: Suppose that I am concerned with objects that are described by N 2D points. These points all have equal importance to me. I begin with a large number R of such objects, each with its corresponding N-element list of 2D points. Let's call this the population. I'm now given a new object B, with the N locations of its 2D points. I'd like to know how confident I can be that B belongs to the population. I'd also be happy with a probability, or any other number with a simple interpretation. I'm willing to assume that everything is uncorrelated, point $i$ in each object is uniformly distributed with respect to point $i$ in the other objects, and other such typical assumptions. When there's just one variable, Student's t-test seems the right tool. That's about the limit of my statistical depth. I've read about multidimensional t-tests, Hotelling tests, and related topics, but understanding how to use these tools requires substantially more sophistication and depth than I have. Frankly, I can't turn the references I've found (such as the Wikipedia page on "Multivariate t-distribution"), into a practical solution. I'm encouraged that these tools are built into Mathematica (the system with which I'm familiar - I don't know R), and that they seem very close to what I need, but figuring out how to use them for this question seems just out of reach. The example uses I've read involve comparing things like the means of two populations and other comparisons, but I don't see how to phrase my question in terms of those sorts of tests. Ideally I'd like to spend the time to learn enough statistical theory to understand these tools well enough to answer this question for myself. But at the moment I'm faced with a practical situation for which I need to get an answer. My goal is a little Mathematica routine that takes two inputs: a matrix of R vectors, each made up of N 2D points, and a second vector of N 2D points, and **it returns to me a number expressing the confidence that the vector belongs to the population**. Or, equally useful, some kind of measure that tells me how "far" the new vector is from belonging to the population. If anyone can help me get to that goal in concrete, practical terms I would be very grateful.


title: How to compare differences in the mean when the data isn't normal

text: I have two sets of mean data that I have done a Kolmogorov-Smirnov test on. The tests came back showing that one data set was normal, but the other was not. The two data sets are independent of each other, and I need to compare the differences in the mean. The fact that one of the data sets isn't normal surely rules out using a t-test or ANOVA, so should I use a Mann Whitney test? The statistical package i'm using is SPSS.


title: Interpreting odds ratios less than 1 with 3-category outcome

text: I have a 3-category ordered outcome (food consumption: 1=no food, 2=less food, 3=more food) and a 3-category ordered predictor (food exposure: 3=no time, 2= less time, 1= more time- whereby 3=no time is taken as reference category in the ordinal regression model). I want to explore hypothesis that more food exposure is associated to more food consumption. I want to know how I can interpret odds ratios less than 1. For example, I have OR= 0.62 for predictor category 2= less time. I have calculated OR as exp(coeff) in Excel, whereby OR of reference category no time is exp(0)=1.


title: kde.test gives different p-values for similar KDEs

text: I have 4 sets of random floats between [0,1]: $d1, d2, d3, d4$. I need to compare these sets by pairs and I make use of the `kde.test` function of R's `ks` package to do so. When I use `kde.test` to obtain the `pvalue` of distributions $d1,d2$ I get a value of ~0.55, but if I do the same for $d3,d4$ I get ~0.00. Here's the code I use in `R` to come up with those values:               d1 <- c(...)     d2 <- c(...)     pv12 <- kde.test(x1=d1, x2=d2)$pvalue          d3 <- c(...)     d4 <- c(...)     pv34 <- kde.test(x1=d3, x4=d2)$pvalue      The complete $d1,d2,d3,d4$ sets are here and their lengths are: `len(d1) = 20, len(d2) = 210, len(d3) = 200, len(d4) = 2100`. Sets $d1,d2$ have fewer elements than $d3,d4$ but the _shapes_ of the KDEs are similar. Here's the KDEs for $d1,d2$ (blue and red respectively): ![d1 red, d2 blue](http://i.stack.imgur.com/f2juZ.png) and here's the KDEs for $d3,d4$ (blue and red respectively): ![d3 red, d4 blue](http://i.stack.imgur.com/d8Zem.png) Because of the similarity of the KDEs I would have expected somewhat similar `pvalues` but the results I get are wildly different. I'm either not understanding what the `kde.test` does or what the resulting `pvalues` means (or both). Could someone explain what is happening and why my expectations are no correct? * * * ## Add I get very similar results if I apply the Kolmogorov-Smirnov test from the scipy.stats.ks_2samp package which also gives a `p-value` as output.


title: Assuming homoscedasticity in linear regression, why are predictions still less precise for marginal predictor values

text: I understand, that the basic model of simple linear regression assumes homoscedasticity, i.e. the variances $\sigma^2$ around the regression line are equal for all predictor levels/values. With this assumption in mind, I wonder why confidence intervals for a predicted $y$ value are not equally large for all predictor levels. The formula for such a confidence interval (e.g. 95%) suggests, that the interval gets larger, when the difference between the predictor value of interest $x$ and the mean predictor value $\overline x$ increases (see the enumerator under the root): $\hat y \pm t_{n-2}(0.975) \cdot s \cdot \sqrt{1 + \frac{1}{n} + \frac{(x- \overline x)^2}{\sum_{i=1}^n \left(x_i - \overline x \right)^2}}$ The increasing size of CIs toward the edges of a distribution can also be seen in standard `geom_smooth()` plots by the `ggplot2`package in R. ![Simple linear regression plot with CI](http://i.stack.imgur.com/0vDD1.png) Now, how can the variances around the regression line be equal at all predictor levels, but the corresponding CIs are larger at the edges?   Is this a problem of sample and population, insofar as the difference in CI range shrinks, when the sample gets larger (When $n \rightarrow \infty$, you can drop the entire root term in the CI formula), while at the same time there is higher evidence for homoscedasticity when $n$ gets larger? Or do I get things completely wrong and the two things are totally independent of each other?


title: How to extract covariance (error) matrix of independent parameters x,y from covariance matrix of dependent parameters x,y,z?

text: Suppose I have a covariance (error) matrix `C_xyz` for system with 3 parameters: `x, y, z`. Then it appears that `x, y` and `z` depends on each other:               det(C_xyz) = 0     ax + by + z = 0      How to find covariance matrix `C_xy` of independent `x,y` parameters knowing covariance `C_xyz` and function `a,b`?


title: Confusion relating to Alpha Cronbach

text: I am still confused in using Cronbach's Alpha. Do we have to use it in sample survey before going for final survey? or is value of Alpha calculated based on the final survey data? Similarly, what does the set of questions measuring a construct mean? Does it mean that the questions are mutually exclusive or do all the questions are to be used together (mutually inclusive) to measure the underlying construct?


title: Error running optim function with STAR from book example

text: I'm running an example of Smooth Transition AR (STAR) Model from the book "Analysis of financial time series, 3rd edition" by Tsay, in section 4.1.3. The script is as follows:               da=read.table("data/m-3m4608.txt",header=T)     rtn=da[,2]     star <- function(par){     f = 0     T1=length(rtn)     h=c(1,1)     at=c(0,0)     for (t in 3:T1){     resi = rtn[t]-par[1]     at=c(at,resi)     sig=par[2]+par[3]*at[t-1]^2+par[4]*at[t-2]^2     sig1=par[5]+par[6]*at[t-1]^2     tt=sqrt(sig+sig1/(1+exp(-1000*at[t-1])))     h=c(h,tt)     x=resi/tt     f=f+log(tt)+0.5*x*x     }     f     }     par=c(.001,.002,.256,.141,.002,-.314)     m2=optim(par,star,method=c("BFGS"),hessian=T)      The data can be found at http://faculty.chicagobooth.edu/ruey.tsay/teaching/fts3/m-3m4608.txt . I get the following error at optim:               > m2=optim(par,star,method=c("BFGS"),hessian=T)     Error in optim(par, star, method = c("BFGS"), hessian = T) :      non-finite finite-difference value [2]     In addition: Warning message:     In sqrt(sig + sig1/(1 + exp(-1000 * at[t - 1]))) : NaNs produced      How can I fix it?


title: Fitting data to gamma distribution to find score which corresponds to pvalue < 0.05?

text: I have data of size 116.667 rows defined as:               iD  Signal     chr17.3620  0.5741552     chr1.7341   0.5680284     chr7.3937   0.5479430     chr17.3890  0.5402434     chr12.3200  0.5298978     chr17.7227  0.5298536      As it is difficult to show results with all the data, because I have no reputation enough to upload images I made my post with the first 200 rows. Summary:   * Min. = 0.3693   * 1st Quartile = 0.3847   * Median = 0.4039   * Mean = 0.4199   * 3rd Quartile = 0.4413   * Max. = 0.5742 Also I tried some exploratory analysis to see which kind of distribution follows this data (i.e. `plot(density(t))` and `qqnorm(t)`; `qqline(t,col=2)`). According to preliminary results, I would say that this data apparently follows a gamma distribution. Here I put the first 200 signal values:               A <- structure(list(V1 = c(0.5741552,0.5680284,0.5479430,0.5402434,0.5298978,0.5298536,0.5282417,     0.5165426,0.5131503,0.5129329,0.5105448,0.5104201,0.5095860,0.5090263,     0.5061467,0.4972821,0.4959428,0.4953381,0.4920510,0.4915160,0.4868505,     0.4843749,0.4825519,0.4823313,0.4809742,0.4788553,0.4775991,0.4770962,     0.4745947,0.4743952,0.4727112,0.4718017,0.4714738,0.4674141,0.4670385,     0.4648104,0.4633502,0.4616054,0.4615068,0.4614247,0.4613338,0.4597812,     0.4551755,0.4535067,0.4528133,0.4508228,0.4494993,0.4494936,0.4442789,     0.4413460,0.4412279,0.4402557,0.4392294,0.4385639,0.4385187,0.4361337,     0.4344499,0.4342413,0.4342331,0.4338879,0.4337806,0.4336820,0.4329372,     0.4325534,0.4323201,0.4312287,0.4292037,0.4281761,0.4279843,0.4279774,     0.4252035,0.4243487,0.4228516,0.4226953,0.4218263,0.4214821,0.4212546,     0.4210894,0.4206089,0.4204235,0.4193896,0.4168915,0.4164699,0.4152126,     0.4127455,0.4126053,0.4113571,0.4105654,0.4099753,0.4088188,0.4085093,     0.4075957,0.4074018,0.4072499,0.4072114,0.4067329,0.4065400,0.4052757,     0.4044982,0.4040699,0.4036509,0.4033471,0.4031712,0.4026698,0.4017872,     0.4011538,0.4011325,0.4011320,0.4008897,0.4006470,0.4003469,0.3996736,     0.3992583,0.3991979,0.3990366,0.3989118,0.3983172,0.3980860,0.3978592,     0.3977522,0.3965371,0.3963045,0.3957640,0.3954328,0.3950159,0.3935825,     0.3934975,0.3932916,0.3931091,0.3929565,0.3922829,0.3919779,0.3919713,     0.3914740,0.3910446,0.3909540,0.3890607,0.3890550,0.3876478,0.3875172,     0.3873815,0.3872299,0.3870533,0.3858995,0.3858361,0.3855984,0.3854444,     0.3852595,0.3849558,0.3847531,0.3844442,0.3842814,0.3831377,0.3822418,     0.3817666,0.3805661,0.3803090,0.3802035,0.3800845,0.3800580,0.3799694,     0.3795814,0.3794039,0.3792874,0.3788970,0.3787295,0.3785160,0.3782523,     0.3782439,0.3779547,0.3778596,0.3777452,0.3770986,0.3767652,0.3767104,     0.3765786,0.3760886,0.3760124,0.3753271,0.3750943,0.3749116,0.3744146,     0.3743998,0.3730250,0.3729932,0.3727007,0.3726170,0.3722539,0.3721743,     0.3721055,0.3720965,0.3718959,0.3714824,0.3711862,0.3709115,0.3708262,     0.3706647,0.3701728,0.3697490,0.3692868)))          t<-as.matrix(A$V1)      My questions are: **1.- How can I fit my data to a gamma distribution?** **2.- Is there any way to determine, according to the model fitting, which signal value corresponds to p-value < 0.05?** Thanks for your help! PD 1: I accept suggestions in order to edit as this is my first post in the community! PD 2 : I've been trying all of this using R


What free tool can I use to do simple Monte Carlo simulations on OS X?


title: Are there any available implementations of density or conditional density tree learning?

text: Greetings, I am working on joint and conditional density trees for approximating clique potentials in Bayesian Belief Networks. A brief introduction to topic is available from this paper: http://www.autonlab.org/autonweb/14653.html in case you'd like to get a better description what I'm talking about. I am looking for an implementation that supports both discrete and continuous variables in a joint probability distribution. The leaf nodes in the learned tree would provide probability distributions for both discrete and continuous variables. I'd be very glad to find an open source implementation, but any implementation in any form would work for me, at least as a test for my own implementation (which I'd have to do anyway) I've seen this approach described in various papers in different contexts, but I could not find any implementations. Kind regards Seref


title: Power Analysis with Existing Data Set

text: I am new to this forum but have found several threads to be highly useful so am posing a question myself. My data was collected ( **fish length = factor** , **fish mercury = response** ) from several rivers over several years for the purpose of environmental (mercury) monitoring. What I would like to do is, using the data I have, perform power analysis to determine how many samples should be collected in the future to get the same results. The purpose is to recommend the least number of samples necessary (thus killing the least number of fish). Is this possible and does anyone have a recommendation on how to go about this? Much thanks in advance!


title: Identifying the experimental unit

text: I have included a diagram of an experiment where cells are pooled from multiple animals, cells from the pool are then allocated to different treatments and a measurement is taken from each individual cell in each treatment. The aim of this experiment is to determine how doses of a drug (treatments) affect cell size. I have previously received different opinions on what constitutes the experimental unit in this case. I would say the experimental unit is the pool so n=1 but when looking at the definition of experimental unit it states: “The smallest division of the experimental material such that any two experimental units can receive different treatments” - surely two individual cells could theoretically receive different treatments since they are separate entities? Most examples I have found are more straight forward, for example had the snails themselves been exposed to the treatment I could see how the cells within each snail would no longer be independent and thus the snail would be the experimental unit. I would appreciate any thoughts or clarification. ![Study design](http://i.stack.imgur.com/m1wol.jpg)


title: R find key influencers

text: In Excel's data mining tools there is a "Key Influencers" tool which will look at a dataset which is perhaps customers and whether or not they converted to a given goal (e.g. a flag that equals 1). It then tells you the most influential factors in reaching that goal,( e.g. Gender=Male and Age=30-45). What would be equivalent algorithm in R to achieve a similar outcome.


title: Using least squares straight line fitting on experimental data with errors

text: I'm trying to fit a line into a graph with points that have errors both in the y and x directions. I found the following document on the subject: https://www.che.udel.edu/pdf/FittingData.pdf My problem is that I don't quite understand what the author means on page 8 with: > To use this method, you must first establish the standard error of each > point; this quantity will be called e_i. If you made repeat measurements, > you should use the standard deviation." How can I compute the standard error for each point exactly? Could someone give an example, how I come up with the standard errors if I have, say, the next set of data:               x = [2, 4, 6], x_error = [0.2, 0.15, 0.1]     y = [1, 3, 5], y_error = [0.1, 0.12, 0.2]      I want to at least take into account the y_error when fitting the line.


title: Is there a way to optimize regression according to a specific criterion?

text: I am curious, is there a way to optimize a regression according to a specific statistic? Let's say I am interested in a model with the best possible AIC statistic (or MSE or whatever measurement I am interested in) - could I somehow direct the regression to give me the top X models that would do this? (Of course I would not ignore the other measures, but would it be possible to ask for this?) What software supports this, or would you write your own code (in say R)? Also in general, when multiple regression results are displayed (let's say an all-possible or best-subset regression is done), is there a ranking and if so, according to what measure/criteria? I'm not saying that this is the best way to evaluate models, but perhaps it would be a way to explore candidate models? (This is really not my main question though.) Thanks .. learning a lot reading this site.


title: How do you decide what your train, validation and test percentages are?

text: When splitting up my labeled data into training, validation and test sets, I have heard everything from 50/25/25 to 85/5/10. I am sure this depends on how you are going to use your model and how prone to overfitting your learning algorithm is. Is there a way to decide or is it all by rule of thumb? Even ELSII seems vague on the subject.


title: How to calculate Beta and coefficient of determination ($R^2$) from unstandardized coefficients in OLS regression?

text: I have a table in which the multiple linear regression results is provided. If I have unstandardized coefficients and standard error for each independent variable, is it possible to calculate standardized coefficient (Beta) and coefficient of determination ($R^2$) from these data. My friend provided $R^2=0.41$ for these data, but I doubt if the results are reliable. I brought the table in the following, would you please compute Beta and R-squared for me to compare with the original table? ![enter image description here](http://i.stack.imgur.com/LJGFb.jpg) the descriptive statistics is the following table: ![enter image description here](http://i.stack.imgur.com/re6g2.jpg)


title: Question on sampling error

text: I have problem understanding following example in parentheses. Maybe you can help me. Ideally, researchers know the full extent of the population they want to study, and they can select a sample from this population at random. Statisticians can calculate the probability that such random samples represent the population; this is usually expressed in terms of sampling error ( **for example, there might be a 95 percent probability that the distribution of responses in a sample will be within 3 percent of the distribution in the population** ).


title: Simulating C4.5 algorithm using Decision Tree node in SAS Enterprise Miner 6.2.

text: In SAS Enterprise Miner 6.2 it's possible to approximate CHAID and CART methods using Decision Tree node, according to SAS Help, but there is nothing about C4.5 algorithm. How can I mimic C4.5 algorithm using Decision Tree node? I would be grateful for any help.


title: How to add odds? (or how I got the wrong sandwich)

text: Jon and Frank ordered 2 footlong sandwiches, but instead got 4 6-inch sandwiches. Jon had ordered spicy italian, and Frank chicken teriyaki. The sandwiches aren't labeled and are wrap, so they can't tell which sandwiches are which. What would be the odds of Jon picking the 2 sandwiches that are are his own? (the 2 spicy italian sandwiches)? I know that there is a 50/50 chance of choosing a correct sandwich, since there are 2 terikaki and 2 spicy, but once 1 is chosen, there is only 1/3 chances to get the second one right.. but how do the the odds of getting the first 50/50 choice AND the second 1/3 choice correct add up? Thanks!


title: Estimating maximum predictive power in noisy data

text: When assessing model performance, one would like to separate prediction errors due to limitations of the model from those errors due to intrinsic noise. For example, in noisy data, AUC in an ROC curve can be less than 1.0 even when using the "true" model. What are some approaches to get a quantitative estimate of a model's maximum predictive power given a certain level of noise?


title: "Central limit theorem" for weighted sum of correlated random variables

text: I'm reading a paper which claims that $$\hat{X}_k=\frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}X_je^{-i2\pi kj/N},$$ (i.e. the Discrete Fourier Transform, DFT) by the C.L.T. tends to a (complex) gaussian random variable. However, I know this isn't true in general. After reading this (fallacious) argument, I searched over the net and found this 2010 paper by Peligrad & Wu, where they prove that for **some** stationary processes, one can find a "C.L.T. theorem". My question is: do you have any other references that try to address the problem of finding the limiting distribution of the DFT of a given indexed sequence (both by simulation or theory)? I'm particularly interested in the convergence rate (i.e. how quickly the DFT converges) given some covariance structure for $X_j$ in the context of time-series analysis, or derivations/applications to non-stationary series.


title: Dealing with poor fit in an Item Response Theory model

text: I'm studying an online course with about 3000 students who each took several quizzes and I'm trying to apply Item Response Theory (using the ltm package in R) to model the questions, determine which items are most or least important to retain and what new items may be required, and to provide an alternate way to score and compare student performance. However, whenever I attempt to model all 16-20 questions on a single quiz together, I discover that the model is a poor fit based on the two- and three- way margins, with many Chi-squared residuals exceeding 10 and some exceeding 20, even with a three-parameter model. I tried two strategies to deal with this. In one I tried partitioning the questions into smaller groups based on similarity of the tested concepts. However I found in order to get acceptable fit (no more than 2-3 margins with Chi-squared over 3.5), I had to limit partition size to 3-5 items at most. After doing this, I have little insight into exactly what latent variable is being measured in each partition and thus how to interpret low discrimination, or how to assemble the scores for each partition into a test score. I also am not sure how to compare items in different partitions regarding importance. The other strategy I tried was to use "ltm" to fit to a two-dimensional latent variable space. This produced reasonable Chi-squared residuals (assuming the usual 3.5 cutoff remains reasonable in the 2D case) but my attempts to qualitatively understand what is being measured by each variable, either by looking at item parameters, 3D ICC graphs, or maximum posterior scores, failed. Lacking this understanding, I'm not sure how to combine the two latent variable scores into a single test score. It may be that the tests I'm working with are just too short and cover too broad of an area to adequately apply IRT to (covering e.g. 2-3 topics over a multiple choice quiz of 16 items). I hope this is not the case. It's also possible I'm being a bit too concerned with fit. I'm new to Item Response Theory and would appreciate any insight on this.


title: Understanding the output of C5.0 classification model using the CARET package

text: The `C5.0` classification model was used in this 4-class problem data with $N_{train}$=165, $P$=11, using `caret` R-package by running the code below. The winnowing option was tuned over in the model, which is a kind of feature selection approach. This excerpt I quote regarding winnowing from the companion book of `caret`, a must-have book in my opinion to realize hidden gems coded in the package:   _Kuhn M, Johnson K. Applied predictive modeling. 1st edition. New York: Springer. 2013_. > C5.0 also has an option to winnow or remove predictors: an initial algorithm > uncovers which predictors have a relationship with the outcome, and the ﬁnal > model is created from only the important predictors. To do this, the > training set is randomly split in half and a tree is created for the purpose > of evaluating the utility of the predictors (call this the “winnowing > tree”). Two procedures characterize the importance of each predictor to the > model: 1. Predictors are considered unimportant if they are not in any split > in the winnowing tree. 2. The half of the training set samples not included > to create the winnowing tree are used to estimate the error rate of the > tree. The error rate is also estimated without each predictor and compared > to the error rate when all the predictors are used. If the error rate > improves without the predictor, it is deemed to be irrelevant and is > provisionally removed.               c50Grid <- expand.grid(.trials = c(1:9, (1:10)*10),                            .model = c("tree", "rules"),                            .winnow = c(TRUE, FALSE))     c50Grid     set.seed(1) # important to have reproducible results     c5Fitvac <- train(Class ~ .,                        data = training,                        method = "C5.0",                        tuneGrid = c50Grid,                        trControl = ctrl,                        metric = "Accuracy", # not needed it is so by default                        importance=TRUE, # not needed                        preProc = c("center", "scale"))       > c5Fitvac$finalModel$tuneValue        trials model winnow     16     70  tree  FALSE        CV tuning output:   ![enter image description here](http://i.stack.imgur.com/yHzDV.png) Excerpt from the C5.0 tree output:               > c5Fitvac$finalModel$tree     [1] "id=\"See5/C5.0 2.07 GPL Edition 2014-01-22\"\nentries=\"70\"\ntype=\"2\" class=\"Q\" freq=\"9,16,60,80\" att=\"IL17A\" forks=\"3\" cut=\"0.92485309\"\ntype=\"0\" class=\"Q\"\ntype=\"2\" class=\"Q\" freq=\"0,4,59,80\" att=\"IL23R\" forks=\"3\" cut=\"0.26331303\"\ntype=\"0\" class=\"Q\"\ntype=\"2\" class=\"Q\" freq=\"0,4,19,80\" att=\"IL12RB2\" forks=\"3\" cut=\"0.41611555\"\ntype=\"0\" class=\"Q\"\ntype=\"2\" class=\"Q\" freq=\"0,4,9,80\" att=\"IL23R\" forks=\         Now importance of predictors:               > predictors(c5Fitvac )      [1] "IL23R"   "IL12RB2" "IL8"     "IL23A"   "IL6ST"   "IL12A"   "IL12RB1"      [8] "IL27RA"  "IL12B"   "IL17A"   "EBI3"      **Questions:**   1. Why is it in the plot, the `accuracy` levels of No-winnowing about two times that of the winnowing? Can you please help interpreting this output when it says winnow = FALSE?    2. How to visualize the tree output, instead of the computed junk text that appeared in my case? is there any way to _behold_ a tree instead of crowded symbols?


title: Does MLE require i.i.d. data? Or just independent parameters?

text: Estimating parameters using maximum likelihood estimation (MLE) involves evaluating the likelihood function, which maps the probability of the sample (X) occurring to values (x) on the parameter space (θ) given a distribution family (P(X=x|θ) over possible values of θ (note: am I right on this?). All examples I've seen involve calculating P(X=x|θ) by taking the product of F(X) where F is the distribution with the local value for θ and X is the sample (a vector). Since we're just multiplying the data, does it follow that the data be independent? E.g. could we not use MLE to fit time-series data? Or do the parameters just have to be independent?


title: Estimation issues for OLS with bounded response variable

text: Suppose we have a first differences regression with a bounded response variable, $$\Delta y_{i} = \beta_1 \Delta X_{i}+\Delta\epsilon_i.$$ For example, suppose $y_{it} \in [0,1]$ where $y_{it}$ could represent participation of a population on a given survey, or the pass rates of students on a standardized test in a given school $i$ for a given year $t$. $X_{it}$ is some variable that changes over time, for example if $y_{it}$ is pass rates in a given school, $X_{it}$ can be average class size. Our goal is to get a unbiased estimate for $\beta_1.$ If we estimate the above model by using OLS, the following problem could occur. For simplicity suppose we have two periods, $t =1, 2.$ Then if $y_{i1} = 1$ then $\Delta y_{i} \in [-1,0],$ whereas if $y_{i1} = 0$ then $\Delta y_{i} \in [0,1].$ Therefore if $y_{i1}$ is close to $1$, it can only increase by little. Whereas if $y_{i1}$ is close to $0$, it can increase by a lot. But realistically this can lead to biased results. This is because if the $y_{it}$ was not bounded, an increase in $X_i$ in the second period may increase $y_i$ in second period, but this is impossible in our case with bounded response when $y_{i1} = 1.$ The issue is similar if increasing $X_i$ leads to a decrease in $y_i$ when $y_{i1} = 0.$ Hence I suspect OLS might not be a appropriate estimator in this case, but under other circumstances it may work? The idea I have to fix this problem is to take into consideration in the model the following. Realistically given a low initially $y_{i1}$ (less than 0.5) it is easier to have a higher $y_{i2}$ due to change in $X_i$, then to increase $y_i$ given initially $y_{i1}$ high (more than 0.5). I am not sure how we would exactly do this, I suspect using a indicator variable might work. Would this type of approach work? Any other ideas to fix this capping issue?


title: Testing whether there is an increase between two regression slopes in a time series

text: I have a data of sea levels from a particular station from 1940 till today. The linear regression on the data from 1940-1980 resulted in a slope of 3. The linear regression on the data from 1980-2010 resulted in a slope of 4. What test can be used to establish whether the increase in the slope is significant?


title: Reconciling notations for mixed models

text: I am familiar with notation such as: \begin{align} y_{ij} &= \beta_0 + \beta_i x_{ij} + u_j + e_{ij}\\ &= \beta_{0j} + \beta_i x_{ij} + e_{ij} \end{align} where $\beta_{0j}=\beta_{0}+u_j$, and \begin{align} y_{ij} &= \beta_0 + \beta_1 x_{ij} + u_{0j} + u_{1j} x_{ij} + e_{ij} \\ &= \beta_{0j} + \beta_{1j} x_{ij} + e_{ij} \end{align} where $\beta_{0j}=\beta_{0}+u_{0j}$ and $\beta_{1j}=\beta_1+u_{1j}$ for a random intercepts model and a random slope + random intercepts model, respectively. I have also come across this matrix/vector notation, which I have been told is "mixed model notation for grown ups" (according to my elder brother): $$ \mathbf{y}=\mathbf{X\beta} + \mathbf{Z b} + \mathbf{e} $$ where $\mathbf{\beta}$ are the fixed effects and $\mathbf{b}$ are the random effects. If I have understood correctly, the latter notation is a more general notation for the former which are specific versions of the latter. I should like to see how the former can be derived from the latter.


title: What is an acceptable proportion of events in logistic regression?

text: In market research I'm building a logistic regression model to estimate the likelihood that clients may change banks. The proportion of events is roughly 10% in my sample. From university I remember that a proportion of events that is too small introduces bias into the estimate. Or is it the standard error that gets biased? (Question 1) As a rule of thumb, what is an appropriate proportion of events? (Question 2) Please give a straight and rough answer, rather than a technical explanation supported by various references.


title: What distribution to choose? Do Kernel estimations help?

text: I recently did some maximum likelihood estimations where the dependent variable justified the use of a normal distribution. Now, however, the dependent variable has a skewness of 0.4 and a kurtosis of 1.5 and also looks quite strange. Now my question is: How to choose/find an appropriate distribution? I read about Kernel Density Estimates but do not know whether this could yield a distribution-formula which could be used in a maximum likelihood estimation.


title: Interpreting summary function for lm model in R

text: What is the meaning of `t value` and `Pr(>|t|)` when using `summary()` function on linear regression model in R?               Coefficients:                                   Estimate Std. Error t value Pr(>|t|)         (Intercept)                    10.1595     1.3603   7.469 1.11e-13 ***     log(var)                        0.3422     0.1597   2.143   0.0322 *


title: How fit a regression model with two time series

text: I want to assess the effect of temperature on fish mortality. Both variables are time series, so the residual will be autocorrelated. Is there a regression method to deal with that? Any advice will be appreciated


title: Error metrics for cross-validating Poisson models

text: I'm cross validating a model that's trying to predict a count. If this was a binary classification problem, I'd calculate out-of-fold AUC, and if this was a regression problem I'd calculate out-of-fold RMSE or MAE. For a poisson model, what error metrics can I use to evaluate the "accuracy" of the out-of-sample predictions? Is there a poisson extension of AUC that looks at how well the predictions order the actual values? It seems that a lot of Kaggle competitions for counts (e.g. number of useful votes a yelp review will get, or number of days a patient will spend in the hospital) use root mean log squared error, or RMLSE. /Edit: One thing I've been doing is calculating deciles of the predicted values, and then looking at the actual counts, binned by decile. If decile 1 is low, decile 10 is high, and the deciles in between are strictly increasing, I've been calling the model "good," but I've been having trouble quantifying this process, and I've convinced there's a better approach. /Edit 2: I'm looking for a formula that takes predicted and actual values and returns some "error" or "accuracy" metric. My plan is to calculate this function on the out-of-fold data during cross-validation, and then use it to compare a wide variety of models (e.g. a poisson regression, a random forest and a GBM). For example, one such function is RMSE = sqrt(mean((predicted-actual)^2)). Another such function would be AUC. Neither function seems to be right for poisson data.


title: Gibbs sampling from full conditionals

text: I have the following joint density: $p(x_1,x_2,y_1,y_2) \propto \exp\left(−\left(x_1^2+x_2^2+c_1(y_2-y_1)^2+c_2(y_2-y_1)^4\right)\right)$ Can I use Gibbs sampling to sample from that? How can I get the full conditional distributions, and how can I draw samples from those conditional distributions which are quite complicated I think?


title: Determine density of min(x,y) and max(x,y) for independently uniform distributed variables

text: Two independent random variables, X and Y, are uniformly distributed on the unit interval (-1,1). Determine the density for U=min(X,Y) and for W=max(X,Y)


title: cokriging and generalised linear spatial models

text: The answer to this question gave me an interesting idea. The part that interested me was a brief description of kriging concluded by > The way, kriging as usually practiced is not quite the same as least squares > estimation, because Σ is estimated in a preliminary procedure (known as > "variography") using the same data. That is contrary to the assumptions of > this derivation, which assumed Σ was known (and a fortiori independent of > the data). Thus, at the very outset, kriging has some conceptual and > statistical flaws built into it. Thoughtful practitioners have always been > aware of this and found various creative ways to (try to) justify the > inconsistencies. (Having lots of data can really help.) Procedures now exist > for simultaneously estimating Σ and predicting a collection of values at > unknown locations. They require slightly stronger assumptions (multivariate > normality) in order to accomplish this feat. The question I had from this was is there a similar method for cokriging like how generalized spatial models are an improvement upon kriging? Basically a technique that doesn't assume anything on Sigma or mu and uses the secondary variate.


title: How to deal with two time variables in panel data modelling using PLM package in R?

text: We have question regarding the analysis of panel data in R (plm) which was observed over two time scales. As way of a background, we have completed an economic experiment with a group of 6-8 students, which is completed over 18 sessions (time 1) with 12 rounds per session (time 2). We have three treatments and each treatment is replicated three times (3 session * 12 rounds = 36 times). In each session there are six players making decisions each round (as some are removed prior due to a bidding process). We have set out our experimental design below in Table 1 for greater clarity. We have also included a table showing what the dataset looks like (Table 2). The treatments test three different scenarios (3 levels) across three different types of fisheries (3 levels). So our question is how should we deal with my two time variables (session and round) given that some participants play in more than one round and session? We understand that we cannot ignore rounds, as a player's decision will be influenced by what happens in previous rounds, in the current and previous sessions. The plm function doesn’t allow you to analyze two time variables; therefore, should we set one time variable as a factor(e.i. `as.factor(sessions)`)? Any help would be appreciated. ![enter image description here](http://i.stack.imgur.com/ag4Hy.jpg)


title: nearest neighbors degrees of freedom

text: For polynomial fitting with a polynomial of degree $n$, we have $n$ degrees of freedom. Is there a similar concept for $k$ nearest neighbors? Is there any way to compare the degrees in general? I come from an EE background so I might not familiar with some of the advanced terminology.


title: Loss sequence for glmnet

text: Is the loss sequence returned by cv.glmnet (the slot "cvm") on the standardized variable scale or the original scale? This only affects the penalty part of the loss, I think. In other words, are models obtained by calling cv.glmnet with standardize=TRUE comparable to those obtained with standardize=FALSE?


title: Two-way anova with interaction term: what is the point of a post-hoc test?

text: I have a statistics interpretation question. I've recently performed a two-way anova to identify an interaction term between my categorical independent variables (genotype + temperature) that influences my continuous dependent variable (speed). My hypothesis is that genotype+temperature strongly interact to reduce speed. In my case, the interaction term is significant - ie., genotype and temperature significantly interact, and I can easily observe that speed is reduced in this case. I am now directed by common practice to perform a post-hoc test, but I don't really understand why: What post-hoc test is appropriate to help confirm the _interaction_ between temp and genotype? I am not interested in anything except the interaction, as I consider all of the other conditions to be controls. That said, is such a post-hoc test necessary, given that the ANOVA itself is designed to (and succeeds) in revealing a significant interaction? Is it necessary to determine the **sign** of the interaction? _Bonus points if you can point me to a nonparametric test for this, as my data is not gaussian (but my large N=100+ helps me be comfortable with using the data in the ANOVA itself)._


title: Help needed with wording: significant effect size but overlapping CIs

text: I have an issue with describing the results of a simple analysis in a self- consistent way. This is for a study of an intervention that took place in 25 different towns, to reduce the number of cot deaths that took place. Data is for the number of cot deaths before the intervention started (the 'baseline' phase), and for the number after the intervention was rolled out. Both baseline and intervention phases went for a year. I've done a Poisson regression for the number of cot deaths, adjusting for clustering at the town level. I get a significant effect (reduction by 20%, with 95% C.I 9 - 30.2), but the issues arise when I try to state the estimated rate of cot death before and after the intervention. Adjusted again for clustering at the town level, I get 113 per 100,000 (95% CI 96 - 130) at baseline, and 92 per 100,000 (95% CI 75 - 120) after the intervention. So although the effect size is significant, the confidence intervals for the rates overlap! I'm assuming that comparing the confidence intervals on the rates would be misleading because they arise from marginal distributions and do not reflect the pairing between towns from the baseline to intervention phase. But how should I describe these results to avoid this misleading conclusion?


title: Given known bivariate normal means and variances, update correlation estimate, $P(\rho)$, with new data?

text: I'm dealing with two correlated random variables which are modeled via a bivariate normal distribution. I have values for the means ($\mu_x, \mu_y$) and individual variances ($\sigma_x, \sigma_y$) of the variables (effectively the marginalize normal distributions). Given new observations of the random variables, could I update my estimate via a Bayesian inference for the correlation coefficient, $P(\rho)$? **Thoughts so far:** Starting with the sampling distribution: $P(x,y|\rho) = \frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_x)^2}{\sigma_x^2} + \frac{(y-\mu_y)^2}{\sigma_y^2} - \frac{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} \right] \right)$ We want to update the prior, $P(\rho)$ (not sure of the prior form) with some observation $(x_i,y_i)$: $P(\rho|x_i,y_i)=P(x_i,y_i)^{-1}*P(x_i,y_i|\rho)*P(\rho)$ What is throwing me off is the form of the prior distribution. I know for my process that the $\rho$ is bounded due to positive correlation. How can I factor this information into the definition of the prior? Additionally, could I implement this procedure in a recursive manner? **Other thoughts:** Perhaps I am overthinking this problem entirely. I could just compute the Pearson correlation coefficient directly given all measured data. However, my expected sample size will be small, and I also know bounds for the correlation ($0<\rho<1$). Given the small sample and a priori known constraints, I still am leaning Bayesian.


title: Model stability in cross-validation of regression models

text: Given multiple cross-validation folds of a logistic regression, and the resulting multiple estimates of each regression coefficient, how should one measure whether or not a predictor (or set of predictors) is/are stable and meaningful based on the regression coefficient(s)? Is this different for linear regression?


title: Correlation vs. causality

text: Correlation does not imply causation. Causation does imply correlation but not necessarily linear correlation. ... So does correlation imply high-order causation? If A and B are correlated, is it always possible to find `K={K1 K2 ... Kn}` variables such that               A~K1, K1~K2, ...,Kn-1~Kn and Kn~B      where ~ denotes causality. In other words, if there is no causality in a correlation, does a common- causal variable(s) always exist?


title: Benchmark dataset for decision tree algorithm

text: I'm implementing a decision tree algorithm, and I'd like to get a feel for how it performs relative to other implementations. Can anyone recommend popular datasets for training and testing decision tree algorithms? I've found some resources like this, but I'm not sure widely used it is.


title: SVAR Model with Short run restrictions

text: I am currently working on implementing `SVAR` model in an economic analysis. I have 10 variables in my analysis and currently struggling to incorporate the short run restrictions. Here is an example with four variables say Y1 (GDP), Y2(Exchange rate), Y3 (interest rate) Y4(Bank risk). I wanted to implement the following short term restrictions. ![enter image description here](http://i.stack.imgur.com/tqXEW.png) I went through the some textbooks and the examples were starting with short run ordering. Unlike here where I do not have an order and hence making the matrix non triangular.


title: Logical reasons for choosing regression through the origin

text: Is it reasonable to choose a regression model with a value of `0` for the intercept when this makes logical sense? For example, I am trying to model a physical geometric relationship, and I know that when `x = 0`, `y = 0`. Yet the consequences of choosing such a model are that the `R^2` value becomes significantly higher (it changes from `0.67` to `0.95`). When I examine the residuals for both models, I can see that they both have roughly the same distribution. The origin option is shown in Figure 1 and non-origin in Figure 2. ![enter image description here](http://i.stack.imgur.com/iFZFo.jpg) ![enter image description here](http://i.stack.imgur.com/XvWbi.jpg) How should I decide which model is more appropriate? I've read through some of the other questions and answers on this topic but I haven't seen any discussion about physical limitations providing the basis for the choice. My dependent variable here is an area calculation, and my independent variable is a measurement of one dimension of the shape. For example, if I had a set of rectangles of length `l`, width `w` and area `A`, I am trying to model the relationship between `l` and `A`. However, as these are not perfectly regular there is some variation in the relationship but it appears to be linear in several cases, and based on some of the comments, not so much in this particular instance.


title: Most probable value given predictions with confidence levels

text: I'm new to http://stats.stackexchange, so please let me know if I'm doing something wrong. I have the following problem: Suppose that there is a car traveling along a flat surface. There are no roads on this surface, so the travel can be thought of as random. There are five people in a helicopter flying above this surface, all of whom are observing this car from the helicopter, through the same window (in the helicopter). After some time, each person is asked to predict where the car will be in one second. Each person gives a prediction, along with how confident they are of that prediction i.e. person `P` says "I am `C%` confident that in one second from now, the car will be at location `L`). Given that there are five such people, how can I determine the most probable location of the car one second from now? I don't think this really matters, but being a stats n00b, I'd rather provide this information than not: the number of people in the helicopter is variable. Five is just an example for this case. It's not fixed. **More Information** : The observers in the helicopter have all been trained for such a prediction task. However, each observer specializes in predicting certain types of behavior. For example, one observer might be trained in predicting all types of left turns, etc. Assume that human error does not exist and that all observers tell the truth. Also assume that they all have the same objective scale of measuring their own confidences. A note on confidences: Confidence does not refer to a confidence interval, which talks about the range of possible/probable values for a variable. A bad example (but one that illustrates the underlying point): if observer A says "I will bet \$100 on my prediction being correct" and observer B says "I will bet \$150 on my prediction being correct", all other factors (such as how much each observer values their money) kept constant, observer B is more confident than observer A. Assume that all observers have been properly trained to correctly use a particular objective scale by which to measure their confidence levels.


title: Replace t-test with modified z-test - anyone ever seen this?

text: I just read about a procedure I'd never heard of, and was wondering if anyone has had any experience with it. In Donald Berry's "Statistics" book, he presents an alternative to doing a t-test by modifying the sample standard deviation, with a factor that increases with lower-N: $s\rightarrow s\left(1+\frac{20}{n^2}\right)$ and then doing a simple z-test. Has anyone ever seen this, or seen a derivation of this? It seems to empirically match the t-distribution test at the 99% level, and overestimates the standard deviation for lower percentiles. I doubt that it would make a lot of practical difference, but if justified, it might help students in intro stats classes. thoughts?


title: Convergence in Probability of Empirical Median

text: I'm stuck with this one. Let $X_1,...X_n$ be an i.i.d. sequence of random variables with CDF F. The empirical CDF of $X_i$ is defined $$ \hat F_n(x) = \frac{1}{n} \sum_{1 \leq i \leq n} I\{X_I \leq x \}. $$ Note that for $x \in \Re$, $\hat F_n(x) \xrightarrow{p} F(x)$. Also, define the smallest median of P as $$ \theta_0 = \inf \{ x \in \Re : F(x) \geq .5 \}. $$ Suppose that the median is unique, i.e., for any $\epsilon > 0$, $P(X_i \leq \theta_0 + \epsilon) > .5$. Define an estimator $\hat \theta_n$ of $\theta_0$ by $$ \hat \theta_0 = \inf \{x \in \Re : \hat F_n(x) \geq .5 \}. $$ Show that $\hat \theta \xrightarrow{p} \theta_0$.


title: How Can I use some variables selected by LASSO?

text: I am very new about statistics. So, please understand if my question is somewhat awkward, and please give me related any advice. I have some data set. X = 500 x 100 (500 observations x 100 predictor variables) Y = 500 x 1 (500 response variables) With that data, I first estimated mean square error (MSE) using Ridge method with 10 fold cross validation. Then, I wanted to find the subset of predictor variables (sparse solution) that might have more impact on the estimation. For this, I used LASSO method with also 10 fold cross validation, and I picked up the best coefficient set that gave me the least MSE value. At this point, can I directly compare the least MSE value estimated by LASSO method and the MSE value estimated by Ridge method? If it is reasonable and if the MSE value estimated by Lasso method is smaller than the MSE value estimated by Ridge method, can I say that the new model constructed by LASSO method would be better than using the one constructed by the whole predictor variables (Ridge method)? I have another option in mind. How about applying Ridge method again to the predictor variables selected LASSO method to compare the performance between the subset and whole predictor variables? I don't know this is reasonable or not. If you know some papers related to my questions, please let me know. Thanks in advance!


title: Ballot problem with ties

text: In an election, candidate $A$ receives $n$ votes and candidate $B$ receives $m$ votes where $n>m$. The probability that $A$ is always ahead in the count of the votes is $\dfrac{n-m}{n+m}$. This is known as Bertrand's ballot theorem. A variation is to allow for ties. So for $n\geq m$ what is $P(A\text{ is never behind})$?


title: Balanced cross validation

text: I understand what is K-fold cross validation, but what does balanced folds mean? I thought k-fold cross validation is balanced. What are the differences between the two?


title: Probability of always more heads than tails

text: A coin that comes up heads with probability $p$ is flipped $n$ consecutive times. What is the probability that starting with the first flip there are always more heads than tails that have appeared? I have seen hints that this is Bertrand's ballot problem in disguise . Here is my attempt: Let $X$ be the number of heads in $n$ tosses with a coin of probability $p$. Let $Y$ be the number of tails. Let $s=\text{floor}(n/2)+1$ $\begin{align*}P(\text{Always more heads}) &= \sum_{k=0}^n P(\text{Always more heads}|X=k)P(X=k)\\ &= \sum_{k=s}^n P(\text{Always more heads}|X=k)P(X=k)\\ &= \sum_{k=s}^n \dfrac{k-(n-k)}{k+n-k}{n \choose k}p^k(1-p)^{n-k}\\ &= \sum_{k=s}^n \dfrac{2k-n}{n}{n \choose k}p^k(1-p)^{n-k} \end{align*}$ That doesn't seem very simple too me. Is there a better way to do this? Thanks.