title: Why are 0.05 < p < 0.95 results called false positives?

text: ~~**Edit:** The basis of my question is flawed, and I need to spend some time figuring out whether it can even be made to make sense.~~ **Edit 2:** Clarifying that I recognize that a p-value isn't a direct measure of the probability of a null hypothesis, but that I assume that the closer a p-value is to 1, the more likely it is that a hypothesis has been chosen for experimental testing whose corresponding null hypothesis is true, while the closer a p-value is to 0, the more likely it is that a hypothesis has been chosen for experimental testing whose corresponding null hypothesis is false. I can't see how this is false unless the set of all hypotheses (or all hypotheses picked for experiments) is somehow pathological. **Edit 3:** I think I'm still not using clear terminology to ask my question. As lottery numbers are read out, and you match them to your ticket one-by-one, something changes. The probability that you have won does not change, but the probability that you can turn the radio off does. There's a similar change that happens when experiments are done, but I have a feeling that the terminology I'm using - "p-values change the likelihood that a true hypothesis has been chosen" - isn't the correct terminology. **Edit 4:** I've received two amazingly detailed and informative answers that contain a wealth of information for me to work through. I'll vote them both up now and then come back to accept one when I've learned enough from both answers to know that they've either answered or invalidated my question. This question opened a much bigger can of worms than the one I was expecting to eat. In papers I've read, I've seen results with p > 0.05 after validation called "false positives". However, isn't it still more likely than not that I've chosen a hypothesis to test with a false corresponding null hypothesis when the experimental data has a p ~~< 0.50~~ which is low but > 0.05, and aren't both the null hypothesis and the research hypothesis statistically uncertain/insignificant (given the conventional statistical significance cutoff) anywhere between 0.05 < p < ~~0.95~~ whatever the inverse of p < 0.05 is, given the asymmetry pointed out in @NickStauner's link? Let's call that number A, and define it as the p-value which says the same thing about the likelihood that you've picked a true null hypothesis for your experiment/analysis that a p-value of 0.05 says about the likelihood that you've picked a true non-null hypothesis for your experiment/analysis. Doesn't 0.05 < p < A just say, "Your sample size wasn't big enough to answer the question, and you won't be able to judge application/real-world significance until you get a bigger sample and get your statistical significance sorted out"? In other words, shouldn't it be correct to call a result definitely false (rather than simply unsupported) if and only if p > A? This seems straightforward to me, but such widespread usage tells me that I might be wrong. Am I: a) misinterpreting the mathematics,   b) complaining about a harmless-if-not-exactly-correct convention,   c) completely correct, or   d) other? I recognize that this sounds like a call for opinions, but this seems like a question with a definite mathematically correct answer (once a significance cutoff is set) that either I or (almost) everybody else is getting wrong. Thanks.


title: Calculating agreement between an ordinal and continuous scale

text: I'm having a bit of trouble wrapping my head around something. I have a data set which contains two columns that essentially attempt to measure the same thing, one on a 1-50 continuous scale and another on a ordinal (5 levels) scale. Let's say they measure the risk of a comet impact (just an example). However the way the two columns are obtained is vastly different, the 1-50 scale is generated by an array of simulations, whereas the ordinal scale is based on expert opinion. Now I want to know how much these two measurements "agree" with each other. So what I did was: For the the continuous scale (call it $u$) I did a Kernel Density Estimation for _each_ of the levels of the export opinion (call it $w \in W$). This estimates Probability Density Functions for each of the classes in $w$ given $u$. By Bayes rule I then computed the posterior probability $P(x|x \sim w_i)$ with $x$ being some sample (row) and $w_i$ some class of the expert opinion. Now this sounded reasonable, to me. Except for choosing the priors (uniform, percentage in the original data, empirical?) and the density estimation (type of kernel, or is spline fitting better...) But here it got iffy. Remember that I want to know agreement, and perhaps more importantly the samples on which the simulation and the expert disagree. So for each of the samples I calculated all the posterior probabilities and then picked the maximum one. This allows me to calculate a Cohens Kappa for agreement and pick the rows that differ for further inspection (who was right?). But this all seems a bit, well, wrong. What am I missing here?


title: How do we relate fuzzy logic and mahanobolis distance?

text: What is the exact relation between fuzzy logic and mahanobolis distance? how can we see them in one application?


title: How will studying "stochastic processes" help me as a statistician?

text: I wish to decide if I should take a course called "INTRODUCTION TO STOCHASTIC PROCESSES" which will be held next semester in my University. I asked the lecturer how studying such a course would help me as a statistician, he said that since he comes from probability, he knows very little of statistics and doesn't know how to answer my question. I can make an un-educated guess that stochastic processes are important in statistics. But I am also curious to know how. That is, **in what fields/methods, will basic understanding in "stochastic processes" will help me do better statistics?**


title: Spearman's rank correlation coefficient over population

text: I have a data set composed by 2.000.000 projects, each of them is defined by the number of active users and the size (KBs) of the project. I would like to check whether a correlation exists between these two features. Does it make sense to use the Spearman's rank correlation coefficent over the entire population? thanks in advance


title: Using R for linear regression to see if a linear regression model appropriate for the data set

text: This is my R code and running result:(See below) How to judge is the linear regression model appropriate for this data set? Except R^2 value, Can the p-value in last row of the running result mean something? What does this p-value mean and can it mean the linear regression model appropriate for this data set? Why? Thanks in advance.               x=c(7,12,10,10,14,25,30,25,18,10,4,6)         y=c(128,213,191,178,205,446,540,547,324,117,75,107)         list(x,y)         reg1 <- lm(y~x)         summary(reg1)         plot(x, y)     abline(reg1)         reg2 <- lm(y~x-1)         reg2         summary(reg2)     #------------------         Call:     lm(formula = y ~ x)          Residuals:         Min      1Q  Median      3Q     Max      -55.805 -21.085   3.139  14.946  80.859           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  -22.753     21.846  -1.041    0.322         x             19.556      1.335  14.652 4.38e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 37.23 on 10 degrees of freedom     Multiple R-squared:  0.9555,    Adjusted R-squared:  0.951      F-statistic: 214.7 on 1 and 10 DF,  p-value: 4.38e-08


title: Op and op Convergence Property Related Question

text: There are some properties of the Op and op notation stated in many statistics textbooks but their proof is given nowhere. I have managed to prove most of them, but I cannot see how the following holds: $$ (1+o_p(1))^{-1}=O_p(1) $$ Can anyone show why this property holds in a thorough manner? Thanks!


title: On nonlinear regression, fits, and transformations

text: I am trying to fit a nonlinear regression model in R using `nls()`. I have a form of the equation I want to fit to: $$y = (a \times x_{1}^c +b \times x_{2}^d) (x_{3}^e)$$ where the coefficients to be found in regression are a,b,c,d, and e. My data is output from a simulation model where $x_{1}$, $x_{2}$, and $x_{3}$ are all integers from $0$ to $10$, with the condition that $x_{1} + x_{2} + x_{3} \le 10$. $y$ is also integer valued and ranges from $0$ to roughly $1000$. The objective is to fit these data to a rate function that will be used in a Markov Chain. When I try to fit this regression model directly using `nls()`, my `nlsResiduals` plot looks like this: ![nls residuals](http://i.stack.imgur.com/6scJ3.png) I know that autocorrelated residuals are problematic, and that non-normal residuals can also be problematic. How can I fix this problem? I was thinking of using transforms on the data like $$\log(y) = \log((a \times x_{1}^c +b \times x_{2}^d) (x_{3}^e))$$ or $$y^{1/n} = ((a \times x_{1}^c +b \times x_{2}^d) (x_{3}^e))^{1/n}$$ where $n > 1$. I've noticed if $n$ increases, my autocorrelation graph and QQ-plot look "better" (i.e., more scattered and more normal, respectively). Both of these seem to correct a lot (but not all) of the autocorrelated residuals, and help to make the residuals more normally distributed. Am I on the right track here, or am I committing some cardinal sin in statistics? Once I settle on a transformation, how can I tell which is best? Any help, suggestions, or comments are very appreciated.


title: Pooled variance

text: I have data for two years on plant characters for a given number of lines. since it involves same parameters, ANOVA is calculated separetely for two years. then how do I test the significance of Pooled variance as there is mention about its calculation but no table value(if exists) and draw conclusion? the main objective is to obtain more precise analysis by pooling two years data.kindly help..


title: What is the best book about generalized linear models for novices?

text: I'm still pretty new to generalized linear models, and I struggle with a lot of the notation in most of the GLM texts I've picked up. Are there extremely popular GLM books that lend themselves better to readability?


title: Why might it be important to specify a family when using GEE if it doesn't make any assumption about the joint distribution?

text: I understand that one of the advantages of GEE is that you make no assumptions about the joint probability distribution, relying instead on the mean, the variance and the associations (corr). Why then do we specify a family in the command for `gee`?               model <- geeglm(y~x, data=data, id=id, family=gaussian('identity'), corstr='ar1')


title: Using low-discrepancy sequence for bernoulli trials in MC sim

text: I need to generate binomial distribution random numbers for my Carlo simulation (I need Bernoulli trials for a parameter). Thus far, I've used R "rbinom" function for that. However, as I understand, I can more effectively use the parameter space by utilizing a low-discrepancy sequence for generating random numbers. Is it possible to use Sobol sequence to generate Binomial random numbers? Or does it even make sense? If yes, then how to convert the sobol sequence to bernoulli or binomial distribution? Thanks!


title: Sampling from the normal-gamma distribution in R

text: Does anyone know of a way to sample from the normal-gamma bivariate distribution or the normal-inverse-gamma bivariate distribution in R? I could create the distribution myself as a function, but then I would not know what to do to sample from it.


title: Understanding the construction of Dirichlet process

text: I'm trying to understand the construction process of DP, however, with little background in measure theory, the original papers are hard to read, but I believe the ideas behind these papers can be followed. Let $y_1, y_2, \ldots, y_n$ be i.i.d. with distribution function $F$, which is characterized by $k$ unknown probabilities $\theta_1, \ldots, \theta_k$, i.e., $y_i$s are discrete, taking a finite number $k$ of possible values. Bayesian inference based on these observations often put a Dirichlet prior on the parameters of $F$. Posterior and predictive distribution can then be calculated. That finite mixture models is not hard to understand, but... When the $y_i$s can take infinitely many possible values ($k \rightarrow \infty$), then the Bayesian nonparametric models can deal with this situation but I don't understand why. Here are my questions:   1. $k \rightarrow \infty$ means that $y_i$ can take any point mass on $\mathcal{R}$, so it is a continuous r.v. isn't it?    2. In the finite case, $F$ is indexed by $1, \ldots, k$, while in the infinite case, $F$ should be indexed by partitions $(T_1, \ldots, T_k)$ (Ferguson, 1972), how can I get the intuition of the latter case from the former?


title: Is time-delay embedding/attractor reconstruction used in some machine learning algorithm?

text: I try to model/forecast blood glucose levels from my diabetes diary, so I have to deal with some 5-7 daily measurements of estimated carbohydrates, physical activity, insulin doses and measured blood glucose. It is an unevenly spaced time series I've got (and it is not so densely sampled that for example it would reveal the peaks in blood glucose in the hour after a meal). As all this data is from some kind of a control cycle, I stumbled upon the topic of Time Delay Embedding or Attractor Reconstruction. I'd like to know if there is some standard machine learning technique that makes use of underlying control cycle structure? Thank you very much for your thoughts. Damian


title: What is standard error used for?

text: I'm using a tutorial I found and plotting mean values along with the standard errors to show my data. But I'm having a problem discussing the results. My plot is as shown below: some of the standard errors (shown as a error bar) vary much and some of them are very close to zero. ![enter image description here](http://i.stack.imgur.com/3sUze.jpg)


title: What's the difference between binomial regression and logistic regression?

text: I've always thought of logistic regression as simply a special case of binomial regression where the link function is the logistic function (instead of, say, a probit function). From reading the answers on another question I had, though, it sounds like I might be confused, and there is a difference between logistic regression and binomial regression with a logistic link. What's the difference?


title: Pairwise compatibility metric

text: I work at a company that sells clothes, and I've had good results with using cosine similarity to determine which products are "similar" to each other simply based on who owns them. I wanted to take product recommendations a step further and recommend a "compatible" product (like a solid tie for a striped shirt) using a new dataset—so not just who owns what. I have a training set with this structure:   * Every product has feature data like color, pattern, etc. (assume that these are ordinal measurements).   * I have a yes/no target on some product **pairs** that says whether two products are compatible. The goal is to be able to say that two new products are compatible, given the feature data for both products. It's the pairwise aspect that's throwing me off in the supervised learning framework. If I were just predicting a binary target based on color, pattern, etc., that's a basic setup with several algorithms to choose from. But to me, the problem seems a little more complicated since the prediction is based on feature data from two products—and the prediction needs to be symmetric so that it is the same no matter which product comes "first." Thus, I'm trying to get a symmetric compatibility metric when I have (1) feature data on products and (2) target data on pairs of products.


title: How to test for outliers in an mlogit model in R

text: I am running a multinomial logistic regression using the `mlogit` package and `mlogit` function in R. Now I need to check for outliers for the model. Is there any approach or function in R for testing outliers in an `mlogit` model?


title: Model suggestion for a Cox regression with time dependent covariates

text: I´m modeling the effect of pregnancy on the outcome of a disease (dead-alive). Approx 40% of the patients did become pregnant after the time of diagnosis-but at different points in time. So far I´ve done KM plots showing a clear protective effect of pregnancy on survival and also a regular Cox model- however these have been modeled using only a dichotomised pregnancy variable and assuming the effect is present from the time of diagnosis which is clearly unrealistic since the median time to pregnancy is 4 years from diagnosis. What kind of model would absorb the effect of multiple pregnancies at different time points after diagnosis? Would it be correct to model the pregnancies interacting with time (which would require some serious data reconstruction-any automated software that could help with this?) or is there another preferred modeling strategy for these problems? Also what is the preferred plotting strategy for these problems?


title: Gibbs sampler on the precision (with a gamma prior) in a hierarchical Bayesian model doesn't converge

text: I am deriving a Gibbs sampler with a model similar to the model in this paper (a graphical model is shown in page 4). To put it simple, my question only concerns $w_i$ (a $K$-dimensional vector drawn from a normal distribution) and its precision $\gamma_w$, thus I will just make these 2 variables unknown and treat the rest as known (hopefully without losing the mathematical rigors): $w_i \sim \mathcal{N}(0, \sigma_w^2I_K)$ $\gamma_w = \frac{1}{\sigma_w^2} \sim \Gamma(\alpha_0, \beta_0)$, here $\alpha_0, \beta_0$ are hyperparameters. And the conditional probability is: $p(w_{ik}|-) \sim \mathcal{N}(\mu_{w_{ik}}, \sigma_{w_{ik}}^2)$, where $\sigma_{w_{ik}}^2 = (\gamma_w + T)^{-1}$. $T$ is actually a pretty complicated term but it's nonnegative. $p(\gamma_w|-) = \Gamma(\alpha_0 + \frac{1}{2}KN, \beta_0 + \frac{1}{2}\sum_{i=1}^N w_i^Tw_i)$ I code it up and generate some simulated data according to the graphical model. I set all the variables other than $w_i$ and $\gamma_w$ to their true value and keep unchanged. The hyperparameters are set to 0.001 which basically forces the prior of $\gamma_w$ to be 0 (I am actually not certain if this is a good choice). The sampler runs with $\gamma_w$ goes to infinity. The log- likelihood of samples goes to much larger than the log-likelihood of the true data. From the conditional distribution, this non-convergence is actually kind of "observable", as each $w_{ik}$ is drawn with a larger precision than the sample from $p(\gamma_w|-)$ (because of the nonnegative term). Then when sampling $\gamma_w$, we are using a smaller $\sum_{i=1}^N w_i^Tw_i$. Thus the samples of $\gamma_w$ get larger. This seems wrong, but I cannot figure out where the problem is. Any comments will be highly appreciated. Thanks! Updated: I think this may be generally true for a hierarchical normal model with unknown variance, then I found this paper which seems to partially agree with my observation. But this paper is a little beyond my crappy math background...


title: Is discriminant analysis supervised learning?

text: Is linear discriminant analysis, specifically Linear Programming Discriminant Analysis (LPDA), supervised learning? Can you provide a valid reference that states so if possible. My study supervisor and I have been disagreeing about it. I'm convinced linear discriminant analysis, whether Fisher LDA or LPDA, is supervised learning. Both techniques use a labelled set of objects to derive a function which can be used to predict class labels for unlabelled objects. My study supervisor does not agree, stating that nothing is "learned" when using discriminant analysis.


title: More details on bootstrap procedure to estimate confidence interval of sample SD

text: I originally asked the question in this post, where I explained that I only have 25 measures at my disposal. Hence, computing the sample SD $\sigma_x$ and assume it's the one from the distribution looks to be very optimistic and I wanted to compute a confidence interval. It was suggested to me to use the bootstrapping method, but I'm not sure how exactly to use it; I believe that I would have to proceed as follows:   * make some assumption over the underlying probability distribution behind my sample, for example, I assume that it is normal, and apply the parameter I estimated : $\mathcal{N}(s, \sigma_x)$   * Generate $n$ samples of size 25 (precisely the same size as my original sample)   * Compute the sample SD of each of the 25   * And then look at the (85%) which are the closest to the mean; they all will be in the 85% interval. I guess I got it right , but I **_have_** to make an assumption on the underlying distribution, which I would really like to avoid. Is there another way to go, or did I misunderstand somehting?


title: Using control variates & antithetic method with Monte Carlo

text: Supposing $g(x)=\sqrt[3]{x}$, I want to calculate the expected value of g, $E(\sqrt[3]{x})$, using Monte Carlo method, by generating $x_i$ from a Weibull distribution with parameters $(1,5)$. After that, I want to use the **control variates method** and the **antithetic method** in order to to reduce the variance of my estimator, which I found with the simple Monte Carlo. **And here is my problem, I do not know how to do these methods.** I would appreciate if someone could help me do that or give any tip/help.Thank you very much for your concern, in advance. ### What I have done so far Supposing $S$ is our estimator, then we know that $S=(\sum \limits_{i=1}^{N} g(x_i))/N$. Using Matlab, I found the expected value $S$ by generating 1000 random numbers from the `weibull(1,5)` distribution and calculate the sum. Here is my algorithm:               N=1000     sum=0;     for i=1:N;       X = wblrnd(1,5);         res(i)=X.^(1/3);        sum=sum+res(i);     end     S=sum/N


title: Meaning of standard deviation of the mean difference

text: I'm trying to understand the meaning of stating the standard deviation (SD) of the mean difference (MD) [or otherwise called the absolute mean difference]. This is for a paper I'm writing up where other examples also quote the SD of the MD as part of the analysis summary. While I understand what the mean difference represents, i.e., being a "measure of statistical dispersion", I'm not sure of the utility of reporting the standard deviation of this value. I can kind of conceptualise what it means, but again the issue is the utility. For reference, the analysis is of repeated measures data for equivalence. I will be reporting the correlation coefficient and also the statistical probability for equivalence (with equivalence range) using a Two One Sided T-test (TOST). Thus my question regarding the utility of the standard deviation of the mean difference. Maybe I'm lacking in statistical knowledge and please feel free to point out a suitable reference for me to consult. Thanks in advance. **Update - example:**   _Machine A_ is the reference measurement system. _Machine B_ is the new measurement system. A certain number of real life measurements are made with each machine ( _N_ ) on a common object to understand how equivalent the machines are. Ideally they're identical - "of course" say the designers, and can even be expected - but this testing will be used to support the equivalence limit of +/- _x_ units. (Given in reality it's not possible to say _A_ is exactly equivalent to _B_.) So a paired measurement comparison. When reporting these results it has been common to report the mean difference and the standard deviation of the mean difference, plus the correlation coefficient (Pearson's) with its p value. **Update - comments:**   1\. This question does indeed refer to the standard deviation of the mean difference, not simply the difference.   2\. I well understand what the standard deviation means in terms of distributions. But my question is really "what is the meaning of the standard deviation of the mean difference in relation to the quantity being measured?"   3\. If it's as simple as a guide to the variability of the difference, then so be it.   4\. If the result was reported as A = B +/- _x_ units with a p < 0.05 then I understand that this is different to quoting a standard deviation type number, but I still get it intuitively. I don't get the SD of a MD in regard to the original measurement entity nearly as well.


title: Difference between "in-sample" and "pseudo out-of-sample" forecasts

text: Probably a very basic question to the forecasters around here, but I was wondering whether there is an explicit difference between _in-sample forecasts_ and _pseudo out-of-sample forecasts_. Both is meant in the context of evaluating and comparing forecasting models.


title: Quick reference for statistics

text: Some Quick References for statistics.   **conditions** :   1- offline, i.e., independent from internet   2- soft-copy, i.e., (workable on a computer) not a hard-copied book or etc   3- intermediate to advanced a comprehensive   4- not an encyclopedia with many links to that page or other term   5- ...   6- free or almost free   7- reliable and professional   8- general, i.e., not dependent to specific software or package   9- ...   10- well classified terms   11- some you add.   By the way, I don't think the list of conditions is too long. That is what I and you need.   **Updates:**   Find a prefect answer which is given bellow!


title: Change detection algorithm - likelihood ratio

text: > Consider a sequence of independent random variables $(y_k)_k$ with a > probability density $p_{\theta}(y)$ depending upon only one scalar > parameter. Before the unknown change time $t_0$, the parameter $\theta$ is > equal to $\theta_0$, and after the change it is equal to $\theta_1 \neq > \theta_0$. The problems are then to detect and estimate this change in the > parameter. _snip... snip..._ > The tools for reaching this goal are as follows. First, our description of > all the algorithms of this chapter is based on a concept that is very > important in mathematical statistics, namely the logarithm of the likelihood > ratio, deﬁned by $$s(y) = \ln\frac{p_{\theta_1}(y)}{p_{\theta_0}(y)}$$ and > referred to as the log-likelihood ratio. The key statistical property of > this ratio is as follows : Let $E_{\theta_0}$ and $E_{\theta_1}$ denote the > expectations of the random variables under the two distributions > $p_{\theta_0}$ and $p_{\theta_1}$, respectively. Then, > > $$E_{\theta_0}(s) < 0\space \text{and} \space E_{\theta_1}(s) > 0 $$ (This appears to be a quotation from Chand & Xiao, _Change-Point Monitoring for Secure In-Network Aggregation in Wireless Sensor Networks_ , 2007.) My question is HOW can we get the 2nd formula? What does $E_{\theta_0}(s)$ mean (I know $E$ is the expectation of $y$, but what about $(s)$ )? How can we get the > 0 or < 0?


Which is a better cost function for a random forest tree: Gini index or entropy? I am trying to implement random forest in Clojure.


title: How to store the standard errors with the lm() function in R?

text: All is in the title... I know how to store the estimates but I don't know how to store their standard errors... Thanks               > x <- runif(100)     > y <- 5 + 3 * x + rnorm(100, 0, 0.15)     > reg <- lm(y~x)     >      > summary(reg)          Call:     lm(formula = y ~ x)          Residuals:          Min       1Q   Median       3Q      Max      -0.32198 -0.12626  0.02584  0.10873  0.31411           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  5.00931    0.03087  162.25   <2e-16 ***     x            2.98162    0.05359   55.64   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.1499 on 98 degrees of freedom     Multiple R-squared: 0.9693,     Adjusted R-squared: 0.969      F-statistic:  3096 on 1 and 98 DF,  p-value: < 2.2e-16           >      > reg$coef     (Intercept)           x         5.009309    2.981617


title: Stopping rule for chi-squared discretization algorithm

text: I developed an algorithm that uses the chi-squared test to perform supervised discretization of a continuous variable. I described it in the paper "ChiD-A Chi-Squared Discretization Algorithm" published in the WUSS 2011 Proceedings available at http://www.wuss.org/proceedings11/ The stopping criterion is not very intelligent, and I would like to know if there are better ways of stopping the discretization process.


title: How to model the size of a trade area?

text: I would like to model the dimensions of the trade areas for a group of stores. I have geocoded sales data for each store, which I used to construct a convex hull using a GIS program. Essentially, I took the nearest 80 percent of sales for each store, and then snapped a rubber band around those points to get a trade area boundary for each store. Then I calculated the area of these regions, and modeled the radii as if the trade areas were circles, using the characteristics of the area right around the store as predictors. This gives terrible predictions since many of the areas are elliptical or have spike-like bulges. Does any one have any suggestions for how to model the shape and size of these regions in a better way?


title: Mapping discrete variables

text: Let's say I have a set of discrete, arbitrary variables sampled at continuous co-ordinates, e.g.               ==================================     = x       = y       = Value      =     ==================================     = 4.1     = 3.4     = A          =     = 4.6     = 7.9     = A          =     = 4.7     = 8.1     = A          =     = 2.1     = 5.0     = B          =     = 3.5     = 9.0     = B          =     = 7.0     = 1.1     = C          =     = 9.8     = 4.2     = C          =     ==================================      I would like to generate a map showing the distribution of the values across a range of co-ordinates (e.g. `0<x<=10, 0<y<=10`), where each value has an arbitrary colour assigned, and every point within those co-ordinates has a calculated interpolated value, where it is reasonable to do so. E.g. we give a colour to each of A, B and C. the point (5.0,5.0) is in the middle of the cluster of As so should be coloured as per A. the point (7.0,4.0) is between the As and Cs so the colour should be calculated via some analysis of distance to each of those surrounding points. Broadly I would expect the top-left quadrant to be a colour assigned to B, the bottom right quadrant to be a colour assigned to C, with a band of colour assigned to A between and the point at which the boundaries between the 3 colours occur calculated via some kind of nearest-neighbour algorithm. Are there any existing algorithms to map the kind of information I wish to display? It should take into account any inliers (if that's the right term).


title: HMM learning from video data?

text: I am having a problem understanding how to learn the parameters for the HMM from observed data. Let's say that my HMM model has one hidden variable for affect(emotion) with three values/states (anger, happiness, confusion) and each hidden variable at t has three observations (voice pitch & volume, facial expression). How do I go about collecting data & training? I was thinking that maybe I could record (get the three observables) from people who are in specific states meaning I will ask my actors to say a few lines using each of the three emotions, separately. Now, after that I don't get what I should do. Let's say I have 3 actors that means I am 9 samples for just training. If each sample is 3 minutes, do I train the model on the length of each of the samples? What if the samples are of slightly varying length? Also, how does the model learn P(anger(t+1)|happy(t)) from the data is separated by specific values (i.e. each sample is a specific emotion only)? Can someone please help me understand this? Thank you.


title: Determining a variable's contribution to the variation in another

text: I have data that look like this, with a tonnage measure and density measure for each carriage in a bunch of trains.               Train Carriage_No Tonnes  Density     A     1           105.5   2.12     A     2           104.9   2.28     A     3           101.2   2.30     A     4           108.7   2.41     B     1           112.3   2.51     B     2           109.7   2.34      etc.. **How do I determine density's contribution to the variation in tonnes? What is doing this called? Some`R` code would be most helpful!** Our goal is to fit more tonnes into each carriage, but the variation around the mean is such that we're overloading an unacceptable percentage of carriages. Therefore we need to tighten the variation to allow us to increase the mean (i.e. more total tonnes) without overloading. There are a few things we can do to improve tonnage variability, but controlling density isn't one of them (reacting to changes in density may be possible). So I'd like to understand how much the density accounts for the variation in tonnes.


title: Maximizing/minimizing product of output and response

text: I need to optimize based on an objective function that is non-standard, as far as I know. If the predictors are $X$, output of the model is $\hat y$, and response is $y$, my objective is essentially: $$j = \sum( \hat y \cdot y ). $$ Minimized it would just be $-j$. Does anyone know of an existing approach to optimization/regression that uses an objective like this? $\hat y\cdot y$ is equivalent to what is called the margin in machine learning, but I cannot find any work that discusses optimizing this directly. I could obviously just cram this into a global optimizer but I am interested in any work that is pre-existing.


title: With Random Forests, is it possible to provide an error term for each variable or each sample?

text: I'm wondering if there are any implementations of random forests that allow one to provide error terms for the inputs. An error term could be based on a per-variable basis, or on a per- sample/variable basis. If an error term was specified on a per variable basis, we would let the random forest algorithm know that "this variable $A$ has a sampling error of $\pm 5\%$, whereas this variable $B$ has a sampling error of $\pm 10\%$". If an error term was specified on a per sample/variable basis, we would specify the error terms on all individual samples. We would let the random forest know that "for sample 1, variable $A$ has a sampling error of $\pm 5\%$, and variable B has a sampling error of $\pm 8\%$, whereas sample 2, variable $A$ has a sampling error of $\pm 6\%$ and variable $B$ has a sampling error of $\pm 12\%$". I am particularly interested in implementations in R, C++, C# or Python.


title: Explaining changes caused by dummy variables quantitatively

text: Consider the following linear model \begin{equation} Y_{j}=\beta_{1}+\beta_{2}E_{j}+\beta_{3}B_{j}+\epsilon_{j} \tag{*} \end{equation} where $Y_{j}$ represents the natural logarithm of the annual salary of an employee currently, $E_{j}$ the number of finished years at an educational institution, and $B_{j}$ is the natural logarithm of the initial salary of the employee. I have regressed $(*)$ on a large data set, and compared the coefficients $\beta_{2}, \beta_{3}$, to another model that includes two additional dummy variables, $D_{1j}$ (gender) and $D_{2j}$ (minority). Although they have not changed drastically, there still is some change. Now, I have no trouble giving a qualitative explanation as to why this is the case. What I am looking for is this: > **Is there some calculation, perhaps a test, that will facilitate a > qualitative explanation of why the coefficients of $(*)$ may change when one > adds two dummy variables?**


title: Word probabilities in a Naive Bayes filter

text: While implementing a Naive Bayes filter, I stumbled across a problem with the calculation of the conditional probabilities $p(w|c)$ of a word $w \in \mathcal{W}$ given a class $c \in \mathcal{C}$. Following this answer on StackOverflow about classifying fruits $\mathcal{C} = \left\\{ \text{banana}, \text{orange}, \text{other} \right\\}$ having the features $\mathcal{W} = \left\\{\text{long}, \text{sweet}, \text{yellow}\right\\}$, this is what my problem is about: Assume we have the three observed sentences (or documents) $\mathcal{D} = \left\\{ D_{1}, \dotsc, D_{3}\right\\}$ in the class $c = \text{banana}$ as $ \begin{align} D_1 = \left\\{\text{long}, \text{sweet}, \text{yellow}\right\\} &\to \left\\{\text{long}, \text{sweet}, \text{yellow}\right\\} \\\ D_2 = \left\\{\text{long} \right\\} &\to \left\\{\text{long}, \neg \text{sweet}, \neg \text{yellow} \right\\} \\\ D_3 = \left\\{ \text{yellow}, \text{sweet} \right\\} &\to \left\\{ \neg \text{long}, \text{yellow}, \text{sweet} \right\\} \end{align} $ then what is the probability $p(w = \text{long}|c = \text{banana})$? According to the linked answer, the probability would be calculated as $ \begin{align} p(w = \text{long}|\text{banana}) &= \frac{\left|w\right|}{\left|\mathcal{D}\right|} = \frac{\left|w\right|}{\left|w\right| + \left|\neg w\right|} = \frac{2}{3} = 0.\overline{6} \end{align} $ because the word appears in two out of three sentences. Regarding the fact that I later on want to classify documents/sentences of **different lengths** (think mail classification), can't we also say, that the probabilities of the word appearing in the individual training sentences are: $ \begin{align} p(w = \text{long}|\text{D}_1) &= \frac{\left|w\right|}{\left|\text{D}_1\right|} = \frac{1}{3} = 0.\overline{3} \\\ p(w = \text{long}|\text{D}_2) &= \frac{\left|w\right|}{\left|\text{D}_2\right|} = \frac{1}{1} = 1.0 \\\ p(w = \text{long}|\text{D}_3) &= \frac{\left|w\right|}{\left|\text{D}_3\right|} = \frac{0}{2} = 0.0 \end{align} $ because the word appears $1$ out of $3$ times in $\text{D}_1$, etc. Then wouldn't the (well, expected) probability be $ \begin{align} p(\text{long}|\text{banana}) &= \frac{1}{N} \sum_{i=1}^{N} p(\text{long}|\text{D}_N) = \frac{1}{3}\left(\frac{1}{3} + \frac{1}{1} + \frac{0}{2}\right) = \frac{4}{9} = 0.\overline{4} \end{align} $ (?) which looks like the complement by accident. Or is the way to think that $ \begin{align} p(\text{long}|D_1) &= p(\text{long}|\text{long}, \text{sweet}, \text{yellow}) = p(\text{long}|\text{long}) = 1 \\\ p(\text{long}|D_2) &= p(\text{long}|\text{long}, \neg \text{sweet}, \neg \text{yellow}) = p(\text{long}|\text{long}) = 1 \\\ p(\text{long}|D_3) &= p(\text{long}|\neg \text{long}, \text{sweet}, \text{yellow}) = p(\text{long}|\neg \text{long}) = 0 \end{align} $ because $\text{long}$ is conditionally independent of all other features (following the Naive assumption)? But then, what about double-counting words (i.e. in mail classification)? The approach of counting over all words would yield $2/6 = 0.\overline{3}$ here, which is another different number. I assume that averaging combination is not the way to go, but how _does_ the combination of the per-sentence probabilities work right?


title: Why isn't the lasso selected variable even not significant?

text: I performed a lasso selection using `lars::lars` for a well normally distributed outcome using a pool of 86 predictors. Here is the plot of the output: ![enter image description here](http://i.stack.imgur.com/21WGb.png) The following covariance test confirmed 'clive' is the best selected:                 Variable Drop_in_covariance P.value     1    clive            13.4357  0.0000     2     pepl             0.2965  0.7436     3     sb5y             0.1028  0.9023     4    radio             0.0656  0.9365     5     roof             0.9746  0.3781      However, when I built a model then using `glm` with 'clive' and some other biological variables (they are all within the 86-variable pool for lasso), I got this results:               Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) -0.41025    0.05518   -7.43  4.2e-13 ***     clive        0.06499    0.04110    1.58   0.1144         cdaysbms2   -0.00798    0.01956   -0.41   0.6836         exp         -0.01194    0.00386   -3.10   0.0021 **      bodr        -0.00052    0.04004   -0.01   0.9896         medu         0.02642    0.01343    1.97   0.0496 *       inco         0.00815    0.00298    2.73   0.0065 **       Here you can see the best selected 'clive' isn't statistically significant while other some unselected are even significant. It is true when the model only includes 'clive', it becomes significant:               Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  -0.4269     0.0304  -14.02  < 2e-16 ***     clive         0.0743     0.0133    5.58  3.9e-08 ***      But I'm still confused how to interpret it. I would think if 'clive' is the best one from a pool, then its p.value should be the lowest one even when have other variables from the pool in the model. Would somebody offer an explanation?


title: Does testing for assumptions affect type I error?

text: I just performed simple simulation. Made two "populations" with different means and the same variance. Since I prepared them I know that they: are normal, differs in location and both have the same scale. Then I wrote simple script in R which 1000 times draws samples from the both populations, test them for normality (Shapiro-Wilk), for equal variances (F) and for difference (t). If the assumptions of t-test are met I run t-test. Well, in fact I run all of them one by one, but combine the results with logical AND operator. So, if **IsNormal** (A) AND **IsNormal** (B) AND **VarsAreEqual** (A, B) AND **MeansAreDifferent** (A, B), then I count the result. Finally I divide number of TRUE values and divide it by the number of iterations. So, every test of assumption may fail. Shapiro-Wilk may falsely (they **ARE** normal) reject H0 as well as F test (they **HAVE** equal variances). And I don't get values near 0.05! So, my question is: is this an example of the "multiple testing" phenomena (like multiple comparisons)? If so, what about testing for assumptions? Such testing is crucial as if the assumptions are not met, the quantiles of sampling distribution of a given test may be completely different than expected and so given results may be completely useless. But testing for the assumptions clearly changes the alpha. Am I correct? Perhaps we should use corrections for alpha, like Bonferroni?               data<- data.frame(populA=rnorm(1000, mean=1), populB=rnorm(1000, mean=3))          iternum <- 1000     alpha   <- 0.05     nsample <- 50     results <- array(FALSE, iternum);          for(i in 1:iternum) {         sampA <- sample(data$populA, 30);         sampB <- sample(data$populB, 30);              result_t <- t.test(sampA, sampB)$p.value;         result_F <- var.test(sampA, sampB)$p.value;         result_normA <- shapiro.test(sampA)$p.value;         result_normB <- shapiro.test(sampB)$p.value;              results[i] <- (result_normA >= alpha) & (result_normB >= alpha) &                        (result_F >= alpha) & (result_t < alpha)     }     > 1-(mean(results))          [1] 0.158


title: Plotting multiple densities on the same graph in R

text: The R code below produces three CDFs using three different approaches. How can I plot these CDFs on one graph in R?               #Method 1     z <- rnorm(n=100000, m=0, sd=1)      z1 <- rnorm(n=100000, m=0, sd=1)      z2 <- rnorm(n=100000, m=0, sd=1)      z_app <- (z1/sqrt(3)) + z*(1+(z2/sqrt(2*3)))          #CDF     F_z_app <- ecdf(z_app)     plot(F_z_app, xlim=c(0,5),ylim=c(0.5,1))               #Method 2     v3 <- rt(100000, 1)/sqrt(3-2)      w3 <- rchisq(100000,2)      z3 <- rnorm(n=100000, m=0, sd=1)      z_eff_3 <- v3 + z3 * sqrt((3*(1+v3*v3))/w3)           #CDF     F_z_eff_3 <- ecdf(z_eff_3)     plot(F_z_eff_3, xlim=c(0,5),ylim=c(0.5,1))                    # Method 3     z <- rnorm(n=100000, m=0, sd=1)          #CDF     F_z <- ecdf(z)     plot(F_z, xlim=c(0,5),ylim=c(0.5,1))


title: How to improve neural network sensitivity with a lopsided binary outcome?

text: I'm working on predicting (not explaining) a 0/1 outcome that generally has only about 10% "1"s (I'm not at liberty to name the variables). N ~40,000. Logistic regression proved unsatisfactory, both when using about 5-10 main effects and after I built in several interaction terms suggested by CHAID procedures. Sensitivity was ultimately only about 25%. I then turned to neural networks (radial basis function networks, in SPSS). I was pretty shocked to see the program fail to classify any cases as "1"s. That is, sensitivity was zero. **First question** : Is this a common or understandable NN result under these conditions? Next I tried randomly excluding a large number of "0" cases in the training set, bringing the fraction of "1"s up to about 40%. Now the program was able to correctly identify a halfway-decent number of cases in the training set, with sensitivity around 30%, but that dropped to 20% when the solution was applied to the test set, which once again consisted of only about 10% "1"s. **2nd question** : How would you get around this problem?


title: Kernel SVM in primal training with Stochastic Gradient Descent

text: In short: I am currently reading Online Learning with Kernels (http://books.nips.cc/papers/files/nips14/AA33.pdf) for fun and I can't figure out how he got to equation 8 from equations 6 and 7. The idea is: We want to minimize a risk function $R_{stoch}[f,t]:=c(x_t,y_t,f(x_t))+\lambda\Omega[f]$. If we want apply the representer theorem on f, writing it as $f(x)=\sum\alpha_i k(x,x_i)$, how can we get to the STOCHASTIC gradient descent update? Say we take the soft margin loss for SVMs. It would be easy to take the gradient w.r.t. to f and loss (well sub-gradient for loss) and do gradient descent. But for online learning with stochastic gradient descent, I'm kinda lost. Thank you! Please do not hesitate to ask further details. Any help would be greatly appreciated.


title: Categorical and dummy variables in multiple factor analysis with FactoMineR

text: I'm interested in using multiple factor analysis to analyze data collected from ceramic sherds. The data is pretty complicated, containing both continuous variables (vessel rim diameter, vessel wall thickness) and grouped categorical variables (presence or absence of base- and over- glazes on the vessel interior and exterior). After scaling and centering the continuous variables, my data looks something like this:               row.names   RimDiameter  SherdThickness2cmBelowRim  BGIBrown    BGIclear...     23          -1.05924698  0.1160736                  0           1...     56          0.15586427   -0.7386500                 1           0...      Using FactoMineR, I tried running MFA as follows:               sherdMFA = MFA(cleanRimGlaze[,-1], group=c(2,17,13,1,1),        type=c("s",rep("n",4)), ncp=5,         name.group=c("RimDiamThickness","BGI","BGE","LI", "LE"))      but it produces this error               Error in MFA(cleanRimGlaze[, -1], group = c(2, 17, 13, 1, 1), type = c("s",  :          not convenient group definition      ~~I can't figure out how to fix this error. Looking at the MFA example on the the FactoMineR website, it seems that MFA is built to handle categorical variables as factors, and converting my dummy variables to factor levels might solve the group definition problem. I have tried this (see my painfully slow learning in the comments), but MFA expects more than two factors per column, so dummy variables seem not to be the desired input. However, because in my data it's possible for multiple glazes to be present, it seems weird to have separate factors for "brown", "clear", and "clear and brown". Am I right in thinking that it's better to use dichotomous variables for presence/absence, and simply mark presence for both brown and clear? Is there a way to run MFA in FactoMineR with categorical variables as dummy variables? Or maybe even a more appropriate statistical analysis?~~ Fixed several errors as follows:   1. Ditched the [,-1] because it was causing the "not convenient group definition" error.   2. Converted qualitative columns to factors with two levels to fix the "factors are not defined in the qualitative groups" error.   3. Removed qualitative columns which contained only one of the two factors to fix the " Error in contrasts<-( _tmp_ , value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels " error that popped up. And now MFA runs. I'm still curious about people's opinions on the appropriateness of multiple factor analysis for such data. My sample (around 300 lustreware sherds from a much larger number excavated in the 1930's) contained 81 rim sherds. For these sherds, rim diameter was estimated, thickness of sherd 2 cm below rim was measured, and presence or absence of difference colored glazes was measured on the interior and exterior of the sherds (for example, opaque blue, transparent blue, white, clear, lustre, turquoise, no glaze visible). Because certain interior-exterior combinations seemed common, like white glaze on the interior and blue glaze on the exterior with lustre overglaze on both, it seemed reasonable to have both "Base Glaze Interior Blue" and "Base Glaze Exterior Blue" variables, as well as "Over Glaze Interior" and "Over Glaze Exterior" variables for lustre. This amounts to 20 variables and a row.names for the sherd identifiers. 20 variables to a sample size of 81 seemed to scrape by as acceptable for such analysis, although obviously a larger sample size would be desirable (I'm looking at Sample size and subject to item ratio in principal components analysis for this, but you might know more recent work it). There is a lot more ceramic material available, so any hypotheses developed based on this smaller sample could be tested with more data collection. I chose MFA because the data mixes continuous and categorical variables. As far as groups go, it seemed like a good idea to combine all of the presence/absence for the base glazes on the interior together, and those on the exterior together. By analogy, the presence of lustre overglaze on interior and exterior each were assigned a group, and finally the highly correlated estimated rim diameter and sherd thickness also got their own group (making 1 continuous group and 4 categorical groups). My first thought was PCA, but the categorical glaze data, while convertible to dummy variables, breaks the normal distribution assumption upon which PCA is predicated. I looked at Latent Class Analysis too, but that didn't seem suited to the continuous data of the rim diameter and sherd thickness variables. MFA seemed capable of handling both, and acknowledging the real groups present in the data collected. Many thanks to @dcarlson for his clarifying questions.


title: Finding criteria for a household financial budget falsification

text: I’m working on a financial problem about budget of households. Households in a state fill a form about their net budget in every year and our insurance company investigate their financial status and find the exact amount of their budget. The net budget can be positive or negative. I’m designing a system with neural network that we can find the households that Falsify their net budget, So my output is a binary [ 0 and 1] which 1 is false net budget and 0 is true net budget. Suppose that households’ net budget is P and our investigated outcome of net budget is R. As you know :               R>=P      Now we have this:               If P=R then  :  Criterion = 0 (true net budget)     If P>0 and R>0   :  Criterion = R-P/R      If P<0 and R<0   :  Criterion = (abs(P)-abs(R)) / abs(P)      If P<0 and R>0    :  ????     (A lot of households have this situation and I dint find any criterion for it)     If P<0 and R=0    :  ????     (Some of households have this situation and I dint find any criteria for it)     If P=0 and R>0    :  ????      (Some of households have this situation and I dint find any criterion for it)      After finding these criterion that must be between **0 and 1 [0 1]** , I set a threshold (like 0.5), If the criterion is higher than 0.5 household is falsifying the net budget (output = 1) and If it is lower than 0.5 we have a true budget and there isn’t any problem in budget declaring of households (output = 0). ** In range of [0 1], if falsification is higher (depends of above functions - P and R) the creation goes toward 1 and if the household is honest it goes toward 0. So every household (sample) has a criterion between 0 and 1. My main problem is finding a good criterion when **(P <0 , R>0)**, After that **(P <0 , R=0)** and finally **(P=0 , R>0)**. IF I have a problem in first, second and third criterion, feel free to say. I think the output of second situation is different from third because structures of criteria is different. I’m going to find a single criterion function If it is possible. **Example :** First conduction : P = 10 , R = 15 : Criteria : (15-10/15) = 33.3% ( falsification rate for this household ) Second conduction : P = -16, R = -10 : Criteria : (abs(-16)-abs(-10))/abs(-16) = 37.5% ( falsification rate for this household ) Now if for example P=-20 and R = 40, How we can calculate falsification rate? . This is main question. after that how we can calculate when (P<0 , R=0) and finally (P=0 , R>0)? **I check this criterion with this procedure:** for example in first criterion when P>0 and R>0 we have P = 5 and R = 10 so E = 50% . Now if P=fix = 5 and R=increase ( for example now R= 20) we have E = 75%. This is true conceptually in our problem. We must have more falsification when P=5 and R=20 Compared with P = 5 and R = 10 . + we must check it when R=fix and P=increase (or decrease). Thank you so much for your helps. PS. I Have 2000 samples with different values of P and R.


title: Feasible GLS of Random effect Panel data

text: When we define $\Theta$ for Feasible GLS in panel data, we have to define $\text{plim} \;\sigma$ in following picture. ![enter image description here](http://i.stack.imgur.com/HoAPv.png) How can I show that $\hat{\sigma}^{*2} \overset{p}{\longrightarrow} \frac{\sigma_\varepsilon^2}{T}+\sigma_u^2$?


title: Why does sign of coefficient of categorical variable changes in logistic regression when running a model on test sample

text: Hi this is my first post on this site, I am building a logistic regression model in SAS using proc logistic, my target is 0 and 1, I have 10-15 predictors 2 are continuous and rest are categorical, I divided total sample in two parts train(80%) test(20%), when I am checking the stability of coefficients on test sample I found that few categorical predictor changes their sign, also few turned insignificant on test data. I am using categorical variables as it is (not creating dummy) in class statement in proc logistics. so in short what would be reason for sign change on test data? Thanks, Abhishek


title: Delta Method for confidence intervals in online A/B testing

text: Scenario: I have two groups (treatment and control) that are matched to each other during a "pre" period. That is, controls were chosen such that their behavior in times before the treatment correlated well with the group that eventually was subjected to treatment. After treatment periods I want to compare the effect of treatment by finding the difference between the expected or organic behavior growth in the treatment group and the actual treatment behavior. The expected behavior is calculated by a linear regression between treatment and control groups during the "pre" period and assuming this functional relationship persists after the treatment period. $$ Y_{expected} = A X_{t>t_{0}} + B $$ Where $t_0$ is the treatment time and $X_{t>t_0}$ is the observed control group values after treatment. The effect is then found by $$ L = \frac{Y_{t>t_0} - Y_{expected}}{Y_{expected}} $$ My question is now: Is it necessary to use the Delta Method to find a confidence interval for $L$? I know $Y_{expected}$ is a functional transformation of the random observation $X_{t>t_0}$ but it is not clear to me why I can't simply find confidence bounds for $Y_{expected}$ by measuring the variation in the residuals $Y_{t<t_0} - Y_{pre-expected}$ during the pre-treatment time. Thanks in advance for any input or pointers.


title: Calculating difference between two means based on linear regression

text: I have a study which takes measurements on 30 subjects at 7 time points. The gradient of the best fit line is then calculated (call it say beta1), which is used to calculate a production rate (beta1*constant). The beta1 in each subject will have a standard error. I want to know the mean beta1 (mean production rate) for the whole group of 30, and the standard error of that measurement. Then, the study is repeated on all 30subjects several months later, and the production rates calculated again. The mean of this production rate is compared with the mean of the initial production rate. Can someone help me out as to the statistics? Many thanks.


title: The derivation of the Mahalanobis distance formula

text: I recently asked about the Mahalanobis distance and I got pretty good answers in this post: Bottom to top explanation of the Mahanalobis distance? I think I got the idea, but what I still felt missing was the derivation of the formula for the Mahalanobis distance. So my question is: "How does one derive the formula for Mahalanobis distance?" Why does the formula have the form: $$D(\textbf{x},\textbf{y})=\sqrt{ (\textbf{x}-\textbf{y})^TC^{-1}(\textbf{x}-\textbf{y})} $$ Could perhaps someone give analogous derivation as user @sjm.majewski gave on Principal component analysis on the link below: Making sense of principal component analysis, eigenvectors & eigenvalues UPDATE: From Wikipedia intuitive explanation was: "The Mahalanobis distance is simply the distance of the test point from the center of mass divided by the width of the ellipsoid in the direction of the test point." So is $C^{-1}$ the width of the ellipsoid in the direction of the test point? I mean this distance I can understand: $$\displaystyle\frac{\textbf{x}-\textbf{u}}{\sigma}$$ But this distance confuses me...:/ $$\sqrt{ (\textbf{x}-\textbf{y})^TC^{-1}(\textbf{x}-\textbf{y})} $$


title: How to test the effect of pre-post, treatment-control intervention with three age groups, two dependent variables and gender effects?

text: The hypothesis of my research project is “a certain training improves the ability in Subject A” . _**Sample details_** Groups: 1A (25), 1B (21), 2A (25), 2B(27), 3A(23), 3B (23). The number within brackets represent the number in each group. Control/experimental: 1A, 2A, 3A are experimental groups and 1B, 2B, 3B are control groups. Group labelled (1A, 1B); (2A,2B); (3A,3B) fall under 3 different age groups respectively. Data collected was as follows:   1. Pre-test and post-test scores for each group taken using subject A ability test    2. Subject A and Subject B achievement scores for each grade (5-6 scores taken at different intervals in an academic year of class unit tests given by the teacher).  _**Statistical analysis_** With limited knowledge in experimental statistics, pre-test and post-test results for all the groups were taken to test Subject A ability. I conducted a paired t-test for each group separately. Subject B unit test scores were also taken just to have a check on the effect of training on one other subject achievement, different from subject A. Questions : Apart from paired t-test what other statistical tool will help in line with my hypothesis?   1\. I wish to compare the Subject A ability scores between control and experiment groups   2\. Compare the Subject A ability scores on the 3 different age groups   3\. Compare the unit test scores conducted in class separately for subject A and Subject B   4\. Compare ability scores on the subject A for male and female in each group and overall PS: I will be using SPSS for statistical analysis.


title: Graphs in regression discontinuity design in "Stata" or "R"

text: Lee and Lemieux (p. 31, 2009) suggest the researcher to present the graphs while doing Regression discontinuity design analysis (RDD). They suggest the following procedure: > "...for some bandwidth $h$, and for some number of bins $K_0$ and $K_1$ to > the left and right of the cutoff value, respectively, the idea is to > construct bins ($b_k$,$b_{k+1}$], for $k = 1, . . . ,K = K_0$+$K_1$, where > $b_k = c−(K_0−k+1) \cdot h.$"               c=cutoff point or threshold value of assignment variable     h=bandwidth or window width.      ...then compare the mean outcomes just to the left and right of the cutoff point..." ..in all cases, we also show the ﬁtted values from a quartic regression model estimated separately on each side of the cutoff point...(p. 34 of the same paper) My question is how do we program that procedure in `Stata` or `R` for plotting the graphs of outcome variable against assignment variable (with confidence intervals) for the sharp RDD.. A sample example in `Stata` is mentioned here and here (replace rd with rd_obs) and a sample example in `R` is here. However, I think both of these didn't implement the step 1. Note, that both have the raw data along with the fitted lines in the plots. Sample graph without confidence variable [Lee and Lemieux,2009] ![enter image description here](http://i.stack.imgur.com/a2KPD.png) Thank you in advance.


title: How to validate a survey when using Likert scale?

text: I have developed a very simple survey to measure the attitudes of students toward science fairs. It makes use of the Likert Scale and Open Ended Questions. Do you have a recommendation for the type of statistics to use to validate this survey?


title: Recoded items as separate factors in factor analysis

text: Hello and thanks in advance for your help! I have conducted and EFA for a scale I made that includes 13 items. Of these items, 2 are recoded. Looking at the rotated component matrix I see that this scale has 2 factors, one with 11 items and one with the 2 recoded items. What does this mean and how do I decribe this? Does anyone know of a journal article or citation that deals with this sort of issue.


title: Quantify the difference between two samples

text: My problem is this: I have a need to do some analysis of student test results, and I am seeking some tools in order to do so. What I have are collections of scores for a number of different tests, and a number of different students. So each student may have zero or more scores, and each test may have zero or more results. For some test (let's say **Test A** ) I have 30 scores, distributed across the possible range (1 to 9). For a second collection point I have a similar result set; the scores are distributed differently, but still across the same range. So, I know I can calculate the mean and look at those. But that doesn't tell me a huge amount. Essentially I want to find (1) some way to compare one collection to another, and (2) a way to quantify it to a numerical value; in this way I could say "the amount improved between collections 1 and 2 was twice as much as that between collections 2 and 3" - or something similar. I've read a bit about effect size, but I'm not sure if it will suit my purpose. Additionally, each test conforms to a different range and scale; I have found an algorithm that works for me to 'normalize' the scores into a common base. Provided the algorithm is perfectly weighted (for argument's sake), would it be safe to use data from multiple different tests all together as one 'collection' of scores? To sum it all up, I suppose I am looking for any tools that would be useful in my endeavor to quantify and compare data sets of test scores.


title: What is model identifiability?

text: I know that with a model that is not identifiable the data can be said to be generated by multiple different assignments to the model parameters. I know that sometimes it's possible to constrain parameters so that all are identifiable, as in the example in Cassella & Berger 2nd ed, section 11.2. Given a particular model, how can I evaluate whether or not it's identifiable?


title: For ungrouped binary data, can the deviance be approximated by a chi-squared distribution?

text: I heard that "for ungrouped binary data, the deviance cannot be approximated by a $\chi^2$-distribution." Is it true? Why? Thanks in advance!


title: Interpretation of a variable in the regression analysis

text: I have a question regarding the interpretation of my regression analysis. Imagine the following situation: There is a very high (and significant) correlation between variable A and the dependent variable in the correlation matrix. In the regression analysis, I have multiple other independent variables and variable A is just one of them. In all of the cases Variable A has a lot of influence (high beta) and this is also significant. However there is an exception. Only when I put variable A in combination with a particular variable (Variable X) the influence of variabdrops and it is not significant anymore. So, only when variable X is not taken in the regression variable A is significant and has a lot of predictive power. When variable X is also taken, it is not significant. There is definitely not multicolliniearity involved. Does that mean that the influence of Variable A on the dependent variable works through variable X?


title: Prediction Interval for Neural Net With Hessian :: nnet in R

text: I would like to construct a confidence interval around prediction from a neural network, without resorting to bootstrapping - given the computational cost. Can I use the Hessian returned in this way to produce a 95% CI? 1) Can you take the negative inverse of the Hessian as the var/covar matrix? I read here that this depends on what is being maximized or minimized. Is this true and how can you know for certain? 2) Is this an accepted routine for producing a confidence interval around a prediction and if so, how would you do it?


title: What does 'more statistically significant' mean?

text: I think have a misconception about statistical significance. For example, I have two treatments A and B. Mean effect of A is -105 with 95%CI:(-168,-42) and for B -88 with 95%CI: (-146, -29). Which of the results is more statistically significant then? When can we conclude in general that result is more statistically significant based on p-values or confidence intervals? However I've read here that > A result does not become "more" statistically significant if the p value is > "a lot smaller" than the alpha value, as opposed to being simply "slightly > smaller". So is it correct to compare statistical significance in such a way? I'm very confused, so I appreciate any hints.


title: How can i add a weight correctly in a classification algorithm created by me of handwritten sign?

text: I have 150 images of handwritten signs (like cross, star, circle, tridents), i extracted from them the number of branched point, ent-point and "close" a s feaetures (Having said that features are not enough): if they are a close figure like a circle or an open figure like a cross. **Classification** :               Cross:   1 branched point and 4 end point, close = 0     Circle:  0 branched point and 0 end point, close = 1     Star:  5 branched point and 0 end point, close = 1        Trident: 1 branched point and 4 end point, close = 0     ...      I have 8 categories of images and for everyone i choose a range of branched points and a range of end points for classifying them in a certain category. So a star could be have not only 5 branched point and 0 end point like above, but a range of 3-8 branched point and a range of 0,4 end point. For classifying any new image empirically, i decided that: position: [cross, trident, unknown ...] with close = 0 (open figure):               with branched point 0  [0,0,100...]     with branched point 1  [100,0,0...]     with branched point 2  [50,50,0...]     ...     with branched point 4  [0,0,100,...]        ...      position: [circle, star, unknown ...] with close = 1 (close figure):               with branched point 0  [100, 0, 0...]     with branched point 1  [0,0, 100...]     ...     with branched point 4  [0,50, 50...] #a star with non precise feature extraction     with branched point 5  [0,100, 0...]     with branched point 6  [0,50, 50...] #a star with non precise feature extraction      position: [cross, trident, unknown ...] with close = 0 (open figure):               with end point 0 [0,0,0,...]     with end point 1 [0,0,100,...]     with end point 2 [0,0,100,...]     with end point 3 [33,33,33...]     with end point 4 [5050,0,...]     with end point 5 [33,33,33,...]      position [circle, star,unknown,... ] with close = 1 (close figure):               with end point 0 [50,50 ,0...]     with end point 1 [0,50,50...]     with end point 2 [0,50,50...]     ...      **After that, i chose to give a weight.** For doing that, i calculate the range of every category (for star: (1) branched point: 6-4 = 2, (2) end point 1-0 = 1. if the figure have perfectly the right feature (star case: 0 end points; i give the maximum weight (30%). If it has 1 end point 15% . If it has 2 end points 0% weight. I add the proper weight to [0,50,0,50,...] and subtract to the others for having always a 100% total. At the end i add Branched points probabilities and end points propabilities and divided for 2.               (Bpoint prob + Epoint prob ) / 2 for every         Example: i get an image with 4 branched point and 1 end point, and it is a close figure: the weight for a star with 1 end point is of 15% so: the weight for a star with 1 branched point is of 20% so:                branched point: [0,50,50...] -->    [0,65,35]      end point:      [0,50,50...] -->    [0,70,30]      calculation: [0,65,35] + [0,70,30] / 2 = [0, 67,5, 47,5] so i have got a probability of 67,5% to be a star and a probability of 47,5% to be unknown. 0 % for other categories. Unfortunately, in many cases i have negative percentages, or percentage over 100, how can i fix this?


title: Explanation of cubic spline interpolation

text: Can someone explain to me what a cubic spline is, and how we could use it to interpolate a function? I have searched on the internet but I would like a simple explanation.


title: Generating realizations of sampling distributions

text: I have a sample of $N$ individuals. I know (assume) that the distribution of species' abundances follows a particular form, e.g., negative binomial. A nested example of these kinds of distributions is below. (This figure assumes non-integer abundances are possible. Mine would not.) ![enter image description here](http://i.stack.imgur.com/vYG9C.jpg) The catch is that I do not know the number of species in my sample. I only know the total number of individuals. I want to calculate the number of species and their abundances given a particular distribution of species abundances. Abundances will have integer values. How can I generate a realization of species and their abundances from the assumed form (e.g., negative binomial) and the total number of individuals? Numeric solutions are fine. _Further thoughts:_ I think there's a very elementary integration step that I'm not seeing. I know what to do for a related problem in which I have a list of species and the abundance of each: For a given sample $N$, I can use a multinomial to draw realizations using relative abundances.


title: Compressed sensing: Optimization in $L_1$ norm and total variation with fourier coefficients

text: I'm reading the article Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information (Candes, Romberg and Tao, 2004). In this article they are talking about recovering the function $f$ whose fourier coefficients are known on some domain $\Omega$, by solving the following optimization problems: $$ \min ||g||_{TV} \space\space\space \text{s.t.} \space\space\space\hat{g}(w)=\hat{f}(w),w \in \Omega$$ and $$ \min ||g||_{L_1} \space\space\space \text{s.t.} \space\space\space\hat{g}(w)=\hat{f}(w),w \in \Omega$$ Can someone please give me a reference that suggests how to actually solve these optimization problems (that combines both $g$ and $\hat{g}$)? A relevant R package would also be nice.


title: How to test over-dispersion for a compositional data in R?

text: I would do an over-dispersion test for a compositional data set (I don't have the original count values) for choosing later an appropriate regression model. Here is an example of my data set:               X<-matrix(sample(1:100,100,T), 20,5)     X<-X/apply(X,1,sum)     colnames(X)<-paste("Comp", 1:5, sep="")     rownames(X)<-paste("Id", 1:20, sep="")      Does anyone know how to do this test in R?


title: How to show $E[|X|]= \sigma$ where $X \sim N(0, \sigma^2)$?

text: Let $X \sim N(0, \sigma^2)$. I want to show that $$E[|X|]= \sigma.$$


title: Why are the fitted probabilities for the linear probability model and the probit model identical?

text: I estimated a linear probability model (LPM) $P(y=1|x_1) = b_0 +b_1x_1 + u $ and a probit model $P(y=1|x_1) = \Phi(b_0 +b_1x_1 + u) $, where $\Phi()$ denotes the cumulative normal distribution. The regressor $x_1$ is the same binary variable. I observed that the predicted probabilities for both models $\hat{y} =\hat{b_0} +\hat{b_1}x_1 $ and $\hat{y} =\Phi(\hat{b_0} + \hat{b_1}x_1)$ are identical across both models. Why is this the case from a theoretical point of view?


title: Model assumptions of partial least squares (PLS)

text: I am trying to find information regarding the assumptions of PLS regression (single Y). I am especially interested in a comparison of the assumptions of PLS with regards to those of OLS regression. I have read/skimmed through a great deal of literature on the topic of PLS; papers by Wold (Svante and Herman), Abdi, and many others but haven't found a satisfactory source. The pub linked below does mention assumptions of PLS, but it only mentions that 1) Xs need not be independent, 2) the system is a function of a few underlying latent variables, 3) the system should exhibit homogeneity throughout the analytical process, and 4) measurement error in X is acceptable. There is no mention of any requirements of the observed data, or model residuals. Does anyone know of a source that addresses any of this? Considering underlying math is analogous to PCA (with goal of maximizing covariance between Y and X) is multivariate normality of (y, X) an assumption? Do model residuals need to exhibit homogeneity of variance? I also believe I read somewhere that the observations need not be independent; what does this mean in terms of repeated measure studies? Wold et al. (2001) PLS-regression: a basic tool of chemometrics. http://www.utdallas.edu/~herve/abdi-wireCS-PLS2010.pdf Thank you for any insight.


title: What's the methodology behind the most-difference-between-groups-tag-cloud?

text: What is the likely stats methodology used in this old OKCupid post?: http://www.economist.com/blogs/johnson/2010/10/sexuality_and_language And this: http://blog.okcupid.com/index.php/the-real-stuff-white-people-like/ Here's the problem setup: You have a corpus of users' self descriptions. Choose the common (issue #1!) words that maximize the difference (issue #2!) between the N groups, here N=4 groups, 2-by-2 by gender and sexuality. So, as I understand this, you're doing a regular NLP thing, constructing a sparse matrix that has frequencies (or some kind of normed frequencies, like tf-idf). Rows are the documents, columns are counts on each of the words used across all documents within that observation. Then you use those word- occurrence-counts to construct a basis on which a set of users, say heterosexual women, has a high score, but a relatively low score on the 3 other bases chosen at the same time. The words/dimensions with the highest loading on that basis /is/ the set of words output on the chart. Issues above that seem tricky to me:   1. What strikes people who have done this as the likely means of deciding on whether words count as "common"? Do you use thresholds? What's the appropriate way to choose most salient words? If you do something like tf-idf, pmi, this may tell you relative topicality of a document versus a corpus, but this seems to be not the same problem as the relative topicality of a set of documents versus the corpus.   2. What is the method for grouping? Since I know the number of groups, k-means seems appropriate, but how does it work mechanically on a giant sparse matrix like the word dictionaries? I more than a little hand-waved above in understanding that the word occurrence as dimensions behave like matrices I'm used to. If there are 10,000 words in the sparse matrix, aren't there an infinite number of bases? What is used to maximize the difference between groups, and is a maximum guaranteed to exist? Besides K-means, this would also naively look to me like Fisher's LDA in the sense that I'm asking for a linear combination of features (a set of words) that maximize the difference between groups.   3. bonus: If you're doing stemming to match the same word across parts of speech, what data management do you do to key the stemmed forms back to the pretty representation? (You can't invert a stem function---so f("jumping") = "jump", f("jumped") = "jump", but there is no 1-to-1 f^-1("jump") = "jumping".) Or in this case, perhaps a filter was applied at the start to isolate noun-phrases? Any advice or especially links to code in R or Python would be helpful


title: Multiple samples vs. bootstrapping

text: I'm somewhat new to understanding statistics meaningfully, and the bootstrap in particular. Let's say you're doing something (polling, running experiment, whatever), and can afford to get 1000 samples. From what I understand, before the bootstrap method existed (> 25-30 years ago), you'd be better off with, say, 10 sets of 100 samples each and calculating your statistics based on sampling distributions. (As an aside, is there a good way to calculate the ratio of sets/samples? I'm know that you get a smaller standard error when n is larger, but don't know if there is a "sweet spot", and from my rather limited understanding, it's preferable to have n > 30 to avoid messing around with t statistics). But now that bootstrapping exists, my understanding is that you get one set of all 1000 samples (since this will be a better proxy for real poplulation than 100 would be), and do the bootstrapping resampling magic based on that. Is this correct? Or is best practice some hybrid of the two? Or does it ultimately depend on what the actual activity is? Thanks!


title: Johansen cointegration test with stationary variables

text: Journals and books about the Johansen cointegration test suggest that only I(1) variables should be put into a cointegration test. However many journals modelling economic growth and inflation use two I(0) variables in the cointegration test. Is this an acceptable procedure? What is then the interpretation of the residual if two I(0) variables are used in the cointegrating regression?


title: Get a representative vector from a large set, and compare it with samples

text: We're a small team of programmers and we're trying to solve a little problem, but we think we need some advices from professional mathematicians. We want to know if a picture of a card is an Identity Card or not, so we've implemented this algorithm (I've simplified a lot so we can be focused in the maths problems).   1. Take a sample of a lot of Id Cards and make a color map of each one of them. This map is a 360 dimensional vector of decimal values (An array of 360 elements).    2. With the previous data, we calculate one representative vector. Actually, we make this representative vector calculating the arithmetic mean for each dimension.    3. We take the image that we want to identify, calculate the color map vector and compare it with the representative vector. (We're using the Pearson correlation coefficient)    4. If the coefficient is near 1, our image is an Id Card.  We're obtaining reasonable results, but this are our questions:   1. Are there any method better than the arithmetic mean to construct the representative vector?    2. Are there any method better than then Pearson correlation coefficient to compare both vectors? **Edited** : Our plan is to have a representative vector of a lot of documents: IdCards, Driver licenses, passports, Residence cards, etc, and try to identify in wich category the card that we're comparing fits better. Each dimension in the vector (that we call color map) represents the percentage of pixels of the image that have that color. The sum of all the values of the vector is allways 100% The main reason because we're using Pearson it's because we thought that with it, the result is not affected by the number of pixels of the image, and the empiric results seems to confirms that fact. (we obtain the exactly the same results using percentages or simply using the count of colors) but we've to admit that we've a limited knowledge of statistics, so we could have take a wrong decision (in this decision, or in others).


title: Flying Bomber aircraft through SAM sites - combining normal distributions

text: I have been puzzling over this for days but I don't think this was covered at school We simultaneously fly a known number of bomber aircraft, K, through three sequential batteries of surface to air missile (SAM) batteries, A B C. Assume the planes are incapable of destroying the SAM batteries. Assume that the SAM batteries are sufficiently spaced out such that only one battery A,B or C can engage the planes at once. I.e. the enemy gets three separate goes, once with each battery to shoot at the planes. E.g if we start with 10 planes, battery A shoots down 3, B then gets a go at the remaining 7 it shoots down 2 etc etc... Assume each SAM battery shoots down a random number of planes. This random number is drawn from a normal distribution unique to each battery i.e. Mu_A, Var_A, Mu_B, Var_B etc. How would one combine these three distributions into a single distribution from which could be drawn the number of planes that make it through all three layers? I need the parameters i.e. Mu_Combined, Var_Combined of the combined distribution   * Is it as simple as adding the distributions together? (This seems a little naive to me)   * I found this question: Combining normal distributions \- Is it simply an average? Thanks very much EDIT: Use of the normal distribution is essential. In the model I am using the number of planes shot down is not particularly dependent on the number there are to shoot as it assumed each battery gets one shot and can only shoot some integer number of missiles. In addition the model uses a continuous distribution with the same mu and var as the binomial. Apparently it is a well established approximation whose accuracy can be improved using a continuity correction.(applied later)


title: Interpretation of Expression in ordinal regression

text: Suppose we are looking at an ordinal variable with $4$ categories. So there are three threshold coefficients $b_1, \dots b_3$ and one probit slope $b_4$. What is the interpretation of the following expressions: $$\frac{b_{1}}{b_{4}^{2}}, \frac{b_{2}}{b_{4}^{2}}, \frac{b_{3}}{b_{4}^2}, \frac{1}{b_{4}}$$ assuming a multivariate normal distribution?


title: How do I do basic hypothesis testing in a Bayesian framework?

text: This seems like it should be a widely-answered, but I've spent a lot of time searching and reading textbooks and have not been able to figure it out. Let's say I give my exciting new treatment pill to a random sample of N people and a placebo pill to a random sample of M people. A of the N people die and B of the M people die. I'd like to know the probability that more people die taking the placebo than the treatment. Currently I do:               > data <- matrix(c(A,B,N-A,M-B), nrow=2)     > chisq.test(data, correct=FALSE)      and look at the p-value (which I realize isn't actually the answer to my question, but that's frequentists for you). And if I want to know the interval there's a 95% chance the actual value is in, I just do:               > pbar = A/N; se = sqrt(pbar * (1-pbar)/N); e = qnorm(.975)*se; pbar + c(-e, e)      But I want the Bayesian equivalents. Now I know I can compute a Beta distribution with a uniform prior and then use Monte Carlo sampling to figure out the answers, but that seems hackish and inefficient. What's the more direct way to compute the answer?


title: Disagreement between normality tests and histogram graphs

text: My data consist of compaction measurements from 3 different cell types (X,Y, and Z). My goal is to know whether there are "significant" differences between these measurements, so I have tested for:   1. Whether my samples are normally distributed     * using the Shapiro–Wilk test     * using the Jarque-Bera test     * plotting qqnorm graphs     * plotting histograms    2. Whether the samples come from the same distribution     * using two-sample Kolmogorov–Smirnov test (K–S test) and comparing X vs Y, X vs Z, and Y vs Z     * using Kruskal–Wallis comparing X, Y, and Z together My data consist of 232 measurements for X, 284 for Y, and 124 for Z. The Shapiro-Wilk and Jarque-Bera tests in R always give me p<0.05, which I accept as not being normally distributed. However, when I plot histograms I get a normal-like distribution. ![http://i49.tinypic.com/x5a91.jpg](http://i.stack.imgur.com/Wi6GJ.png) The qqnorm plots also don't look that skewed, but maybe this is just my inexperience in interpreting qqnorm graphs (this is my first time making them). ![http://i49.tinypic.com/2nqqohl.jpg](http://i.stack.imgur.com/6SVjG.png) ![http://i47.tinypic.com/207c4.png](http://i.stack.imgur.com/NWOae.png) ![http://i48.tinypic.com/o01ptv.png](http://i.stack.imgur.com/aMUyU.png) Because of the supposedly non-normal distribution, I compared my data using KS test and Kruskal-Wallis, which always give me the result that my population Z is drawn from a different distribution compared to X and Y. However, I do not know if this is true, as R always reports for my two-sample Kolmogorov–Smirnov tests:               Warning message:     In ks.test(dataX, dataY) : cannot compute correct p-values with ties          Warning message:     In ks.test(dataX, dataZ) : cannot compute correct p-values with ties          Warning message:     In ks.test(dataY, dataZ) : cannot compute correct p-values with ties      probably because the samples have different sizes. I'd like to know what you think about it, and whether I should consider using more parametric tests rather than the non-parametric ones I've used, or whether the tests I've used are valid regardless of the normality of the data. Also, my measurements seem to differ very little among themselves, for example:               ---Summary stats for WT cells     Min.    1st Qu. Median  Mean    3rd Qu. Max.      0.1450  0.3720  0.5000  0.5598  0.7102  1.9290           ---Summary stats for Df cells     Min.    1st Qu. Median  Mean    3rd Qu. Max.     0.0550  0.4030  0.5445  0.5857  0.7210  1.5350           ---Summary stats for Dp cells     Min.    1st Qu. Median  Mean    3rd Qu. Max.     0.0670  0.4790  0.6255  0.6782  0.7897  2.0160      Here's a boxplot of the data: ![http://i49.tinypic.com/6qvgio.png](http://i.stack.imgur.com/kHWJV.png) Green=X, Blue=Y, Red=Z So I'm unsure about the conclusions I may derive from them.


title: What is "symmetric property" for stationary distribution

text: I have the one step transition matrix $$\pmatrix{0 & \alpha & 0 & \beta \\\ \alpha & 0 & \beta & 0 \\\ 0 & \beta & 0 & \alpha \\\ \beta & 0 & \alpha & 0 \\\\}$$ I want to work out the stationary distribuion. So I end up with $$ \pi_1 = \alpha \pi_2 + \beta \pi_4$$ $$ \pi_2 = \alpha \pi_1 + \beta \pi_3$$ $$ \pi_3 = \beta \pi_2 + \alpha \pi_4$$ $$ \pi_4 = \beta \pi_1 + \alpha \pi_3$$ and it says that by noticing the symmetric property of all four states, we can deduce that all the $\pi$'s are equal to each other. What is the symmetric property? Is it because we can see they are all in a similar form to $A = \alpha A + \beta B$ and $B = \alpha B + \beta A$ and so we get that they're symmetric?


title: equivalence testing for a chi-squared test

text: I'm familiar with the concept of equivalence testing for testing whether two means are the same Is there an _equivalent_ equivalence test as the counterpart to the chi- squared test for contingency tables (specifically demonstrating that two or more variables are independent). As a starting point, taking the two one-sided t-test as a guide, I imagine performing two one-sided chi-squared test and adjusting the chi-squared value +/- the meaningful difference range I have previously chosen. two problems arise 1) whether a one-sided chi-squared test is the correct way to think about things given that it is already one sided. 2) choice of my meaningful difference range. For the means example the range is in meaningful units (that of the original variable). I'm not sure how to go about choosing the range in the units of the chi-squared value. any thoughts are welcome, thanks


title: Confidence intervals based on asymptotic normality

text: I wrote a report about `confidence intervals based on asymptotic normality`. However, my supervisor said that the definition in the second paragraph is wrong, but I do not know what is wrong. I have attached the report here. Can you please take a look and let me know what is wrong? ![enter image description here](http://i.stack.imgur.com/vJ5BT.png)


title: The gradient of a bivariate probit model

text: I am working on a bivariate probit model and I want to calculate the gradient. In order to do that I have to take the devirative w.r.t. $\beta$ in the following expression: $$\int_{-\infty}^{x_{g1}\beta} \int_{-\infty}^{x_{g2}\beta}\phi_2(z_1, z_2, \rho)dz_2dz_1$$ Could somebody give me a hint how to do it? I tried to use differentiation under the integral sign, however, I could not get the required result.


title: Location test of the ratios of random variables

text: I've been struggling with this problem for a while now, with not much progress. Any suggestions would be helpful! Suppose I have $X_1, .., X_5, Y_1, ..., Y_5, X'_1, ..., X'_5, Y'_1, ..., Y'_5$, i.e., I have 4 groups of 5 random variables each. The r.v.s in each group are i.i.d. and follow the unknown distributions $X, Y, X', Y'$, respectively. I would like to test the null hypothesis $H_0: \mathbb{E}\dfrac{X}{Y} = \mathbb{E}\dfrac{X'}{Y'}$. I've been fairly stuck on how to do this. If it's any easier, I'm also fine with testing the null hypothesis $H'_0: \dfrac{\mathbb{E}X}{\mathbb{E}Y} = \dfrac{\mathbb{E}X'}{\mathbb{E}Y'}$. Here's the progress I've made so far:   1. I think that non-parametric tests don't give me enough power here, due to the small number of samples, so I'm trying to go with a parametric approach.   2. Based on domain-specific knowledge, I know that each of $X, Y, X', Y'$ can be modeled as a negative binomial distribution (I'm looking at count data). I've already been able to fit MLEs of the mean and shape parameters for each distribution.   3. I'm pretty sure that I want to use some simulation-based method to solve this problem, but I haven't figured out what is the right quantity to simulate.   4. If it helps, I have about 30,000 sets of these 20 observations, of which I believe ~80% would not reject the null hypothesis. I tried to construct an empirical null distribution (using the ratio of ratios as the test statistic) from these 80%, but given the need for FDR correction, my p-value threshold becomes very low and this empirical test has not given me enough power. Any ideas, pointers, etc. would be greatly appreciated! Thank you in advance!


title: Correlation between a categorical predictor and a continuous outcome variable

text: How do I perform correlation between a categorical antecedent variable and a continuous outcome variable? Like for example, correlating EACH attachment style with another variable, i.e, social support?


title: Modeling a non-stationary bounded series

text: I'm trying to model a time series variable that represents a percentage, strictly bounded between 0 and 1, that is also non-stationary about the mean. **Is there a model form that is able to account for non-stationary behavior while producing forecasts within the [0,1] interval?** I came across the following work while researching the topic which I believe could be extended to answer this particular question (just a hunch based on a quick reading): > Park, Joon Y., and Peter CB Phillips. "Nonlinear regressions with integrated > time series." Econometrica 69.1 (2001): 117-161. > > Park, Joon Y., and Peter CB Phillips. "Nonstationary binary choice." > Econometrica 68.5 (2000): 1249-1280. > > Hu, Ling, and Robert de Jong. Nonstationary Censored Regression. mimeo, Ohio > State Uni, 2006. Unfortunately, I'm unaware of any implementation (preferably in `R`) of the ML estimators discussed. * * * Per Alecos' requests, below is a plot of the series that led to this question ![graph](http://imgur.com/uQPNeEG.jpeg) and the de-trended series ![detrend](http://imgur.com/W6PYAiK.jpeg)


title: Analysis of Deviance in R - Which test?

text: I have two generalized linear models `mod1` and `mod2`.               mod1 <- glm(Score ~ Height + Gene, data=mydata, family=binomial)          mod2 <- glm(Score ~ Height * Gene, data=mydata, family=binomial)      I want to perform an analysis of deviance to test the significance of the interaction term. At first I did `anova(mod1,mod2)`, and I used the function `1 - pchisq()` to obtain a p-value for the deviance result I got from the anova table. I also did another test: `anova(mod2, test="Chisq")`. This gave a table for all the terms added sequentially first to last, and I obtained a very different p-value value for the interaction term `Height:Gene`... **Which one should I use?**


title: Orwin's fail safe N test

text: In Orwin's fail safe N test how to decide the values of criterion for a trivial log odd's ratio and mean log odds ratio in missing studies. I am a medical doctor. Please tell me in simple english. The data are               1. Classic fail-safe N      Z-value for observed studies                                    27.97543     P-value for observed studies                                    0.00000     Alpha                                                          0.05000     Tails                                                           2.00000      Z for alpha                                                    1.95996     Number of observed studies                                    5.00000     Number of missing studies that wouldbring p-value to >alpha  1014.0000          2. Orwin's fail-safe N     Odds ration in observed studies                              5.7339      Criterian for a ‘trivial’ odds ratio                           ?     Mean odds ratio in missing studies                           ?


title: How to make inferences on group SD and and the SD of the group SD in a hierarchical Bayesian model?

text: The hierarchical model specified below is quite "standard" and easy to implement in for example JAGS/BUGS. It has an hierarchical gamma prior on the subject's precisions ($\tau_j$) which in turn has gamma priors on its mean ($m$) and SD ($d$). After having fitted this model I have access to estimates of the mean group precision ($m$) and SD of the group precision ($d$) but I rather make inferences on the natural scale. That is, the mean group SD and the SD of the group SD. **How can I make inferences on the mean group SD and the SD of the group SD from such a model? Is it at all possible?** Here is a Kruschke style diagram of the model taken from the book Doing Bayesian Data Analysis ![enter image description here](http://i.stack.imgur.com/WITDz.png) And here is the corresponding BUGS code:               model {       for( i in 1 : Ndata ) {         y[i]  ̃ dnorm( mu[ subj[i] ] , tau[ subj[i] ] )       }       for ( j in 1 : Nsubj ) {         mu[j]  ̃ dnorm( muG , tauG )         tau[j]  ̃ dgamma( sG , rG )       }       muG  ̃ dnorm( 2.3 , 0.1 )       tauG  ̃ dgamma( 1 , .5 )       sG <- pow(m,2) / pow(d,2)       rG <- m / pow(d,2)       m  ̃ dgamma( 1 , .25 )       d  ̃ dgamma( 1 , .5 )     }      One solution would perhaps be to not put the hyperpriors on the precision ($\tau_j$) but rather on the SD ($\sigma_j$) though I'm not at all sure what priors I should use then.


title: The bias of Zellner estimators in dynamic SUR models

text: I have been playing around with a seemingly unrelated regression (SUR) estimation. However, for dynamic SUR models it is known that -- analogous to the ARIMA case -- an OLS/GLS estimate is biased. For example this article provides a correction. So here's my question: Is there a _R_ package or some other implementation that provides this correction? Thanks in advance!


title: Is z-score meaningful in classification or regression?

text: I am a beginner in statistics, therefore I hope I can state my problem in a correct manner. I have a some instances or samples and I can collect below statistical parameters for **classification and regression** problem:   * Sample Size   * Minimum value   * Maximum value   * Standard deviation   * Variance   * Mean And, I want to use **z-scores** to compare or classify samples, my question is: does using z-score **make sense** or what can I use instead of z-score to obtain meaningful classification parameter?


title: Sum of Binomial and Poisson random variables

text: If we have two independent random variables $X_1 \sim \mathrm{Binom}(n,p)$ and $X_2 \sim \mathrm{Pois}(\lambda)$, what is the probability mass function of $X_1 + X_2$? **NB** This is not homework for me.


title: Sequential Monte Carlo for hierarchical models

text: Does anybody know, can Sequential Monte Carlo be applied for multi-dimensional problems i.e. simulating more than 1 distribution like in hierarchical models? Maybe you know some following literature


title: How to create a density plot for data from multiple years with each year represented by a different colour?

text: I have data on a variable spanning multiple years. I want to see how the distribution has evolved over the years.   * Is there an easy way to produce a density plot for each year in single plot with different colours for each year?