title: When estimating variance, why do unbiased estimators divide by n-1 yet maximum likelihood estimates divide by n?

text: I am totally confused: On the one hand you can read all kinds of explanations why you have to divide by n-1 to get an unbiased estimator for the (unknown) population variance (degrees of freedom, not defined for sample size 1 etc.) - see e.g. here or here. On the other hand when it comes to variance estimation of a supposed normal distribution all of this doesn't seem to be true anymore. There it is said that the maximum likelihood estimator for variance includes only a division by n - see e.g. here. Now, can anyone please enlighten me why it is true here but not there? I mean normality is what most models boil down to (not least due to the CLT). So is the choice "division by n" yet the best choice for finding the best estimation for the true population variance after all?!? Thank you!


title: Panel Data: In a fixed effects model, does auto-correlation introduce bias?

text: Given a panel of countries over time, a fixed effects estimator makes sense to control for country-specific effects. My intuition tells me that if the dependent variable is correlated with lags of the independent variables, then bias will be introduced into the estimator. However, I'm having difficulty rigorously understanding why this would be the case. Additionally, is there any easy way to tell if the bias would be towards or away from zero?


title: Making a heatmap with a precomputed distance matrix and data matrix in R

text: I have made a heatmap based upon a regular data matrix in R, the package I use is `pheatmap`. Regular clustering of my samples is performed by the `distfun` function within the package. Now I want to attach a precomputed distance matrix (generated by Unifrac) to my previously generated matrix/heatmap. Is this possible?


title: Is it problematic if one predictor in a set accounts for almost all the prediction?

text: I am running a logistic regression with customer event data with multiple predictors. However, one variable is extremely important, alone predicting 60% of the customers for the event. When this main predictor is included in the model, other predictors add very little to prediction over and above this main predictor. This main predictor is not a post event variable. The variable has full business support to be in the model.   * Given this, is it still okay to retain this main predictor variable in the model?   * Does this suggest that there is anything wrong with the model?


title: Likert Item, Regression analysis

text: HI there: I have a one likert item I would like to analyze. I know that most scales require multiple likert items, but this is just an exploratory exercise of an existing survey. It's out of my control. When I treat the likert categories (strongly support - support - oppose -strongly oppose) as a numeric variable (coding 0 to 1) and conduct a multiple regression, I get some predictors in my models as statistically s ignificant. When I try to create some predicted values based on different configurations of my indendent variables, however, some combinations of values create predicted values larger than 1, which is impossible given the original scale. Incidentally, when I recode the four categories to two (support - oppose) and rerun it with a binomial logistic regression, the same predictors do not appear as statistically significant. I have three questions: in the first scenario, is it meaningful and / or appropriate to rescale the predicted values so that the maximum predicted value is 1 and the minimum is 0? Or is that not appropriate? Second question, can someone explain why the predictors would not be showing up as statistically significant when turning to logistic regression? Is there important variance at the poles of the DV (between strongly support and support) that is being lost when squishing them into binomial? Third question, should I just analyze this as ordinal? Thank you very much for your time. I am conducting this on R 3.0.2. Simon


title: Importing tables from MS access into R

text: I am a complete beginner at using R but have recently been running some code that a colleague gave me to connect R to a Microsoft Access database and pull out specific records. However a few days ago it stopped running. The error appears when trying to connect to the database (using the RODBC package) and import the table. It says 'Error in as.POSIXlt.character(x, tz, ...)'. From googling I can see that this means something such as a date etc. is in an incorrect format. However nothing on the database has been changed. So I'm completely baffled as to why it has suddenly stopped running! Does anyone have any ideas or possible solutions? Thank you!


title: how to impute missing values on numpy array created by train_test_split from pandas.DataFrame?

text: I'm working on the dataset with lots of NA values with sklearn and pandas.DataFrame. I implemented different imputation strategies for different columns of the dataFrame based column names. For example NAs predictor 'var1' I impute with 0's and for 'var2' with mean. When I try to cross validate my model using train_test_split it returns me a nparray which does not have column names. How can I impute missing values in this nparray? P.S. I do not impute missing values in the original data set before splitting on purpose so I keep test and validation sets separately.


title: Visualising changes in many-to-one (directional) relationships (votes) between a set of objects over time

text: I have a set of objects, each of which can be assigned to another object within the set in a many-to-one, directional assignment, like a vote. Objects cannot be assigned to themselves reflexively, but they can be in an unassigned state ("not voting"). So for example, the set `{A, ..., Z}` plus relation could be in the following state:               A -> E     E -> C     B -> C     C -> B     D -> Y     {F, ..., Z} are not voting.      This state will change over time, in discrete steps, one vote at a time. For example, at time `t = 0` the set could be as above. Then at `t = 1`, `Z` votes for `C` (`Z -> C`), resulting in the state:               A -> E     E -> C     B -> C     C -> B     D -> Y     Z -> C     {F, ..., Y} are unassigned.      Can anyone come up with a neat way to illustrate graphically this set and the state changes over time? Ideally, it should be evident at a glance   * How many votes there are for one object at a given time.   * Which objects are voting for another given object. Thanks!


title: Computing the difference between means and resolving a new/composite sample size

text: When computing the difference between means I can resolve a new mean, variance, standard error of the mean, and margin of error, but is there a way to compute a new/composite sample size value? I ask because I use the difference between means as a way to calibrate the mean and later when comparing other calibrated means sample size comes into the formulas. (http://stattrek.com/AP-Statistics-4/Difference-Means.aspx) ### Update In the code example below I subtract the calibrations mean from the benchmarks mean time to execute a test. The calibration is an empty test to clock the cost of the benchmarks overhead. The only part I am missing is how to create a composite sample size (if that is possible).                      mean -= cs.mean;            variance += cs.variance;            sd = sqrt(variance);            moe = sd * getCriticalValue(getDegreesOfFreedom(me, cal));            rme = (moe / mean) * 100;            extend(me.stats, {              'ME': moe,              'RME': rme,              'deviation': sd,              'mean': mean,              'size': ???,              'variance': variance            });


title: How can I treat blocks in a split plot design?

text: In a field experiment involving crops, what is the difference in considering block as random or otherwise as fixed factor? As far as I understood, random means that conclusion can be extended to other levels not included in the study; fixed factor on the contrary restricts the analysis only to a specific experiment.


title: Postscript in R: How to remove a large left margin (on y axis)

text: I have a problem. I am creating a wide barplot (50x grouped 4 columns) and printing it to ps. However the figure does not fit to printed ps - the x axis goes out of a page. It would fit if only I could remove a large left margin on yaxis. But I have no idea how to move picture left and reduce the margin. I have spend some time on asking uncle google but I have not found anything. Could somebody help me? Here is my R script:               mx <- matrix( c(5,3,8,9,5,3,8,9),nr=4)          postscript(file="<fileName>.ps");          barplot(mx, beside=T,  col=gray.colors(4), cex.axis = 1.4, cex.names=1.2, xlim=c(1,60), width=0.264,  xaxs = "r", yaxs = "r")          legend("topleft", c( "A1",  "A2","B1","B2"), pch=15,              cex=1.1,            col=gray.colors(4),             bty="n")


title: Integral of a function with uniform kernel

text: I am trying to understand question 9-1 on p.334 in Cameron & Trivedi (link) where I have to calculate the bias of a kernel density estimate at $x=1$ and $n=100$, where we assume that the underlying density is standard normally distributed, $N(0,1)$. The kernel is uniform, that is we have $$K(v)=0.5 \times 1(|v|<1)$$ where $v$ is defined as $(x_i-x)/h$ and $\mathbf1$ is the indicator function, that is either $1$ if its argument is true and $0$ otherwise; $h$ is the bandwidth, given to be $1.0$ in this example. The bias of the kernel is then "simply" obtained as $$\text{bias} = 0.5h^2f''(x)\int v^2 K(v)dv$$ This doesn't seem to difficult at first. I know $h$, I can calculate $f''(1)$ from the normal distribution analytically. However, I am stuck on finding $\int v^2 K(v)dv$. Is there a standard trick, a table or anything that can help me calculating it? I could do it numerically but I believe it should be possible to do it analytically too. ps this is not a homework question, just curiosity how to apply Kernel density estimates without relying on numerical computations.


title: Best way to specify a mixed ANCOVA in R?

text: After using `ezANOVA` as my primary way of specifying mixed ANOVAs, I've hit a stumbling block when it come to adding a covariate to the model. I am using an ANCOVA in order to determine if there is a developmental trajectory in my data; namely, I need to be able to see the F-statistic and p-values for interactions with the covariate (see p.466 onwards here if you want an example). Using `ezANOVA`, I can include covariates but the output does not show the F-statistic and p-values for interactions with the covariate - the main effect of the covariate is also not tested using this method. My `ezANOVA` model is as follows:               aov.model<-ezANOVA(           data=textureView.child.outliersRemoved           , dv=.(x)           , wid=.(ID)           , within=.(Texture,View)           , between=.(TNOGroup)           , between_covariates=.(Age)           , type=3           , return_aov=TRUE           )      Another option is to use `lm` or `Anova`, but I don't know how to specify the error terms properly for either and I'm limited because I want to use Type-III sums of squares (`drop1` doesn't work in the cases where I've tried to use the `aov` wrapper for `lm`; it fails while reporting `'Error in formula.default(object, env = baseenv()) : invalid formula'`). Finally, I've heard about using the `nlme` package to specify my ANCOVA as a mixed model instead, but I don't know where to begin here (despite spending a while reading about it). To give a summary, I'm trying to do a 2 (between; TNOGroup) x 2 (within, Texture) x2 (within, View) mixed ANCOVA, with age as a covariate. I want to use Type-III sums of squares, and see the F-statistic and p-values for interactions with the covariate, as well as for the main effect of the covariate. Any advice on the best way to do this would be much appreciated.


title: Sample size needed for Fleiss' Kappa?

text: A group of raters (about 20) will be watching a series of videos and will be classifying them into 4 categories. I will be running a Fleiss' kappa to measure the agreement. How does one compute for the sample size to arrive at 0.8 power, 0.05 alpha? Also, will that sample size be the number of videos to be evaluated?


title: How do I weight population statistics in SPSS?

text: I'm struggling with analysing my survey results in SPSS and really hoping you can help. I surveyed occupational therapists around aspects of their practice, on a national level and the seniority of my sample were made up of the following: Band 5's- 8 Band 6's- 25 Band 7's-35 Band 8's- 6 This distribution compares nationally within the NHS to: Band 5- 24% (3,214) Band 6- 45% (6130) Band 7- 25% (3336) Band 8 -6% (838) So, how do I translate this into SPSS to ensure that my responses are weighted according to seniority, ensuring that my findings have generalizability? I am currently analysing my categorical data via chi square tests. Thanks so much for any help you can offer-I've trawled through the net and nothing seems to offer straight forward advice! Amy


title: Distribution of an euclidean norm

text: If $X = \begin{pmatrix} X_{1} \\\ X_{2} \end{pmatrix}$ (with $X_{1},X_{2} \in \mathbb{R}^{n}$) follows a gaussian distribution $\mathcal{N}(\mu,\Sigma)$. What is the distribution of $\Vert X_{1} \Vert^{2}$ (where $\Vert \cdot \Vert$ denotes the euclidean norm on $\mathbb{R}^{n}$). If $\Pi=\left( \begin{array}{c|c} I_{n} &\ 0 \\\ \hline 0 & 0 \end{array}\right)$, I know that $X'=\Pi X = \begin{pmatrix} X_{1} \\\ 0 \end{pmatrix} \sim \mathcal{N}(\Pi \mu, \Pi \Sigma {}^t \Pi)$ and $\Vert X_{1} \Vert^{2} = \Vert \Pi X \Vert^{2}$. I had the following idea : if $\mu' = \Pi \mu$ and $\Sigma' = \Pi \Sigma {}^t \Pi$, let $\Gamma {}^t \Gamma = \Sigma'$ be the Cholesky decomposition of $\Sigma'$. Then, $\Gamma^{-1} X' \sim \mathcal{N}(\Gamma^{-1}\mu',I_{n})$. So, I know that $\Vert \Gamma^{-1} X'\Vert^{2}$ follows a $\chi^{2}(n,\Vert \Gamma^{-1} \mu' \Vert^{2})$ but I don't know if that helps.


title: When is the coefficient of variation considered good

text: I collect some data out of a survey. And now I want to analyse the answers. The question consists of needed features in a Tool. These features can be rated from 1 to 5. Table of the features an their rating:               F.1 F.2 F.3 F.4 F.5 F.6     5   4   5   5   5   4     4   4   5   5   3   5     5   3   4   4   2   4     5   5   5   5   5   5     4   4   4   3   4   3     2   3   5   5   4   4     3   5   5   5   4   4     5   5   5   4   4   4     3   2   4   5   4   3     ...      with a corresponding mean of:               4   3,8888888889    4,6666666667    4,5555555556    3,8888888889    4      with a Standard deviation of:               1,1180339887    1,0540925534    0,5 0,7264831573    0,9279607271    0,7071067812      and with a coefficient of variation of:               27,95084971 27,10523708 10,71428571 15,94719125 23,86184726 17,67766952      What can I conclude out of these data ? I mean, Is there a known value (alpha) for a good coefficient of variation ? So that I could say : _The coefficient of variation for feature 1 is 27,9% > alpha, therefore ..._ I have about 58 datasets (Answers) like in the table listed above. But I don't how could I conclude something meaningful out of it. I calculated the mean, the standard deviation, the coefficient of variation and the standard error. What Information can I get out of these values ? Can I say something like: _For feature 1, I have a coefficient of variation of 27,9% this means that ... (but/and also)a standard error of ... this means ... therefore ..._ I m kinda stuck here, I would be thankful for an suggestion!


title: Bonferroni Multiple comparisons -- testing for any true

text: I am doing multiple comparisons with several thousand independent alternatives. None are statistically significant when Bonferroni corrected; however, many more are significant when uncorrected than expected by chance, so I compute this probability to show that most of my low p-valued tests are true. Can anyone point me to a citable paper (for a referee) on this approach. I found it discussed (with citation missing!) in wikipedia: http://en.wikipedia.org/wiki/Multiple_comparisons#Assessing_whether_any_alternative_hypotheses_are_true thanks!


title: Cor.test and statistical test values for spearman correlation, how to make sense of them?

text: I am trying to make sense of the spearman test I am getting as a result of using cor.test on R package. This page seemed to me to inform that I should care not only for p-values and rho values from the output but as well its relation with the t tests. Sometimes I get as result p value of 0 and rho value of 1 and the value of the statistical test is zero as well. I saw a similar question here but seems too focused on some performed spearman correlation. Since this is not really R related but involves conceptual understanding of spearman test statistic as mentioned on R, I thought this would be more appropriated than stack overflow. Thank you.


title: Comparing the variability of two different large data sets

text: I have two data sets, one data set contains value in the range of 10,000- 1,000,000 and other data set contains the value between 1,000-10,000. I want to compare the variability means,which data has more variability. Do I need to first standardized / normalized the data and then compare the variance, but standardization will give me mean 0 and variance 1. Also, which is better Coefficient of variation or Interquartile Range Thanks in advance


title: post-hoc multiple comparisons with chi squared using spss v. 22

text: I am carrying out an analysis using a 4$\times$2 crosstab. I found an overall significant difference but I would like to find if there are significant differences among the 4 groups. Is there a way to carry out these multiple comparisons?


title: How to assess prevalence of grossly inaccurate data in recorded body weights in medical records?

text: I have around 5 million observations of weights of patients over the course of a year. Some patients have one measurement during that year and others have 30-40 (or more). I've noticed a great deal of variability. For example, on one day a patient will have a recorded weight of 120 lbs and the next a recorded weight of 220 pounds - obviously, one of those is incorrect. In this data set, we also have the locations (multiple) that are entering the weights, so it would be helpful to know which of the locations are entering the inaccurate data. Eventually, we'd like to use this weight data in various studies but would first like to get an idea of how often these bogus weights appear in this collection. Any suggestions?


title: Averages of averages (of averages, of averages...)

text: Consider the following cell biology experiment. We are comparing $T$ different treatments of cultured cells. Each treatment $t$ is replicated in several _(microtiter) wells_ , indexed by the variable $w \in \\{1, 2, \cdots, W\\}$. To measure the response to the treatment in well $w$, a total of $F_w$ non- overlapping micrographs, or _fields_ , are recorded. Then, for each field $f$ in well $w$, a total of $C_{wf}$ cells are identified computationally, whereby each cell $c$ (in well $w$, field $f$) is represented by a set of $P_{wfc}$ pixels. Finally, associated with each pixel $p$ is a measurement $x_{wfcp}$ (derived from the intensities of various fluorescence signals recorded at that pixel). > The problem is to aggregate all the pixel measurements $x_{wfcp}$ to produce > a "reasonable measure" $X_t$ of the effect of the treatment $t$ on the cells > treated with it, as well as some measure of the "spread" of $X_t$. The standard approach to such problems is to use the mean as "the measure" and the variance (or standard deviation) as "the spread". In this case, however, there are multiple, non-equivalent ways in which means and variances could be computed. Focusing for now on the means, at one extreme, one could simply add up the $x_{wfcp}$ over all the pixels (disregarding their distribution over cells, fields, and wells), and divide this sum by the total number of pixels $P$ (for treatment $t$): $$ \frac{1}{P}\sum_{w=1}^W\sum_{f=1}^{F_w}\sum_{c=1}^{C_{wf}}\sum_{p=1}^{P_{wfc}} x_{wfcp} $$ At the opposite extreme, we could average at each level: first compute the average $x_{wfc}$ of $x_{wfcp}$ for each cell, then compute average $x_{wf}$ of the $x_{wfc}$ for each field, and so on: $$ \frac{1}{W}\sum_{w=1}^W \left[\frac{1}{F_w} \sum_{f=1}^{F_w} \left[\frac{1}{C_{wf}}\sum_{c=1}^{C_{wf}} \left[\frac{1}{P_{wfc}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right]\right]\right] $$ In general, these two expressions will not be equal. Plus there are several variations in-between. By my count, there are 8 ways of doing this (including the two above); I've listed all in their full glory at the end of this post. For example, one could compute this (number 6 in the list below): $$ \frac{1}{W}\sum_{w=1}^W \left[\frac{1}{C_w} \sum_{f=1}^{F_w} \sum_{c=1}^{C_{wf}} \left[\frac{1}{P_{wfc}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right]\right] $$ ...where $C_w = \sum_f \sum_c \; 1$ is the total number of cells (summed over all fields of) well $w$. (The recipe encoded by this expression says: "compute the average value of $x_{wfcp}$ for each cell, namely $x_{wfc} = \left[\sum_p x_{wfcp}\right]/P_{wfcp}$; then, for each well $w$, compute the average of these $x_{wfc}$ averages over all the $C_w$ cells in well $w$—disregarding their distribution over fields—, namely $x_w = \left[ \sum_f \sum_c x_{wfc}\right]/C_w$; and, finally, average the $x_w$ over all $W$ wells, $\left[\sum_w x_w\right]/W$.") Faced with all these different ways to "use averages" to measure the effect of treatment $t$, the immediate question is, of course, which one to choose? A sharper version of the question would be: how can I determine under which scenarios a given variant would be appropriate/informative/useful? And, more generally: are there any pitfalls in computing averages of averages (of averages... )? Thanks! * * * (corrections welcome) $$ \small \begin{array}{lrl} 1\. \;\; \frac{1}{P}\sum_{w=1}^W\sum_{f=1}^{F_w}\sum_{c=1}^{C_{wf}}\sum_{p=1}^{P_{wfc}} x_{wfcp} && && \\\ 2\. \;\; \frac{1}{W}\sum_{w=1}^W \left[\frac{1}{P_w} \sum_{f=1}^{F_w} \sum_{c=1}^{C_{wf}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right] && \mathrm{where} && P_w = \sum_{f=1}^{F_w}\sum_{c=1}^{C_{wf}}\sum_{p=1}^{P_{wfc}} \; 1 \\\ 3\. \;\; \frac{1}{F}\sum_{w=1}^W \sum_{f=1}^{F_w} \left[\frac{1}{P_{wf}}\sum_{c=1}^{C_{wf}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right] && \mathrm{where} && F = \sum_{w=1}^W \sum_{f=1}^{F_w} \; 1 \, , \, P_{wf} = \sum_{c=1}^{C_{wf}}\sum_{p=1}^{P_{wfc}} \; 1 \\\ 4\. \;\; \frac{1}{C}\sum_{w=1}^W \sum_{f=1}^{F_w} \sum_{c=1}^{C_{wf}} \left[\frac{1}{P_{wfc}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right] && \mathrm{where} && C = \sum_{w=1}^W \sum_{f=1}^{F_w} \sum_{c=1}^{C_{wf}} \; 1 \\\ 5\. \;\; \frac{1}{W}\sum_{w=1}^W \left[\frac{1}{F_w} \sum_{f=1}^{F_w} \left[\frac{1}{P_{wf}}\sum_{c=1}^{C_{wf}}\sum_{p=1}^{P_{wfc}} x_{wfcp}\right]\right] && && \\\ 6\. \;\; \frac{1}{W}\sum_{w=1}^W \left[\frac{1}{C_w} \sum_{f=1}^{F_w} \sum_{c=1}^{C_{wf}} \left[\frac{1}{P_{wfc}} \sum_{p=1}^{P_{wfc}} x_{wfcp} \right]\right] && \mathrm{where} && C_w = \sum_{f=1}^{F_w} \sum_{c=1}^{C_{wf}} \; 1 \\\ 7\. \;\; \frac{1}{F}\sum_{w=1}^W \sum_{f=1}^{F_w} \left[\frac{1}{C_{wf}}\sum_{c=1}^{C_{wf}} \left[\frac{1}{P_{wfc}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right]\right] && && \\\ 8\. \;\; \frac{1}{W}\sum_{w=1}^W \left[\frac{1}{F_w} \sum_{f=1}^{F_w} \left[\frac{1}{C_{wf}}\sum_{c=1}^{C_{wf}} \left[\frac{1}{P_{wfc}} \sum_{p=1}^{P_{wfc}} x_{wfcp}\right]\right]\right] && && \hspace{3in} \end{array} $$


title: Proper test for time series

text: I have monthly precipitation for one site over a 35 year period. I would like to test if the total precipitation for a specific year is equal to the mean total precipitation for the rest of the period of record. In addition, I would also like to split the data into wet season [Nov-April] vs dry season [May- Oct]. I realize there are likely issues to take into account regarding auto- correlation, but I am yet analyze these types of data personally. Would someone please be kind enough to point me in the right direction? Thanks ahead of time. EDIT: IrishStat, The data is below, I use R, so below was created using dput(). I could figure this out, but it would take me a while and I would feel more confident if I were shown once.               PPT <- structure(list(YEAR = 1965:2010, JAN = c(3.15, 4.81, 4.94, 2.49, 4.26, 4.31, 3.96, 5.06, 3, 3.58, 6.27, 1.79, 6.14, 14.43, 5.04, 6.22, 0.94, 3.23, 6.6, 5.74, 5.08, 3.8, 6.35, 4.38, 0.87, 8.66, 15.48, 13.08, 7.52, 5.02, 5.45, 5.23, 5.88, 16.33, 5.18, 2.95, 3.33, 4.91, 2.49, 5.75, 2.08, 2.62, 4.38, 7.02, 2.16, 6.73), FEB = c(5.95, 9.37, 6.09, 2.84, 3.87, 4.27, 6.72, 3.78, 7.09, 6.56, 3.97, 2.39, 1.87, 4.48, 7.73, 1.46, 11.59, 9.58, 11.61, 5.95, 7.03, 4.66, 8.96, 6.4, 1.64, 10.3, 2.4, 8.06, 4.29, 1.88, 2.84, 6.36, 6.52, 6.19, 0.82, 0.87, 3.56, 2.65, 7.03, 8.12, 5.6, 4.22, 2.07, 4.88, 4.86, 5.43), MAR = c(4.89, 2.58, 0.69, 1.98, 8.94, 7.17, 4.9, 7.12, 8.23, 6.88, 6.83, 9.95, 4.73, 5.79, 9.29, 11.6, 2.59, 7.84, 9.38, 5.58, 5.06, 3.91, 5.43, 6.36, 4.32, 10.41, 4.89, 2.9, 7.31, 7.83, 10.03, 10.2, 3.75, 6.11, 5.79, 4.12, 8.73, 4.18, 5.2, 0.68, 4.28, 0.39, 0.54, 4.32, 14.43, 5.16), APR = c(0.61, 2.45, 1.68, 2.25, 4.55, 1.94, 0.48, 2.29, 9.22, 5.08, 6.84, 1.96, 1.83, 4.76, 5.77, 13.53, 1.24, 2.9, 12.97, 3.19, 1.01, 3.33, 0.63, 3.71, 2.88, 2.47, 9.1, 2.72, 2.55, 4.74, 6.71, 11.68, 6.32, 4.54, 0.11, 1.22, 0.3, 3.12, 3.06, 2.32, 20.5, 6.09, 3.43, 5.49, 2.07, 2.04), MAY = c(0.42, 5.66, 2.86, 2.7, 9.89, 2.86, 4.53, 7.8, 2.63, 3.4, 11.11, 10.7, 2.8, 10.71, 5.5, 9.61, 11.88, 1.64, 1.54, 2.49, 7.85, 5.35, 12.35, 0.31, 7, 4.94, 13.75, 2.4, 5.81, 2.89, 6.04, 0.48, 7.96, 0.81, 3.15, 0.65, 0.55, 2.34, 5.84, 2.03, 7.08, 3.44, 1.86, 9.26, 7.3, 7.03), JUN = c(15.43, 2.34, 4.33, 3.79, 4.27, 9.77, 2.82, 4.42, 6.91, 3.01, 3.52, 4.19, 0.71, 12.32, 3.37, 3.81, 2.51, 7.53, 6.6, 2.8, 8.08, 5.54, 15.58, 5.91, 18.52, 6.22, 5.85, 4.62, 3.38, 7.32, 4.53, 7.42, 5.88, 2.23, 8.32, 4.21, 13.77, 3.4, 9.54, 10.8, 10.44, 1.33, 6.35, 3.32, 3.69, 5.04), JUL = c(13.72, 7.6, 3.69, 6.95, 10.88, 1.93, 6.88, 6.56, 7.67, 5.29, 10.03, 8.03, 5.04, 23.67, 8.4, 4.06, 5.05, 11.94, 0.82, 3.9, 8, 4.29, 5.16, 4.8, 8.92, 5.84, 8.64, 5.03, 8.08, 10.94, 8.65, 6.55, 28.58, 6.17, 9.71, 3.21, 13.69, 10.76, 18.44, 4.73, 11.43, 5.34, 7.09, 5.4, 5.18, 2.18), AUG = c(10.1, 8.51, 7.88, 3.57, 12.91, 11.84, 6.4, 4.14, 4.3, 6.75, 8.69, 2.08, 9.33, 6.24, 5.41, 1.31, 6.28, 6.03, 6.08, 14.23, 5.78, 5.01, 12.48, 13.19, 2.17, 2, 6.74, 6.2, 7.57, 6.49, 9.61, 5.99, 1.25, 5.9, 5.99, 3.02, 11.02, 5.83, 5.2, 8.28, 11.4, 7.13, 5.96, 14.14, 6.24, 10.33), SEP = c(10.74, 3.43, 6.37, 4.19, 5.52, 5.53, 11.07, 1.01, 10.13, 7.6, 12.88, 6.48, 7.53, 3.97, 9.45, 5.28, 1.11, 2.31, 7.02, 0.57, 6.01, 3.07, 2.4, 16.12, 4.49, 1.64, 3.13, 1.44, 4.99, 1.46, 3.05, 7.72, 1.29, 24.11, 2.27, 9.53, 4.66, 14.77, 3.66, 12.63, 4.7, 5.25, 6.64, 7.67, 8.16, 6.33), OCT = c(2.57, 3.84, 7, 0.38, 1.75, 7.1, 0, 2.77, 1.28, 0.56, 5.52, 6.81, 4.13, 0, 0.69, 3.77, 1.9, 2.22, 5.53, 1.78, 13.08, 5.89, 0.23, 1.97, 2.28, 2.88, 2.28, 3.85, 5.3, 5.88, 12.13, 2.93, 4.88, 1.72, 4.42, 0.39, 3.43, 8.44, 2.68, 2.44, 0.25, 3.76, 8.78, 4.14, 7.67, 1.68), NOV = c(0.89, 4.88, 0.58, 3.83, 0.8, 2.69, 2.29, 4.64, 3.59, 3.05, 8.83, 5.57, 6.73, 3.76, 7.55, 5.43, 0.81, 4.78, 6.84, 2.26, 2.98, 6.06, 5.07, 4.44, 7.76, 1.84, 3.8, 11.7, 1.8, 3.71, 11.63, 2.13, 12.03, 6.27, 1.02, 7.34, 1.78, 5.69, 3.84, 11.64, 3.22, 3.7, 3.91, 3.47, 3.82, 4.73), DEC = c(3.87, 4.92, 5.77, 6.2, 8.27, 5.21, 5.93, 7.51, 7.07, 3.58, 3.63, 4.95, 4.76, 3.93, 1.96, 1.67, 5.88, 8.04, 5.5, 1.7, 3.29, 4.62, 5.3, 2.47, 5.21, 2.37, 2.87, 5.28, 4.73, 2.05, 4.84, 6.14, 2.52, 5.35, 5.4, 3.49, 2.74, 4.87, 4.13, 5.58, 3.9, 4.51, 5.77, 3.34, 13.52, 2.72)), .Names = c("YEAR", "JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC"), class = "data.frame", row.names = c(NA, -46L))


title: How are we defining 'reproducible research'?

text: This has come up in a few questions now, and I've been wondering about something. Has the field as a whole moved toward "reproducibility" focusing on the availability of the original data, and the code in question? I was always taught that the core of reproducibility was not necessarily, as I've referred to it, the ability to click Run and get the same results. The data-and-code approach seems to assume that the data are correct - that there isn't a flaw in the collection of the data itself (often demonstrably false in the case of scientific fraud). It also focuses on a single sample of the target population, rather than the replicability of the finding over multiple independent samples. Why is the emphasis then on being able to re-run analysis, rather than duplicate the study from the ground up? The article mentioned in the comments below is available here.


title: How to interpret this output in Dimension Reduction?

text: I'm running Linear Discriminant Analysis on a dataset and then performing clustering on it. I'm reducing it to dimensions 2,6,10. On comparing metrics like Accuracy and Normalized Mutual Information, the values at dimension 6 reduce by almost half. While at dimension 10 it is close to dimension 2. I'm not able to comprehend why the accuracy might go down and then up again. The point is why do we see for 2D and 10D a better answer and for 6D a worse answer. I mean u go from 2D to 6D, the accuracy goes worse and then goes back up again for 10D


title: Factor Analysis - principal axis analysis and MAP-Test

text: I have got an important question concerning principal axis factor analysis: Is it possible to extract the factors according to the MAP-criteria, or is this criteria only valid when using principal component analysis?


title: Correlation of handedness to sexual orientation?

text: There have been studies done that report that heterosexual individuals are somewhat more likely to be right-handed than homosexual individuals. Does it necessarily follow that left-handed individuals are more likely to be homosexual than right-handed individuals? I feel like this is a mathematical concept I can't place my finger on right now, and could be solved with a formal proof. Wikipedia article on this topic Edit: I think this is might be the method. Assume individuals can only be either heterosexual (group $A$) or homosexual (group $B$), and can be either right-handed (group $R$) or left-handed (group $L$). Given, $P(A,R) > P(B,R)$ $\because P(A,R)+P(A,L)=1, P(B,R)+P(B,L)=1$ Combine the two equations, $P(A,R)+P(A,L)=P(B,R)+P(B,L)$ move terms, $P(A,R)-P(B,R)=P(B,L)-P(A,L)>0$ $\therefore P(B,L)>P(A,L)$ Therefore, a person who is left-handed is more likely to be homosexual than (fixed) a right-handed person is. Is this correct, based on the (obviously) simplified assumptions?


title: Chi square goodness of fit test

text: The table below displays the number of accidents recorded at a particular intersection during each of the four seasons last year:                   Season           Spring Summer Fall Winter         no. of accidents 13      24     18  25      We would like to conduct a chi-square goodness-of-fit test to determine whether accidents are uniformly distributed over the four seasons.The value of the test statistic for the appropriate test of significance is? I'm not quite sure how to proceed with this one.. Binomial distribution would not work, I think, since it's asking for a uniform distribution. How would I compute p-hat?


title: Any suggestions for making R code use multiple processors?

text: I have R-scripts for reading large amounts of csv data from different files and then perform machine learning tasks such as svm for classification.   Are there any libraries for making use of multiple cores on the server for R.   or   What is most suitable way to achieve that?


title: How can I link item by relevance?

text: I am building a web application for used book trading and I am adding a feature to propose other book that would be interesting when they view an offer. Currently the data that I store are the following (they are updated each time someone visit an offer) : ISBN (book of the offer)   SessId (a unique id that everyone has when visiting the website)   NumberOfVisit (The number of time someone has view an offer of that book) I also have access to some user update data which categorize the book by subject and course. It isn't necessarily up-to-date and precise, but it's nonetheless data. What are the approach to list the most interesting books for a book ?


title: How to present box plot with an extreme outlier?

text: I could use some guidance about presenting some data. This first plot is a case-control comparison for the cytokine IL-10. I've manually set the y axis to include 99% of the data. ![IL-10 With Manual Y axis](http://i.stack.imgur.com/mm8Tk.jpg) The reason I set this manually is because the case group has an extreme outlier. ![With outlier](http://i.stack.imgur.com/IVDNe.jpg) My collaborators are hesitant to perform an outlier removal to our dataset. I'm OK with it, but they'd rather not. That'd be the obvious solution. But if I'm going to keep all the data and not remove this outlier, how can I present this boxplot optimally? Split axis? Is it acceptable to use just the first graph and note that it was constructed to include all the data? (This option feels dishonest to me). Any advice would be great.


title: Common words that have particular statistical meanings

text: I am not a statistician but my research work involves statistics (analyzing data, reading literature, etc.). I was again reminded from a comment on one of my questions posted here that there are some common words that have particularly specific meanings or connotations for those who are well- practiced in the field of statistics. It will be helpful to have a list of such words and may be phrases along with some comments.


title: Is it necessary to run 10 replications of boosted regression model if using dismo in R?

text: I am trying to use the dismo package to build up a BRT model with presence- only data. I have randomly sampled 10 000 background data. My questions is: do I need to run 10 (or 5) replications of the brt model on the dataset like a conventional cv procedure? For example, for each replication, 60% of the presence data and 60% of the background data will be used for training the model, and the rest presence data for test. The average of the prediction maps will be the result. However, in the dismo package, the gbm.step method already uses cv to derive the optimal number of trees, so I am wondering if I can just put all the presence data and the background data to generate the optimal model? Thanks a lot!


title: Convert nonlinear mixed model to log-linear mixed model in R (nlme)

text: I want to analyze fish aquaculture data comparing growth rates by state and controlling for tank effects as a mixed effect. I will use the equation $$ w = aL^b $$ but I want `b` to be a function of state of origin. I can do this using the `nlme` function in the `nlme` package. Here is a toy data set similar to mine.               tl <- seq(100, 300, 1/5)     state <- c(rep(x = c("DE", "TX", "FL", "VA", "SC"), times = length(tl)/5), "DE")          ai <- 0.000004     bi <- 3.3          lw <- NULL     set.seed(1234)     for(i in 1:length(tl)) {       lw[i] <- log(ai) + bi * log(tl[i]) + rnorm(n = 1, 0, 0.1)     }          w <- exp(lw)         df <- data.frame(cbind(w, tl))     df$state <- as.factor(state)         df$tank <- as.factor(c(rep(1:15, length(tl)/15), 1:11))      This is the code I would use:               nl <- nlme(w ~ a * tl ^ b, fixed = list(a ~ 1, b ~ state),                 random = b ~ 1|tank, data = df, start = c(a = rep(0.0000001, times = 1),                 b = rep(3.1, times = 5)))     summary(nl)      However, a plot of the data reveals that the error is likely to be multiplicative on the arithmetic scale (variability in `w` (weight) increases with increasing `tl` (total length)):                   plot(tl, exp(lw))      My current model using nlme is $$ w = aL^b + \epsilon \quad \epsilon \sim N(0, \sigma^2) $$ where _b_ is a function of state. However, I want to make the equivalent log- linear model: $$ w = aL^b\epsilon \quad \epsilon \sim N(0, \sigma^2) $$ $$ \log(w) = \log(a) + b \cdot \log(L) + \log(\epsilon) $$ If it weren't for the fact that I want _b_ to be a function of state I could use the `lme` function in `nlme`:               lm1 <- lme(log(mass) ~ log(tl), data = df, random = ~ 1|tank)     summary(lm1)      How would I make the equivalent log-linear model to my nonlinear model except for the error distribution? Xiao et al. 2010 discuss this issue but it's for non-mixed models and for functions when _a_ or _b_ are a function of a factor. http://www.esajournals.org/doi/pdf/10.1890/11-0538.1 The final part of the question would be whether I could use AIC to compare the two models or if I should just use the log-linear if the data clearly look to have that error structure?


title: Express answers in terms of original units, in Box-Cox transformed data

text: For some measurements, the results of an analysis are appropriately presented on the transformed scale. In most of the cases, however, it's desirable to present the results on the original scale of measurement (otherwise your work is more or less worthless). For example, in case of log-transformed data, a problem with interpretation on the original scale arises because the mean of the logged values is not the log of the mean. Taking the antilogarithm of the estimate of the mean on the log scale does not give an estimate of the mean on the original scale. If, however, the log-transformed data have symmetric distributions, the following relationships hold (since the log preserves ordering): $$\text{Mean}[\log (Y)] = \text{Median}[\log (Y)] = \log[\text{Median} (Y)]$$ (the antilogarithm of the mean of the log values is the median on the original scale of measurements). So I only can make inferences about the difference (or the ratio) of the medians on the original scale of measurement. Two-sample t-tests and confidence intervals are most reliable if the populations are roughly normal with approximately standard deviations, so we may be tempted to use the `Box-Cox` transformation for the normality assumption to hold (I also think that it is a variance stabilizing transformation too). However, if we apply t-tools to `Box-Cox` transformed data , we will get inferences about the difference in means of the transformed data. How can we interpret those on the original scale of measurement? (The mean of the transformed values is not the transformed mean). In other words, taking the inverse transform of the estimate of the mean, on the transformed scale, does not give an estimate of the mean on the original scale. Can I also make inferences only about the medians in this case? Is there a transformation that will allow me to go back to the means (on the original scale) ? _This question was initially posted as a commenthere_


title: How do I get the number of rows of a data.frame in R?

text: After reading a dataset:               dataset <- read.csv("forR.csv")        * How can I get R to give me the number of cases it contains?   * Also, will the returned value include of exclude cases omitted with `na.omit(dataset)`?


title: Multiple comparisons on a mixed effects model

text: I am trying to analyse some data using a mixed effect model. The data I collected represent the weight of some young animals of different genotype over time. I am using the approach proposed here: http://blog.gribblelab.org/2009/03/09/repeated-measures-anova-using-r/ In particular I'm using solution #2 So I have something like               require(nlme)     model <- lme(weight ~ time * Genotype, random = ~1|Animal/time,               data=weights)         av <- anova(model)      Now, I would like to have some multiple comparisons. Using `multcomp` I can do:               require(multcomp)     comp.geno <- glht(model, linfct=mcp(Genotype="Tukey"))     print(summary(comp.geno))      And, of course, I could do the same with time. I have two questions:   1. How do I use `mcp` to see the interaction between Time and Genotype?   2. When I run `glht` I get this warning: `covariate interactions found -- default contrast might be inappropriate` What does it mean? Can I safely ignore it? Or what should I do to avoid it? **EDIT:** I found this PDF that says: > Because it is impossible to determine the parameters of interest > automatically in this case, mcp() in multcomp will by default generate > comparisons for the main effects only, **ignoring covariates and > interactions**. Since version 1.1-2, one can specify to average over > interaction terms and covariates using arguments interaction_average = TRUE > and covariate_average = TRUE respectively, whereas versions older than 1.0-0 > automatically averaged over interaction terms. **We suggest to the users, > however, that they write out, manually, the set of contrasts they want.** > One should do this whenever there is doubt about what the default contrasts > measure, which typically happens in models with higher order interaction > terms. We refer to Hsu (1996), Chapter~7, and Searle (1971), Chapter~7.3, > for further discussions and examples on this issue. I do not have access to those books, but maybe someone here has?


title: Video/Audio online material for getting into Bayesian analysis and logistic-regressions

text: Well, I'm an engineer by day. Although most of my work revolves around modeling, we generally do pretty basic stuff. An "Advanced" model would be a monte carlo simulation validated using R2 tests. Currently, in my field, there is a lot of research using Logistic and bayesian analysis. My question is, which courses would you recommend someone to take from MIT's open course site or any other sites, for someone who learns best by video/audio first, and reading second? What i'd like to learn are the following:   * Be able to understand the models and when to employ them   * able to take in field data (which is generated once and cannot be regenerated) and design and perform experiments   * Able to understand the results, look at them, and figure out if something is off, "show stopper" or "outliers", or if everything is fine and dandy   * Be able to validate and calibrate the model, to actual "As-built" results   * Be able to forecast the results using appropriate sensitivity analysis   * be able to forecast / "plug" missing data   * be able to write journal papers related to my field my field in a nutshell is: transportation demand modeling for passenger vehicles, using either the generic four step model, or socio economic activity/tour based models such as PECAS or urbansim


title: how to test if a proportion of alleles is different between conditions?

text: I am in doubt about which type of test I could performed on my data. I have a data frame of the following type: contig, pos, A1, A2, Ctrl_A1, Ctrl_A2, Treat_A1, Treat_A2 comp1293557_c0_seq1, 284, T, C, 38, 20, 43, 17, comp1478014_c0_seq1, 252, A, G, 27, 14, 30, 16, comp15049_c0_seq1, 209, A, G, 56, 11, 28, 26 contig = the name of the "gene" which has a variant, pos = the position of the variant (mutation) A1 corresponds to Allele 1, A2 to Allele 2. The numbers in the table are the number of occurrences of allele 1 or 2 in either the control or the treatment samples. What I am interested in, is whether the treatment sample uses significantly different proportions of the alleles A1 and A2 when compared to the control. And that for each row. I thought about applying a Chi-squared test or an Exact test of Fisher to each row but I seriously doubt this is the way to answer my question. Do you have any insight or ideas to solve this problem? Thanks, Sophie


title: Crosscorrelation of stochastic process

text: Let $Z_1,Z_2 $ i.i.d. standard normal $$ X(t) = \begin{cases} 0, & \text{if } t<Z_1, t<Z_2\\\ 1, & \text{if } Z_1\le t <Z_2 \text{ or } Z_2\le t <Z_1 \\\ 2, & \text{if } t\ge Z_1, t\ge Z_2 \end{cases} $$ ## Asked $E[X(t)X(t')]$ ## Approach $E[X(t)X(t')] = P(X(t)X(t')=1)+2 P(X(t)X(t')=2)+ 4 P(X(t)X(t')=4)$ I work out $P(X(t)X(t')=1)$ as follows: $P(X(t)X(t')=1)=2 P(Z_1<t<Z_2,Z_1<t'<Z_2)$ (x2 because switching Z1,Z2 gives same probability) If $t'>t$ $P(Z_1<t<Z_2,Z_1<t'<Z_2)= 2*P(Z_1<t'<Z_2 |Z_1<t<Z_2)P(Z_1<t<Z_2)$ $=2P(t'<Z_2)P(Z_1<t<Z_2)$ $=2(1-\Phi(t')) * \Phi(t)(1-\Phi(t))$ This can be generalized to $P(X(t)X(t')=1)=2 (1-\Phi(\max(t,t'))\Phi(\min(t,t'))(1-\Phi(\min(t,t'))$ However, according to my solutions this is wrong. ## Correct answer $P(X(t)X(t')=1)=2 (1-\Phi(\max(t,t'))\Phi(\min(t,t'))$ ## My Question What am I doing wrong?


title: Training sample division vs extra binary/nominal factor

text: Suppose, there is some classification/regression problem. It seems for me, that sometimes division of training sample by some feature (or maybe by some other reasonable method, say feature generation and then devision by it.) and then fitting two models is better then fitting one model. For example, you anylize life duration of people of two contries. There are lots of medical, biological, etc. features of people. you can add one more binary feature that takes value 0 if a man from country_1, and 1 if he is from country_2. What if these contries are "very different"? Intuitively I feel, that fitting models separately for two countries can even improve forecasting quality, despite the fact of lower traing sample. The data coming from one country seems to be more specific and "purer". So, under what conditions devision of training sample is reasonable? What rule of thumbs or explanations in plain English could you advice?


title: How to interpret variance ratio test results?

text: I should do Variance Ratio Test on two vectors in R. I found the a builtin function: **var.test** but I don't know how to interpret its results (p-value). What result do I have to wait if these two vectors have constant variance? I think I also could use it to understand if the ratio is mean reverting or not, right? Thank you


title: Markowitz portfolio mean variance optimization in R

text: I have 5 emerging market foreign exchange total return series, for which I am forecasting single period future returns (1 year). I would like to construct a Markowitz mean variance optimized portfolio of the 5 series, using historical variances and covariances (1) and my own forecast expected returns. Does R have an (easy) way/library to do this? In addition how would I go about calculating (1) is there a built in function? For interest sake my currencies are USDTRY, USDZAR, USDRUB, USDHUF and USDPLN.


title: Comparing mclust() and k-means centroids

text: I have some code that looks for clusters in x,y data. To check the number of clusters I use, I want to get the BIC. This is not possible (easily) using `kmeans()`, and so I've switched to the mclust package. Specifically, I'm trying to replace `kmeans()` from the R stats package, with `Mclust()` from the mclust package. Using `Mclust()` requires me to specify which model should be used for the clustering. According to `?Mclust`, the following models can be used in `Mclust()`:               univariate mixture           "E"  =   equal variance (one-dimensional)     "V"  =   variable variance (one-dimensional)     multivariate mixture             "EII"    =   spherical, equal volume     "VII"    =   spherical, unequal volume     "EEI"    =   diagonal, equal volume and shape     "VEI"    =   diagonal, varying volume, equal shape     "EVI"    =   diagonal, equal volume, varying shape     "VVI"    =   diagonal, varying volume and shape     "EEE"    =   ellipsoidal, equal volume, shape, and orientation     "EEV"    =   ellipsoidal, equal volume and equal shape     "VEV"    =   ellipsoidal, equal shape     "VVV"    =   ellipsoidal, varying volume, shape, and orientation     single component             "X"  =   univariate normal     "XII"    =   spherical multivariate normal     "XXI"    =   diagonal multivariate normal     "XXX"    =   ellipsoidal multivariate normal      I'm presuming that k-means in stats is a "spherical, unequal volume" model, ie. to get `k-means(x = data, centers = 6)` to match `mclust()`, I should use `mclust(data, G = 6, modelNames = c("VII"))`. However, in the limited tests I've done, this gives different cluster centroids. The example below uses 6 clusters with some test data. The centroids obtained through each method are shown. ![enter image description here](http://i.stack.imgur.com/NJrQ0.png) Can anyone confirm which `mclust()` model is equivalent to `kmeans()`?


title: GLM for comparing 2 populations with binomial data?

text: I've been using R to run GLMs with the logit link to compare clutch sizes (a binomial data set) and proportional data among several years. I now need to compare 2 groups of years (good years and poor years) between 2 sites for the same data. I am not interested in any interactions between year-type (good or poor) and site, but really am just interested in comparing good vs. good and poor vs. poor between sites. If my data were continuous, I'd use a 2-way ANOVA followed up with t-tests for each comparisons, using the Holm–Bonferroni method for these multiple (2) comparisons. However, I'm not entirely sure how to approach these comparisons with binomial data. To start, I'm considering clutch sizes (proportional data will come next). I'm thinking that this code (below) is akin to the 2-way ANOVA for binomial data?               model<-glm(clutchsize~Site+Year_type, binomial)      Then, since I got a significant result here, I can go forward and make direct comparisons between the good years and the poor 'independently' of each other (using the Holm–Bonferroni method). But how can I make those comparisons? Can I use the `glm` model for that as well, like this:               model<-glm(clsz_gd~Site_gd,binomial)      where `clsz_gd` is binomial and `Site_gd` is a factor with only 2 levels? I presume I can't use t-tests! Is this the way to approach these data? Thanks!


title: Efficient way to merge multiple dataframes in R

text: I am building a big data frame by merging the content of a few files together. These files share the same columns layout.               c = read.delim('bigfile1.txt')     c1 = read.delim('bigfile2.txt')     c2 = read.delim('bigfile3.txt')          ctmp1 = merge(c, c1, all=TRUE)     ctmp2 = merge(ctmp1, c2, all=TRUE)      Is the above code efficient? Should I reuse the same variable name instead, e.g.               tmp = merge(c, c1, all=TRUE)     tmp = merge(tmp, c2, all=TRUE)


title: Sampling from a posterior distribution in SAS

text: Lets say I have a dataset where I want to estimate the relative risk of outcome X based on a binary treatment level Y, using PROC GENMOD to fit a logistic regression model. I can use the BAYES statement to have SAS 9.2 produce nice, MCMC-based prior distributions for my data. My question is, what if I then want to sample _from_ that posterior distribution? For example, if I estimate a OR and its 95% credible interval for the effect of Y, and subsequently want to use that estimate in another analysis. If I wanted to perform a sensitivity analysis based on the possibility of error in my estimate of the effect of Y, it makes sense to just run a Monte Carlo simulation based on the distribution of Y. But in this case, rather than being a really straightforward random number generation, I need to pull from that prior. Is the best way to do this using the OUTPOST option to get the data set of generated posterior samples and then just randomly sample from that new data set? Or is there some other, more clever mechanism I'm not seeing?


title: Multivariate multiple regression in R

text: I have 2 dependent variables (DVs) each of whose score may be influenced by the set of 7 independent variables (IVs). DVs are continuous, while the set of IVs consists of a mix of continuous and binary coded variables. (In code below continuous variables are written in upper case letters and binary variables in lower case letters.) The aim of the study is to uncover how these DVs are influenced by IVs variables. I proposed the following multivariate multiple regression (MMR) model:               my.model <- lm(cbind(A, B) ~ c + d + e + f + g + H + I)      To interpret the results I call two statements:   1. `summary(manova(my.model))`   2. `Manova(my.model)` Outputs from both calls are pasted below and are significantly different. Can somebody please explain which statement among the two should be picked to properly summarize the results of MMR, and why? Any suggestion would be greatly appreciated. Output using `summary(manova(my.model))` statement:               > summary(manova(my.model))                Df   Pillai approx F num Df den Df    Pr(>F)         c           1 0.105295   5.8255      2     99  0.004057 **      d           1 0.085131   4.6061      2     99  0.012225 *       e           1 0.007886   0.3935      2     99  0.675773         f           1 0.036121   1.8550      2     99  0.161854         g           1 0.002103   0.1043      2     99  0.901049         H           1 0.228766  14.6828      2     99 2.605e-06 ***     I           1 0.011752   0.5887      2     99  0.556999         Residuals 100                                                   ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1      Output using `Manova(my.model)` statement:               > library(car)     > Manova(my.model)          Type II MANOVA Tests: Pillai test statistic       Df test stat approx F num Df den Df    Pr(>F)         c  1  0.030928   1.5798      2     99   0.21117         d  1  0.079422   4.2706      2     99   0.01663 *       e  1  0.003067   0.1523      2     99   0.85893         f  1  0.029812   1.5210      2     99   0.22355         g  1  0.004331   0.2153      2     99   0.80668         H  1  0.229303  14.7276      2     99 2.516e-06 ***     I  1  0.011752   0.5887      2     99   0.55700         ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


title: Generate two sequences that have a pre-specified rank correlation

text: I have two samples $ X_1,\cdots,X_n $ and $ Y_1, \cdots, Y_n $, and would like to re-order them so that they have a pre-specified value of rank correlation such as Kendall's tau. I do know that there are many ways to generate two samples that have a specific Pearson's correlation coefficient, but don't know whether there is a way to do this about Kendall's tau.


title: Standard error of mean vs. standard error of coefficient in Gaussian glm

text: I'm having trouble understanding why calculating the standard error of the mean: $$ \widehat{SE} = \frac{\widehat{SD}}{\sqrt{n}} $$ gives me a different SE estimate than the SE estimate of a coefficient from a Gaussian GLM (identity link function) in R. Can someone explain why the maximum-likelihood estimators of the standard error of the model coefficients are different from the simple formula for the standard error of the mean? Here's the code I'm using to play around with this:               set.seed(55)       # make up normally-distributed data     x    <- c(rep('Treatment1', 20), rep('Treatment2', 20))     y    <- c(rnorm(n=20, mean=5, sd=5), rnorm(n=20, mean=8, sd=5))     dat1 <- data.frame(X=x, Y=y)     dat1       # calculate summary stats     library(psych)     sum1 <- describeBy(x=dat1, group=dat1$X)            # run glm with normally-distributed residuals     m1 <- glm(y~x, family='gaussian', data=dat1)     summary(m1)            # create summary table     table1 <- data.frame(ID=c('Mean', 'Coefficient', 'SE', 'Coefficient_SE',                                'Upper 95% CI     (mean)', 'Upper 95% CI (Coefficient)',                                'Lower 95% CI (Mean)', 'Lower 95% CI (Coefficient)'),                          Treatment1=c(  # mean                                       sum1$Treatment1$mean[2],                                          # coefficient                                       m1$coefficients[1],                                              # SE                                       sum1$Treatment1$se[2],                                           # Coefficient SE                                       summary(m1)$coefficients[1,2],                                          # Upper 95% CI (mean)                                       sum1$Treatment1$mean[2] + 1.96*sum1$Treatment1$se[2],                                          # Upper 95% CI (coefficient)                                       m1$coefficients[1] +1.96*summary(m1)$coefficients[1,2],                                         # Lower 95% CI (mean)                                       sum1$Treatment1$mean[2] - 1.96*sum1$Treatment1$se[2],                                         # Lower 95% CI (coefficient)                                       m1$coefficients[1]-1.96*summary(m1)$coefficients[1,2]),                          Treatment2=c(  # mean                                       sum1$Treatment2$mean[2],                                         # coefficient                                       m1$coefficients[2] + m1$coefficients[1],                                          # SE                                       sum1$Treatment2$se[2],                                           # Coefficient SE                                       summary(m1)$coefficients[2,2],                                           # Upper 95% CI (mean)                                       sum1$Treatment2$mean[2] + 1.96*sum1$Treatment2$se[2],                                          # Upper 95% CI (coefficient)                                       m1$coefficients[2] + m1$coefficients[1] +                                                 1.96*summary(m1)$coefficients[2,2],                                           # Lower 95% CI (mean)                                       sum1$Treatment2$mean[2] -                                             1.96*sum1$Treatment2$se[2],                                          # Lower 95% CI (coefficient)                                       m1$coefficients[2] + m1$coefficients[1] -                                             1.96*summary(m1)$coefficients[2,2])                          )


title: Trial with two sites, two years - Different number blocks - Problem? (Q1)

text: I have got two questions on an agricultural field trial that was conducted at two sites in two conscutive years. Virtually everything was the same in all trials (crop variety, planting density...). The data permits that all data can be analysed in a joined model. As there are two questions I open two threads and put the link to the other question here Trial with two sites, two years - lme model building and simplification (Q2)               Site 1, Year 1: 4 repetitions      Site 1, Year 2: 5 repetitions     Site 2, Year 1: 6 repetitions     Site 2, Year 2: 6 repetitions      I'm using a mixed model for data analysis (`lme` of `nlme` package in R - would also be open to lmer), modelling block (repetitions) as a random factor, nested inside **Site * Year**. The model seems to run fine (details about modelling problems in other question). Are there any problems resulting from the lack of repetitions or is the analysis valid this way?


title: how to represent geography or zip code in machine learning model or recommender system?

text: I am building a model and I think that geographic location is likely to be very good at predicting my target variable. I have the zip code of each of my users. I am not entirely sure about the best way to include zip code as a predictor feature in my model though. Although zip code is a number, it doesn't mean anything if the number goes up or down. I could binarize all 30,000 zip codes and then include them as features or new columns (e.g., {user_1: {61822: 1, 62118: 0, 62444: 0, etc.}}. However, this seems like it would add a ton of features to my model. Any thoughts on the best way to handle this situation? Thanks for any help.


title: Trial with two sites, two years - lme model building and simplification (Q2)

text: I have got two questions on an agricultural field trial that was conducted at two sites in two conscutive years. Virtually everything was the same in all trials (variety, spacing, planting density...). The data permits that all data can be analyised in a joined model. There is another question about the trial analysis. Here's the other question: Trial with two sites, two years - Different number blocks - Problem? (Q1) Factors in my model area.               Site (2 sites)     Year (2 years)     Block (4-6, see other question)       treat (6 treatments)      If toy data is needed than I will provide... Site should be treated as a fixed factor. It is known that the vegetation period is longer at Site 2 and that generally higher yields result. Year should - at first - be a random factor. As modelling random factors with less than 4-6 levels does apparently not make sense I also decide to treat it as a fixed factor. As blocks are nested inside trials (and treatments inside blocks) I create an interaction variable               data$trial <- with(data, interaction(Year, Site))       Then I setup the model like so:               object.lme0 <- lme(dm ~ (Site + Year) * treat, random = ~ 1 | trial / block, data)          > anova(object.lme0 )                  numDF denDF F-value p-value     (Intercept)     1    90    4270  <.0001     Site            1     1     185  0.0467     Year            1     1       0  0.9444     treat           5    90       1  0.5235     Site:treat      5    90       2  0.1204     Year:treat      5    90       0  0.8447      Am I choosing the right model? Is this a sensible approach to do the analysis over all sites and years or should trials be examined seperately or do I have to do both (seperate for each site and year (= trial) and joined analysis (like above))? Furthermore: I am confused about model simplification. My trial design is dictating the above model and some statisticians here told me that simplifying the model is bad statistical practice. What about the Year:treat interaction in the above model. Can / should it be kicked out?


title: OLS regression from first principles

text: I'm trying to work out how you do this from first principles. The Wikipedia page on linear regression gives me enough to solve it with matrix operations through the origin but I can't find much literature on implementing an algorithm for an OLS fit returning coefficients, a t-statistic and an $r^2$ value for the fit Can anyone point me to good reference? Thanks


title: Odds ratio terminology

text: Suppose you have two groups $A$ and $B$ and you want to find the odds ratio of some outcome. How do you denote this odds ratio? Is it just $OR(\text{A} \ \text{and} \ B)$?


title: Course of Action for 2x2 tables with 0's in cell and low cell counts

text: I'm comparing two groups A, B after cross tabulating them with cases and controls. I get a table as such:                   Control Cases       A   8       0       B  14       0      Obviously, I can't do an odds ratio because of the 0's. Doing a Fisher's exact test gives me a p-value of 0.5152 but the confidence interval goes from 0.16 to Infinity. Adding .5 to each cell doesn't change Fisher's exact test (although I'm not sure why the ODD's ratio is still infinite in R) and I'm hesitant to use a normality approximation due to the small cell counts. However, the signal seems to be quite strong that A and B both act the same... Neither of them are push an observation into the "case" category. (Here A & B are toxic substances). What would be some other ways I can tackle this to get this point across? Thank you!


title: Regression of family of marginal density functions

text: From a large sample of triples $(X, Y, U)$ I need to estimate a function $(x, y, u) \mapsto f(x, y, u)$, such that for each fixed $x, y$, the function $u \mapsto f(x, y, u)$ is a density function; namely, the marginal density function of $U$, given $X = x$ and $Y=y$. Is there a standard way of doing this? My initial thought was to do something along the lines of 1) For a collection of points in the $(x, y)$-grid, estimate (e.g. by some standard kernel density technique) $u \mapsto \widehat{f}_{x, y}(u)$ based on data points $(\tilde{x}, \tilde{y}, \tilde{u})$ where $\tilde{x}, \tilde{y}$ resides in some suitable neighbourhood about $(x, y)$ and $\tilde{u}$ varies freely. 2) Bilinearly interpolate between these points to obtain $f$. Is this a reasonable approach?


title: Should the distribution of composite scores be normally distributed?

text: **Should the distribution of composite scores be normal (bell curve)?** The composite scores that I have developed is from a weighting of their separate components. The component scores are in response to a normal survey (i.e. proper sampling proceedures were used). I have then weighted the scores (according to a basic formula that factors in their level of importance).


title: Why is my probit analysis resulting in numbers outside 0 and 1?

text: I am doing a probit on a binary dependent variable. The binary dependent variable has observed outcomes of 227 successes and 1704 failures. If I run predict after the probit I get probabilities that make sense for each set of observations, but I want to be able to make sense of the beta coefficients that I get from a straight probit analysis. The same thing happens with logit.


title: Can it be valid to standardize non-normal data going into a PCA?

text: I have a collection of highly non-normal behavioral variables measured for each animal in a given test. I would like to reduce to 2 or 3 scores to characterize the individual's test using PCA. I see in the literature where others have referenced Multivariate Statistics books to justify using non- normal data for this kind of use of PCA. The measures are in very different scales of units (e.g. one that runs to 30, others in the tens of thousands). If I use [R] prcomp with centering and scaling on my non-normal raw data, the resulting PC loadings make biological sense. But can one center and scale non- normal data like this?


title: Standard deviation of error from leave one out

text: I have a classification process which I test using leave-one-out. Precisely, how do I compute the standard deviation of the result? If I compute it normally, I get a very large deviation which makes sense because each test is either 100% or 0%. What is the accepted convention in this case? Thanks.


title: Colinearity and scaling when using k-means

text: I'm trying to gain a better understanding of kmeans clustering and am still unclear about colinearity and scaling of data. To explore colinearity, I made a plot of all five variables that I am considering shown in the figure below, along with a correlation calculation. ![colinearity](http://i.stack.imgur.com/W5MZJ.jpg) I started off with a larger number of parameters, and excluded any that had a correlation higher than 0.6 (an assumption I made). The five I choose to include are shown in this diagram. Then, I scaled the date using the `R` function `scale(x)` before applying the `kmeans()` function. However, I'm not sure whether `center = TRUE` and `scale = TRUE` should also be included as I don't understand the differences that these arguments make. (The `scale()` description is given as `scale(x, center = TRUE, scale = TRUE)`). Is the process that I describe an appropriate way of identifying clusters?


title: Is there a stats tool for this analysis I run in excel?

text: I am trying to find a statistical or machine learning tool that replicates this analysis I am doing manually in excel. Each row in my data set is a user.               User | sites visits | pages viewed | items added to cart | Purchased (1 or 0)      I set up a pivot table so Purchased is in the columns field. I have it set to % of row.               Purchased   Didn't Purchase         Total      80%              20%                100%      I then add all the other variables as filters in the Pivot Table. I start by filtering out the people who didn't view 1 page and see how that impacts that % of people who purchased. I then filter out the people who didn't view 2 pages, etc. I then start filtering out other variables to see how the interactions work together. I am looking for a combination of inputs that leads to an 80% purchase rate. So maybe it's 2 visits, 6 pageviews, and 2 add to carts. Is there a tool that automates this type of analysis?


title: Understanding the Canonical Correlation Analysis

text: In canonical correlation analysis where max { a'COV(X,Y)b / sqrt[a'COV(X,X)a * b'COV(Y,Y)b] } where U = a'X, and V = b'Y the standard tools to solve the canonical correlation analysis (MATLAB, R, etc.) yield solutions (as many as min rank of both variables) where CORR(U1,U2) = CORR(V1,V2) = CORR(U1,V2) = CORR(U2,V1) = 0 for all Ur, Uk, Vr, Vk where r is different from k. Correlations are not transitive so this seems too restrictive to me. I would like for example only Vs to be orthogonal among each other, and Us to not necessarily be orthogonal. Are these strict restrictions a must or just there to use the eigen decomposition techniques these tools are using. If not how can I solve this optimization problem? Thanks


title: Statsmodels OLS and MSE

text: So, I have data set and I calculate the model parameters and errors using statsmodels:               result = sm.OLS(y, X).fit()     result.summary()      Now, `result.mse_resid`, `result.mse_total` provide MSE of the residuals and total mean squared error. Reading the statistics tutorials I see that $MSE=\frac{RSS}{DFE}$, where DFE is degrees of freedom for error. I should be able to calculate MSE as follows:               reg=LinearRegression()     reg.fit(X,y)     yp=reg.predict(X)     resid=y-yp     rss=np.sum(resid**2)     MSE=rss/(result.nobs-2)      What is the MSE calculated using OLS and why is it different from this one (or what am I not understanding correctly)?


title: Confusion between empirical Bayes and Hierarchical Bayes Model for beta-binomial model

text: I am confuse with difference between empirical Bayes method and hierarchical Bayes method. Take an example (ref : http://www.cs.cmu.edu/~xuerui/papers/ctr.pdf). Suppose we have :   * $C$ is the number of clicks we observed on an advertising ad   * $r$ which represent the probability that an ad will be clicked (basically it is the ratio of $y$ and $n$ in frequentist way)   * $I$ is the number of impressions of an advertising (which mean the number of internet user that have seen an advertising)    * Natural hierarchy : advertising ads are grouped into an account ![Graphical Model Representation of the r estimation with empirical Bayes method](http://i.stack.imgur.com/C3GEL.png) Figure above :   * For each account that has $N$ ads we have a distribution Beta$(\alpha,\beta)$   * The observed clicks $C_{i}$ are assumed to be generated from an underlying -unobserved- $r_{i}$ according to a Binomial distribution Bin($I_{i}$,$r_{i}$), where $I_{i}$ are the number of observed impressions of the ad $i$    * The parameters $r_i$ are considered as samples drawn from a population charachterized by hyperparameters according to a Beta($\alpha,\beta$) Well, for me it is a beta-binomial model with a two stage empirical Bayes method (ML estimation) to estimate $r_i$ for each ad $i$. But I don't see any difference with the literature called "Hierarchical Bayes models". We can also add an extra "level" by assuming that the parameters $\alpha,\beta$ are drawn from a distribution? Somebody can explain it to me? Furthermore what is the difference in this model to observe or not the $r_{i}$? For each distribution that we have (binomial, beta) it is "conditionnal" distribution isn't it? For instance $C|I,r \sim$ Bin(I,r) ? I'm very confuse with all of that but I think I am quite close from understanding something that is new to me. Thank by advance for any advice! BimBimBap


title: Question about standardizing in ridge regression

text: Hey guys I found one or two papers which use ridge regression (for basketball data). I was always told to standardize my variables if I ran a ridge regression, but I was simply told to do this because ridge was scale variant (ridge regression wasn't really part of our course, so our lecturer skimmed through it). These papers I read didn't standardize their variables, which I found a bit surprising. They also ended up with large values of lambda (around the 2000-4000 level) via cross-validation, and I was told that this is due to not standardizing the variables. How exactly does leaving the variable(s) unstandardised lead to high lambda values and also, what are the consequences of not standardizing the variables in general? Is it really such a big deal? Any help is much appreciated.


title: Conditional kernel density plot with R's np package

text:
I tried to use the Kernel Density plot method from Hayfield and Racine (2008) `np` package for my own data, but somehow ended up with different type of plots and I have no idea what the difference is between my data and the example data provided by the package. I used the Italy GDP example that is provided in the np package and described in Racine's primer. The example itself worked well and I ended up with a nice 3D plot. When I used it with my own data I just got two 2D plots instead of one 3D plot. Why is that? Do I miss something content wise (that's what I think and it's the reason being why I guess this is rather a CV than a SO question) or is it a syntax problem?               # get a bandwith object       bw2 <- npcdensbw(formula=mydf$values~ordered(mydf$datefield), tol=.1, ftol=.1)        str(mydf)     # returns       data.frame':  780 obs. of  2 variables:       $ datefield: Ord.factor w/ 104 levels "1984-04-01"<"1984-07-01"<..: 1 1 1 1 1 1 1 1 1     1 ...
           $ values : num  50.7 58.5 56.1 55.5 62.7 ...     tables(mydf$datefield)     # shows that I have 13 entries per quarter.          quartz()     plot(bw2)     npplot(bw2)     # same 2D plot for both :(      all of these properties are exactly in line with the example data. Still I get only 2D plots instead of a 3D plot. I do have three dimensions: quarters,values and conditional density. What's wrong?


title: Exploring data attributes

text: I have a database with many attributes. I would like to know which attributes has the minimum variation in the data. Is there some standard technique? It should be like clustering without split records in clusters. I would like to know what the records in particular cluster have in common. I was going to compute the mean ($\bar{x}$) and st.d. ($s$) for each continuous attribute $x$. After computing the coefficient of variation $CV=\frac{s}{\bar{x}}$ I would say that attributes with $CV\leq0.1$ are the similar ones. For categorical ones I would choose attributes with more than $90\%$ relative frequency for the mode. Is there some standard technique?


title: Understanding detailed balance condition

text: I guess I understand the equation of the detailed balance condition, which states that for transition probability $q$ and stationary distribution $\pi$, a Markov Chain satisfies detailed balance if $$q(x|y)\pi(y)=q(y|x)\pi(x),$$ this makes more sense to me if I restate it as: $$\frac{q(x|y)}{q(y|x)}= \frac{\pi(x)}{\pi(y)}. $$ Basically, the probability of transition from state $x$ to state $y$ should be proportional to the ratio of their probability densities. My question is, how does MCMC fulfilling detailed balance yield the stationary distribution?


title: weak exogeneity in VAR analysis

text: I have a problem interpreting the notion of "weak exogeneity in a VAR process". Assuming we have the following structural form: $y_t = b_{12}z_t + \gamma_{11}y_{t-1} + \gamma_{12}z_{t-1} + \epsilon_{yt}$ $z_t = b_{21}y_t + \gamma_{21}y_{t-1} + \gamma_{22}z_{t-1} + \epsilon_{zt}$ Does weak exogeneity imply that either $b_{12}=0$ or $b_{21}=0$, but not both?


title: Can you justify this very common (but seldom discussed) treatment of $p$-values?

text: This is a question about a seemingly illogical approach to hypothesis testing and _p_ -value interpretation that I find almost universally endorsed. I’ll cite the answers at Can I trust ANOVA results for a non-normally distributed DV? as an example. Wanting to conduct, say, an ANOVA but faced with non-normal distributions of residuals, people make choices and give rationales for them as if there were just 3 options:   1. Decide that the residuals are virtually normal, and treat the _p_ -values of the resulting tests as if they were exact. Cite those values using as many decimal places as desired.   2. Decide that normality fits only roughly, but conduct a parametric test as if these data were normal, and again treat the _p_ -values as if they were exact—up to as many decimal places as desired.   3. Decide that the normality assumption is not tenable, and abandon the parametric test. Under this logic, we would have to believe that, as departures from normality become more and more severe, _nothing_ happens to the accuracy of _p_ -values --until some critical threshold, at which point they become completely invalid. ![-image-](http://i.stack.imgur.com/TRc59.jpg) Under a more nuanced view, there whould be another alternative, which is to conduct the parametric test despite _some_ non-normality and treat the resulting _p_ -values as useful approximations, good to perhaps one decimal place. We accept uncertainty and settle for estimates in so many other situations: why not treat _p_ -values the same way?


title: Top-coded survey data

text: Top-coded income data are often modeled with a Pareto distribution, but that is controversial. What would be wrong with declaring those values as missing and then using multiple imputation (MI) to generate values for those cases using a normal model while assigning a lower bound as is allowed in IVEware (as long as the number of missing cases is small, say <5%)? IVEware uses sequential regression multiple regression imputation. From p.2 in the article below the algortihm is: > The basic strategy creates imputations through a sequence of multiple > regressions, varying the type of regression model by the type of variable > being imputed. Covariates include all other variables observed or imputed > for that individual. The imputations are defined as draws from the posterior > predictive distribution specified by the regression model with a flat or > noninformative prior distribution for the parameters in the regression > model. The sequence of imputing missing values can be continued in a > cyclical manner, each time overwriting previously drawn values, building > interdependence among imputed values and exploiting the correlational > structure among covariates. To generate multiple imputations, the same > procedure can be applied with different random starting seeds or taking > every P th imputed set of values in the cycles mentioned above. From p. 3 the conditional regressions are based on: Each conditional regression is based on one of the following models: 1\. A normal linear regression model on a suitable scale (for example, a BoxCox power transformation may be used to achieve normality) if Y j is continuous; 2\. A logistic regression model if Y j is binary; 3\. A polytomous or generalized logit regression model if Y j categorical; 4\. A Poisson loglinear model if Y j is a count variable; and 5\. A twostage model where zeronon zero status is imputed using logistic regression, and conditional on nonzero status, a normal linear regression model is used to impute nonzero values, if Y j is mixed.


title: What happens to the variance of missing covariates in poisson models

text: In normal multivariate regression you can leave out covariates (assuming they aren't correlated with any covariates you leave in) and the variance of covariates left out gets subsumed into the error term. In survival analysis poisson models are used. In these models the variance is a function of the mean. What then happens to the variance of covariates left out of the model since there is no free parameter for estimating the error? My, maybe naive, thought is that the variance estimate for poisson models is only correct if the model is fully specified. But that would mean that the variance for any model leaving out covariates is underestimated. Which, in turn, would mean that the p-values for any such would be very optimistic. But such models are very common in epidemiology and I can't believe they would all be this flawed. What am I missing?


title: Confidence intervals for bootstrap and FPC

text: I have a population that I am technically sampling without replacement using a stratified design. The resultant ratio estimator of the sample mean in any stratum $i$ is $\hat{R}_{i} = (\sum y_{i})/(\sum x_{i})$. I also know that $x_{i} \in \theta_{x,i}$ in any stratum, such that I can obtain $\hat{\theta}_{y,i} = \hat{R}_{i}\theta_{x,i}$. My question, then, is whether I can treat this situation as simple random sampling with replacement, as I'm applying the ratio estimate to the known, total value of $\theta_{x,i}$. In other words, $\theta_{x,i} = N$ and not $\theta_{x,i} = N-n$. I'm using a non-parametric bootstrap routine to come up with confidence intervals. The $y_{i}$ are rare events, so I cannot assume that the bootstrap distribution can be normally approximated, so I'm just using the 95% highest density region of the bootstrap replicates. I'm curious if this is appropriate, or if I need to include the finite population correction factor somehow, as some strata are sampled pretty well (i.e. $n/N \geq 0.1$). Thanks @StasK for your suggestions, although I need a little clarification. I did not use the percentile method to obtain confidence intervals. Instead I found the HDR. In other words, if $f(x)$ is the density function of random variable $X$, then I found the region in the sample space of $X$ such that $HDR(f_{\alpha}) = \\{x: f(x) \geq f_{\alpha}\\}$, where $f_{\alpha}$ is the largest constant such that the $Pr(X \in HDR) \geq 1-\alpha$. I did this to somewhat account for the skewed nature of the bootstrapped distribution. Also, the distribution of the bootstrapped replicates is polymodal - yet another reason why I used the HDR. Do you still think this is inappropriate? If so, then what you're recommending is finding the following interval $(\hat{\theta}-\hat{t}^{(1-\alpha)}\cdot \hat{se},\hat{\theta}-\hat{t}^{(\alpha)}\cdot \hat{se})$. Is that right?


title: ROC50 score - how is it calculated?

text: I have never found a good explanation of what the ROC50 score is and how it is calculated. To take an example definition found on the net: > The ROC50 score is the area under the ROC curve, up to the first 50 false > positives So if I understand correctly, you look at how many correct positive classifications you have, until you get the first fifty wrong ones. But this makes the ROC50 score dependent on which positive classifications you look at first- if they are true or false. What ordering does one use? Does one look at the objects most confidently classified as positive first? I.e. for an SVM-classified list of objects you begin with the object with the highest score from the decision function and work your way downwards until you get to fifty false positives? If so, is it correct to say that the ROC50 score gives you information about how well the most confidently classified objects are classified?


title: Internal reliability for ordinal scale using Spearman correlation matrix

text: I have a question about estimating the reliability (coefficient alpha) of a pilot scale by assessing its internal consistency. The pilot questionnaire generates ordinal data, 9 items each with 9 ordered catagories of response. The data set has 183 complete readings of the pilot scale. My question is: Is it legitimate to calcualte alpha using cronbach's formula from a correlation matrix consisting of the spearmans's rank order correlation coefficients between the 9 items?


title: Model for probability of song reaching top 10 ranking, over time?

text: I'm trying to model the probability of a song reaching Billboards top 10 over time. My data has the columns "Day since release", "If reached top 10". For example, [12,1] means the song hit top 10 on the 12th day since release, and [350,0] means the song has been out for 350 days and has not become a top 10 hit yet. Any suggestions for the best type of technique to approach this problem? Thanks for any help!


title: Pattern recognition training base division

text: 1) A dataset (100 samples, 10 classes, 8 chars each) was multiclassed obtaining a good accuracy with an external dataset. 2) The same dataset was divided in half (50 samples, 5 classes, 8 chars each) and trained as 2 classification models (like, model 1 with classes 1-5, and model 2 with classes 6-10). Why would the results of 1 be better then 2? Any ideas? Theoretically, my assumption is that 2 has "a local" view of the whole dataset, and should get better results.


title: Reasons for transforming multiple class classification problem into a set of binary sub-problems?

text: Does anyone know of a good reference that list the reasons for transforming multiple class classification problem into a set of binary sub-problems? In response to comment: One reason to transform a multiple class classification problem into a set of binary sub-problems are that a binary network is usually simpler and therefore better suited to fit small data sets (Less prone to over fitting). Furnkranz argues that a set of simpler binary classifiers can also provide better performance than one complex multiple class classifier because it is often easier to learn how to distinguish between two classes, than it is to learn how to distinguish between multiple classes. Simpler models reduce training time also. Milgram shows that a set of binary SVM can preform better than a multi-class MLP. Though this might be due to choice of classifier rather than binarization. If we focus on a MLP classifier, training a set of binary one-against-all MLPs and combining using softmax function would be very much equal the the multi- class MLP however this would allow us to use different hyper-paramater settings for each class.


title: The biological significance of ANOVA 3-way interaction

text: Dear statistics experts, I have trouble to find a sensible statistical approach to back up some very obvious (at least to my eyes) interpretation of a dataset (see descriptive plot below). I measure a response to two different stimulus (Stimulus 1 and 2) at 3 different intensities (weak medium high), in 2 different populations (population 1 and 2). What I see and would like to support with a statistical test, is that overall the response exhibit a dynamic range along the intensity range, but more importantly, that this range is very steep for stimulus 1 in the population 2 but, not in population 1. In terms of design, I have   * a continuous dependent variable (my response).   * two categorical variables with two levels each (stimulus and population)   * one ordinal variable with 3 levels (stimulus strength) Samples are not independent, as response to each stimulus/strength pair is assessed on the same subjects in the two populations. I thought of using a mixed effect ANOVA with subject as a random variable to account for the repeated measure. First problem : How to measure the departure from the homoscedasticity assumption ? A Levene's test reject the null hypothesis of equality of variance between groups. Anova is said to be robust to violations of its assomptions, but is there a way of judging the severity of this violation, especially given my unbalanced design? Second, would my biological statements be reinforced if I was to find a significant 3-way interaction population x stimulus x strength ? I have two concerns :   * Does a 3-way interaction p-value have any meaning on such a design and especially given my first point (disparity in variance) ?   * even if this effect was real, is it allowed (not to say feasible), to conclude anything biologically relevant from a three-way significant interaction, which essentially mean, (correct me if I am wrong), that at least one of the two-way interaction between a pair of factors is influenced by the third one ? Specifically here, I would like to test if the stimulus x strength interaction is influenced by the population factor. Should I then test for significance of this two-way interaction beforehand ?  My meditations have carried me far beyond what my very basic understanding of statistics allow me to comprehend. Thanks in advance for your input, Benjamin ![enter image description here](http://i.stack.imgur.com/7lc9g.jpg)


title: Outliers and the mean

text: I would like to know what the following example is called in mathematics: In a gymnastics competition the judges scored a competitor as 10, 8, 3, 7, 7, 9, and 8. I recall that the ending score was always calculated after they discarded the highest and lowest numbers in the group and used what was remaining to get the resulting score. What is this kind of "math" called? Would this be an example of discarding the outliers? And is the result a mean of the remaining numbers?


title: Finding a correlation or statistically significant relationship in data

text: I have a data set that is represented below, but much longer. How can I find if there is a statistical relationship between business category and business score?               Business Name     Business Category        Business Score     Business1         Shopping                      1     Business2         Eating                        1     Business3         Shopping                      2     Business4         Shopping                      2     Business5         Eating                        1     Business6         Manicure                      3     Business7         Haircut                       2     Business8         Haircut                       3     Business9         Eating                        3     Business10        Manicure                      3     Business11        Haircut                       3     Business12        Eating                        1     Business13        Shopping                      2     Business14        Eating                        3     Business15        Manicure                      1     Business16        Haircut                       3     Business17        Haircut                       2     Business18        Shopping                      3     Business19        Eating                        2     Business20        Manicure                      3


Suppose we have a time series $X_t$ s.t. $X_t \sim^{iid} (0,1)$. How do you prove that if $ X_t \sim^{iid} (0,1) $, then $ E(X_t^{2}X_{t-j}^{2}) = E(X_t^{2})E(X_{t-j}^{2})$? Or, I guess, if $X,Y\sim^{iid} (0,1)$ (that is, $E(XY)=E(X)E(Y)$), why then is it that $E(X^2Y^2)=E(X^2)E(Y^2)$? This spins-off from another question wherein apparently "If the squares were dependent, there's a form of dependence among the unsquared values." This makes sense, but how does one prove this exactly? My attempt: Instead of dependence => dependence (which I think would involve probability distributions), I try to prove uncorrelatedness => uncorrelatedness as follows: $E(X^2Y^2) \neq E(X^2)E(Y^2)$ $\implies E(X^2Y^2) \neq (Var(X)+E(X)^2)(Var(Y)+E(X)^2)$ $\implies Var(XY)+E(XY)^2 \neq (Var(X)+E(X)^2)(Var(Y)+E(Y)^2)$ $\implies Var(XY)+(E(X)E(Y))^2 \neq (Var(X)+E(X)^2)(Var(Y)+E(Y)^2)$ $\implies Var(XY)+(E(X)E(Y))^2 \neq Var(X)Var(Y)+Var(X)E(Y)^2+Var(Y)E(X)^2+(E(X)E(Y))^2$ $\implies Var(XY) \neq Var(X)Var(Y)+Var(X)E(Y)^2+Var(Y)E(X)^2$ $\implies ...$ $\implies E(XY) \neq E(X)E(Y) \ QED$ Ugh...


title: Do these residual plots violate the linearity and homogeneity assumptions for linear regression?

text: There seem to be too many points clustered around negative values for all the plots And while 3 & 4 seem to have random enough patterns, 1 & 2 seems to have negatively sloped trend. If these were to violate linearity and homogeneity assumption, I should stop using the regression model, correct? ![Residual plot variable 1](http://i.stack.imgur.com/oVR4u.png) ![Residual plot variable 2](http://i.stack.imgur.com/jY8yb.png) ![Residual plot variable 3](http://i.stack.imgur.com/dpqmX.png) ![Residual Plot variable 4](http://i.stack.imgur.com/EKpXn.png)


title: Rank histogram (or Talagrand Diagram) of precipitation data

text: Given a dataset of lines containing 6 wind **forecast** values plus 1 observed (actual) value in each like:                FCT1   FCT2     FCT3   FCT4   FCT5     FCT6   OBSERVED     -3.17   3.51    -5.71   1.37   -0.22   -0.65    -2.38        -2.7    2.21    -0.71   2.73   -0.33   -2.62    -1.38        -1.2    3.15    -4.17   3.33   -0.48   -1.65    -2.30        ...     -3.0    3.50    -1.79   3.37   -0.18   -0.62    -2.32         To make a rank histogram (or Talagrand Diagram), I understand that I need to loop through the lines, sort the **forecast** values for each one and assume that the ordered values (six, in this case) are the inner limits of each bin on the diagram. 6 limits generates 7 bins. Then, I need to take the corresponding **observed** value and increase the bin it fits (the bin that its range contains the observed value). I need to do it for every row, so each row has its limits. It has to do with What PDF should be fit to a rank histogram?. I think it is not a simple histogram built via `hist()` in R. Am I wrong? > ![enter image description here](http://i.stack.imgur.com/KlyGa.png) How about a precipitation **forecast** data? Like:                FCT1   FCT2   FCT3   FCT4   FCT5   FCT6   OBSERVED       0       0     0.1    0     0        0        0       0       0     0      0     0.02     0        0       0       0.1   0      0     0        0        3      What am I supposed to do to know the bin that has the range that fits the **observed** value like 0, for example? How can I build 7 bins from this **forecast** data?