title: Estimating values of a sequence from observed differences

text: I have a sequence of random variables $S_1, S_2 \dots S_N$ that is guaranteed to satisfy $$S_1 + S_2 + \cdots + S_N = 0$$ I can't observe any of these random variables directly, however I _can_ observe their differences: $$S_{ij} = S_i - S_j$$ I typically have many observations of each difference $S_{ij}$ - call the observations $s_{ij}^{(1)}\dots s_{ij}^{(n_{ij})}$. Each observation has some error - I'm happy to assume that the error is the same for each observation, and is independent of $i$ and $j$. I would like to estimate the means of the $S_i$ (and possibly get a confidence interval for the means). I can form a matrix $\hat{S}_{ij}$ of estimated values for the mean of the differences in the following way: $$\hat{S}_{ij} = \frac{1}{n_{ij}} \sum_{t=1}^{n_{ij}} s_{ij}^{(t)}$$ The catch is that I don't observe every possible difference, so the matrix has some NA values (put another way, some of the $n_{ij}$ are zero). I can assume that the adjacency matrix of non-NA values has exactly one connected component. How can I come up with an estimates $\hat{S}_j$ for the mean of the original, non-differenced values?


title: Ordinal regression models vs. nominal regression models

text: Suppose we are trying to decide whether to use an ordinal regression model or a nominal regression model. If we know that the outcome variable is ordered should we still use a nominal regression model to compare with an ordinal model?


title: Trying to Refine SARIMA models

text: I have a SARIMA forecast from statewide Real Estate Sales data.. but I'm not happy with it. The SARIMA parameters are confusing to say the least. I am finding that the current model is not forecasting high enough, although the month by month fluctuations look reasonable. Currently the projected year ahead is 4.9% above this year. but the year over year gain for the current year is about 10%+, so in short the model is not weighted heavily enough to the recent year over year change. Between the ARIMA and Seasonal components I can't see how to fix/specify this? The Parameter descriptions remain a mystery no matter how many versions I read. Looking at the chart it is clear that the growth rate over the last three years is higher than the current forecast, Surely there is a simple way to deal with this sort of specification in SARIMA, Bonus if it is understandable by mere mortals!? R code follows (source data is public) I'm a real newb at R as well.               require(RCurl)     myCsv <- getURL("https://docs.google.com/spreadsheet/pub?key=0Ak_wF7ZGeMmHdFZtQjI1a1hhUWR2UExCa2E4MFhiWWc&single=true&gid=1&output=csv")     fullmatrix <- read.csv(textConnection(myCsv), stringsAsFactors = FALSE)     fullmatrix$date <- as.Date(fullmatrix$date, format = "%m/%d/%y")     fullmatrix$Unit.Sales <- as.numeric(gsub(",", "", fullmatrix$Unit.Sales))     fullmatrix$Average.Selling.Price <- as.numeric(gsub(",", "",  fullmatrix$Average.Selling.Price))     fullmatrix$Median.Selling.Price  <- as.numeric(gsub(",", "", fullmatrix$Median.Selling.Price))     fullmatrix$Total.For.Sale <- as.numeric(gsub(",", "", fullmatrix$Total.For.Sale))     order.date <- order(fullmatrix$date )         fullmatrix <-fullmatrix[order.date, ]         library(xts)         require(forecast)         require(astsa)         sarima(fullmatrix$Unit.Sales,12, 0, 1, 1, 0,1,0) # (slightly higher fore)     fore <- sarima.for(fullmatrix$Unit.Sales,12, 12, 0, 1, 1,0,1,12)#specifies same model          fullmatrix.xts <- as.xts(x=fullmatrix[,-1],order.by=fullmatrix$date)     unitsales.xts <- as.xts(as.numeric(fullmatrix.xts$Unit.Sales),order.by=index(fullmatrix.xts))         forecast.xts <- as.xts(as.numeric(fore$pred), order.by=as.Date(c("2013-11-01","2013-12-01","2014-01-01","2014-02-01","2014-03-01","2014-04-01","2014-05-01","2014-06-01","2014-07-01","2014-08-01","2014-09-01","2014-10-01")))     orig_fore <-rbind(unitsales.xts,forecast.xts)     colnames(orig_fore) <- "Unit.Sales"     plot(orig_fore)      ![Actual with Forecast](http://i.stack.imgur.com/v26XT.png)


title: Interpretation of R's lm() output

text: the help pages in R assume I know what those numbers mean. I don't :) I'm trying to really intuitively understand every number here. I will just post the output and comment on what I found out. There might (will) be mistakes, as I'll just write what I assume. Please correct me, and I will edit the wrong parts.   Mainly I'd like to know what the t-value in the coefficients mean, and why they print the residual standard error. I hope someone can clarify that.               Call:     lm(formula = iris$Sepal.Width ~ iris$Petal.Width)          Residuals:          Min       1Q   Median       3Q      Max      -1.09907 -0.23626 -0.01064  0.23345  1.17532       A 5-point-summary of the residuals (Their mean is always 0, right?). The numbers can be used (I'm guessing here) to quickly see if there are any big outliers. Also you can already see it here if the residuals are far from normal distributed (which they should be).               Coefficients:                      Estimate Std. Error t value Pr(>|t|)         (Intercept)       3.30843    0.06210  53.278  < 2e-16 ***     iris$Petal.Width -0.20936    0.04374  -4.786 4.07e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       Estimates $\hat{\beta_i}$ , computed by least squares regression. Also, the standard error $\sigma_{\beta_i}$ . I'd like to know how this is calculated.   Also, no idea where the t value and the corresponding p come from. I know $\hat{\beta}$ should be normal distributed, but how is the t value calculated?               Residual standard error: 0.407 on 148 degrees of freedom      $\sqrt{ \frac{1}{n-p} \epsilon^T\epsilon }$ , I guess. But why do we calculate that, and what does it say us?               Multiple R-squared: 0.134,  Adjusted R-squared: 0.1282       $ R^2 = \frac{s_\hat{y}^2}{s_y^2} $ , which is $ \frac{\sum_{i=1}^n (\hat{y_i}-\bar{y})^2}{\sum_{i=1}^n (y_i-\bar{y})^2} $ . The ratio is close to 1 if the points lie on a straight line, and 0 if they are random.   What is the adjusted R-squared?               F-statistic: 22.91 on 1 and 148 DF,  p-value: 4.073e-06       F and p for the **whole** model, not only for single $\beta_i$s as previous. The F value is $ \frac{s^2_{\hat{y}}}{\sum\epsilon_i} $ . The bigger it grows, the more unlikely it is that the $\beta$'s do not have any effect at all.


title: How to interpret the output of the summary method for an lm object in R?

text: I am using sample algae data to understand data mining a bit more. I have used the following commands:               data(algae)     algae <- algae[-manyNAs(algae),]     clean.algae <-knnImputation(algae, k = 10)     lm.a1 <- lm(a1 ~ ., data = clean.algae[, 1:12])     summary(lm.a1)      Subsequently I received the results below. However I can not find any good documentation which explains what most of this means, especially Std. Error,t value and Pr. Can someone please be kind enough to shed some light please? Most importantly, which variables should I look at to ascertain on whether a model is giving me good prediction data?               Call:     lm(formula = a1 ~ ., data = clean.algae[, 1:12])          Residuals:       Min      1Q  Median      3Q     Max        -37.679 -11.893  -2.567   7.410  62.190             Coefficients:                     Estimate Std. Error t value Pr(>|t|)          (Intercept)  42.942055  24.010879   1.788  0.07537 .        seasonspring  3.726978   4.137741   0.901  0.36892          seasonsummer  0.747597   4.020711   0.186  0.85270          seasonwinter  3.692955   3.865391   0.955  0.34065          sizemedium    3.263728   3.802051   0.858  0.39179          sizesmall     9.682140   4.179971   2.316  0.02166 *        speedlow      3.922084   4.706315   0.833  0.40573          speedmedium   0.246764   3.241874   0.076  0.93941          mxPH         -3.589118   2.703528  -1.328  0.18598          mnO2          1.052636   0.705018   1.493  0.13715          Cl           -0.040172   0.033661  -1.193  0.23426          NO3          -1.511235   0.551339  -2.741  0.00674 **       NH4           0.001634   0.001003   1.628  0.10516          oPO4         -0.005435   0.039884  -0.136  0.89177          PO4          -0.052241   0.030755  -1.699  0.09109 .        Chla         -0.088022   0.079998  -1.100  0.27265          ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1            Residual standard error: 17.65 on 182 degrees of freedom       Multiple R-squared:  0.3731,    Adjusted R-squared:  0.3215        F-statistic: 7.223 on 15 and 182 DF,  p-value: 2.444e-12


title: Interpreting output from anova() when using lm() as input

text: I am learning about building linear regression models by looking over someone elses R code. Here is the example data I am using:               v1  v2  v3  response     0.417655013 -0.012026453    -0.528416414    48.55555556     -0.018445979    -0.460809371    0.054017873 47.76666667     -0.246110341    0.092230159 0.057435968 49.14444444     -0.521980295    -0.428499038    0.119640369 51.08888889     0.633310578 -0.224215856    -0.153917427    48.97777778     0.41522316  0.050609412 -0.642394965    48.5     -0.07349941 0.547128578 -0.539018121    53.95555556     -0.313950353    0.207853678 0.713903994 48.16666667     0.404643796 -0.326782199    -0.785848428    47.7     0.028246796 -0.424323318    0.289313911 49.34444444     0.720822953 -0.166712488    0.323246062 50.78888889     -0.430825851    -0.308119827    0.543823856 52.65555556     -0.964175294    0.661700584 -0.11905972 51.03333333     -0.178955757    -0.11148414 -0.151179885    48.28888889     0.488388035 0.515903257 -0.087738159    48.68888889     -0.097527627    0.188292773 0.207321867 49.86666667     0.481853599 0.21142728  -0.226700254    48.38888889     1.139561277 -0.293574756    0.574855693 54.55555556     0.104077762 0.16075114  -0.131124443    48.61111111      I read in the data and use a call to `lm()` to build a model:               > my_data<- read.table("data.csv", header = T, sep = ",")     > my_lm <- lm(response~v1 + v2 + v3 + v1:v2 + v1:v3 + v2:v3, data=my_data)     > summary(my_lm)          Call:     lm(formula = response ~ v1 + v2 + v3 + v1:v2 + v1:v3 + v2:v3,      data = my_data)          Residuals:         Min      1Q  Median      3Q     Max      -2.0603 -0.6615 -0.1891  1.0395  1.8280           Coefficients:              Estimate Std. Error t value Pr(>|t|)         (Intercept)  49.33944    0.42089 117.226  < 2e-16 ***     v1            0.06611    0.82320   0.080  0.93732         v2           -0.36725    1.06359  -0.345  0.73585         v3            0.72741    1.00973   0.720  0.48508         v1:v2        -2.54544    2.21663  -1.148  0.27321         v1:v3         0.80641    2.77603   0.290  0.77640         v2:v3       -12.16017    3.62473  -3.355  0.00573 **      ---      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 1.375 on 12 degrees of freedom     Multiple R-squared:  0.697, Adjusted R-squared:  0.5455      F-statistic:   4.6 on 6 and 12 DF,  p-value: 0.01191      Following along with their code I then use a call to `anova()`:               > my_lm_anova <- anova(my_lm)     > my_lm_anova     Analysis of Variance Table          Response: response               Df  Sum Sq Mean Sq F value   Pr(>F)        v1         1  0.0010  0.0010  0.0005 0.982400        v2         1  0.2842  0.2842  0.1503 0.705036        v3         1  9.8059  9.8059  5.1856 0.041891 *      v1:v2      1  4.3653  4.3653  2.3084 0.154573        v1:v3      1 16.4582 16.4582  8.7034 0.012141 *      v2:v3      1 21.2824 21.2824 11.2545 0.005729 **     Residuals 12 22.6921  1.8910                         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      However, I am not sure:   1. Why I would use the call to ANOVA in this situation, and   2. What the ANOVA table is telling me about the predictor variables. From the code they appear to use the ANOVA table as follows. For predictor variable v1, the result of   * Adding the 'Sum Sq' entry for v1 together with half of the 'Sum Sq' entry for v1:v2 and half of the 'Sum Sq' entry for v1:v3,    * Dividing by the sum of the entire 'Sum Sq' column, and   * Multiplying by 100 gives the percent of variance of the response variable that is explained by predictor variable v1 in the `lm()` model. I don't see why this is nor why half of the 'Sum Sq' entry for v1:v2 is attributed to v1 and half to v2. Is this just convenience?


title: R statistics interpretation

text: I am learning exploratory statistical testing methods in R. Specifically I am trying to learn single variable linear regression. I was hoping that someone would be able to give an explanation of the info that is given as output. Disclaimer: I am a beginner in statistics and not that well versed so I apologize for the general nature of this question. My code (I am using a data set with 177 points): `plot(b$Months.On.Program, b$leads)` Output                Coefficients:         (Intercept)  b$Months.on.Program                0.4420               0.2537      I understand what the intercept means but I do not know what .2537 is..(slope?) Code: `fit<-lm(b$Leads~ b$Months.on.Program)` `summary(fit` Output:                   Residuals:         Min      1Q  Median      3Q     Max         -10.844  -3.740  -0.949   0.036  59.171              Coefficients:                     Estimate Std.      Error   t value   Pr(>|t|)           (Intercept)           0.4420     1.0368   0.426     0.67           b$Months.on.Program   0.2537     0.0483   5.253     4.31e-07 ***            Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1            Residual standard error: 8.616 on 175 degrees of freedom       Multiple R-squared:  0.1362,  Adjusted R-squared:  0.1313        F-statistic: 27.59 on 1 and 175 DF,  p-value: 4.313e-07      I realize that my lack of knowledge in basic statistics is hurting me here but I was hoping that someone would break down the information that I am receiving here. Any links to solid statistics resources would also be greatly appreciated. Things I am trying to understand: Is my linear regression accurate? (I suspect it isn't b/c of outliers) How do I test to see if it's accurate? If not how do I compensate for outliers? What other info can I glean from the summary statistics? Thanks in advance!


title: Explaining R linear regression output - in lay terms

text: Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) -57.6004     9.2337  -6.238 3.84e-09 ***     InMichelin    1.9931     2.6357   0.756    0.451         Food          0.2006     0.6683   0.300    0.764         Decor         2.2049     0.3930   5.610 8.76e-08 ***     Service       3.0598     0.5705   5.363 2.84e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       Above is the output of linear regression in R. How do I understand the `Std. Error`, `t value` and `Pr(>|t|)` values in lay man language (someone with very elementary knowledge of statistics)? More importantly, how do I interpret the 'statistical significance' for the weights.


title: R: Explanation of a multiple linear regression summary

text: I am quite new with R and while i am able to perform the basics i am not yet able to understand the output results. For example: `summary(lmodel)` generates the following:               Call:     lm(formula = temperature ~ altitude + sea.distance, data = meteodata)          Residuals:          Min       1Q   Median       3Q      Max      -1.40156 -0.46083 -0.04078  0.50961  1.49844           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  19.736000   0.257222  76.728  < 2e-16 ***     altitude     -0.016509   0.004006  -4.121 0.000714 ***     sea.distance  0.531111   0.414166   1.282 0.216927         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.7425 on 17 degrees of freedom     Multiple R-squared:  0.864, Adjusted R-squared:  0.848      F-statistic:    54 on 2 and 17 DF,  p-value: 4.317e-08      I know about the `?summary`, `help(summary)` and i am familiar with the basic concepts, but i would kindly ask for some thorough explanation of the summary in regards to my model (`residuals`, `coefficients`, `F-statistic`, `Residual standard error`, etc.). It seems that there are several similar answered questions in this forum, however none of them macthes exactly my case. Any help is welcome. Thanks in advance


title: Interpreting the output of linear regression

text: I used scipy.stats.linregress function for my dataset. It returns 5 output variables. > Slope = 0.0453 > > Intercept = 14.905 > > r-value = 0.91 > > p-value = 0 > > standard error of estimate = 0.00015 This is the documentation link if anybody would like to see Linear Regression Scipy Can somebody help in understanding the output? Questions I would like to ask are the two variables highly correlated. The r-value is something that I think points to the correlation, but since it is not 1.0, what would be ways to find out which ones do not satisfy this constraint, can I use any of the parameters here? In general when you see such values how to grasp them?


title: interpreting R lm(..) output when a variable is used as a factor

text: Below is a print screen of a summary(lm(..)). I called the response variable `response` explained by a continuous variable `X` and a factor `Y`. ![enter image description here](http://i.stack.imgur.com/mOEAG.png) I'm trying to understand (not to the point of understanding the algorithms used to compute these results) the meaning of this output. My questions might be duplicates, but I haven't found the answer anywhere or if I did, I couldn't understand it!   * The model is `response = -4.461010 + blabla..`. What does the p.value for this intercept mean? Does it mean that the predicted line of my data significantly does not pass by the Origin?   * `factor(Y)2` is significant. Does it mean that for the subset of `data2` where `Y` equals 2, the values of `response` are significantly different than in the subset where `Y`equals 1 (which is taken as a reference)?   * Imagine a graph where the lines for `subset(data2,Y==1)` and `subset(data2,Y==2)` look alike but the line for `subset(data2,Y==3)` is totally different. If `factor(Y)3` is taken as a reference, we'll get the two others to be significant and if `factor(Y)1` is taken as a reference, `factor(Y)3` will be significant but not `factor(Y)2`. Is it correct?   * For the interactions `X:Y2` and `X:Y3`. Does the low p.values mean that for both the subset `Y==2` and `Y==3`, the variable `X` acts differently than it does for the subset `Y==1`? Thanks a lot for your help!


title: Basic questions concerning the interpretation of results from summary(lm(...~...)) in R

text: set.seed(11)     a = runif (12)     b = rep(c(1,2,3),4)     summary(lm(a~b))$coeff     summary(lm(a~b-1))$coeff      What does a p.value for the intercept means ? What differences makes the second formula ? Why is the p.value lower and the R squared higher in the second formula ? What does the first column ("Estimate") means ? Thanks a lot for your help.


title: Book/Website recommentation on probability distributions usage

text: Is there a book(s) / Website(s) that list the different type of distributions and real world usage or application for each one of them (beyond the typical know ones such as normal, binomial, poisson, etc). Looking for practical application information that would allow to pick the right family for the problem that one would be addressing.


title: Can within- and between-subjects data ever differ systematically?

text: I've recently asked whether it is valid to treat repeated measures data as both within- and between-subject data and compare the analysis of both to see if there are differences. My question was motivated by a particular paper (PDF) about presence questionnaires that, _in my opinion_ , implied that such questionnaires are essentially not applicable in between-subject designs, as people would only be able to reliably rate their sense of presence in a particular environment relative to another environment. This assumption was questioned by some of the answers, and hence I was wondering whether treating the same data as between-subjects and within- subjects data can ever lead to a systematic difference. That is, is it possible that an effect will show up if I analyse data as within-subjects but not show up when analysed as between-subjects? By way of a very contrived example, assume I get subjects to rate on a scale of 1 to 10 how much they like an offer of free beer. In one condition, people are offered one litre of beer for free, in the second condition, people are offered two litres of beer for free. If I ran this experiment with a between- subjects design, I would personally assume that there isn't any discernible difference between the two conditions, because - hey, free beer! But as a within-subjects design, I think I stand a reasonable chance of seeing an effect, because more free beer is better than less free beer.


title: Rules of Thumb to choose an initial number of class intervals and refine that choice (potentially automatically)

text: I was wondering if there are established rules of thumb (or algorithms) that, given a set of observations can help:   1. choose an initial number of class intervals.   2. refine that choice to a better number. I could find talk of using square-root(N), where N is the number of observations as an initial guess of the number of class intervals. Thanks in advance.


title: Making sense of a factor analysis application

text: I'm a biochem guy so the stats in this paper are giving me a real hard time. > Kliman, RM, Irving, N, and Santiago, M. Selection Conflicts, Gene > Expression, and Codon Usage Trends in Yeast. J Mol Evol (2003) 57:98–109 > (PDF copy) In this paper I'm having a hard time figuring out exactly what they did to find factors 1&2\. On page 3, paragraph 4 they give only the briefest of explanations. They use the frequencies of synonymous codons for each AA class and mysteriously enter these values into factor analysis. What exactly does this statement mean?


title: Repeated measures: aggregate, fixed effects, or random effects?

text: I have a repeated measures design with data collected across 8 treatments, with 5-6 repeats for each treatment in closely spaced time intervals (time is not important in the study). I only have 4 subjects in the study, and one dropped out two-thirds of the way through. As far as I can see, I have three options to analyze these data: 1) use a mixed model with subject as a random effect. The main problem here is that there are only 4 subjects and so estimates of the variance of the random effects will be imprecise. 2) treat subject as a fixed effect. Someone suggested this to me, but I'm unclear how this accounts for the fact that the repeated measurements are not independent. Surely it doesn't? 3) collapse the data across the repeated measures so that the response variable is either the sum or mean of all repeated values for a given subject/treatment combination. This obviates the need for a repeated measures design, but it seems I lose some power. With only 4 subjects, I'm inclined to go with option 3, but I'd love to hear people's thoughts.


title: Comparing predictive validity of separate IVs on single DV

text: I'm attempting to understand some data and would greatly appreciate help in picking the appropriate analysis measure. **Context:**   * I am conducting psychometric analysis of the predictive validity of two separate self-report measures on a dependent variable (also self-report, also continuous).    * All participants completed both IV's and the DV.   * There is reason to believe that the subscales of each measure should contribute to the variance explained in the DV in a particular way. **For example:**   * Measure 1 (MA) has 5 subscales MA1, MA2, MA3, MA4, MA5.   * Measure 2 (MB) has 5 subscales MB1, MB2, MB3, MB4, MB5. If I only wanted to know the predictive ability of one measure I would use a hierarchical multiple regression as follows...   1. Control for gender   2. MA4 and MA5   3. MA3   4. MA2 **QUESTIONS:**   * What is the best way to assess the best predicting measure?    * Can I run separate heirachical MRs and compare beta values?    * Should I put the competing measures constructs in each step of the hierachical MR and look at what is a significant predictor?


title: When is regressing out factors serially, equivalent to a simultaneous model?

text: I have two factors $a,b \in R^M$ which I wish to regress out of my variables $Y \in R^{M \times N}$. $b$ is a function of the variables $Y$ and is defined as: $b = \frac{1}{N}Y1$ Regressing out the variables serially, the formula I have is: $$\begin{eqnarray} Y_{\perp a} &=& Y - m(m^Tm)^{-1}m^TY \\\ Y_{\perp a,b} &=& Y_{\perp a} - b_{\perp a}(b_{\perp a}^Tb_{\perp a})^{-1}b_{\perp a}^TY_{\perp a} \\\ \text{is equivalent to:} \\\ Y_{\perp a,b} &=& Y - [a,b]([a,b]^T[a,b])^{-1}[a,b]^TY \end{eqnarray}$$ Is there a formal way to prove that regressing out $[a,b]$ simultaneously to get $Y_{\perp a,b}$ is equivalent to regressing out $a$, then calculating $b_{\perp a}$ on the residual, $Y_{\perp a}$ and regressing that out to get the same $Y_{\perp a,b}$ . I've tested this empirically in matlab and this appears to be true, but I can't come up with a formal proof.


title: LETOR data sets, which contain the queries and documents?

text: While researching information retrieval evaluation, I found LETOR. This is a test collection used for evaluation of ranking algorithms. Since the data contains relevance judgements for queries, it can also be used as a test collection for information retrieval evaluation (precision, recall...). I checked the whole site, but didn't find the queries and the documents, which are the foundation for the relevance judgements. **Where can I find the actual data?** (I am afraid it is not freely available)


title: What is the method for calculating stdres in {MASS}

text: The documentation for the `MASS` package does not detail the calculation method for the standardized residuals using `stdres()`. I want to be sure that it matches the method I intend to use:               residuals/(stdevres*sqrt(1-hat))      Does someone know if this is the calculation applied in this function?


title: Non-nested model comparisons - how to interpret output of coxtest in R

text: I'm trying to compare the outputs of two non-nested models in my dataset and am using the `coxtest` in the R package `lmtest` and don't know how to interpret the output. I have two predictors which are just two ways of categorising the same categorical scale: one is more general than the other and collapses the factor levels into two categories. My general investigation is to look at whether 2 categories can account for variation as well as 5. Maybe the model with 5 gets a slightly better fit, but is this significantly better given that we need five factor groups to account for the behaviour, when 2 does a decent job? I'll use the `diamonds` dataset from `ggplot2` to show what I'm trying to do.               library(ggplot2)     diamonds$cut2 <- diamonds$cut #creating a new column to classify 'cut' in a different way     levels(diamonds$cut2) <- c("B", "B", "B", "A", "A")          mod1<- lm(price~cut+carat, data=diamonds)     mod2<- lm(price~cut2+carat, data=diamonds)      Then I run the `coxtest`:               library(lmtest)     coxtest(mod1, mod2)      And I get this output:               Cox test          Model 1: price ~ cut + carat     Model 2: price ~ cut2 + carat                     Estimate Std. Error     z value Pr(>|z|)         fitted(M1) ~ M2        0 1.1825e-10 -8.8620e-01   0.3755         fitted(M2) ~ M1     -947 1.3691e-10 -6.9171e+12   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      I've read the documentation of `lmtest` here, but still don't really know what these numbers are saying. Mainly, I'd like to know whether it is favouring one model or another, and if I can really argue that one is better. Also, when both rows are significant, does this mean that each model contributes different things?


title: How to simulate survival times using true base line hazard function

text: I want to simulate survival times from the model $\lambda(t|X_1,X_2) = \lambda_0(t) \exp(\gamma_1X_1+\gamma_2X_2)$ where covariates are given by $X_1\sim$Binomial(1,0.5), $X_2\sim$Uniform(0,1), the true parameter values are $\gamma_1=-0.5$, $\gamma_2=1.5$ and the true baseline hazard function is as follows $ \lambda_0(t)= \left\\{ \begin{array}{ll} \exp(-0.3t); & 0<t\leq 1 \\\ \exp(-0.3); & 1<t\leq2.5 \\\ \exp(0.3(t-3.5)); & t>2.5 \end{array} \right. $ The censoring time is from the exponential distribution with mean 2.5. For a simple baseline hazard function(e.g. $\lambda_0(t)=\lambda$), i used the following algorithm to generate survival times $$T=H_0^-\\{-\log(u) \times \exp(\gamma_1X_1+\gamma_2X_2)\\} =\lambda^-\\{-\log(u) \times \exp(\gamma_1X_1+\gamma_2X_2))\\}$$ where $u\sim$Uniform(0,1). The R codes are as follows               gamma <- c(-0.5, 1.5)     X1 <- rbinom(n=200,size=1,prob=.5)     X2 <- round(runif(n=200,0,1),2)     time <- (1/.3)*(-log(runif(200))/exp(gamma[1]*X1+gamma[2]*X2))     cen <- rexp(200,rate=1/2.5)     event <- as.numeric(time <= cen)      But i'm struggling to generate for the piecewise baseline hazard function that i mentioned above. Would anyone please help me to generate survival times for the above settings using R.


title: The usage of data mining in pharmaceutical companies?

text: I know that data mining applications are being used in pharmaceutical companies, but my question is: what do they use them for? Sometimes I read: "drug discovery", but how? How is it used for drug discovery?!


title: How to analyze survival data with possible temporal dependencies

text: I want to analyze some survival data. I have measurements of a biomarker (real-valued variable) before a first treatment, after that first treatment, and then after a second treatment (different to the first). I also have patient age and disease stage (ordinal variable) at enrollment. I am particularly interested in how values of the biomarker affect survival. All else being equal, I guess I'd run a Cox analysis. However, we have reason to believe that treatment efficacy may depend on the value of the biomarker, and that the first treatment may change the biomarker in some patients. For example, if the biomarker is high before treatment, the first treatment may be effective, but may decrease the biomarker, reducing the scope for the the second treatment to be effective. I could run a Cox analysis, for example using as explanatory variables the raw pre-treatment biomarker and relative changes from this value at the subsequent two time points. However, I've not been able to convince myself that this adequately models the possible temporal dependencies (does it?). How should I analyze this data? Our sample size is about 45, so I doubt there's scope for a complex model. Thanks in advance.


title: Identifying Instrumental Variables to use in tsls function in sem R package

text: I want to reproduce the Example from Introduction to Econometrics by G. S. Maddala in `R` using `tsls` function from `sem` package. I'm confused how to identify the instrumental variables for the equations. Any help will be highly appreciated. Thanks               ## ---- fm42 ----     fm42 <-       tsls(           formula = Q~R+RD+X         , instruments = ~R+RD+RS #Not Sure         , data    = Data2SLS         )          summary(fm42)     summary(fm42)$coeff      ![enter image description here](http://i.stack.imgur.com/8ua46.png)


title: gsem goodness of fit

text: does anyone have any advise how to obtain goodness of fit information after a gsem command (for a logit path model) in Stata verion 13? Or would you advise to use BIC and AIC instead? Thank you for any help! Kind regards.


title: Identify features corresponding to high singular/principal component values

text: In MATLAB, while SVD gives a diagonal matrix S of decreasing singular values, PCA gives a column vector LATENT of decreasing principal component values.   1. How can S be used to obtain a subset of features that best describe the data?   2. How can LATENT be used to obtain a subset of features that best describe the data? NOTE: QR decomposition in MATLAB gives a permutation vector E which enables the above.


title: How do I test to see if the demographic proportions between two unequally sized groups are statistically different?

text: Sample of 102 broken into two groups based on a factor of interest Group 1, $n=77$ Group 2, $n=25$ Trying to figure out how to test if the demographic (age, race/ethnicity, gender, education level, have kids, etc) break down is significantly different between groups. I've calculated proportions and it seems like most of it will be fine, but there are a few that look like they could be significantly different. I'm hoping to be able to run this in SPSS.


title: Can you combine z scores from similar but non-identical populations?

text: Suppose you have a group of people estimating the heights of trees. Suppose you have many different types of trees, some trees with a very low height, and others that are very tall. Suppose each person in the group estimates the height for each type of tree. Something like the following:                           |          col         Guesser |     Tree  1    Tree 2     Tree 3      -----------+----------------------------------               1 |        7m         39m        65m                2 |        5m         33m        64m                3 |       10m         45m        70m               4 |        9m         44m        70m      -----------+----------------------------------      Clearly, we may compute z scores for the guesses for each tree, and construct composite scores using covariance, as in the following: http://www.psychassessment.com.au/pdf/mod5.pdf But the above construction uses paired data, and the covariance of the pairs is used in the construction (for example, guesser 2 is consistently low, and guesser 3 is consistently high). What happens if you have one group of people guessing for tree 1, a second group of people guessing for tree 2, and a third group of people guessing for tree 3, as in the following:                           |          col         Guesser |     Tree  1    Tree 2     Tree 3      -----------+----------------------------------               1 |        7m                                       2 |        5m                                       3 |       10m                                      4 |        9m                                      5 |                   40m                          6 |                   35m                          7 |                   43m               8 |                            62m               9 |                            68m              10 |                            75m              11 |                            72m              12 |                            68m      -----------+----------------------------------      In this case, the data is not paired in a meaningful way. Is there any way to group the z scores you get for the different guesses into a single collection of z scores?


title: Post-hoc Tests for Fisher's Exact Test

text: I am running a series of Fisher's test to analyse some 2x2 contingency tables with small values. I have found significance. Is there anything equivalent to Cramer's V or the Odds Ratio or some sort of post-hoc test that I can use to tell me more about my results?


title: glm inflated error...why?

text: I'm pretty new to stats, so this may be dumb. I've been running a bunch of models on randomly generated data to try and develop my understanding of type 1 error. I've noticed that using `glm(family=binomial)` I get more type 1 errors than I should when giving a binomial input (a two-column matrix of success and failure). In the code below, the first loop generates a thousand logistic regression models for random data, but I get type 1 errors (p < .05) about 15% of the time! The second loop runs the a similar thing as a Bernoulli test (just single vector of zeros and ones in the y). Here I get what I want to see, about 50 type 1 errors per 1000 models. Can anyone explain this to me? I see that if I change the possible values for the random numbers in my success/failure y-matrix (the variable I call range here) I can lower the type 1 errors, but I don't understand why. This gives me about triple the number of type 1 errors that I expect.               #Binomial glm     fit.p=c()     for(i in 1:1000){         range=0:10         y=matrix(sample(range,2000,replace=T),ncol=2,nrow=1000)         x=rnorm(1000,100,50)               fit=glm(y~x,family=binomial(link='logit'))                        fit.p[i]=anova(fit,test='Chisq')[2,5]                                 }     print(length(which(fit.p<.05)))      This works fine. About 50 errors per 1000               #Bernoulli glm     fit2.p=c()     for(i in 1:1000){         y=sample(0:1,1000,replace=T)         x=rnorm(1000,100,50)               fit2=glm(y~x,family=binomial(link='logit'))                        fit2.p[i]=anova(fit2,test='Chisq')[2,5]                                 }     print(length(which(fit2.p<.05)))


title: What is the diff between singly censored and progressive censored data in survival analysis?

text: I have a question regarding survival analysis . To my understanding, the singly censored data are those if there is one point in time, i.e, say, if the patient died (bulb is still working?) after certain time or not. So, what is progressive censored data? When does it arise? I did a quick google search and the links are either advanced (theoretical) or unsatisfactory. Can anyone give an example in a simple language? The issue came about as Cox F-test says it cannot be used for "progressive censored" data (Lee 1992).


title: Channel capacity

text: > A binary channel is described with following conditional probabilities: > > $p(y_1|x_1)=1, \quad p(y_1|x_2)=p(y_2|x_2)=\frac{1}{2}$ > > Additionally, we know that "1" appears 13 times more frequently than "0" in > the output. What is a capacity of the channel? My question is if the additional piece of information is really relevant as I get two different solutions with and without using it ( I put $p(y=1)=13 \cdot p(y=0)$ )? P.S. Correct answer is when omit the tip.


title: Test for change over time in a time series

text: I have a time series of anomalies in some data, with 69 observations in total. I want to statistically test if the number of anomalies in the time series have increased over time? Can I use a t-test for the purpose, possibly by dividing my data into 2 or 3 groups, and then comparing these groups? In case, I'm completely off, please suggest a better test/approach. Statistics, unfortunately, isn't my strong suit.


title: Determining polynomial model coefficients forcing slope = 1 and intercept =0

text: I have two observational variables, $\text{C}_{\text{obs}}$ and $R$, both subject to measurement error. I believe that a model of the form $$\text{C}_{\text{mod}} = a_0 + a_1R + a_2R^2 + a_3R^3$$ would be a reasonable representation of the relationship between $\text{C}_{\text{obs}}$ and $R$. I want to determine the coefficient vector $[a_0, a_1, a_2, a_3]$ in this model. An often-cited reference in my field describes finding these parameters using a procedure that forces the slope and intercept of the linear relationship between $\text{C}_{\text{mod}}$ and $\text{C}_{\text{obs}}$ to be 1 and 0, but does not provide any details on the method. I have attempted to contact the author directly without success, and I have been unable to find references to such a procedure (but I may not be using the correct search terms). I am an IDL user and am familiar with R; I would appreciate it if someone could explain how to solve this problem or point me to IDL or R routines that might be applied.


title: A question on Monte Carlo

text: We simulate $E(\exp(z))$ by Monte Carlo method,where $z\sim{}N(0,1)$. For sample sizes $2^{16}$ and $2^{17}$, the variance errors are 0.00531 and 0.00364, respectively. The ratio of these two errors is approximately $1.5$. My questions are: (1) What information can we learn from this ratio? (2) Why should the ratio approach $\sqrt 2$ as the sample size approaches infinity?


title: Tool for calculating AUC Value

text: Is there any tool that can calculate the AUC value from a ROC curve if I already have true positive, true negative, false positive, false negative values.


title: Comparing scatterplots with lots of points

text: I have two sets of data of protein-protein interactions in a matrices entitled: s1m and s2m. Each DB and AD pair make an interaction and the one matrix looks like:               > head(s1m)          DB_num AD_num     [1,]      2   8153     [2,]      7   3553     [3,]      8   4812     [4,]     13   7838     [5,]     24   3315     [6,]     24   6012      I can then plot the density of the points basically showing where the points are the most concentrated: s1m: ![enter image description here](http://i.stack.imgur.com/ZX50C.png) s2m: ![enter image description here](http://i.stack.imgur.com/FDWBJ.png) The code I used in R to make these plots was:               z <- kde2d(s1m[,1], s1m[,2], n=50)     filled.contour(z)     z <- kde2d(s2m[,1], s2m[,2], n=50)     filled.contour(z)      I want to be able to somehow compare how simiarly these plots are rather than just looking at them by eye. Is there someway to do this? By the way, I know very little about statistics. These are very large datasets also, something like 10,000 points among a matrix of 15k by 15k.


title: Fitting HLM models with Heavy-Tailed distribution: robust aproach to lmer?

text: I've been using package lmrob in R to estimate models with heavy-tailed distribution of the residuals, which I cannot correct even with transformations of the dependent variable. I was now planing to step up the models to a HLM framework, but I would like to know if there is a procedure to estimate robust HLM (maybe in lmer package?), even if it is necessary to code something manually. Any insights are welcome here.


title: MLE estimation together with volatility estimation

text: When I am estimating a linear time series model (i.e. ARMA) with Maximum Likelihood Estimation, if I also estimate volatility with a model like GARCH (or some other model) and use that variance value in MLE procedure, do I improve the estimation of the mean in any way? Is the procedure an iterative one or one-shot, meaning that I fit an ARMA then run GARCH over residuals then plug that in then refit ARMA and so on until it converges, or is it just run GARCH first then run ARMA and end the procedure? Thanks


title: Effective way to speed-up my mass-univariate problem (i.e. a large set of per-point optimization problems)

text: I am new to numerical methods, and I have to solve a problem of medical imaging. My background is computer science. I have a naive, general question. **Problem Statement:** I have an extremly large set of points, called ' _voxels_ ' subsequently. The set is denoted by $V$. For each voxel $v\in V$, I have to solve a standard optimization problem related to the voxel $v$, _i.e._ , > Find the parameter $p_{min}$ s.t. $d(o^v,p_{min}) \leq d(o^v,p)$ holds for > all possible parameters $p$ Above, $o^v$ is an observation at the voxel $v$, $d$ is a problem-specific (non-linear in my case) objectiver function depending on observations and model parameters. Solving this optimization problem voxel by voxel gives a correct result. However, this naive approach is not feasible as the number of voxels in V quickly becomes large. **My question is as follows:** Is there a general approach to deal with such difficulty? I am trying to find some "global method" that gives approximate solution in reasonable time. However, whenever I type "global optimization problem", I will fall into articles talking about solving global minimum (v.s. local minimum) which is different from what I need. Any ideas? Thank you for your help.


title: Multivariate exponential smoothing and Kalman filter equivalence

text: Suppose the time-series $X$ is hidden state Gaussian random walk and we observe $Y = X + e$, where $e$ is gaussian white noise independent of $X$. The Kalman estimator of $X$ in this case has a steady state closed form solution and corresponds to an exponential moving average smoother with constant smoothing parameter. The optimal smoothing parameter looks like $\lambda = \frac{p}{1-p}$ where $p$ is a quadratic formula of signal to noise ratio between $e$ and $\Delta{x}$. See formula for closed form (by searching ' _kalman solution random walk noise_ '). If we have instead _multiple independent_ observations $Y_1, Y_2,...$, i.e. independent $e_1, e_2,...$ is there a closed form solution for the optimal estimator of $X$? What would it look like? If the $e_1, e_2$ were independent _and_ identically distributed (same standard deviation), I can imagine the best estimator for $X$ is to simply average the estimators obtained treating the problem as univariate for each $Y$. If the $e_1,e_2,...$ have different variance, then some $Y$ series should have low weightings in the overall estimator as their signal to noise ratio is poorer. Perhaps the closed form solution is linear with coefficients proportional to each $Y$ signal to noise ratio. Is there a known closed form solution? Google seems not to be too helpful on this problem.


title: Coefficient of Determination ($r^2$): I have never fully grasped the interpretation

text: I want to fully grasp the notion of $r^2$ describing the amount of variation between variables. Every web explanation is a bit mechanical and obtuse. I want to "get" the concept, not just mechanically use the numbers. E.g.: Hours studied vs. test score $r$ = .8 $r^2$ = .64   * So, what does this mean?   * 64% of the variability of test scores can be explained by hours?   * How do we know that just by squaring?


title: How to adjust for a continious variable when the value 0 is distinctly different from the others?

text: Lets say I want to regress a variable on a covariate which has distinct value zero but the all other values are following some smooth function? It could be e.g. Years of beeing a mom (once you are a mom you are very different from 0 years, i.e. the effect of going from 0 to 0.0001 year should be a jump, but from 0.0001 to 0.0002 sould be smooth). Or maybe exposre is 'Time of exposure to HPV' when we study cervix cancer? One simple thing would be to add two covariates 'Years of being mom' and 'Are you a mom (yes/no)', but then all values would be 'No' exactly when the other variable is '0'. Any ideas? Hierarchical model? If so, how?


title: What is the best mathematical transformation for a variable with many zero values?

text: > **Possible Duplicate:**   >  How should I transform non-negative data including zeros? I have a continuous variable that has many zeros values and is NOT normal, so I can't use parametric statistics on it. I tried using `ln` and root-squared transformations, but they didn't work.


title: How to preserve points near zero when taking logs?

text: > **Possible Duplicate:**   >  How should I transform non-negative data including zeros? I have a question: > Suppose we have a data point $(0,10)$. We then convert to a log scale. Excel > gets rid of this data point since the minimum x-value is $1$. How do we keep > this data point? **Edit.** Can't you change the x-coordinate to $1$? So instead we have $(1,10)$ instead of $(0,10)$. The rest of the values are not changed.


title: How to perform boxcox transformation on data in R tool

text: > **Possible Duplicate:**   >  How should I transform non-negative data including zeros?   >  Error in boxcox.default(y ~ x) : response variable must be positive I want to perform box cox transformation on data using R tool .If anybody can help me with examples or commands to show how it is done in R tool. Thank you in advance. My main problem is I want to take log of data but if data entries contain 0 then log 0 becomes infinite .To handle this situation I am searching for the past 3 days .I found one solution to use boxcox transformation..but again it is giving this error.


title: Regarding transformation of a variable to follow normality

text: When I plotted histogram of a variable to see whether it follows normality, I got a highly positive skewed graph, more like that of a Poisson/exponential distribution. Now, what kind of transformation should I do on the variable to make it follow normality?


title: How to handle data which contains 0 in a log transformation regression using R tool

text: > **Possible Duplicate:**   >  How should I transform non-negative data including zeros? I want to perform log transformation regression in R tool but the problem is that I don't know how to handle situations when data contains 0 since log 0 is infinite... Please help..


title: Data-transformation of data with some values = 0

text: > **Possible Duplicate:**   >  How should I transform non-negative data including zeros? I want to log-transform some of my data because the Levene's homogeneity of variances test rejected the null hypothesis for those sets, before performing a two-way balanced ANOVA. Now, some of the data I collected have null values, because they are related to my plants' tillers, and some of them didn't have any. So for my log-transformation, should I simply use the formula `=LOG(x+1)` (in LibreOffice Calc) so those null values go back to 0? Or the value of the constant I add depends on how the datasets look? For example, I have weights that range from 0 to 12 grams, other weight that range from 0 to 1.5 grams, and leaf area ranging from 0 to 565 cm². Would `=LOG(x+1)` be ok for data ranging from 0 to 565, but not for data ranging from 0 to 1.5? Here are some histograms to show how my data is distributed: ![enter image description here](http://i.stack.imgur.com/3DNEN.png) What do you think? Thanks in advance!


Error in boxcox.default(y ~ x) : response variable must be positive      I am getting this error in `R` when I am performing a Box-Cox transformation on data. Why is this error happening? Here is my data. This is a time series data and I have to perform logarithmic regression of the form: $$y=a+b(\log x_1)+c(\log x_2)$$ I need to find a, b, c and then, check if any such type of relation exists or not.


title: Machine learning to catch fraud

text: I work for a company that ships material (a couple of thousand shipments per day) around the world. In order to ship anything a customer has to declare the weight of a shipment and declare the contents To circumvent customs they under declare material for example T-Shirts might be the declared contents but the actual contents might be a Hugo boss t-shirts which are worth more. Normally they are caught by random inspections and black listed (I have this data set and it grows daily) if it is a re-occurring shipper but this is costly to honest users of the service and impacts the time taken for a shipment to arrive. Repeated offenders are blocked manually based on addresses but can and do sign up again with a different account varying the address enough so as not to be recognized. For example 1 fake street could be applied again to fkae street 1. A human can see the difference but we dont have the staff to go through it. I think machine learning could be the way forward I was wondering if anyone would have any ideas how to classify these shipments. Research papers, ideas, brain storming are all welcome. What i would like to do is use the caught shipments to try an identify customers who have tweaked the address slightly so as to direct customs appropriately and make our black listing more effective Thanks for your time


title: Example of a stochastic process that is 1st and 2nd order stationary, but not strictly stationary (Round 2)

text: This question follows from my previous question , where Robin answered the question in the case of weak stationary processes. Here, I am asking a similar question for (strong?) stationary processes. I'll define what this means (the definition can also be found here). Let $X(t)$ be a stochastic process. We say that $X(t)$ is Nth-order stationary if, for every $t_1, t_2, \dots, t_N$ we have that the joint cumulative density functions $$F_{X(t_1),X(t_2),\dots,X(t_N)} = F_{X(t_1 + \tau),X(t_2 + \tau),\dots,X(t_N + \tau)}$$ for all $\tau$. This is quite a strong condition, it says that the joint statistics don't change at all as time shifts. For example, a 1st order stationary process is such that $F_{X(t_1)} = F_{X(t_2)}$ for all $t_1$ and $t_2$. That is, the $X(t)$ are all identically distributed. It is quite easy to see that a 1st order stationary process need not be 2nd order stationary. Simply assign a correlation structure to say $X(t)$, $X(t+1)$, $X(t+2)$ that _does not_ correspond to a (symmetric) Toeplitz matrix. That is, in vector form, the covariance matrix of $[ X(t), X(t+1), X(t+3)]$ could be given as $$\left[\begin{array}{cc} \sigma^2 & a & b \newline a & \sigma^2 & c \newline b & c& \sigma^2 \end{array}\right]$$ for $a,b,c$ distinct. This is now not 2nd order stationary because $E[X(t)X(t+1)] = a$ and, time shifting by 1 we have $E[X(t+1)X(t+2)] = c \neq a$. In a similar way (presumably) a process that is 1st and 2nd order stationary need not be 3rd order stationary and this leads to my question: > Does somebody have a nice example of a stochastic process that is both 1st > and 2nd order stationary, but not 3rd order stationary?


title: Bandwidth value for vector of equal values in kernel density estimation

text: How to define the value of bandwidth when we have a vector where all the values is the same? In this case, IQR and sample variation are both equal zero and Silverman's rule also result in zero. Using R, we have for $\vec{v} = (1,1)$ the bandwidth of 0.7835               density(c(1,1))$bw      and for $\vec{v} = (10,10)$, the value 7.834955               density(c(10,10))$bw


title: Fastest SVM implementation

text: More of a general question. I'm running a rbf SVM for predictive modeling. I think my current program definitely needs a bit of a speed up. I use scikits learn with a coarse to fine grid search + cross validation. Each SVM run takes around a minute, but with all the iterations, I'm still finding it too slow. Assuming I eventually multi thread the cross validation part over multiple cores, any recommendations on speeding up my program? Any faster implementations of SVMs? I've heard of some GPU SVMs but haven't digged into it much. Any users and is it faster? Thanks.


title: Estimating same model over multiple time series

text: I have a novice background in time series (some ARIMA estimation/forecasting) and am facing a problem I don't fully understand. Any help would be greatly appreciated. I am analyzing multiple time series, all over the same time interval and all of the same frequency, all describing a similar type of data. Each series is just one variable, there are no other corresponding predictors that I'm looking at. I have been asked to estimate a single model that describes ALL the series - for example, imagine I could find one ARIMA (p,d,q) with the same orders, coefficients, etc. that could fit all the series. My supervisor does not want me to separately estimate each series, nor does he want me to do some sort of VAR model with dependencies between the series. My question is: what would I even call such a model, and how might I go about estimating / forecasting it? If it's easier for you to use code examples, I speak both SAS and R.


title: Mackey-glass regression experiments

text: So I'm running some experiments using long-short term memory networks for temporal sequence regression, and one of the sequences I have to use is generated using the mackey-glass differential equation. To make it comparable to another temporal sequence prediction algorithm, we have to use the same protocol as [this other paper] for training and testing, as always. Thing is, mackey glass seems to be deterministic since there are no random variables, and [this other paper] trains their algorithm on the first 6000 points (with the initial transient cut off). They then test on 50 new sequences of 250 points each. But I'm hard pressed to understand how you have 50 "new" test sequences, other than just generating >250*50 new samples (presumably with t>training set samples), and then just chopping that up whichever way. Any ideas? Context: ![enter image description here](http://i.stack.imgur.com/mPc6m.png)


title: Similarities between different size matrices, rescaling problem

text: Given a series of matrices {$M_i$($m_i\times n_i$),`i=1...k`,$m_j,n_j \in$`random`} if we rescale (resize) all matrices into a certain size e.g., `pxq` then the following questions appear:   1. Is the summation of all matrices element-wise statistically legal for any type of variable would be in the cells? Think about additive and non-additive variables and consider they are associated with different supports in terms of size.   2. Is finding correlation matrix legal between those matrices?   3. Does correlation care about additivity?   4. Any other statistics / method to investigate similarities between those matrices? **Updates:**   The following figure shows the procedure.   ![enter image description here](http://i.stack.imgur.com/TohDA.png) **Note 1:** that the variable of interest is non-additive which means the value measured for a specific size of sample is not linearly relevant to what could be measured for different size. **Note 2:** one example for non-additive variables could be considered the variable of energy consumption for softening a block of rock as an initial stage in mineral processing. It is experimentally proved that the energy used for softening a small block has no linear relationship with larger one. While for grade (the concentration of an element) it is correct to rescale the measurement from small block for larger ones, this is not valid for energy consumption since it is non-additive.


title: How do we make a distinction between standardized and unstandardized coefficient?

text: In `SPSS`, regression analysis yields two types of coefficients (standardized and unstandardized). What is a standardized coefficient in statistical sense?


title: Is chi-squared always a one-sided test?

text: A published article (pdf) contains these 2 sentences: > Moreover, misreporting may be caused by the application of incorrect rules > or by a lack of knowledge of the statistical test. For example, the total df > in an ANOVA may be taken to be the error df in the reporting of an $F$ test, > or the researcher may divide the reported p value of a $\chi^2$ or $F$ test > by two, in order to obtain a one-sided $p$ value, whereas the $p$ value of a > $\chi^2$ or $F$ test is already a one-sided test. Why might they have said that? The chi-squared test is a two-sided test. (I have asked one of the authors, but gotten no response.) Am I overlooking something?


title: Which performance measure to use when using SVM: MSE or MAE?

text: It is a common practice to measure an SVM model's performance by calculating its **MSE** (Mean Square Error). Why not use **Mean Absolute Error** (averaging errors' absolute values instead of squared values)?


title: Compare two nonparametric curve estimation approaches: kernel and orthogonal basis

text: I know there're at least two approaches for nonparametric curve estimation: kernel and orthogonal basis. What are their advantages and disadvantages over each other ? And what are the typical application scenarios of each approach ? Could anyone give a hopefully not-too-long explaination ?


title: Statistical comparison of 2 independent Cohen's ds

text: Is it possible to statistically compare the Cohen's $d$ of two studies in order to determine of one of the two Cohen's $d$s is statistically significantly larger (or smaller) than the other one? I have calculated Cohen's $d$ for two virtually identical studies. In one study the cohen's $d$ for control group vs. treatment group is $1.05$; in the other study it is $1.31$. In the second study Cohen's $d$ is larger than in the first study. But can I perform a test to determine if the difference is statistically significant? If a statistical test is not possible is there a method / guidelines that would allow me to say if a difference of $0.26$ between the two Cohen's $d$s is small, medium or large (in other words relevant)?


title: Calculate the variance and the 1st order autocorrelation

text: The stochastic difference equation is $$ {{y}_{t+1}}=\frac{1}{a}{{y}_{t}}-\frac{b}{a}{{x}_{t}}+{{\varsigma }_{t+1}} $$ where   $ {a}>1; $   $ {{x}_{t}}=(1-\rho )\overline{x}+\rho {{x}_{t-1}}+{{\varepsilon }_{t}}, 0<\rho<1; $   $ {\varsigma }_{t+1} $ is a forecast error, that is, $ {{y}_{t+1}}={{E}_{t}}{{y}_{t+1}}+{{\varsigma }_{t+1}} $. $ {{\varsigma }_{t}} $ is serially uncorrelated, and uncorrelated with ${\varepsilon }_{t}$. $ {{E}_{t}} {{\varsigma }_{t+1}}=0 $.      How to calculate the variance and the first order autocorrelation of $ {y}_{t} $?   The reference answers are,   $ \begin{align} & \sigma _{y}^{2}=\frac{{{b}^{2}}(a+\rho )}{({{a}^{2}}-1)(a-\rho )}\sigma _{x}^{2}+\frac{{{a}^{2}}}{({{a}^{2}}-1)}\sigma _{\varsigma }^{2} \\\ & {{\rho }_{y}}(1)=\frac{1}{a}\left[ 1+\frac{{{b}^{2}}\rho ({{a}^{2}}-1)\sigma _{x}^{2}}{{{b}^{2}}(a+\rho )\sigma _{x}^{2}+{{a}^{2}}(a-\rho )\sigma _{\varsigma }^{2}} \right] \\\ \end{align} $      My answers are different from the above ones, could someone give me a detailed derivation? Thanks a lot.


title: Bivariate Normal distribution and correlation

text: Is the CDF of a bivariate normal distribution with mean $(0,0)$ and $\Sigma = ((1,\rho),(\rho,1))$ monotone in the correlation coefficient $\rho$?


title: Calculate p values based on differences between bootstrapped 95% confidence intervals

text: I have bootstrapped model fits to some data from different three groups, using the $2.5$ and $97.5$ percentiles to generate a $95\%$ confidence interval. I know that if $95\%$ confidence intervals do not overlap then there there is a significant difference between values of at least $p<0.05$. I would like to calculate the exact $p$ values for pairwise comparisons between the following groups:                          Mean     LowerCI   UpperCI      Group A:  0.05113   0.02279   0.07926          Group B:  0.10979   0.07954   0.14352          Group C:  0.10937   0.09400   0.13040      Any help would be greatly appreciated!


title: Appropriate statistical test

text: I know little about statistics. Help me, StackExchange! The software I'm working on does what I believe may be a dubious test. Suppose we measure three values (not real data): $$A) Mean = 10.5, sigma = .2, 15 \ samples $$ $$B) Mean = 12.3, sigma = .1, 17 \ samples $$ $$C) Mean = 22.5, sigma = .3, 10 \ samples$$ The hypothesis is that $A+B=C$. Let's assume a normal distribution for now. What is the appropriate test here, and could you work through the numbers in the example? What the code currently does is compare $ABS(C-A-B)$ with $z*\sqrt(Sa*Sa+Sb*Sb+Sc*Sc)$ (where Sx are the sigmas above and z is typically 3). We also apply the same pattern when testing $A+B+C+D=E$, for example (summing the variance as above and comparing against the residual). It feels to me that simply summing the variances like this is not appropriate -- surely the sample count should come in somewhere? (BTW, the sample count can range anywhere from 3 to 1000, depending on the situation, if that makes a difference). Is this test still appropriate as the number of summed values increases? In a similar vein, here's another test our code performs:   $$A) Mean = 10.5, sigma = .2, 15 \ samples$$   $$B) Mean = 52.3, sigma = .3, 4 \ samples$$ The hypothesis is that $5*A=B$. A similar test is applied, where $|B-5*A|$ is compared to $z*\sqrt(5*Sa*Sa+Sb*Sb)$ (the theory being that since A is used 5 times, the sigma should be included 5 times as well?). Anyway, all of this seems suspect to me, but I don't know what I'm doing.


title: Take the residual of a logistic regression model and use it as a response variable in another model

text: In linear regression we can do something like this $$ y' = resid(y \sim \beta_0 + \beta_1 c1 + \beta_2 c2 ) $$ where $c1, c2$ are covariates, and then fit $y'$ in another model: $$ y' \sim b_0 + b_1 x_1 $$ The advantage of this approach is when you have a lot of independent variables like $x_1, x_2, \cdots x_n$ (for example, a large number of DNA mutation sites), and they all share the same covariates (e.g. age, sex, weight etc), you can simplify things and reduce computation. But it looks like this will only work with continous variables, like height, blood pressure, etc, what if I have a binary variable, such as tumor/no-tumor, infected/uninfected etc, is it still possible to do something similar?


title: How can I interpret a regression model with a variable that appears as an inverse?

text: I have an equation where $$Employ = \beta_1 + \beta_2 * \frac{1}{GDP} + e $$ Then how can I explain the relationship between Employ and GDP? For example when GDP increase a certain unit or percentage, measure the change in Employ? Is there any other way to interpret this type of regression model? @AdamO: Thanks. I just want to know if there existed a way to explain the relationship between Y and X in the regression model : $$ Y=\beta_1 +\beta_2\frac{1}{X}+\epsilon $$. (such as something similar to the interpretation of level-level model (or log-log one): $$ Y=\beta_1 +\beta_2X+\epsilon $$: if X increases by 1 unit then Y rise by $$\beta_2$$ unit @all: thank you!


title: Transforming a relationship from quadratic to linear

text: Assume that I'm running a path analysis and I have discovered that certain relationships, empirically, are quadratic rather than linear. In order to model the data as such, I want to transform the relationship so that it becomes linear. Then, do I transform both variables involved through taking their square roots and then running the path analysis again (in an SEM software such as AMOS), or do I just transform one of the variables? For example, Imagine I have three variables arranged like this: A -> B -> C The relationship between A and B is linear, and the relationship between B and C is quadratic. if I want to model the data as such, do I then only transform variable C through taking it's square root? What is the general principle for which variables to transform, when what I actually want to transform is the relationship between variables?


title: Classes of distributions closed under maximum

text: Let $Q_p$ be a class of probability distributions on non-negative reals parametrized by $p$, so that $$ Q_p([0,\infty)) = 1. $$ I wonder which known classes of distributions are closed under taking the maximum and, i.e. if $X_1\sim Q_{p_1}$ and $X_2\sim Q_{p_2}$ are independent then $\max(X_1,X_2)\sim Q_{p_3}$.


title: What are the dangers of using a (log)normal distribution for a discrete response?

text: I have seen in the engineering field some papers (one example) using normal or lognormal distributions to model discrete outcomes. Typically, the explanatory variable is binned (into equal intervals) to allow for each point to represent a probability to belong to a given outcome (i.e. number of points in the bin/total, for each bin). Then the set of probabilities obtained is assumed to belong to a normal (lognormal) distribution as follows: $$ P(Y) = \phi[((\log)X-\mu) / \sigma], $$ so $$ X = \sigma\phi^{-1} + \mu. $$ Which allows them to evaluate the distribution parameters using OLS regression. Because of the behaviour of the inverse of the normal distribution function at 0 or 1 they also have to dismiss all points that have such probabilities, which obviously leads to much data not being used, and this is one issue. Plus discrete outcomes should be modelled using appropriate distributions, i.e. binomial and multinomial. But I cannot get my head around the other shortcomings of using such a method for discrete outcomes, I feel there must be many. Can anyone provide more insight?


title: Wondering what this bean plot analysis chart means

text: How does one interpret and generate bean plot charts. Here is one example taken from Walkes et al. 2010. What kind of data is it most useful for? ![bean plot](http://www.biomedcentral.com/content/figures/1471-2407-10-503-1.jpg)


title: Joint and Conditional probability

text: I'm trying to prove a result for Gibbs Sampling with multiple latent variables. I am not sure what the expansion of both the joint and conditional probability would be. In particular, let's say that I have two latent variables x and y (joint probability) with conditional probability to all the variables except x and y and the Dirichlet parameters a and b. Therefore I would have something like: Pr(x,y|w,a,b) How would you expand this? Would the probability chain rule apply here?


title: JAGS: unobserved parents and intialization

text: I have a very simple multinomial logistic regression problem; it's a simplified version of the alligator lake/food/size BUGS example. In my case I have (effectively) lake/food and have just stripped down the alli.bug example to handle two dimensional information: worker and IP address. There are 18 categories of 'worker' and 23 categories of 'IP address'. Unfortunately I am getting an error                 Error in node X[1,1:23] Observed node inconsistent with                           unobserved parents at initialization      Granted, there is a strong possibility I've more than one bug, but I understand that this is typically associated with an initialization problem. Unfortunately, my initialization matrices are rather large. A shotgun approach to finding appropriate values seems unlikely to succeed within a reasonable period of time. Any hints or suggestions about possible next steps would be greatly appreciated. EDIT: Here's the code. As I mentioned, it's just a stripped down version of the Alligator example. I've also fiddled with the subscripts on X[i,] a bit so the error message may be slightly different (will refer to node X[1,1] rather than X[1,1:23]).               var      X[I,K],      # observations I=no of workers=18, K=no of jobs=23      n[I],        # total for each covariate pattern      E[I,K],      # fitted values      OlogOE[I,K], # O log O/E      G2,            # goodness-of-fit statistic      phi[I,K],    # exp (beta[k] ' x[i,j])      p[I,K],      # fitted probabilities      alpha[K],      # factor for jobs = 2,...,23      beta[I,K],     # factor for workers = 2,...,18  for each job      b[I,K];         # factor for workers = 2,..,,18, relative to job 1, centred     model {     # PRIORS        alpha[1] <- 0;       # zero contrast for baseline workers        for (k in 2:K){           alpha[k] ~ dnorm(0,0.001);  # vague priors        }        # Loop around jobs:        for (k in 1:K){           beta[1,k] <- 0;   # corner-point contrast with first worker         }         for (i in 2:I) {                beta[i,1] <- 0;   # zero contrast for baseline worker           for (k in 2:K) {              beta[i,k] ~ dnorm(0,0.001); # vague priors           }        }          # LIKELIHOOD             for (i in 1:I) {                 # loop around workers           for (k in 1:K) {           # loop around jobs        # Multinomial response              X[i,k] ~ dmulti(p[i,k] , n[i]);                 p[i,k]        <- phi[i,k] / sum(phi[i,]);                 log(phi[i,k]) <- alpha[k] + beta[i,k];              }           }            # FITTED VALUES        for (i in 1:I) {     # loop around worker              for (k in 1:K) {     # loop around jobs                 E[i,k] <- p[i,k] * n[i];                 OlogOE[i,k] <- X[i,k] * log( X[i,k] / E[i,k] );              }         }        G2 <- 2 * sum( OlogOE[,] );     }


title: Constructing a Bayesian network from the begining

text: I am trying to construct a bayesian network which detects fraud.I have got a huge data set, with elements such as country, top spend merchant just to name a couple. The first part of my problem is to use a test set of data which is about 11,000 rows in excel and use this to firstly use Bayes formula to see the $P(\text{Fraud}|\text{each variable})$. After I have done this I have to perform this on the rest of the data and see if it is a good algorithm. This is the first part that I am confused with. How do I do this? I know $P(\text{Fraud}|X_1)=P(X_1|\text{Fraud})\frac{P(\text{Fraud})}{P(X_1)}$ and I dont understand why I have been told to do this on a test set of data.


title: Converting a spearman correlation to a euclidian dissimilarity

text: I am applying ward hierarchical clustering on a data set for which I have pairwise similarities. Since hierarchical clustering need a dissimilarity matrix, I am trying to convert my similarity matrix into a dissimilarity one. Besides, ward algorithm needs a euclidian dissimilarity, so I tried several conversions like the one proposed here `Warning: ward's linkage specified with non-Euclidean dissimilarity matrix` And, when I compute the cophenet, I get a poor value (around 0.39, smaller than all the values I previously got). How can I convert my spearman similarity matrix into a euclidian dissimilarity one?


title: Survival Analysis -- Probability Definition and Distribution application

text: edit: I think really the crux of my question is more meta, really how should I think of survival data and how do I know what kind of distribution I should fit into it? Are there any resources on how real world data can be applied with probability distributions beyond the normal? Like if I wanted to use a weibull distribution for this data, how does that actually happen? \-- I am trying to create a survival curve for a subscription company, where a user signs up, receives a weekly product, and then can cancel at any point in the future. What is the best way to set up the data? I currently aggregate everyone by weeks in existence, and not calendar time, e.g. everyone has had a week[1], everyone that has been around longer than a year has a week[52], everyone that canceled three weeks after signing up, whether they signed up in december or august, canceled at week[3]. The data that I am feeding into R's K-M survival package is structured as customer_id, weeks of survival, is_cancel (as a 1 for cancel or 0 for not) I know that Kaplan Meier looks at, for every week t, 1-(deaths[t]/survivals[t]). The company currently looks at churn in terms of survivals[t]/total[t], where total[t] would equal survivals[t]+deaths[t]. These curves have very different shapes, with the Kaplan Meier version looking logarithmic, while the current company model looks like much more like a*x^-b. Does anyone know if there's a standard way to look at customer churn? Is one way better than another, especially if I were to fit a curve so as to predict future survival probabilities to eventually model LTV/CLV? I would appreciate any insight. Thanks!


title: Please provide an example of when bootstrap has less bias than classically approximated estimates?

text: The recent question "Why does my bootstrap interval have terrible coverage?" has got me wondering if anybody has some really good examples of distributions in which bootstrapping standard errors systematically outperforms classic estimators (I am not sure what the right terminology for the classic set of estimators... perhaps moment approximation estimators?). In order to respond to the comment that bootstrapped estimators might be less biased by Ben Ogorek, I have also modified the code (originally posted by Flounderer) to include an estimate of bias and mean squared error. After repeating the simulation 10,000 times, it appears to me that the mean squared errors are statistically identical. Thanks for your consideration in this matter!               tCI.total <- 0     bootCI.total <- 0     m <- 10 # sample size          Trep <- 10000 # number of repetitions of the proceedure     Brep <- 1000 # number of repetitions of the bootrap          sampv <- mbootv <- rep(0,Trep)          true.mean <- exp(2) + 1          for (i in 1:Trep){       samp <- exp(rnorm(m,0,2)) + 1       sampv[i] <- mean(samp)            tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(m)            boot.means <- rep(0,Brep)       for (j in 1:Brep) boot.means[j] <- mean(sample(samp,m,replace=T))            mbootv[i] <- mean(boot.means)            bootCI <- sort(boot.means)[c(0.025*length(boot.means), 0.975*length(boot.means))]            if (true.mean > min(tCI) & true.mean < max(tCI)) tCI.total <- tCI.total + 1       if (true.mean > min(bootCI) & true.mean < max(bootCI)) bootCI.total <- bootCI.total + 1      }     tCI.total/Trep     # estimate of t interval coverage probability     # 0.5634     bootCI.total/Trep  # estimate of bootstrap interval coverage probability     # 0.5416          # Let's look at bias esimate for the sample mean and the bootrapped population mean estimate     (true.mean - mean(mbootv)) # bias estimate of bootstrapped means     # 0.170623     (true.mean - mean(mbootv))^2 + sd(mbootv)^2 # mean squared error of bootstrapped means     # 198.5914          (true.mean - mean(sampv))  # bias estimate of sample means     # 0.170475     (true.mean - mean(sampv))^2 + sd(sampv)^2 # mean squared error of sample means     # 198.4912


title: Why discriminative models are preferred to generative models for sequence labeling tasks?

text: I understand that discriminative models, such as CRF(Conditional Random Fields), model conditional probabilities $P(y|x)$, while generative models, such as HMM(Hidden Markov Model), model joint probabilities $P(y,x)$. Take CRF and HMM for example. I know that CRF can have a larger range of possible features. Apart from that, what else makes CRF (discriminative models) preferable to HMM(generative models) in sequence labeling tasks such as Part-of-Speech tagging and NER(Name Entity Recognition)? **Edit:**   I found out that HMMs will have to model $P(x)$, while CRFs don't. Why would it make a big difference in sequence labeling tasks?


title: Which arma model is best one?

text: I am studying ARMA models. ARMA(25,25) or ARMA(1,1) which one better model? Why? I think that the reason is the ommission of irrelevant variable


title: How do I "nest" data in a mixed-effects model? and MANY related questions about mixed effects models!

text: I treated 80 people with drug X and 80 with drug Y. I presented the drugs to them in groups of 10 (meow groups). The drug was consumed for 15 weeks. They reported the severity of their headaches (from 1-100) once a week. All individuals are only in one meow group and everyone within a meow group received the same drug. The meow group that a person is in (not only the drug, but the actual meow group of 10) will influence treatment outcomes because the groups meet weekly. In essence, I have:   * 2 drugs (treatment)    * 16 meow groups   * 160 cases   * 15 time points I want to know whether people who received drug X reported less headache intensity over the 15 weeks. Basically, I want to know the headache severity slope for drug X, and whether that differs from drug Y. Please note that I am not particularly interested in the slopes in general, how quickly/when they change, or how the meow groups differed. I am mostly interested in comparing the efficacy of the drugs. I want to do mixed-modeling. I think treatment condition should be a fixed- factor, time points should be a covariate fixed-factor so that I can estimate a slope over time, and that time points need to be estimated as repeated- measures estimates. I also appreciate that I should run different models and test them and examine fit. Overall though, **what is the most appropriate way to test my hypothesis using a mixed-model design in SPSS (given what information you have)?** This is the SPSS syntax I have so far: ` MIXED pain BY drug WITH time /CRITERIA=.... /FIXED=drug time drug*time | SSTYPE(3) /METHOD=REML /RANDOM=INTERCEPT | SUBJECT(id) COVTYPE(ID) /REPEATED=time | SUBJECT(id) COVTYPE(AR1).` (meow group is missing from this syntax because I don't know where it is supposed to go). I have a general understanding of what a random and a fixed factor is, but I don't quite understand how these impact each other in practice (please note that I can't read most formulas, so plain English will be needed if you want to help on that front!). For example, **how would adding drug (or time) as a random effect influence the estimate of the slopes I am interested in?** Similarly:   * Should time points also be a random effect?   * How should I handle the "meow group" variable? Fix or random, or both?   * How does specifying a subject variable impact the random effect’s association with the dependent variable?   * Is there a minimum number of levels necessary for a variable to be a random effect? I've read that only having two levels does not make it possible.    * What are the "levels" in mixed-effects modelling referring to? Is student within class within school within school board within state within country a "6-level" model, or are the levels something else?  Lastly, how can I “nest” pain estimates within times, time within cases, cases within groups, and groups within treatments? Are pain ratings automatically nested within cases and time automatically nested within cases and cases within groups? I cannot "nest" group factors in the subject factors using SPSS it seems. Sorry for all the questions... I figured putting this all in one related post would be better than spamming the front page with 10 questions. I have read previous related answers on this website without finding answers to these questions (but answers to other questions of course!). Answers to any of these questions would be appreciated (even info on a single question would be great)!


title: Prove the following are independent

text: I have been unable to figure out the following question $f_{Y|X}(y|x)=N(x,x^2)$ $f_X(x)=U(0,1)$ Prove that $\frac{Y}{X}$ and $X$ are independent


title: Calculate $t$ value for expanded sample using just $t$ and $n$

text: Imagine a sample of observations $x_1 \dots x_n$ for which a $t$-value (i.e. $\bar x/SEM$) is known. Say an additional observation of zero is added to the sample so that we have $x_1 \dots x_{n+1}$ where $x_{n+1}=0$. Can a new $t$-value be calculated from just the known original $t$ and $n$? (I suspect that it can, but my mathematics is inadequate for me to work out exactly how.)


title: Correlations in count data

text: I have two count data variables X and Y that contain many zero values (90% in X, 60% in Y). I would like to check if a correlation exists between these variables, but I'm not sure how to proceed due to presence of zero values. Any suggestion? thanks in advance


title: Different Correlation/dependecy measures

text: I'm trying to figure what kind of correlation/depency measure to use on a multivariate data set I'm working on. I've tried the (pearson's)correlation, Kendall's Tau and Spearman's Rho. However, all 3 seem to give me similar correlation matrices. What does that imply about my data? Also, are there any other copula based dependency measures (better?) ?


title: Standardize/normalize power law distribution for machine learning

text: If my data follows a normal distribution I can standardize it for a machine learning algorithm, e.g. logistic regression, by subtracting the mean and dividing the result by the standard deviation. But what should I do with data that is power law/Pareto distributed, e.g. words in a corpus? How should data in this case be standardized? Or should I prefer simple scale normalizatiom?


title: Change detection in a Bernoulli sequence

text: Suppose $X_1,X_2,\ldots, X_n$ is a sequence of independent Bernoulli random variables with (say) $$X_i\sim \begin{cases} \mathrm{Ber}(0.1) & 1\leq i\leq an \\\ \mathrm{Ber}(0.2) & an+1\leq i\leq n\end{cases},$$ for some $a$ that is unknown to us satisfying, say $0.1\leq a\leq 0.9.$ How accurately can we detect the change point $a$ using knowledge of a realization of the sequence $X_1,X_2,\ldots, X_n$? Any references will be useful. Thanks.


title: In ANOVA, if the null hypothesis is true is the expected value of F guaranteed to be 1?

text: By null hypothesis is true I mean that the means of each group is the same at population level. I am wondering whether the expected value would depart from 1 in the case that the variances of the groups are different from one another at population level, or if the groups are non-normally distributed at population level, or in any other case. If it does, why does this change the expected value? I'm thinking of a One-way independent groups ANOVA with three groups, but wonder if the answer to this question would generalise beyond that.


title: How to reproduce Google Content Experiments' Probability of Outperforming Original?

text: I'd like to get a deeper understanding of how Google Content Experiments calculates the Probability of Outperforming Original in practice. The example below gives a case where the reported value disagrees with both the documented method (reproduced in R below) and intuition. The reason I'd like to do this is I am trying to ensure I understand the properties of experiments (e.g., false positive rate) run via Google Content Experiments. The available documentation is insufficient at the moment - for example, no A/A (not a typo) test results are provided. The most straightforward way to do so is via simulation, which requires reproducing the calculation of Probability of Outperforming Original. At the moment, I have an experiment with:   * Original: 1 session, 0 conversions.   * Variation 1: 2 sessions, 1 conversion, with a Probability of Outperforming Original of 50.7%. According to Google's page on how Content Experiments works, they are using Bayes' theorem to update the conversion rate estimate for each variant. The appendix implies that they are using a uniform beta prior (i.e., alpha = 1, beta = 1) and calculate the probability of each variant beating the original by simulating from the posterior distribution. In this case, that is (see http://en.wikipedia.org/wiki/Conjugate_prior#Example):   * Original: beta(alpha = 1, beta = 2)   * Variation 1: beta(alpha = 2, beta = 2) Simulating and calculating the Probability of Outperforming Original is straightforward in R:               > r0 = rbeta(100000,1,2)     > r1 = rbeta(100000,2,2)     > sum(r0 < r1)/100000     [1] 0.69889      So I'm getting ~70%, while Google Content Experiments gives ~50%. ~50% is also the Probability of Outperforming Original when you have no information, which makes sense. It doesn't make a lot of sense when we have information that Variation 1 is better, though.


title: How to bootstrap the best fit distribution to a sample?

text: If I have a sample:               set.seed(0)     x <- rlnorm(500)      Then I can use the fit.distr function to find the best fit among two candidate distributions, e.g.               library(MASS)     find.bestfit <- function(x){        logN <- fitdistr(x, "lognormal")        gam  <- fitdistr(x, "gamma")        ans <- ifelse(AIC(logN) < AIC(gam), "logN", "gam")        return(ans)     }          find.bestfit(x)     [1] "logN"      However, there is some probability that I will not recover the "true" distribution that was sampled (in this case "lognormal" was used to simulate `x`). How can I calculate this probability? I have only gotten so far as to consider using a bootstrap approach, but I am not familiar with this technique and am not sure exactly where to start:               ## create an empty vector     fit.samps <- rep(NA, 100)     ## determine fit to subsamples from original distribution     for(i in 1:100){       fit.samps[i] <- find.bestfit(sample(x, 10))     }      I suspect that my approach is wrong, because the sample size is arbitrary, and ultimately, the best fit distribution based on the `fitdistr` function should be selected most of the time. I would appreciate some pointers on how I might apply the bootstrap approach to answer this question.


title: Estimate a proportion using priority sampling (I just made that up)

text: I have this idea in my head that is either bunk or has a name I don't know. (I'm not naive enough to think I'm breaking new ground here!) Here's my scenario: I would like to know the proportion of a population that has disease x. Instead of taking a random sample of the population, I'm taking what I'm calling a priority sample, which is I start with the people who are most likely to have the disease, then after finding them I go to the next priority or the next most likely. After I test so many, I have a positivity proportion. The more I test, the further down the priority list I get, and the lower my positivity proportion. At some point I stop and then taking my positivity proportion and the number I've tested, I apply it to some distribution (skewed I would guess) and get an estimate of what the actual proportion of the population has the disease with some kind of confidence interval. Already with this I could imagine I could guess some kind of upper bound after some 10-15 positive tests since the true proportion in the population is unlikely to be higher than my positivity proportion since I started with the highest priority. Does this method exist? If so, can you steer me in the right direction to learn how to do it properly? **EDIT** After mulling it over a bit, I realize there's one more detail I haven't included here, which is that I know the size of the total population and therefore can tell you the testing coverage. In lieu of some fantastic mathematical proofs, I made a few simulations. First I take a population of 1000 and divide them up into 100 groups of 10. Each of these have some proportion of the disease normally distributed around the population mean which I want to estimate. (I'm setting it at 0.2 here.) ![enter image description here](http://i.stack.imgur.com/NBL0x.png) I then order these groups into the highest proportion (priority) to lowest. ![enter image description here](http://i.stack.imgur.com/BW49x.png) So if I were keeping track of my 10-person groups, when I reached the median group (50) I would have the population proportion. But the more common way is a cumulative collection. So if I graph the proportion cumulatively I get this. ![enter image description here](http://i.stack.imgur.com/wLJUe.png) This is surprisingly linear after about 20% coverage. And the true proportion seems to lie about 10% below the coverage at 20% for all but the extreme population proportions. So if this slope is consistent, then I could take the population coverage and extend it linearly to 100%.


title: Test to determine if one strategy has better results in a game than another?

text: I have an engine that runs a game, and bots that perform strategies in it. I want to do different tests of versions of a single bot, to see how variations influence success. For instance, if I have Bot A1 and Bot A2, I might run the following trials: ### Trial 0   * Bot A1   * Bot B   * Bot C    * Bot D ### Trial 1   * Bot A2   * Bot B   * Bot C    * Bot D Each trial consists of many games, so the end result is something like this:               Bot \ Place    1st  2nd  3rd  4th          A1             3    8    2    1     B              1    4    7    2     C              0    2    3    9     D              10   0    2    2      So, for the above trial, Bot A1 won the 3 games, was 2nd 8 times, came in last once, etc. **I'm looking for a way to compare two trials to see if the difference between them is statistically significant.** What method is best for doing this? (It's been a few years since I've done stats, so please forgive me if I'm missing information in this question. Also, if you think I'm completely off- base and there's another approach to determine which strategies are most effective, please let me know.) **Update** : I will arbitrarily decide that the way to rank bots is by looking at who comes in most in each place. (First, take whichever bot came in first the most times, then take whichever bot came in second the most times, etc.) For the above example, that would look like `D, A1, B, C`. Also, I'm using F#, for what it's worth. Here is a chart generated from a single trial (actual run-through of the program; not related to the sample data above). The bots are named alphabetically by how well they did in the ranking I described in the previous paragraph (so `A` was better than `B`, etc). ![chart showing the ranking of bots is A, B, C, D, E, F](http://i.stack.imgur.com/wjBop.png)


title: Are $\sigma_{70\%} = cte * \sigma /\sqrt{N}$, & $\sigma_{\bar x_{70\%}} = cte * \sigma /\sqrt{N}$ valid with non-normal distributions?

text: My working population is non normally distributed and from time to time I need to take a sample of it. Sometimes my sample exists of 100 units, at other times of 150 units. Of my sample I calculate the mean, $\mu$, and with bootstrapping I determine the SE, $\sigma_{\bar x}$. Since $\sigma_{\bar x} = \sigma /\sqrt{N}$, I calculate $\sigma$. So far so good, by using control charts ($\mu,\sigma$) I can now see if there is anything happening to worry about. By use of deviation calculation rules you can prove quite easily that the formula $\sigma_{\bar x} = \sigma /\sqrt{N}$ also can be used for non normal distributions. I play with the idea to make control charts using the trimmed mean and the 70%-value. Now I am not so sure anymore. (Frankly, I have my doubts whether you can do this for normal distributions as well.) Can I use the same simple trick to remove the $n$ from my equation? Let me put it like this: Is what I am writing valid? $$\sigma_{70\%} = cte * \sigma / \sqrt{N}$$ with a $cte$ independent of the value of N? $$\sigma_{\bar x_{70\%}} = cte * \sigma /\sqrt{N}$$ with a $cte$ independent of the value of $N$? To use my control charts I don’t need to know the value of my $cte$. For my s-chart I can just plot out $cte* \sigma$ in function of my sample number. But I really need to get rid of the $N$! I hope you understand what I am getting at. I would already be happy with a yes or a no, you can not do that anymore. If you can direct me to a place where I can find more information about this problem I would be even more grateful.


title: Random effect on scale parameter

text: I am trying to use random effects on both the location and scale parameters in a lognormal regression. As you can see in the JAGS code below, I am essentially fitting a accelerated failure time model with interval censored log-normal failure times. I am having problems coming up with decent non- informative priors for the random effects of the scale parameter (for which I use disp in the code). The way I have it set up now is really sensitive to priors specification. Anyone have any suggestions on how to proceed? With no random effect on the scale parameter and using a $U(0,big)$ on the inverse of the srqt of the scale, the model seems to converge just fine with estimates that make sense. Eventually I want to add fixed effects to this regression but am simply trying to fit this simpler model at the moment.               model{           for (i in 1:study){               for (j in offset[i] : (offset[i+1] - 1) ){                   y[j] ~ dinterval(t[j], lim.s[j,1:2])                   t[j] ~ dlnorm(mu[i],disp[i])                }             }           for (j in 1:study){            beta1[j] ~ dnorm(0,taup)            mu[j] <- beta0 + beta1[j]            alpha1[j] ~ dnorm(0,taup.d)            inv.disp[j] ~ dlnorm(log.inv.disp.hat[j],disp.inv.disp)            log.inv.disp.hat[j] <- alpha0 + alpha1[j]            disp[i] <- 1/inv.disp[j]^2          }              ## random effects         alpha0 ~ dnorm(0,0.001)         beta0 ~ dnorm(0,0.001)              ## hyperprior on  sqrt(inverse dispersion) for prior on dispersion         inv.disp.inv.disp ~ dunif(0.01,100)         disp.inv.disp <- 1/inv.disp.inv.disp^2              ## hyper prior for beta1         sigmap ~ dunif(0.01,100)          taup <- 1/(sigmap*sigmap)              ## hyper prior for alpha1          sigmap.d ~ dunif(0.01,100)         taup.d <- 1/(sigmap.d*sigmap.d)     }


title: Sparsity by discarding least squares' coefficient

text: Suppose I wish to regress $Y$ against a normalized $X$, but I would like a sparse solution. After regression, why is discarding the coefficients with smallest magnitude not allowed? For the record, I have heard of, and often use, LARS and LASSO methods. I'm just curious why the above approach is not applicable.


title: What test - close to zero

text: I have data grouped into category A and category B. My hypothesis: the data in category A is closer to zero than category B. What statistical tests can be used? Thanks


title: How do you call a situation or a point at wich statistics data stops changing?

text: When I calculate a prediction, for instance I am trying to find out who is going to win elections and I do that by asking people who they voted on. After a certain number of answers my data will stop changing. The chance percentage of being elected will stop changing. Is there a special analytical or statistical definition of such situation where no more data is necessary to give proper results? Thanks you.


title: How to learn Data Analysis?

text: I completed the ml-class offered by Prof Andrew Ng last fall I started trying one of the problems in Kaggle When I fit a logistic Regression, the algorithm gave me poor results. I saw that people do lot of data analysis before applying any machine learning algorithm I completed a book in Statistics and want to learn about Data Analysis Please help me identifying resources/courses that I can take to learn Data Analysis