title: Estimate density from neighbor distances

text: Assuming I have a data set of known size, and there is one object that I want to test for being in an approximately uniform distributed region of the data set. For the query object, I know the $k$ nearest neighbors and the associated distances. If the data set is uniform there, the distances should be $c\cdot\chi(d)$ distributed, I believe (with $c$ a constant, representing some kind of data set diameter). So I have a spherical subset of the total data set that contains $k$ objects and has a known radius of $d_k$, the distance to the $k$th object, or equally the maximum distance within the $k$ nearest neighbors. Note that the $d_k$s - if I would be looking at more than one such subset - likely are not $\chi(d)$ distributed, as they are derived from the distances. Is there any assumption that I can make on the distribution of the $d_k$ for a dataset? So what is the best way to test if the object is in a somewhat uniformly distributed area of the data set (should I try a Kolmogorov-Smirnov goodness of fit test to $\chi(d)$ of the known distances? But I'm probably not seeing the complete distribution!), and what is the best guess at the density? My key idea is that the object might be an outlier if is in a low-density and approximately uniform area of the data set. Methods such as LOF fail here, because the outlier is too far out. LOF works good for detecting outliers nearby the actual data, but not within a uniform sparse area.


Repeatedly rolling a six sided die four times and summing the highest three results gives you a distribution with what mean and standard deviation? I've only taken AP statistics, but I would like to learn how to do this.


title: Test if one predictor's lack of a relationship with the DV is because of a negative relationship with another predictor correlated with the DV

text: I recently collected some survey data where people made a series of judgments about others' actions and their willingness to make some choice about that action. Looking at the zero-order correlations, I see that some of the judgments are significantly correlated with the choices of the participant, whereas some others are not. I am concerned that actual relationships may be hiding between the uncorrelated judgments and choices. Let me give a generalized example. Imagine Predictor A and B are correlated with DV X, but Predictor C is not correlated with DV X. I was thinking that there could actually be a relationship between C and X, but that C may be negatively correlated with A or B resulting in a zero-order correlation between X and C coming out as a wash (and therefore not looking significant). I imagine I should be able to test for this using regression, but I am not quite sure on the exact procedure. Would simply putting all A, B, and C predicting X in a regression showing that C remains insignificant be sufficient so show that it is not due to this inverse correlation? Thank you very much for your suggestions.


title: Deriving the conditional distributions of a multivariate normal distribution

text: We have a multivariate normal vector ${\boldsymbol Y} \sim \mathcal{N}(\boldsymbol\mu, \Sigma)$. Consider partitioning $\boldsymbol\mu$ and ${\boldsymbol Y}$ into $$\boldsymbol\mu = \begin{bmatrix} \boldsymbol\mu_1 \\ \boldsymbol\mu_2 \end{bmatrix} $$ $${\boldsymbol Y}=\begin{bmatrix}{\boldsymbol y}_1 \\ {\boldsymbol y}_2 \end{bmatrix}$$ with a similar partition of $\Sigma$ into $$ \begin{bmatrix} \Sigma_{11} & \Sigma_{12}\\ \Sigma_{21} & \Sigma_{22} \end{bmatrix} $$ Then, $({\boldsymbol y}_1|{\boldsymbol y}_2={\boldsymbol a})$, the conditional distribution of the first partition given the second, is $\mathcal{N}(\overline{\boldsymbol\mu},\overline{\Sigma})$, with mean   $$ \overline{\boldsymbol\mu}=\boldsymbol\mu_1+\Sigma_{12}{\Sigma_{22}}^{-1}({\boldsymbol a}-\boldsymbol\mu_2) $$ and covariance matrix $$ \overline{\Sigma}=\Sigma_{11}-\Sigma_{12}{\Sigma_{22}}^{-1}\Sigma_{21}$$ Actually these results are provided in Wikipedia too, but I have no idea how the $\overline{\boldsymbol\mu}$ and $\overline{\Sigma}$ is derived. These results are crucial, since they are important statistical formula for deriving **Kalman filters**. Would anyone provide me a derivation steps of deriving $\overline{\boldsymbol\mu}$ and $\overline{\Sigma}$ ? Thank you very much!


title: Validating an existing questionnaire into another language

text: I am doing research and trying to use a questionnaire that is available in English in another language.   1. Is there a systematic approach for validating a questionnaire in another language?    2. What statistical test should I be undertake and why?


I am interested in applying Bayesian additive regression trees (BART) for classification analysis of gene expression data. I am relatively new to R (and Bioconductor packages) and I am unable to find some code or vignette that I can use to learn from. I will be thankful if someone can point me in a good direction.


title: Non-linear regression fails to converge, but fit appears good

text: I'm fitting a 4 parameter nonlinear regression model to multiple datasets, some of which fail to converge, however, the parameters output after a failure provide a fit that looks good, if not exceptional to my (and other's) eyes. I've explored convergence criteria and they do converge eventually but the visual fit is terrible. Is there any precedent for taking the visual fit ad ignoring convergence or are there some other things I can try? I'm fitting a model of the form $\sqrt{c_1 x+c_2\exp(x)^y}/\sqrt{\exp(x)^y}+c_3$ where $y$ is a known constant, using nlminb in R.


title: Compression theory, practice, for time series with values in a space of distributions (say of a real random variable)

text: **Example of problem:** Part of our research team is working on providing operationally wind power forecast. Usually, since there are different time scalse that interest forecast user, a forecast is issued every 15 min (it has even happened that 5 seconds was requirer) for every 15 minutes ahead up to serveral days. Obviously when you contract with a wind farm owner he has many wind farms, can be 10 or even 1000 wind farms depending of who you are working with. It is likelly that this number is going to be bigger in the near future. Now if the forecast user wants not only the point forecast but the forecast of 100 intesreting quantile (i.e. the whole marginal distribution) and if he wants you to store everything in a data base year after year, then you start tinking ..... ok that makes a lot of data. I'm not an international meteorological institute and cannot afford that. **Edit (according to @whuber's comment):** to make it critical it can be tough of with a 5 second time frame. But anyway this is not meant to be the question, just an example (people coule provide other examples) **My question is** : is there a compression theory/best practice for time series of distributions (i.e function of time with values in a space of distribution). I agree that you could work quantile by quantile and apply simple compression algorithm dedicated to signals (could be based on wavelet) but I am searhing for something more dedicated. For example if you know everything is well approximated by gaussians, the quantile by quantile approach would be stupid ! **Edit (according to @whuber's comment):** The question is really how to integrate the temporal dimension in a compression scheme, hence nI don't want to make a choice of a subset of quantiles that would be the same for ALL time. The remark with gaussian is not dedicated to the example, it is an example to illustrate a case when very good compression is required. An other trivial remark (from the theoretical side) would be: imagine you have a very complicated distribution that cannot be summarized with few quantiles and that you observe it for 1000 successive times. It would be sad to store it for each of these times.


title: Kurtosis of made up distribution

text: Take a look at the image below. Blue line indicates standard normal pdf. The red zone is supposed to be equal to the sum of grey areas (sorry for awful drawing). I wonder can we create a new distribution with higher peak by shifting grey zones to the top (red zone) of the normal pdf? ![new distribution with higher peak](http://i.stack.imgur.com/uTOL3.png) If such transformation can be made, than what do you think about the kurtosis of this new distribution? Leptokurtic? But it has the same tails as the normal distribution does! Undefined?


title: How to report Kruskal-Wallis test?

text: As part of the analysis of the data collected from a survey, I was carrying out Kruskal-Wallis Tests between some Likert scale questions and demographics such as firm size and job position. I am using this test instead of One Way Anova since the data is non-normal. For example, testing the responses of the following statement vs Job position (manager or junior): I arrive late at work. Rated on a 5-point scale from never to always. When reporting the results of the test (if result is significant) can I state something like: _Managers are more likely to arrive late than juniors (H=14.338, p <.01)_ or is this the only way one can report the result: _A statistically significant difference (H=14.338, p <.01) exists between late arrivals at work by managers and juniors._


title: Mixture of normals, dependent on past

text: I have the following probability model: $(X_k|\text{PastHistory}_{k-1}, \theta_0,\theta_1,\theta_2) \sim (\pi\cdot N(\theta_1+\theta_0\cdot X_{k-1},1)+(1-\pi)\cdot N(\theta_2+\theta_0\cdot X_{k-1},1))$ and $X_k=0.5X_{k-1}+\epsilon_k$ with $\epsilon_k \sim N(0,1)$ I want to build a function that gives a sample of size n, from the mixture. I've tried something like this in Mathematica:                 Sample[n_] := Module[{lst},         lst = {RandomVariate[NormalDistribution[0, 1/(1 - 0.5^2)]]};         element = 1;         While[element <= n,           AppendTo[lst,              RandomVariate[               TransformedDistribution[                p*X1 + (1 - p) X2, {                X1 \[Distributed]                 NormalDistribution[ θ1 + θ0*lst[[element]], 1],                 X2 \[Distributed]                 NormalDistribution[ θ2 + θ0*lst[[element]], 1]             }]]];          element = element + 1;         ];        lst        ];      Is this correct? Many thanks for any help. ;)


title: Confusion about the term "stochastic drift"

text: After reading lots of material about the subject, I believe that **the term "stochastic drift" is defined in a two different ways.** These two different definitions make the term unambiguous and I assume that this is the reason you find so many questions regarding unit root testing and ADF. **Definition one** (quoted of "Basic Econometrics" by Gujarati): $Y_t = \beta_1 + Y_{t-1} + \epsilon_t$ which is a random walk with drift and therefore non stationary. If we write it as $ \Delta Y_t = \beta_1 + \epsilon_t$ this means that $Y_t$ will exhibit a positive ($\beta_1 > 0$) or negative ($\beta_1 < 0$) trend. Such a trend is a stochastic trend. **Definition two** (found in our professor his slides, wikipedia and this forum): The random walk model with drift is given by $Y_t = \beta_1 + Y_{t-1} + \epsilon_t$ Backward iterating the process yields: $Y_t = \beta_1 + Y_{t-1} + \epsilon_t$ $=\beta_1 + \beta_1 + Y_{t-2} + \epsilon_t + \epsilon_{t-1}$ $=\beta_1 + \beta_1 + \beta_1 + Y_{t-3} + \epsilon_t + \epsilon_{t-1} + \epsilon_{t-2}$ ... $=y_0 + \beta_1 t + \sum_{i=0}^{t-1} \epsilon_{t-i}$ The intercept in the non-stationary AR model adds a deterministic trend, i.e. $\beta_1 t$ to the stochastic trend $\sum_{i=0}^{t-1} \epsilon_{i-1}$. **Wikipedia definition** (http://en.wikipedia.org/wiki/Dickey%E2%80%93Fuller_test) Wikipedia derives exactly the same notation as in the slides of our professor, being: $y_t=y_0 + \beta_1 t + \sum_{i=0}^{t-1} \epsilon_{t-i}$ But adds "with a deterministic trend coming from $\beta_1 t$ and a stochastic intercept term coming from $y_{0}+\sum_{i=0}^{t-1} \epsilon_{t-i}$, resulting in what is referred to as a stochastic trend." **Conclusion and final question** What terminology do you use? I think I like trend stationary and difference stationary the best. Furthermore, I also think the terms deterministic trend and stochastic trend are now clear to me. I don't like the way Gujarati and Wikipedia phrases it, I prefer definition two of the terms.


title: Comparing variability of two groups of multivariate data points

text: I have a number of samples (data points) from two experimental groups, A and B (500 altogether). Each sample corresponds to a large number of variables (genes in my case); together, they form a matrix of $500 \times 40000$. I would like to answer the question, are samples in group A on the general more variable than samples in group B? The biology of the problem suggests that this might be the case (read, this is the hypothesis that I want to test). I thought I might do the following: for each group, calculate the multivariate center of the group (by taking sample-wise means for each variable, for example), then calculate the euclidean (or some other) distance to the centre for each sample in the group and then the variance of the distance. Then, do the same for the other group. I am well aware that this is a very naive approach. What would be the advisable approach?


title: Framework for Symbolic Classification

text: Please advise good framework for symbolic classification. I am currently using GPTIPS for that but I belive there are better options. What I am trying to do is following:   1. for set of features of my data I am solving standard symbolic multigene regression and as output for each individual in my generation I receive set of trees, where each tree represents some formula.   2. Each tree in my model encodes new feature in some new space. In this new feature space I am trying to solve classification problem. Currently I use kNN for that, but any other procedure should also apply. So basically I am solving standard classification problem, with additional step of feature space transformation using symbolic regression. Framework can be written on following languages: Java, Matlab, R, Python. I personally would prefer Java or Matlab, to Python or R. Thanks.


title: Show series of exponential decay distributions in one chart

text: Every week we have new customers sign up, and some of them submit support requests. I'd like to graph the change in the number of support tickets per customer over time. I'm doing this by bucketing customers into one-week cohorts, and then graphing each distribution as a Tufte box plot, as in this image. The graph I have now looks _something_ like this: ![tufte graph](http://i.stack.imgur.com/jk82x.jpg) The problem is that the distribution of real data, instead of the fake data I'm using here, resembles an **exponential decay** with only integer values, and is thus skewed toward "zero tickets". So we have 52 cohort distributions that look like this: ![exponential decay](https://img.skitch.com/20111122-f47p8di6f4agyh5pmn8ph2hxhp.jpg) (where the X axis is the number of tickets, and the Y axis is the number of customers who sent in that many tickets), making it tricky to get the actual distribution from the box plot. In general, the median doesn't tell you a whole lot about the distribution. How could I better show the series of exponential decay distributions, to get a grasp of how the shape is changing over time? Thanks, Kevin


title: Group comparison: joint model vs. separate model, contradictory result

text: Any hint on this issue would be really welcome: The goal is to employ a two- country comparison analysis. When running two separate OLS regressions, one for each country, the coefficient of variable x is statistically insignificant in country1 regression (Ho:b1=0 is not rejected), while the coefficient of the same variable x is statistically significant in country2 regression (Ho: b2=0 is rejected). Afterward, I run a joint model by including a country dummy (country1=1 for country1, and country1=0 for country2) and an interaction term: inter=x*country1. The coefficient of the interaction term appears to be positive and statistically significant, implying that the effect of variable x is greater in country1 than in country2. So, I am wondering whether there is any statistical explanation for this contradictory result; the separate regression analysis showed that the effect of variable x is insignificant in country 1 and significant in country2, but the joint model analysis showed that the effect of variable x is greater in country 1 than country2... Thanks a lot!


title: Calculating and comparing histograms of complex numbers

text: If I have a 1D sequence of complex numbers with real and imaginary parts, how can I compute the histogram of this sequence? Do I separately compute the histogram of both the real and imaginary parts? Moreover, if I get a 2D histogram of a sequence of complex numbers, how do I apply goodness of fit tests (such as the $\chi^2$ test) to compare two histograms of complex numbers? Is there a good reference book on the computation of histograms of complex numbers? How would such a histogram of complex numbers be visually represented?


title: Basic easy rules for statistics

text: In a binomial experiment, if we observe $x=0$ positive individual among $n$ individuals, then the proportion of positive individuals is significantly lower than $3/n$ with a type 1 error less than and very close to $5\%$. This fact, sometimes called the "rule of three", is a consequence of the inequalities $$\exp\left(-\frac{np}{1-p}\right) \leq \Pr(X=0) \leq \exp(-np).$$ Do you know other such basic easy rules for statistics? I find them very interesting and useful. This principle is not really a "rule of thumb" because it has a reliable theoretical foundation, but I don't see another tag for this question (I hope it is not off-topic)


title: Modeling consumer demand for a multi-tiered product that is subject to cannibalization?

text: I'm looking for research/white papers on demand estimation for goods/services that are offered on a tiered-basis with varying prices (i.e. more complex the good/service, higher the price.) I'm not sure if this kind of good/service would fall under "bundled goods/services". A peculiarity of the marketing strategy for this good is that the different tiers of goods cannibalize each other during sales (e.g. demand for full price cheaper model drops when a higher priced model goes on sale.) I'm interested in researching if such a tiered good is treated as a single good or multiple goods in demand estimation, and also the specific econometric methods to approach it (e.g. multivariate regression with demand and supply equations.) Thank you for your time!


title: Expected value, correlation, and independence

text: I need help with a problem. Supposed x, y, and z are events in F (algebra of sets) in a probability space (universal set, F (algebra of sets), P). Define two random variables: a(omega) =                  {1 if omega is an element of set X        {0 if omega is not an element of set X      b(omega) =                  {2 if omega is an element of Y        {1 if omega is an element of Z        {0 if omega is not an element of Y or Z      a) Find E[A], E[B], and E[AB] b) Are A and B uncorrelated? Under What conditions on (X,Y,Z) is sufficient for zero correlation between A and B? c) Are A and B independent? What conditions on (X,Y,Z) is sufficient for zero independence of A and B? Attempts: a) Find E[A], E[B], and E[AB] E[A] = sum(xi*p(i)) = 1*p(X)+0*(1-p(X) = p(x) E[B] = sum(xi*p(i)) = 2*p(Y)+1*P(Z)+0*P(YUZ) = 2p(Y)+p(Z) E[AB] = E[X]*E[Y] if independent and E[AB] = E(A)E(B) + Cov(A,B) if not independent How can I find the Cov(A,B) of this problem? b) Are X and Y uncorrelated? What conditions? X and Y are uncorrelated when cov(A,B) = 0 which is when E[AB] - E[A]-E[B] = 0. c) Are X and Y independent? What conditions? X and Y are independent when P(X intersect YUZ) = P(A)*P(YUZ)


title: How to test Simultaneity

text: I am researching the influence of ΔENVIRONMENTAL TAXES+ΔENV.EXPENDITURE =>ΔCO2 EMISSIONS. and further, the influence of ΔCO2 EMISSIONS =>ΔENV.TAXES after a certain period of time. I have a time series of 16 years for 19 countries. I have done OLS, FE, RE and Hausman test (to check for endogeneity; and if Hausmant test is significant I choose RE, and when RE did not work, I used an INSTRUMENT (e.g. government public debt (in %). Is this procedure appropriate? However, this is not the end of the story. Furhter, I want to test both equations in a system together. I wonder if I can use lagged variables? For example in 1. equation I use ΔENV.TAXES_(t-1) with lag 1, and in 2. equation I use ΔENV.TAX_(t) without lag. May I still count ΔENV.TAXES_(t-1) and ΔENV.TAXES_(t-1) as same variable? For example in every equation I use the most significant lag: ΔENV.TAXES_(t-1) + ΔENV.EXPENDITURES_(t-1) Izdatki_(t-2)=> ΔEmisije_t ΔEmisije_(t-2) => ΔDajatve_t In case, lagged variables are not appropriate solution, than I suppose I should use only Δvariables in time t. However, if I do this, than this is not the same substantive issue I want to address. e.g. that taxes and expenditure influence on the change in emission after a certain period of time. Is there a test I could use to check whether I have simultaneity in my model at all? Which statistical program is the most appropriate for tackling this issue? Which test should I do to test both equations in a system together?


title: Can ANOVA be significant when none of the pairwise t-tests is?

text: Is it possible for one-way (with $N>2$ groups, or "levels") ANOVA to give a significant difference when none of the $N(N-1)/2$ pairwise t-tests does? Apparently it is possible, but I do not understand how. When does it happen and what the intuition behind such a case would be? Maybe somebody can provide a simple toy example of such a situation? Two further remarks:   1. The opposite is clearly possible: overall ANOVA can be non-significant while some of the pairwise t-tests erroneously report significant differences (i.e. those would be false positives).   2. My question is about standard, non-adjusted for multiple comparisons t-tests. If adjusted tests are used (like e.g. Tukey's HSD procedure), then it is possible that none of them turns out to be significant even though the overall ANOVA is. This is covered here in several questions, e.g. How can I get a significant overall ANOVA but no significant pairwise differences with Tukey's procedure? and Significant ANOVA interaction but non-significant pairwise comparisons. **Update 1:**   * I do not actually have an example of such a situation! My question is purely theoretical; I am interested whether such a situation can happen in principle.   * I thought that it is not possible. But then I saw an answer by @whuber in this thread, saying that "It is well known that a global ANOVA F test can detect a difference of means even in cases where no individual [unadjusted pairwise] t-test of any of the pairs of means will yield a significant result". I got confused and asked this question. **Update 2:**   * My question originally referred to the usual _two-sample_ pairwise t-tests. However, as @whuber pointed out in the comments, in the ANOVA context, t-tests are usually understood as _post hoc_ contrasts using the ANOVA estimates (which is not the same thing as a two-sample t-test). So there are actually two different versions of my question, and the answer to both of them turns out to be positive. See below.


title: Why is variance calculated as $p(1-p) / (n -1)$

text: I had to translate several given statistics equations into code, and I came across this formula: Variance of a simple random sample $= \frac{p(1-p)}{n-1}$ The sample in question are test letters sent to measure the efficiency of a postal operator. Each letter has a value of 1 if it's delivered on time and a value of 0 otherwise. The parameter $n$ is the total number of sent letters. The parameter $p$ is the "number of letters delivered on time" and an estimator for the true efficiency. It's defined as the sum of all letter values from 1 to $n$. I'm not an expert in statistics, but as far as I know the variance is the sum of squared errors from the sample mean divided by $(n-1)$. I couldn't find any explanation for this online. Can anyone explain this formula?


title: Linear regression when you only know $X^t Y$, not $Y$ directly

text: Suppose $X\beta =Y$. We don't know $Y$ exactly, only its correlation with each predictor, $X^\mathrm{t}Y$. The ordinary least-squares (OLS) solution is $\beta=(X^\mathrm{t} X)^{-1} X^\mathrm{t}Y$ and there isn't a problem. But suppose $X^\mathrm{t}X$ is near singular (multicollinearity), and you need to estimate the optimal ridge parameter. All the methods seems to need the exact values of $Y$. Is there an alternative method when only $X^\mathrm{t}Y$ is known?


title: Variance-gamma distribution: parameter estimation

text: I have quite a general question about Variance-gamma distribution. I am interested in how to estimate it's parameters given a set of training points? I tried to find the answer in the internet, but surprisingly managed to find only a couple of relative links:   1. VarianceGamma R package (and it's manual)   2. The paper which is referenced in R package manual. I am not an expert in R, but what I saw in R package was (as I understood) maximum likelihood estimation. They perform some iterative optimisation with different methods, starting with Skew Laplace to initialise the estimator. In the paper they first estimate mean, variance, skewness, kurtosis through moments. Then they assume that asymmetry parameter $\theta$ is small, and set $\theta^2 = \theta^3 = 0$. And after all they solve combined equations to get all the parameter estimators. That's what I understood. So my questions are:   1. How exactly do they perform optimisation in R package? Maybe you know the place where it is described?   2. What is the difference in the two methods? Which one is better to use in real life (well, maximum likelihood is better, but it is much more difficult to implement and the performance is not as good)? Are the assumptions about small $\theta$ strict in the first method?   3. Where can I read more about Variance-Gamma parameter estimation? I would be really grateful for every relative replies, papers and links. Thank you!


title: Chi square for survival data

text: Can the chi-squared test be used for comparing survival curves between or among treatment groups? In the particular case of interest, there are 4 treatment groups.


title: How would you visualize the difference between Cox/Weibull regression?

text: I'm trying to figure a way of properly displaying the difference\resemblance between various regression values on the same data set, using cox ph, weibull regression and log-normal regression. Weibull and log normal are similar in their result (predicted values) but cox ph returns the survival function, so I'm not sure how to compare them. What would be a proper graph\visualization method for this? Thanks!


title: What can be concluded about coefficients?

text: I run multiple variable regression and got the following result. What can i conclude regarding the significance of coefficients? And why? Thanks               SUMMARY OUTPUT                                          ANOVA                                                    Coefficients    t Stat   P-value         Intercept   0.00191         6.27072  0.00000         RPt         0.98782         56.85940 0.00000         SMB         0.87223         31.19128 0.00000         HML         0.68983         31.76355 0.00000


title: Entering data on SPSS - Pre-Test Post-test Questionnaires on Emotion (with negative and positive values)

text: I'm fairly new here in this site. I'm currently taking Bach. Psychology, and Research Methods is one of my core subjects. I have conducted an experiment titled "Gender Differences and Effect of Types of Music on Emotion". So before presenting the problem I'm having, let me briefly describe my experimental design. So basically, it's a 2 (Male/Female) x 2 (Sad/Happy Music) Between Subjects design with 24 participants. For the first session, 12 participants (6 males and 6 females) are exposed to Sad song. For the second session, another 12 participants (6 males and 6 females) are exposed to Happy Song. I used the Genova Emotional Music Scale (GEMS-45) for both pre-test and post- test to measure if there are any changes in Emotion after the participants listen to a music (either sad or happy song). The pre-test serves as the baseline of their current emotion. Some confusions I'm having are:   1. From the first session (with sad song), I added the scores from 3 items (sad, tearful, and sorrowful). While for the second session (with happy song), I added the scores from 3 items as well (joyful, happy, and feel like dancing). Or should I take into account all 6 items (sad, tearful, sorrowful,joyful,happy, feel like dancing) for both sessions?   2. In the questionnaire, the Likert Scale ranges from 1-5. Now, for the sad session, I added negative values on the scores? This is to indicate that negative scores = negative emotion which is sad. So for example, a total score of a participant (in the sad session) for the pretest is 3, and 7 for the postest. The difference between the scores (postest - pretest) would be 4. It means that the participant actually felt sadder after listening to the song. So if the participant scored -3 on the pretest, and 7 on the post-test, the difference would be -4 (SAD) Now, same goes to the happy session. If the pretest is 4, and the post-test score is 9, the difference would be 5 (or +5, which indicates positive emotion = happy) correct? The reason why i came up with the negative value for sad is to avoid SPSS being confused. But I'm not too sure about this. Or, would I have to create a new scale, which consists of very sad to very happy on a 7-Point likert scale (1 being very sad, 7 being very happy)? I hope to get some feedbacks about this as I'm having trouble with it. I know I can run ANCOVA (pretest as covariate, and post-test as DV) or Independence Factorial (Emotion, which is the differences between post-test and pre-test as the DV, and the gender and music as IVs), but I'd hope to enter the data correctly before running any stat-test.


title: Non-parametric correlation for continuous and dichotomous variables

text: I have two variables I want to test with correlation, one is continuous and the other dichotomous. My data are non-normally distributed, plus the variance is heterogeneous, so I have to apply a non-parametric test of correlation. The problem is, I don't know of any test that fits my requirements. Spearman's correlation coefficient requires ordinal data for both variables. Point-biserial correlation coefficient fits my data type, but it's a parametric test. I'm unaware if there's some test else out there. What should I do?


title: Determine number of runs is significant at 5% level of significance

text: I have the following question to answer: > The total thickness X of the four pads of 36 half-ring (H-R) mounts for > aircraft engines, taken periodically from the production line, were found to > be as shown below. Determine whether the total number of runs above and > below the median of this sequence of values of X is significantly different > from its expected value at the 5% level of significance, under the > hypothesis that X is under statistical control. ![enter image description > here](http://i.stack.imgur.com/We6N2.png) I believe this involves Chi Square but I honestly don't know where to begin.


title: Why must these regressions on a constant equal the mean?

text: Show that the OLS, MoR and RoM estimators for the regression of a variable on a constant (ie a variable which always takes the value one) are all equal to just the mean of that variable.


title: Question about long vectors and covariance

text: I have a dataset with triplets of vectors $\{v_1,v_2,v_3\}_{i=1\dots n}$. Each of them can be understood as a time series (or as a window of a time series). In which way can I give a number "similar" to the covariance of these vectors? I mean, I have a covariance matrix of each of them, and I could get a cross- covariance matrix of the pairs ($\Sigma_{11},\dots,\Sigma_{33}$)... But these matrices are not very useful. Can I use any metric of the matrices? (for example the biggest eigenvalue), or is there another, more standard way of getting a number out of them? Thank you very much.


title: How can I estimate the probability of a random member from one population being "better" than a random member from a different population?

text: Suppose I have samplings from two distinct populations. If I measure how long it takes each member to do a task, I can easily estimate the mean and variance of each population. If I now hypothesise a random pairing with one individual from each population, can I estimate the probability that the first is faster than the second? I do have a concrete example in mind: the measurements are timings for me cycling from A to B and the populations represent different routes I could take; I'm trying to work out what the probability is that picking route A for my next cycle will be faster than picking route B. When I actually do the cycle, I've got another data point for my sample set :). I'm aware that this is a horribly simplistic way to try to work this out, not least because on any given day the wind is more likely to affect my time than anything else, so please let me know if you think I'm asking the wrong question...


title: MCMC sampling of decision tree space vs. random forest

text: A random forest is a collection of decision trees formed by randomly selecting only certain features to build each tree with (and sometimes bagging the training data). Apparently they learn and generalize well. Has anybody done MCMC sampling of the decision tree space or compared them to random forests? I know it might be computationally more expensive to run the MCMC and save all the sampled trees, but I am interested in the theoretical features of this model, not the computational costs. What I mean is something like this:   1. Construct a random decision tree (It would probably perform horribly)   2. Compute likelihood of the tree with something like $P(Tree|Data) \propto P(Data|Tree)$, or perhaps add a $P_{prior}(Tree)$ term.   3. Choose a random step to change the tree and select based on the likelihood $P(Tree|Data)$.   4. Every N steps, save a copy of the current tree   5. Go back to 3 for some large N*M times   6. Use the collection of M saved trees to do prediction Would this give a similar performance to Random Forests? Note that here we are not throwing away good data or features at any step unlike random forests.


title: Running time dependent covariate model in R, getting warning

text: I am trying to fit this model, but it's giving a warning message I don't understand:               cfit<- coxph(Surv(final$start,final$end,final$event1)~imp+sec+cluster(final$compid))          Warning messages:     1: In fitter(X, Y, strats, offset, init, control, weights = weights,  :       Ran out of iterations and did not converge     2: In coxph(Surv(final$start, final$end, final$event1) ~ imp + sec +  :       X matrix deemed to be singular; variable 21 22 24 28 29 31 32 33 34      Further to this, does `coxph` allows continuous predictors?


title: How to run JAGS as fast as possible on Mac?

text: I am running lots of simulations in JAGS 3.1.0 , which I have compiled from source. I am using a mid-2009 Macbook Pro 13 2.53GHz. Is there any step by step tutorial (from scratch) to compile ATLAS or GotoBLAS for JAGS and R (where I use the rjags package)? There are information on ATLAS, GotoBLAS, JAGS, and R compilation, but for one reason or another I was not successful in getting them to work together. And lastly, does it even help the speed of simulations with the machine that I am using? Note: I am not looking for parallel processing since I am running two separate simulations simultaneously which use the 2 cores of CPU.


title: What are the assumptions for quantile regression?

text: What assumptions must be fulfilled in quantile regression?


title: Covariance of transformed random variables

text: I have two random variables $X > 0$ and $Y > 0$. Given that I can estimate $$\text{Cov}(X, Y),$$ how can I estimate $$\text{Cov}(\log(X), \log(Y))?$$


title: Sample size and power detection

text: We are currently doing a research on physiological changes of platelet aggregation during three trimesters of pregnancy and postnatal period, then we compare these groups to each other and the control arm as well. We have recruited 10 patients per each trimester, 10 postnatal and 6 as control arm. The total sample size is 46. We wonder how we can calculate the study power? FYI we used Kruskal-Wallis Test to compare mean of the five groups, then if the resulted P-value was significant or close, we tried to perform post hoc Mann-Whitney Test. We also used same agonist (interference) for all groups but the participants are different in demographic point of view.


title: Test for linear separability

text: Is there a way to test linear separability of a two-class dataset in high dimensions? My feature vectors are 40-long. I know I can always run logistic regression experiments and determine hitrate vs false alarm rate to conclude whether the two classes are linearly separable or not but it would be good to know if there already exists a standard procedure to do that.


title: In convolutional neural networks, how to prevent the overfitting?

text: Given certain amount of labeled data, we define the net structure, such as number of layers, types of layers, the number of convolutional layers, the number of pooling layers, etc. And train the parameters using back propagation, meanwhile we show the loss in training procedure and view the testing accuracy in validating data set. But, the loss in training set is nearly zero, and the testing accuracy is kept unchanged no matter how to decrease the learning rate. In this circumstance, is it overfitting? Should we change the net structure? More layers for more parameters?   Could you please recommend some suggestions or references?


title: I am getting a number below zero when calculating out two standard deviations from the mean. Is this ok?

text: I am not a math wizard, so please keep your response simple enough. I need to complete a statistics screening exam for a methods course later on today, and I am hung up on one topic that came up during the practice test. The data set I got was in reference to the number of homicides that have occurred in a number of cities. The range of this data is 0-5. When I am putting together confidence intervals and calculating out as much as two standard deviations from the mean I am getting low values that are negative. Obviously you cannot have a negative number of homicides. When calculating the confidence intervals out to two standard deviations from the mean should I present the low value at ZERO or should I actually present the negative number? For example, if a 95% CI caused the calculation to be -1.5 to 3, would I present that or would I present 0 to 3? Thanks.


title: interpreting the equivalent of one SD in percentage

text: Assuming normal distribution, one SD is equivalent to an effect size of 1.0, which is associated with a gain of 34 percentage points from the average. From a simple regression model, if past achievement, a score from E to A grade, predicted a change of one SD in current achievement, a standardized score (mean of 0), can I infer that an increase of one letter grade in past achievement (e.g., C to B) can result in a 34 percentage gain in current achievement? Or would it be more correct to say "can explain about 34 percentage of the variance"? Any comments appreciated. Thanks!


title: Fitting a model to correlated data

text: Trying to fit a model (bloch equations from MRI) to experimental data (12 parameters, 12 equation). Monte Carlo Multistart ('lsqcurvefit' in MATLAB) gives expected values for all parameters except for two, which are correlated: ![http://imgur.com/1PRayGj](http://i.stack.imgur.com/Z2Per.jpg) Briefly about the model. It is a system of differential equations (solved numerically by Cramer's rule). The fig. contains 9 diff. eq. My model is a bit more complicated and contain 12 diff. eq. http://imgur.com/E5LiYdq (the model) I am passing initial guess, lower and upper bounds for all parameters to the solver. 1st parameter always takes the value of its upper limit and 2nd one the value of its lower limit (please, see the fig.). Bad luck, because I am only after those two parameters. Are there any good fitting algorithm for correlated data?


title: Posterior distribution for true population variance

text: I have a process which may be assumed Gaussian - for now I'm assuming it to be zero-mean. The variance is unknown. I will estimate variance in the usual unbiased way to give the sample variance and thus the ratio of the sample variance to the true variance follows the Chi Squared distribution. In parallel i have developed an equation for a quantity of interest (a risk) which is a function of the pdf of the true population variance. What inference approach do/should i take to give me a pdf of this unknown population variance? Thanks in advance.


title: Poisson process and uniform distribution

text: **Question:** > A single-pump petrol station is running low on petrol. The total volume of > petrol remaining for sale is $100$ litres. > > Suppose cars arrive to the station according to a Poisson process with rate > $\lambda$, and that each car fills independently of all other cars and of > the arrival process, an amount of petrol that is distributed as a uniform > random variable over $(0, 50)$ - assume for example that all car tanks have > a capacity of $50$ litres and drivers decide "at random" when to refill. > > We assume that service is instantaneous so that there are no queues at the > station. > > **(a)** On average, how many cars will the petrol station fully service > (sell the full amount requested) before it runs out of petrol (and before > any refilling occurs)? > > **(b)** How much time will it take on average before the station runs out of > petrol (and before any refilling occurs)? * * * Attempt at part **(a)** : I don't think I need to use the assumption that cars arrive according to a Poisson process in this part. So, define $U_k$ as the random variable that denotes the amount of petrol that car $k$ fills, $k = 1, 2, 3, \cdots$. Thus, $U_k, k = 1, 2, 3, \cdots$ are independently and identically distributed as a uniform random variable over $(0,50)$. Let $N$ be the random variable that denotes the number of cars that the petrol station can _fully_ service. Now I don't quite get the question. Say we have the following scenario: Car 1 comes with 10L remaining in its tank, so it will fill up 40L, hence the amount of petrol left in the station is now 100-40 = 60L. Car 2 comes with 10L remaining in its tank, so it will fill up 40L, hence the amount of petrol left in the station is now 60-40 = 20L. Car 3 comes with 20L remaining in its tank, so it will fill up 30L, but the petrol station only has 20L left, so does this mean Car 3 just leaves the petrol station filling 0L? My gut feeling is that this cannot happen because each car can only fill an amount BETWEEN 0 and 50, ie, (0,50) [note that the end points are not included]. Hence in this scenario, the petrol station runs "out" of petrol at N=3 because it does not have enough to FULLY service Car 3, even though it still has 20L left in the pump. Thus, the petrol station can only service N=2 cars. Is this interpretation correct? If so, how do I find $E[N]$? * * * EDIT: working on further... Define $G_N = \sum_{k=1}^N U_k$, then $E[N] = \sum_{n=1}^{\infty} n P(N=n) = \sum_{n=1}^{\infty} n P(G_n \le 100)$ This equivalence comes from the fact that the event $\{N=n\}$ will only happen if $100 - G_n \ge 0$. However, two questions remain, what is the distribution of $G_n$? (How do I derive the distribution of the sum of $n$ iid uniform random variables and second, how do I compute the infinite sum?


title: Fit between two curves

text: I am trying to determine the fit between two curves (example shown below). Both curves are plotted here as a bunch of XY data points (not functions). The curves represent the trajectories of two planes taxiing, and I'm trying to figure out how far apart they are on average. (axes are in x/y position in meters) I've done correlation and RMSE, but neither of those is really what I want-- the two planes are travelling at different speeds, and thus one curve has many more data points than the other. So RMSE is sensitive to time, whereas I want deviation independent of time. My next thought is to do some sort of integration, but I'm not sure how that would work as the curves crisscross frequently. As you can see, the curves are pretty similar so I'm looking to get a small number as a result... probably around 0-5 meters. What sort of analysis should I be doing here? ![enter image description here](http://i.stack.imgur.com/okkz1.jpg) A close up of the turn: ![enter image description here](http://i.stack.imgur.com/8Em4x.jpg)


title: How many copies of shredded text are needed for getting 99.9% reconstruction of mutated text?

text: This is a question posed by NASA astrobiologist Chris McKay (who wants to put life together back from its pieces), who really wants an answer to this question since it could lead us to reconstruct ancient genomes that are billions of years old (but mutated through years of internal radioactivity)   * A text of length N composed of an alphabet of m distinct units (eg Hamlet has 29551 words and 60 distinct units - 52 letters, space, plus 7 punctuation marks)   * Shred the text in pieces of average length X   * How many copies C of the shredded text are needed to give 99.9% reconstruction?


title: Validation of mixed-effect models

text: I want to use linear mixed effect model for a set of data. After using lme4 package and `lmer()` function and fitting model, I want to validate my model for other period time. But I don't know how to validate my model and generate predicted values. Do you have any ideas about useful validation method? Is there any command about validation in R software? * * * Thanks for your useful comments. I fitted model to one set of data and now I want to check its performance for an other data set. I want to plot figure by "predicted value" and "observed value" and calculate R^2 to evaluate model performance. Thanks


title: How to use the extracted components from varimax rotation as independent variables?

text: I want to use body measurements to predict body weight in broiler chickens. There were seven body parameters measured. I did a principal components analysis and three components were extracted after varimax rotation. How do I use these components as independent variables to predict body weight (the dependent variable)?


title: Why do people use p-values?

text: Roughly speaking a p-value gives a probability of the observed outcome of an experiment given the hypothesis (model). Having this probability (p-value) we want to judge our hypothesis (how likely it is). But wouldn't it be more natural to calculate the probability of the hypothesis given the observed outcome? In more details. We have a coin. We through it 20 times and we get 14 heads (14 out of 20 is what I call "outcome of experiment"). Now, our hypothesis is that the coin is fair (probabilities of head and tail are equal to each other). Now we calculate the p-value, that is equal to the probability to get 14 or more heads in 20 flips of coin. OK, now we have this probability (0.058) and we want to use this probability to judge our model (how is it likely that we have a fair coin). But if we want to estimate the probability of the model, why don't we calculate the probability of the model given the experiment? Why do we calculate the probability of the experiment given the model (p-value)?


title: Robust test for time series count data

text: I'm going to analyse suicide rates for a time series, and I'd like to use robust tests, but I don't know which would be a good one. My purpose is to compare the variation of the suicide rates through the months and years and see if there are significant changes in respect to GDP. Then I also want to check differences among sex, age, etc. What would you recommend to me? I've heard about robust tests in _Introduction to Robust Estimation and Hypothesis Testing_ by Rand Wilcox. They are the best choice when normality and other assumptions are violated. Thanks.


title: What are the properties for independence

text: For 2 variables to be independent of each other, should the correlation = 0 or mutual information = 0 or covariance = 0. I have seen different conditions and all these are really confusing.


title: Statistics collaboration

text: As a biologist, many of the research projects I work on at some point involve collaboration with a statistician, whether it be for simple advice or for implementing and testing a model for my data. My statistics colleagues admit that they do a significant amount of collaboration, insomuch that the tenure review process only considers papers on which they are the first or last author. What would make me (or any other scientist) a better collaborator? What would make it easier for you (as a statistician) to work with me? Specifically, what is one statistics concept you wish all of your scientist collaborators already understood?


title: Can you convert three-way ANOVA to one-way ANOVA?

text: I'm analyzing some data from an experiment. There is one response variable (response) and 3 factor variables, each with 2 levels (profile, drug, disease). I want compare specific combinations of (profile, drug,disease), e.g. (profile 1, drug 1, no disease) vs. (profile 1, drug 2, no disease). I've been approaching this problem as if it were an ANOVA design but I've hit some trouble when trying to make the specific comparisons I want. I initially tried this analysis in R but couldn't find any information on how to do contrasts with a 3-way ANOVA design. I'm now using Stata, which is only slightly more user friendly, but even that is giving me some trouble. So what I've tried doing is writing some code to convert my different factor variables into one factor variable with more levels. So instead of (profile, drug, disease) we have profile:drug:disease E.g., profile1:drug1:disease1 and profile2:drug1:disease1 are different levels of the same factor variable (call it newFactor). So can I run my analysis as a one-way ANOVA on newFactor (in R, this would be aov(response ~ newFactor) and then use contrasts to make comparisons between, e.g., profile1:drug1:disease1 and profile1:drug2:disease1? Is there any thing wrong with this approach, statistically speaking? Sorry for the long question, and thanks in advance for your help. tl;dr version: I want to convert a 3-way ANOVA to a 1-way ANOVA by lumping all my factor variables into one new factor so that I can easily use contrasts.


title: What fraction of a parent’s genes exist in total across n offspring?

text: This problem has perplexed me for a long time, and I would love to learn how to solve it.   1. Each parent contributes, on average, half of the genetic material of each child. If a parent has $n$ children, what percentage of the parent’s genes exist in total across all $n$ offspring?   2. How can one generalize the solution to #1, to make it applicable to more distant relationships like uncles & aunts, grandparents, and cousins, where the fraction of shared genes is smaller? For the purposes of these questions, disregard the complexities of chromosomes, meiosis, and fertilization, and just assume that genes are chosen from each parent independently and at random.


title: Simulate predicted probabilities after a logistic regression

text: I'm generating predicted probabilities (propensity scores) in Stata using the predict command. The basic structure of my code looks like this:               logit depvar indepvars     predict ps               // Gives Pr(depvar = 1 | indepvars)     predict s, stdp          // Gives the standard error of the prediction      The predicted probabilities/propensity scores are obviously estimates with an error distribution around them. What I would like to do is use the results obtained from the `predict` command to simulate some plausible values for the predicted probabilities/propensity scores. Any advice on how I might proceed? Thanks in advance.


title: Which variables to retain in order to preserve the same clustering pattern?

text: Suppose I have 50 scale parameters, these are all genes measured for one sample from a subject at the clinic, after data reduction by PCA, two meaningful components were extracted. This was followed by cluster analysis and turned out to be 4 meaningful clusters of subjects based on the two components of these 50 genes. Since investigating 50 genes for one subject would be costy, one would like to reduce that number so that the same clustering pattern can still be obtained but with minimal costs possible ( there should be some measures here to say acceptable clustering or not, I wonder what kind of measures would fit this case though). Of course, the more genes investigated, the more information gained, but there should be some measure to tell when to stop wasting more money when the same result is _satisfactorily_ achievable will less number of genes. Is there any R package that already implemented this approach? what would be the statistical approach in this case to select the most important genes that would preserve the clustering pattern? what criteria to be used in order to reach the minimum clustering pattern?


title: Neural Networks output range in simulation

text: I am learning some model based on examples ${((x_{i1},x_{i2},....,x_{ip}),y_i)}_{i=1...N}$ using a neural network of Feed Forward Multilayer Perceptron (newff) (using python library neurolab). I expect the output of the NN to be positive for any further simulation of the NN. How can I make sure that the results of simulation (testing) of my learned NN are always positive? (how I do it in neurolab?)


title: "Repeated measures" when the measures are not independent?

text: I have a set of data consisting of measurements in ten subjects of the same physiological parameters. The measurement and processing of the data have been carried out at two sites, A and B. For each measure and each subject, I thus have four versions of the same parameter: Measured at A, processed at A (AA);   Measured at B, processed at B (BB);   Measured at A, processed at B (AB);   Measured at B, processed at A (BA). The main hypothesis is that results AA and BB are the same, which I can easily compare using a paired t-test. But then the secondary question is, given the t-test shows a difference, how can I tease out whether it is to do with the measurement or the analysis? As I understand it, there are two main options:     1. A repeated measures mixed linear model with the subject as a fixed factor and measurement and processing as random factors.     2. A set of paired t-tests (Holm-corrected?) as follows:   a) AA v. BB (for overall difference between the sites)   b) AA v. AB and BA v. BB (testing for difference in processing)   c) AA v. BA and AB v. BB (testing for difference in measurement)   My questions are as follows:     1. Is using approach 1 (MLM) valid in this case, specifically when the four 'repeated measures' are not independent of each other?    2. If no, is approach 2 the best way to deal with the data?    3. How would one combine the p-values (or t-values) from the two tests in 2b (and 2c)? I am also contemplating Holm correction to control my Type I error-rate, but I'm not convinced it's worth doing with only three tests (and I'm not sure it's even valid when they're not really independent of each other), so if anyone has a strong opinion I would be interested to hear it. Thanks for any advice, and I hope I've asked this question comprehensibly. I am a physicist who has strayed into the muddy world of medicine so this is all new territory!


title: How to do a regression analysis of bank systemic risk and liquidity ratio

text: I am looking to determine the effect the Liquidity Coverage Ratio (From Basel III) may have towards the issues of systemic risk or on bank lending rates. How would I go about doing a regression analysis for this? Thanks! Please help if you can :)


title: Cumulants in chinese restaurant process?

text: I already wrote a similar question on StackOverflow but was not welcomed there! So I decided to ask it from the folks here. I have wrote a code in Python for CRP problem. I think everybody here is familar with the subject but nevertheless: Short description of it: Suppose we want to assign people entering to a restaurants to potentially infinite number of tables. If $z_i$ represents the random variable assigned for the $i$'th person entering the restaurant the following should hold: With probability $p(z_i=a|z_1,...,z_{i-1})=\frac{n_a}{i-1+\alpha}$ for $n_a>0$, $i$'th person will sit in table $a$ and with probability $p(z_i=a|z_1,...,z_{i-1})=\frac{\alpha}{i-1+\alpha}$ $i$'th person will sit around a new table. I am not quite sure if my code is correct cause I am surprised how small the final number of tables are. I would be happy if somebody could give me cumulants for the distribution associated with this process.               import numpy as np     def CRP(alpha,N):         """Chinese Restaurant Process with alpha as concentration parameter and N          the number of sample"""         #Array which will save for each i, the number of people people sitting         #until table i         summed=np.ones(1) #first person assigned to the first table         for i in range(1,N):             #A loop that assigns the people to tables                  #randind represent the random number from the interval [1,i-1+alpha]             randind=(float(i)+alpha)*np.random.uniform(low=0.0, high=1.0, size=1)             #update is the index for the table that the person should be placed which             #if greater than the total number, will be placed in a new table             update=np.searchsorted(summed,randind,side='left')             if randind>i:                 summed=np.append(summed,i+1)             else:                 zerovec=np.zeros(update)                 onevec=np.ones(summed.size-update)                 summed+=np.append(zerovec,onevec)         #This part converts summed array to tables array which indicates the number         #of persons assigned to that table         tables=np.zeros(summed.size)         tables[0]=summed[0]         for i in range(1,summed.size):             tables[i]=summed[i]-summed[i-1]         return tables     a=CRP(5,1000)     print a


title: Is it possible that grid search would fail in two dimensional feature space?

text: Grid search suffers from the curse of dimensionality. But is there any case(any hypothetical distribution of data) in a two dimensional feature space where the data’s binary classification using Grid Search would yield very poor results?


title: Iteration of the $\alpha$ parameter in a beta distribution yields a beta distribution?

text: I'd like to know whether starting from a beta distribution and iterating it in the way described below, I will get a stationary (beta) distribution again. More specifically this is the problem I am facing: Starting with a $\mathrm{Beta}(2,1.5)$, for each element $\phi$ drawn from this pdf, draw from an associated distribution; $\mathrm{Beta}(\phi+1,1.5)$. Iterating $n$ times, as $n$ goes to $\infty$, with this procedure do we get back a stationary distribution? Is it a beta distribution?


title: MULTIPLE TESTING: Probability of specific outcomes from MULTIPLE experiments using randomized data

text: This is a follow-up question to this one: Significance test across multiple simulated experiments There's one answer I'm leaning towards accepting, but I wanted to make sure I understood how significance can be calculated across multiple experiments OR estimated from a set of simulated experiments. I have 6 datasets representing millions of coin-flip type experiments over thousands of samples, in which some samples MAY be non-randomly distributed; i.e. cumulative binomial prob < 0.05 in 1,2,3,4,5 or all 6 datasets. Because I'm concerned about multiple testing in this setting, I want to know: How often do these samples score a binomial p value < 0.05 in 1,2,3,4,5 or all 6 simulated datasets (where I randomly flip the coin 100 times for each sample). Simple question: How would you approach this question? Also, please let me know if I should be more specific or ask the question in a different way. Thank you! UPDATE: Here's a specific example:               heads/tails:     exp. 1: 88/11, p < 0.05     exp. 2: 38/12, p < 0.05     exp. 3: 115/3, p < 0.05     exp. 4: 39/47, p > 0.05     exp. 5: 70/13, p < 0.05      exp. 6: 33/30, p > 0.05      4 out of 6 experiments show a binomial prob < 0.05, the other 2 above. Note that the total number of coin tosses differs between experiments. Although I could multiply the six individual p-values to calculate an overall probability of observing these 6 results, I want for each experiment to count EQUALLY, independently of the total number of coin tosses. That's important, because the number of "coin tosses" in the actual data can differ by orders of magnitude! Equally important, I'm concerned about multiple testing. I have > 30,000 samples in each experiment. If I have an overall p-value cutoff of 0.01, i will make 300 incorrect observations! That's why I wanted to simulate each of the 6 experiments 100 times over, with their number of coin tosses = the original data:               heads/tails with fair coin (Pr=0.5):     exp. 1: 99 tosses, observed p < 0.05 in 100 simulations = 12     exp. 2: 50 tosses, observed p < 0.05 in 100 simulations = 13     exp. 3: 118 tosses, observed p < 0.05 in 100 simulations = 9     exp. 4: 86 tosses, observed p < 0.05 in 100 simulations = 10     exp. 5: 83 tosses, observed p < 0.05 in 100 simulations = 7      exp. 6: 63 tosses, observed p < 0.05 in 100 simulations = 11      So, how would you calculate or use the simulated data to estimate the likelihood of observing the original 6 results by chance?


title: Best way to visualize predictions from a linear model

text: Let's say I'm doing some predictive analytics and am trying to predict US GDP per month using a two or three month lag. After every month, I generate new predictions and am able to compare my prediction to the actual gap. What would be the best way to represent this via plot. Line chart? In all the research I've done, there doesn't seem to be a 'standard' way to visually represent time series predictions compared to the actual values over time.


What are the assumptions of the Games-Howell multiple comparisons procedure? Please back up any answers with references.


title: Problems plotting GLM data of binomial proportional data

text: This might be a question for the programmers, but I thought I would ask here first. Im comparing browsing pressure on plants between sites. Ive a value to indicate Browsing Pressure and count data of trees that have been damaged between locations. Ive been using Crawleys R example (page 574) regarding sex ratios and my issue is with the lines command to plot the information. I can't get a fitted regression line to the points to show up for my data. At least not using the commands Crawley suggested. Ive included Crawley's script for quick reference:               density = c(1,4,10,22,55,121,210,444)     females=c(1,3,7,18,22,41,52,79)     males=c(0,1,3,4,33,80,158,365)     numbers=as.data.frame(cbind(density,females,males))     attach(numbers)     par(mfrow=c(1,2))     p<-males/(males+females)     plot(log(density),p,ylab="Proportion male")     y<-cbind(males,females)     model<-glm(y~log(density),binomial)     xv<-seq(0,6,0.1)     plot(log(density),p,ylab="Proportion male")     lines(xv,predict(model,list(density=exp(xv)),type="response"))      My Script:               data<-####see the dput() data below to insert here     attach(data)     p1<-MS1/(MS1+M1)     y<-cbind(MS1,M1)     GM1<-glm(y~BPT,family=binomial (logit),data=data)     summary(GM1)     plot(BPT,p1,col="black",pch=1,main="Relationship a",xlab="Browsing pressure",     ylab="Moose Damage Survey")      This is where I am having difficulties interpreting Crawley's `lines` command. I can create a sequence for my data like his xv values - mine go from 0 to 0.7.   `xv<-seq(0,0.7,0.01)` But as Im not sure what is happening in the command line, blindly inserting my data results in no line showing on my graph. I have some success with `regLine` from the `car` library, but its not a line fit to my data. Any help explaining why this works for Crawleys but not my data would be appreciated My DATA:               structure(list(S = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,      14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,      30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,      46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60),          BPT = c(0.195778643105884, 0.0651216427326909, 0.0199432997190648,          0.255717971810655, 0.21031730503132, 0.0380060076389563,          0.0592940237882193, 0.00425756881206537, 0.107351336677244,          0.353626952501855, 0.134358354834345, 0.0774003112699534,          0.0571275803730281, 0.102480049705623, 0.506048974065245,          0.637795327349053, 0, 0.00825369912299262, 0.123063816485164,          0.244214078077056, 0, 0.231982659809539, 0.0265824936284919,          0, 0, 0, 0.197470272555566, 0, 0.0481572367104686, 0.339072491473947,          0.0640004726042278, 0.099465609734403, 0, 0.153951519788959,          0.103791449528225, 0.22066571012681, 0.101286659375451, 0.0210055739071239,          0, 0, 0.00811999120846242, 0.0576334384710714, 0.0514027186783361,          0.00907843357263548, 0.0785740169740806, 0.0801612249152947,          0.0439884643248279, 0.189608233745472, 0.434349352638054,          0.188190501082836, 0.016381369054353, 0.215189432724754,          0.0948563822493123, 0, 0.0111745795351048, 0.0289751701436048,          0.0962515835377456, 0.199970627051051, 0.0306277369169669,          0), M1 = c(51, 123, 137, 23, 69, 98, 80, 59, 84, 87, 63,          88, 64, 75, 30, 19, 86, 152, 55, 22, 115, 66, 75, 124, 87,          179, 39, 66, 117, 37, 59, 103, 57, 65, 83, 68, 87, 104, 134,          89, 52, 61, 65, 75, 71, 97, 45, 37, 86, 82, 36, 139, 40,          56, 63, 64, 37, 79, 70, 121), M2 = c(108, 195, 255, 31, 104,          157, 135, 88, 117, 152, 96, 124, 102, 118, 37, 25, 135, 277,          97, 47, 174, 121, 134, 196, 163, 270, 75, 123, 214, 46, 109,          156, 92, 118, 127, 112, 134, 152, 194, 132, 107, 101, 105,          140, 107, 163, 87, 60, 122, 122, 66, 237, 66, 95, 108, 95,          50, 128, 128, 189)), .Names = c("S", "BPT", "M1", "M2"), row.names = c(NA,      -60L), class = "data.frame")


title: PCA with covariance matrix calculated using random matrix theory in R

text: I would like to perform a PCA and use the covariance matrix obtained by the random matrix theory. Is there an implementation of this in R? I am currently using the standard _prcomp_ function from **stats**.


title: Effects on R^2 of adding interactions to a linear regression model

text: I am trying to model a linear regression to predict the number of interested students in a class based on a couple of properties of that class-   1. whether that class has an undergraduate curriculum (assigned 0 or 1)   2. whether that class has mandatory homework (assigned 0 or 1)   3. whether that class is held for multiple weeks (assigned 0 or 1)   4. the difficulty of that class (assigned a number from 1.0 to 4.0) I have a little over 80 data points, but the R^2 values I get vary upon whether I take the products of certain variables.               summary(lm(AA.Int~undergrad + multi + HW + Spiciness, data=realdata))      returns an R^2 value of 0.06419 whereas               summary(lm(AA.Int~undergrad * HW * multi * Spiciness, data=realdata))      returns R^2 = 0.4231 And similarly, all summaries involving these "products" seem to give me higher R values. Is this a good thing? How should I interpret my data? Finally, should I model my data using something other than a simple linear model? In case it will help, here is my data:                  undergrad multi HW AA.Int Spiciness      1          0     0  1   6.00       4.0      2          0     1  1  27.50       2.5     3          1     0  1  54.50       2.0      4          1     1  0   7.50       2.0      5          1     1  1  16.50       2.0      6          1     0  0  15.00       3.0      7          0     0  0  18.00       3.0      8          0     0  0  29.00       1.5      9          0     0  0   9.50       3.0      10         0     1  1  33.83       3.0      11         0     0  0  22.93       2.0      12         0     0  1  12.93       1.0      13         0     0  1  13.30       3.0      14         1     1  1  13.20       2.0      15         0     1  1  16.00       3.5      16         0     0  0  21.50       1.0      17         0     0  0  21.10       3.0      18         0     1  0  19.50       1.5      19         0     0  0  12.20       3.0      20        NA    NA NA     NA        NA      21         0     1  1  25.30       2.5      22         0     0  0  26.20       3.0      23         0     0  0   3.50       4.0      24         0     0  0  21.40       3.5      25         0     1  0  10.00       2.0      26         1     1  1  12.00       2.0      27         0     1  0  20.80       3.0      28         1     1  0  20.00       3.0      29         1     0  0  20.50       3.0      30         0     0  0  16.50       2.0      31         0     0  0  15.90       2.0      32         0     1  1  37.80       3.0      33         0     0  0   4.00       3.0      34         0     0  0  17.60       1.5      35         0     0  1  11.80       3.5      36         0     0  0   4.60       1.0      37         0     0  0  18.40       3.0      38         0     1  1  13.05       3.0      39         1     1  1  12.50       2.0      40         0     1  1   6.50       3.5      41         0     0  1  15.95       3.5      42         0     1  0  16.00       1.5      43         0     1  1   9.50       3.0      44         0     0  0  13.30       3.0      45        NA    NA NA     NA        NA      46         0     0  0  18.90       3.5      47         0     0  0  21.80       2.0      48         0     1  0  24.70       2.0      49         0     1  0  19.60       3.0      50         0     0  0  20.30       1.0      51         0     1  0   5.80       4.0      52         0     0  0  20.50       2.0      53         0     0  0  13.90       4.0      54         0     0  0  16.00       2.0      55         0     0  0   9.80       2.5      56         0     1  1  29.00       3.0      57         0     0  0  14.30       3.0      58         0     0  0  17.00       2.0      59         0     0  0  12.70       3.5      60        NA    NA NA     NA        NA      61         0     1  0   6.70       3.0      62         0     1  0  15.00       2.5      63        NA    NA NA     NA        NA      64         0     0  0   7.83       1.0      65        NA    NA NA     NA        NA      66         0     0  0  30.83       2.0      67         0     0  0  14.80       4.0      68         0     0  0  18.70       2.5      69         0     1  0  16.30       3.5      70         0     0  0  11.60       2.0      71         0     0  0  37.80       2.0      72         0     1  0   9.40       4.0      73         0     0  0  12.80       1.0      74         0     1  0   4.00       4.0      75         0     0  0  11.80       3.0      76         0     0  0  20.50       2.5      77         0     0  0  23.90       2.0      78         0     0  0  11.00       3.5      79         0     1  0  12.50       4.0      80         0     0  1  25.60       3.0      81         0     0  0  16.20       1.5      82         0     0  0  17.80       2.5      83         0     0  0   4.50       2.5      84         0     0  1   7.90       1.0      85         0     0  1   5.50       4.0      86         0     0  0  15.85       3.0      87         0     0  0  11.50       1.5


title: VECM , VAR or Granger causality

text: In my study, the variables are cointegrated. Now I want to test causality so should I perform VAR, VECM or the Granger causality test?


title: How to compute the probability associated with absurdly large Z-scores?

text: Software packages for network motif detection can return enormously high Z-scores (the highest I've seen is 600,000+, but Z-scores of more than 100 are quite common). I plan to show that these Z-scores are bogus. Huge Z-scores correspond to extremely low associated probabilities. The values of the associated probabilities are given on e.g. the normal distribution wikipedia page (and probably every stats textbook) for Z-scores of up to 6. So... > **Question** : How does one compute the error function > $1-\mathrm{erf}(n/\sqrt{2})$ for n up to 1,000,000, say? I'm particularly after an already implemented package for this (if possible). The best I've found so far is WolframAlpha, which manages to compute it for n=150 (here).


title: Weighted loss function for non-random sample

text: When comparing a regression estimation method (Y vs X) I currently use a weighted squared loss function: $$ \int_{-\infty}^{\infty}(\hat{f}(x)-f(x))^2 \, \hat{p}(x) \, dx $$ Where $\hat{f}(x)$ is the estimated function, $f(x)$ is an estimation method I'm comparing against, and $\hat{p}(x)$ is the estimated X density. This works well when I assume the X observations are randomly sampled from the true density. However if X is not a random sample then weighting by the estimated X density does not make sense. If I have some knowledge that some intervals of X are oversampled and some intervals are undersampled how should I adjust the weights of my loss function?


title: Installing R on Ubuntu

text: I don't know if this is the right StackExchange for this question, but I am using Ubuntu 12.04 and I am trying to install R, but all the guides I get are really painful and involve long and complicated descriptions of guides I can't find my way around in. Does anyone know an easy guide to install R that I can just follow the steps and get things done?


title: Generating sorted pseudo-random numbers in Stata

text: Today I opened two STATA windows and ran the following command in both:               set obs 100     gen x = rnormal()     sort x      (the difference is that on the second window I generated a variable called y). Summing up: I asked STATA to give me 100 pseudo-random numbers taken from a standard normal distribution, then I sorted it. To my surprise, the numbers of the x and y vectors are the same! I did this at home, and then at work, and my impression is that all of these vectors are the same. Is there an explanation for this, to me, strange behavior? If this is a problem in STATA, does R have a better pseudo-random number generator procedure? A side question. I came up to this "problem" because I was trying to generate two pseudo-random columns in Stata (x and y, say), and then sort then separately. But the two commands I know for sorting (sort and gsort) sort the whole database, not separate columns. Would you know of a Stata command that allows me to sort a column while keeping the other columns fixed?


title: Should we compare the individual monthly forecasts with actual values?

text: Hi I am using Linear and exponential forecasting models to do sales forecasting. In the model itself, we use the forecasts of period t to get next forecast and so on. While analyzing the accuracy of the forecast using Mean Absolute Percentage Error, I get good results. But when I compare the intermediate forecast values of the model against the actual time series data, I can see some big Percentage Errors like MAPE might be 15%, but some of the intermediate percentage errors might be as high as 80% and some will be pretty low. So the MAPE averages out to low. I wanted to know 1) Is it wise to compare these intermediate forecasts with actual values? 2) Can we use these high intermediate forecast values to say that forecasts for so and so month might be unreliable?


title: Can the forecasts using exponential smoothing be negative in value?

text: Hello I am trying to forecast using different exponential smoothing methods(Linear and Winter's). For the optimal parameters, I am getting negative values of the forecasats. I am assuming it means that the values will be zero, since it is a sales forecast. I wanted to know if negative values denote something wrong with the model or is it possible to have negative values in forecast. Since it is just a model, I think we can get negative values for certain type of time series.


title: Assessing differential probability of an event across a large number of plant observations

text: I have 4 files, each containing the observations for a different plant. The observations are the probability of occurrence of a particular "event" over many intervals over the plant's genome. That is, the observation file looks like this:               file/plant 1:       ch01:a1-b1       p = 0.41       ch01:a2-b2       p = 0.005       ...      And in the same format, over the same intervals, I have the probability of occurrence of this event for other 3 plants/files. There are a total of 15,000 observations in each file over identical intervals of chromosomal locations. I am interested in finding those intervals where the event is "significant" in 1 or more plants, but NOT ALL; i.e., events that stand out (or different among the plants). If the event is equally likely in all 4, then it of course is not interesting. So, it would be great if I could get some insight on how to go about this. The things I have managed to find out so far are this: ### 1 two-sample t-test I could perform a two-sample t-test on all pairs, but maybe it will be biased. So I decided to skip this. ### 2. F-test I could perform a F-test, but this is (also) a parametric test and so I will have to find out whether the distribution of residuals is normal (I do not know why the residuals have to be normally distributed, but was told to do this. So it would be nice if someone could explain the logic behind this or share a link that explains it). So, I decided to find out the distribution. But I am not sure how. I use R. If I understand it right, for each of the chromosomal interval (or observation) over all 4 plants, lets say, for example, for observation ch01:a1-b1, p = (0.004, 0.33, 0.19, 0.021). Then my understanding is that, the residuals of these values must be normally distributed. If not, a parametric test is not possible. I am not used to R. I was told to use the "lm" syntax followed by "anova". I would like to know if my understanding and this approach is right? If so, how do I find the distribution of the residuals? And if this is the way, then I have to find the distribution for all 15000 observations?? ### 3. Non-parametric method If F-test is not the method to proceed, then I guess a non-parametric method is the way to go. In this case, which test should I adopt to test the significance of each event over all 4-plants. I guess its multiple (hypothesis) testing? ### Question Recap   * How do I check that the residuals are normally distributed in R?   * Which test should I use?   * If I use a non-parametric test, which test should I use?


title: Deriving $\beta$ from expected loss function (Hastie, Tibshirani)

text: I am looking at equation (2.16) from Elements of Statistical Learning and can't seem to be able to derive it. I used $f(x) = x^T\beta$ as the linear model, and tried calculating $\beta$ by minimizing the expected prediction error (EPE). $$ EPE(f) = E[(Y-f(X))^2] $$ $$ EPE(f) = E[(Y-X\beta)^T(Y-X\beta)] $$ $$ \partial EPE(f)/ \partial \beta = E[X^T(Y-X\beta)] = 0 $$ $$ \partial EPE(f)/ \partial \beta = E[X^TY]-E[X^TX\beta)] = 0 $$ $$ E[X^TY] = \beta E[X^TX] \text{ (as beta is constant)} $$ $$ \beta = [E(X^TX)]^{-1}E(X^TY) $$ Now my result is different from equation (2.16) in the book and I am clueless why is that the case. I am guessing there is a linear algebra rule I am missing. Why is my $X^T$ their $X$ and vice-versa? Any word will help!


title: The effect of the number of replicates in different cells on the results of ANOVA

text: I have written 3way-anova in c++. I have 3 factors, lets say a,b and c and my aim to check the strength of all possible interactions and main effects.The result of my code is the same as MATLAB when I use type I sum of squares.   But when I change the data so that the number of replicates is so high in some cells and too low in others: I don't get the same results as MATLAB. (to be precise, only SSt, SSe and SSa are the same as MATLAB) . My question is that is it possible that since I have a large difference in the number of replicates, I should use type III sum of squares? Or is there is special way that matlab treats the data in such cases so its results differs from mine? Thanks for your help. Pegah


title: Statistics for gambling machine validation

text: Problem is that government wants to close electronic roulette and they claim that roulette failed at statistical test. Sorry for my language but this is translated from Slovenian law as good as possible Official (by law) requirements are:   * frequency of each event should not differ from expected frequency by more than 3 sigma   * chi square test of normal distribution has to be within risk level of 0.025   * test of consecutive correlation has to pass 3 sigma test and chi squared test I have tested first 2 requirements and they pass tests, but I have problems whith understanding 3rd requirement. (keep in mind that this is translated and "consecutive correlation" can be something else) How should I test 3rd requirement? Data if somebody is interested:   http://pastebin.com/ffbSKpr1 EDIT: chi squared fails 2% of the time (what I expect that is expected due to fact that alpha is 0.025) and sigma3 test fails 5% where i expect 9% failure for 3sigma (it looks like that frequencies are not distributed according to normal distribution even for random numbers) I might not understand this law correctly, but it is almost 0% probability to pass 3sigma test for all autocorrelation vectors, since it is 9% probability to fail in single run and 2.5 for chi squared test. Python code:               from math import sqrt     from itertools import *     import random          #uncoment for python 2.x          #zip = izip          #range = xrange               #with open("rng.txt","r") as wr:     #   n = [int(i) for i in wr]     n = [random.randint(0,36) for i in range(44000)]                def get_freq(n):         r=[0 for i in range(37)]         for i in n:             r[i] += 1         return r          def trisigmatest(freq):         Ef = 1.0*sum(freq)/37         sigma = sqrt(sum(i**2 for i in freq)/37-Ef**2)         return all((abs(i - Ef )< sigma*3) for i in freq)               def chiquaretest(freq):         Ef = 1.0*sum(freq)/37         chi2 = sum((i-Ef)**2 / Ef for i in freq)         # values are from http://itl.nist.gov/div898/handbook/eda/section3/eda3674.htm         # (EDIT) I recaluclated these valuse from inverse cdf chi2          # distribution for interval (0.025/2,1-0.025/2) (alpha = 0.025)         return 20.4441 < chi2 < 58.8954                  #whitout autocorelation       gf = get_freq(n)     if not trisigmatest(gf):         print("failed")         raise     if not chiquaretest(gf):         print("failed")         raise               actests = 1000     trifailed = 0;     chifailed = 0;     for i in range(1,actests + 1):         f=((b-a+37) % 37 for (a,b) in zip(n,n[i:]))         gf = get_freq(f)         if not trisigmatest(gf):             trifailed += 1;         if not chiquaretest(gf):             chifailed += 1;     print("trisigmatest failed ", 1.0 * trifailed / actests )     print("chiquaretest failed ", 1.0 * chifailed / actests )


title: Confusion related to a derivation in a paper

text: I was reading this paper: > Ristovski, K., Das, D., Ouzienko, V., Guo, Y., Obradovic, Z. Regression > Learning with Multiple Noisy Oracles. This paper is basically about regression. The idea is to learn a regression model using predictions made by several noisy predictors and learn the best regression model combining the predictions from the predictors. I am actually confused a bit in the equations. In the paper in equation 24, they have given the regression model as $$\begin{eqnarray} f(x^{(i)},w) & = & w^{T}\phi(x^{(i)}) \\ & = & a^{T}\phi\phi(x^{(i)}) \\ & = & k(x^{(i)})^{T}(K+LI_N)^{-1}t. \end{eqnarray}$$ I didn't get how the last line in the above derivation came up. Also it has mentioned that for a new input $x$ $$y(x) = k(x)^{T}(K+LI_N)^{-1}t,$$ where $k(x)$ is the vector with elements $k_i(x) = k(x^{(i)},x)$. The kernel function used in the Gaussian kernel is $$k(x^{(i)},x^{(j)}) = e^{-\frac{(x^{(i)}-x^{(j)})}{2{\sigma}^2}}.$$ I didn't get it. I mean lets say I have a new input $x=(x_1,x_2,\dots,x_m)$ where $m$ is the number of features in an input. How am I supposed to calculate $k(x)^{T}$? I didn't get this $k$ thing. Can anyone please explain? I will be really grateful.


title: Calculating moments of a binomial distribution

text: I have the following problem: If $X$ satisfies a binomial distribution with parameters $n$ and $p$, its expected value is $np$ and variance is $np(1-p)$. Prove these using formulas given in pages 8 and 10. On page 8 we have been given that $E(X)=\sum_{i=1}^k P(X=x_i)x_i$, $\operatorname{Var}(X)=\sigma^2=\sum_{i=1}^kp_i(x_i-\mu)^2=\sum_{i=1}^kP(X=x_i)(x_i-\mu)^2$ and on page 10 we have formulas for $E(a+bX), \operatorname{Var}(a+bX), E(X+Y), \operatorname{Var}(X\pm Y)$. I already managed to find a proof using the identity $x^2\binom{n}{x}=xn\binom{n-1}{x-1}$ which does not use the formulas above but as we haven't given that identity in the course, I was wondering whether there is an alternative proof.


title: How to compute for Bivariate Logistic Distribution

text: This is the logistic distribution of single random variable (taken from Wikipedia). $x$ = random variable $\mu$ = mean of all random variables $s$ = variance. ![enter image description here](http://i.stack.imgur.com/qUHkD.png) Now, I want to do a Bivariate logistic distribution (having two random variables $x_1$ and $x_2$). My dataset is going to be image pixel values! When I find the covariance of two random variables, it turns out to be a 2 x 2 square matrix but I need a single number! How to actually compute a bivariate model?


title: Regression Puzzle GLS var/covar matrix of B given summary statistics

text: There is a problem posted back in 2009 that I found that has me puzzled. I believe the estimates are from a generalized linear regression because after some thought it seems $X'X$ is impossible to obtain from the given info. So it must be something to do with the Cholesky decomposition of $\Sigma$, but I do not see the answer. The problem is posted at http://www.unc.edu/~jjharden/methods/vcv_week3.pdf, but a friend of mine had trouble with the link so i provide it here.. ![enter image description here](http://i.stack.imgur.com/3zNmP.png) Here is $R^2$, but the var/cov matrix eludes me.               ########Given info     n=15     Ybar = 68.7333     b   = c(5.5027, 4.1535, 3.3787, 8.2903); dim(b) = c(4,1)     yty = 259651      xty = c(1031, 11370, 35981, 10176); dim(xty) = c(4,1)               ########Derived info     TSS = yty - n*Ybar^2     EXPSS = t(b) %*% xty - n*Ybar^2     RESSS = yty - t(b) %*% xty - n*Ybar^2     R2 = EXPSS / TSS     sig2 = RESSS / (n-4)


title: In spatial regression, what is a spherical autocorrelation structure?

text: I have a large gridded dataset for the globe (i.e a spherical, wraparound surface) that I'm applying spatial regression to (using a CAR model). I've been using the default autocorrelation function, however as my data is global (point 0,0 lies next to max(x),0 and 0,max(y)), I wondered if a spherical autocorrelation function would work better. Though I've alot of references on its use, I've been unable to find a simple answer or diagram as to what a spherical autocorrelation function actually is! Is it as I presumed, an autocorrelation where data wraps around a sphere, or a specific mathematical function of the variogram shape, or something else entirely?


title: How to report  Games-Howell post hoc tests following ANOVA?

text: I have run a 1 way ANOVA and carried out some post hocs to see which direction the results go in. If I report the direction of the findings - e.g. this group was higher than the other, which statistic do I quote? I have carried out Games-Howell post hoc and all it seems to come with is an associated p value, is this all that I report?


title: Calculate the quantile for a mixed pdf

text: I am trying to draw random numbers that follow a two-part pdf a) $|x|< x_0: f(x)=\text{constant}\quad \to \quad F(p)=a+(b-a)p$ b) $|x|> x_0: f(x)=\exp(-|x|)\quad \to \quad F(p)=-\ln(1-p)$ that is the probability is constant between $(-x_0,x_0)$ and falls exponentially outside it. Knowing the quantile functions $F(p)$ ($0\leq p\leq 1$) of the pdf, how can I construct the quantile that describes the above pdf comprised of two parts?


title: Rate changes over time

text: I need two sets of data collected over 3 years. The first data has 12 waves where each wave represents a quarter. The other has 24 waves where each wave include some part of the previous wave (e.g. 1st wave 1/1/07-4/30/07; 2nd wave 3/1/07-6/30/07 ...) The variables measured were rates (specific lab tests taken from certain patient populations). The data were collected to examine the effect of new reporting system. So finding whether there is a significant increase of test rates over time or not is my interests. I don't have 0 rates in the data.   * Would negative binomial regression appropriate for this type of analysis?   * Or if I take only the first and last waves to compare the differences in rates, what would be the best way?


title: Random Forest regression model in R and data overfitting

text: I have trained my random forest model on a 74,000 training examples where each example consists of two proteins Amino Acids sequence (20 characters) and some numeric values representing the similarity between each individual pair of sequences, and finally a numeric value representing the overall similarity between the two proteins, this is a regression model, so I wish for testing I can use just protein sequence and use my trained model to predict the distance between my test case and each of my training protein sequences. A sample of my test case is:                  test= G,Y,L,P,P,S, A,N,L,F,S,N, 1,-2,16,-4,-1,11, 21      where "G,Y,L,P,P,S," represent a 6 character fragment of the first protein (my testing) and "A,N,L,F,S,N," represent a 6 character fragment of the second protein (my training database) and the numbers "1,-2,16,-4,-1,11," each number represent the similarity between individual pairs of the 6 Amino Acids, e.g., the similarity between "G and A" is 1 and the similarity between "Y and N" is -2 and the similarity between "L and L" is 16, and so on offcourse the higher the number means the higher the similarity between the pairs of characters, finally the last number "21" represent the sum of the previous 6 numbers which represents the overall similarity between the two sequences. when I trained my model on 74,000 of such training datasets the correlation between the predicted distance and actual distance was as high as 0.86, however, when I used the trained model in for testing the correlation was very low 0.17, I strongly believe that this over-fitting problem, however, I'm not sure is it due to the may be not good training datasets or that my features aren't strong enough to give a good prediction especially since the ranking of the features according to their importance was very high for all the features? the following is my features importance according to each node purity. any help on how to recognize the source of the overfitting is highly appreciated:                      IncNodePurity      V3      24564.326      V4      22503.744      V5      25030.450      V6      24583.235      V7      24661.309      V8      20757.662      V9      22985.824      V10     22189.759      V11     23875.170      V12     23674.853      V13     23339.595      V14     19576.762      V15     10169.309      V16     19527.972      V17      5430.600      V18      4415.307      V19     12897.114      V20      3963.717      V21     62614.692


title: I am looking for a methodology for categorizing keywords

text: I.e. I want to put in Baseball and get out Sports & Outdoors.


title: Evaluation of methods under the null by p-value histogram

text: I am using **real data** for a comparison of two groups: there are 3 replicates in one group and 3 replicates in the other group. All 6 replicates come from the control group, so that their responses should, ideally, be very similar, and a test should report very few significant results. Now I got the results of raw p-values for the two-group comparison, and plot the histogram. Unfortunately, the histogram of raw p-values are not **uniform- like** – there are more small raw p-values there than expected. However, when I apply the Benjamini-Hochberg FDR control and check the q-values, there is 0 significant result! I am confused here as I am **evaluating the method under the null** – from the histogram of raw p-values, the method seems to report many false positives, but from the adjusted p-value (q-value) results, it seems good as I shouldn't report any significant result… I admit that this is based on real data, not simulation data, so I cannot guarantee that what I am comparing are in reality "strictly under the null".


title: What's the easiest way to separate two populations in a scatterplot?

text: I have to separate two populations by a line in a scatterplot: ![Dot plot](http://i.stack.imgur.com/B7NpS.png) I would like find a threshold that separates the two populations. In @Waynes words, I would like to cluster the points into two categories, then calculate the line which best separates them. What is the easiest way to accomplish this in `R`? I've tried the Otsu algorithm on the raw data, but I don't know if this one is robust enough for other datasets. **Edit:** I don't know which dot belongs to which population. The data is not multivariate, I only have an X and a Y value, as seen in the dotplot.


title: Comparable traing and test cross-entropies result in very different accuracies

text: # Premises I'm training a _convolutional neural network_ (ConvNet) on 51 subclasses in the _ImageNet_ dataset. In order to keep an eye on _overfitting_ , I have been suggested to plot training and testing **loss function** values (using a _negative log-likelihood_ criterion) and **accuracy** (correct guesses over total number of sample). So far, I have obtained consistent results, similar to the example below. ![Clear overfitting](http://i.stack.imgur.com/eC3mx.png) Then, I killed the overfitting by using a _dropout_ technique, producing the following chart. ![Killed overfitting with dropout](http://i.stack.imgur.com/FoiiZ.png) # Question Now I'm puzzled. I have comparable cross-entropy errors, but quite different accuracies. Shouldn't a similar loss function value, for train and test, give me similar precisions? Rewording the question: why $e_{\mathrm{test}}\simeq e_{\mathrm{train}}\not\Rightarrow a_{\mathrm{test}}\simeq a_{\mathrm{train}}$ where $e$ stays for _cross-entropy error_ (defined below) and $a$ stays for _accuracy_? # Definitions I've been using the terms _cross-entropy errors_ and _loss function values_ interchangeably. What I am referring at is the average prediction error which is given by $$e = - \frac{1}{m}\sum_{i=1}^{m}{\ln[p(y_i|x_i)]}$$ where $m$ is the number of images in the training or testing sets, $y_i$ is the true label and $p(y_i|x_i)$ is the model's output probability associated to the input $x_i$ to be classified as $y_i$.


title: Over dispersed Poisson model fit comparison

text: I am comparing various GLMs fit to a series of data (loss development actuarial triangles) and I am wondering how to rank Over Dispersed Poisson Models (ODP) against other GLMs (e.g. Negative Binomial or Gamma Glms) since AIC is not defined for ODP and I'm not sure deviances are comparable.


title: Validating correctness of ranking algorithm

text: Stats noob here working on an algorithm for ranking players in a tournament. I have a test set of players, each with a "skill" value. My algorithm simulates a "tournament", in which players play against each other (1v1) in some order (e.g. Swiss pairing). At the end of the tournament, I have each player's resulting rank. For example, a sample data set could be something like:               Name | Tournament Rank | "True" Rank     Alex |        5        |      7     John |        1        |      2     Mike |        3        |      1     ...      So John won the tournament even though he was the second best player there, and Mike is actually the best player, but happened to finish third due to luck or whatever. Given each player's "true" skill and their rank in the tournament, how can I best quantify how well my tournament algorithm performed in ranking everyone? e.g. I'm looking to be able to say things like "With this set of input parameters, my tournament simulator places people 10% more accurately than with this other set of input parameters" (Bonus question:) Also, arguably more important than the exact ranking to me is that players are in the correct "bucket" of players. For example, if I split the results into 5 sections (top 20%, 60th-80th percentile, etc.), I need my tournament algorithm to reliably place people into the bucket they deserve to have ended up in. Would you do anything differently than above to check the correctness of how people end up in buckets?


title: Differences between hierarchical regression and stepwise regression analysis

text: I have to run an analysis using different steps, to include first the control variables, second the independent variables, and finally the interaction terms (moderation). However, I don't know if I should use hierarchical regression or stepwise. Could anyone explain me what are the differences between both methods? When do we use each one of them? Thanks a lot!!