title: Mixed/hierarchical/multilevel model for n-of-1/single-subject design: incorporating baseline data?

text: So I often do little self-experiments where I blind & randomize things; these can be formulated as your normal _t_ -tests, but sometimes the measured metrics have extensive baselines which seem like they could be used for more accurate answers. A bunch of reading upon _n_ -of-1 and single-subject designs suggested that people have been moving to mixed/hierarchical/multilevel models for analyzing such setups (eg. Nelson 2012 "Hierarchical linear modeling versus visual analysis of single subject design data" or "The Aggregation of Single-Case Results using Hierarchical Linear Models"). As I understand it, the idea is to split the subject's data into experiment vs baseline, and treat those as the groups. I'm trying to understand how sensible this is with a recent experiment, so hopefully someone can point out if I go wrong in using `lmer` here. * * * We start with a regular linear model which examines purely the experimental data (the numeric `Response` vs the binary `Intervention` variables) and ignores the extensive baseline phase before, during, and after the experiment:               R> experiment <- read.csv("http://dl.dropboxusercontent.com/u/85192141/data.csv")     R> summary(lm(Response ~ Intervention, data=experiment))          ...     Residuals:         Min      1Q  Median      3Q     Max     -1.0156 -0.8889 -0.0156  0.1111  1.1111          Coefficients:                  Estimate Std. Error t value Pr(>|t|)     (Intercept)    3.0156     0.0889    33.9   <2e-16     Intervention  -0.1267     0.1262    -1.0     0.32          Residual standard error: 0.711 on 125 degrees of freedom       (145 observations deleted due to missingness)     Multiple R-squared:  0.008, Adjusted R-squared:  6.73e-05     F-statistic: 1.01 on 1 and 125 DF,  p-value: 0.317          R> confint(lm(Response ~ Intervention, data=experiment))                    2.5 % 97.5 %     (Intercept)   2.8397  3.192     Intervention -0.3765  0.123      The estimated coefficient is not statistically-significant: -0.38-0.12. But it's definitely slanted towards being negative. So this is the 'conservative' case, where we ignore the baseline entirely. What's the optimistic case? Well, it seems to me that the optimistic case is when we take the entire baseline and assume it is exactly the same as the 'off'/0 intervention in the experiment, in which case we get a narrower CI (because our estimate of the intercept has halved its standard error):               R> experiment$Intervention[is.na(experiment$Intervention)] <- 0     R> summary(lm(Response ~ Intervention, data=experiment))          ...     Residuals:         Min      1Q  Median      3Q     Max      -1.9924 -0.8889  0.0076  1.0076  1.1111           Coefficients:                  Estimate Std. Error t value Pr(>|t|)     (Intercept)    2.9924     0.0375   79.88   <2e-16     Intervention  -0.1036     0.1012   -1.02     0.31          Residual standard error: 0.746 on 458 degrees of freedom     Multiple R-squared:  0.00228,   Adjusted R-squared:  0.000101      F-statistic: 1.05 on 1 and 458 DF,  p-value: 0.307          R> confint(lm(Response ~ Intervention, data=experiment))                    2.5 %  97.5 %     (Intercept)   2.9188 3.06607     Intervention -0.3025 0.09538      It's narrowed to -0.30-0.10; still not statistically-significant, but closer. It seems to me that a hierarchical model ought to produce a CI intermediate between the pessimistic and optimistic cases: it loses some power because it's estimating how different the two phases are before it does any combining. Here is my multilevel model, split between baseline and experimental phases:               library(lme4)     experiment <- read.csv("http://dl.dropboxusercontent.com/u/85192141/data.csv")     experiment$Phase <- ifelse(is.na(experiment$Intervention), TRUE, FALSE)     model <- lmer(Response ~ Intervention + (1|Phase), data=experiment); summary(model)          ...      AIC BIC logLik deviance REMLdev      286 297   -139      273     278     Random effects:      Groups   Name        Variance Std.Dev.      Phase    (Intercept) 0.0106   0.103         Residual             0.5057   0.711        Number of obs: 127, groups: Phase, 1          Fixed effects:                  Estimate Std. Error t value     (Intercept)     3.016      0.136    22.2     Intervention   -0.127      0.126    -1.0          Correlation of Fixed Effects:                 (Intr)     Interventin -0.461          m <- mcmcsamp((lmer(Response ~ Intervention + (1|Phase), data=experiment)), n = 100000)     HPDinterval(m, prob=0.95)$fixef                          lower   upper     (Intercept)  -45.3107 56.6558     Intervention  -0.3742  0.1191      The estimated CI comes out exactly in the middle, as expected:   1. pessimistic -0.38 0.12   2. hierarchical -0.37 0.11   3. optimistic -0.30 0.10 So, my basic question is: is this a sane approach to take? It's spitting out answers that seem intuitively correct, but that might just be a coincidence. * * * Incidentally, one might be worried about time trends. The randomization/blocking would fix that in the experimental period but not the baseline. Fortunately, that doesn't seem to be an issue:               experiment$Time <- 1:nrow(experiment)     summary(lmer(Response ~ Intervention + Time + (1|Phase), data=experiment))          ...      AIC BIC logLik deviance REMLdev      298 312   -144      272     288     Random effects:      Groups   Name        Variance Std.Dev.      Phase    (Intercept) 0.0106   0.103         Residual             0.5055   0.711        Number of obs: 127, groups: Phase, 1          Fixed effects:                  Estimate Std. Error t value     (Intercept)   3.42517    0.42325    8.09     Intervention -0.12398    0.12621   -0.98     Time         -0.00132    0.00129   -1.02          Correlation of Fixed Effects:                 (Intr) Intrvn     Interventin -0.128            Time        -0.947 -0.021


title: Is this a stationary time series?

text: I have a wait-time time series for 10 weekdays(2 weeks) with 10 minutes intervals. I'm having hard time to interpret this? Is this stationary? I also applied Philips-Perron Unit root test and I got the following values. Any help is appreciated. Dickey-Fuller = -8.4726, Truncation lag parameter = 5, p-value = 0.01 ![enter image description here](http://i.stack.imgur.com/mTuNi.png)


title: GLMM for SNA and non-independency data

text: I contact you because my case is particular and I don’t know much about GLMM. I have data of social networks (network metrics) of a nonhuman primate species. These data are by nature non independent (behaviour between A and B can influence those of C). I also have individual attributes (gender, age, kinship and hierarchical rank). I would like to test the effect of individual attributes on the network metrics. For this purpose I have to take into count the non-independency of data (I could make link filtering or permutation test but I am trying something new). I was thinking about doing a General Linear Mixed Model (GLMM) including individuals as random factors. I see that this analysis can take into count the non-independency of data, but I only see this test for multiple measures on same individuals (which is another kind of non- independency) and not for a single measure on individuals. Do you think even if this is not the same type of non-independency data GLMM can take into count my type of data? Regards


title: Inner-line correlation (consistency)

text: In this paper, I found a quite interesting approach concerning measurement of consistency in survey data (provided that each line=case represents a participant). Still, I cannot figure out exactly what is meant by "within- person correlation". Could anybody help, please? > Consistency indices can be formed by examining the differences in two items > that are highly similar in content. Goldberg, 2000 (cited in Johnson et al., > 2005) suggested a method called Psychometric Antonyms in which correlations > among all survey items are computed post-hoc and 30 item pairs with the > largest negative correlation are identified. The Psychometric Antonyms index > is then computed **as the within-person correlation across these thirty > pairs of items**. My interpretation goes as follows: First, I identify which variables in a sample are highly (negatively) correlated (based on data from all participants). Second, for each single participant (n): I take the particular values of these pairs from that single participant, note them down in a separate table (pairwise) and then I compute a correlation over these two columns. This would mean that I would create n new tables and n new correlations. Is that correct, or am I doing it wrong? Thanks so much for help.


title: Need help comparing graphs

text: I have data from Google Trends and Center for Disease control that I want to compare. I used regression analysis from Excel to find a curve of best fit. The original data is in this form: ![enter image description here](http://i.stack.imgur.com/xEX7p.jpg) It was impossible to entirely model this using excel so I broke it down into two parts and found a curve of best fit for those two. ![enter image description here](http://i.stack.imgur.com/8m2v4.jpg) The goal is to compare this data with another graph I made with Google Trends (same process, break it down and regression) and assess to what extent is it accurate to use search engine queries to predict flu seasons. If I want to find the graph of the derivative of the original data to look at rates of change, can I take the derivatives of the two smaller functions and compile them together?


title: Non-linear regression, plot mean, SE in R

text: I need to create a graph like the one below.The idea is that I have three data groups with 5 replications on the first two groups and 3 replications on the last one. I must obtain a regression curve for each group and link it to the mean and SD for each group, like the image below. I chose to do the analysis in R; I can get the models but have no idea about how to link it to the mean and SD per group. If any of you has some tutorial, example or code which can help it would be really appreciated. ![enter image description here](http://i.stack.imgur.com/zFicp.jpg) Thanks in advance.


title: How to simulate rare events with extremely low probability

text: Suppose I have a time series with $2\times 10^8$ data points. One kind of rare event will happen at each time point with an equal and extremely low probability of $1\times 10^{-8}$.   Can anyone tell me how I can do this?


title: Finding conditional pdf of $Z$ given $X=x$ where $Z=X+Y$

text: Finding conditional probability density function (pdf) of $Z$ given $X=x$ where $Z=X+Y$ and $X$ and $Y$ are independent normal random variables. Trivially $Z$ is normal. But $Z$ and $X$ are not independent. When means and variances of $X$ and $Y$ are given, how could I get the conditional pdf of $Z$ given $X=x$? I have tried to calculate using the definition of pdf of normal random variable, but it was so messy. Any hint for it?


title: How to understand output from R's polr function (ordered logistic regression)?

text: I am new to R, ordered logistic regression, and `polr`. The "Examples" section at the bottom of the help page for polr (that fits a logistic or probit regression model to an ordered factor response) shows               options(contrasts = c("contr.treatment", "contr.poly"))     house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)     pr <- profile(house.plr)     plot(pr)     pairs(pr)        * What information does `pr` contain? The help page on profile is generic, and gives no guidance for polr.   * What is `plot(pr)` showing? I see six graphs. Each has an X axis that is numeric, although the label is an indicator variable (looks like an input variable that is an indicator for an ordinal value). Then the Y axis is "tau" which is completely unexplained.   * What is `pairs(pr)` showing? It looks like a plot for each pair of input variables, but again I see no explanation of the X or Y axes.   * How can one understand if the model gave a good fit? `summary(house.plr)` shows Residual Deviance 3479.149 and AIC (Akaike Information Criterion?) of 3495.149. Is that good? In the case those are only useful as relative measures (i.e. to compare to another model fit), what is a good absolute measure? Is the residual deviance approximately chi-squared distributed? Can one use "% correctly predicted" on the original data or some cross-validation? What is the easiest way to do that?   * How does one apply and interpret `anova` on this model? The docs say "There are methods for the standard model-fitting functions, including predict, summary, vcov, anova." However, running `anova(house.plr)` results in `anova is not implemented for a single "polr" object`   * How does one interpret the t values for each coefficient? Unlike some model fits, there are no P values here. I realize this is a lot of questions, but it makes sense to me to ask as one bundle ("how do I use this thing?") rather than 7 different questions. Any information appreciated.


title: Scalable multinomial regression implementation

text: Greetings. I need to do high dimensional biological data analysis, where my data can easily exist in hundreds of thousands of dimensions. I am looking for an implementation of multinomial logistic regression that will scale well to data of this size. Ideally, it should allow me to also do Ridge and Lasso regressions also. Which software should I be using? Thanks. Andy


title: Confusion regarding random walk model

text: I was referring to this book where it is given that If we assume equally spaced nodes $i$ for $i=1,...,n$. The first order random walk is constructed using independent increments $$ \Delta{x_i} \sim N(0,k^{-1}), \ \ \ \ \ \ i=1,...,n-1 $$ The density for $${\bf{x}}=(x_1,x_2...x_n)'$$ is derived from the $n-1$ increments as $$ \pi({\bf{x}}|k) = k^{(n-1)/2} \cdot \exp \left(-k/2 \cdot \sum_{i=1}^{n-1}\Delta{(x_i)^2}\right) $$ I didn't get how the density of $x$ was represented in the form of increments. I mean how we can use increments to show the density how are they related or can be used in the density of $x$.


title: Problem with multiple group CFA with censored variables using Mplus

text: I'm having problems running a multigroup CFA in Mplus for my dissertation, based on the established 4-factor model of the Hamilton Depression Scale (21 items). My dataset includes 4 left censored variables which I think may be the problem. After the original CFA failed to output normally, I've rerun the analysis increasing the number of iterations, as advised by Mplus, but at `Miterations = 2000` the output is stating;               THE LOGLIKELIHOOD DECREASED IN THE LAST EM ITERATION.  CHANGE YOUR     MODEL, STARTING VALUES AND/OR THE NUMBER OF INTEGRATION POINTS.          THE MODEL ESTIMATION DID NOT TERMINATE NORMALLY DUE TO AN ERROR IN     THE COMPUTATION.  CHANGE YOUR MODEL AND/OR STARTING VALUES.      I have no idea what I can do next, any advice would be gratefully received.


title: Theory behind non-linear mixed effect models

text: I'm fairly new to this site. I'm currently doing some research work that involves the use of non-linear mixed effect model on longitudinal data (R). I'm presently using Pinheiro and Bates' book, _Mixed-Effects Models in S and S-PLUS_ , with the nlme package. By non-linear, I mean non-linear functions of independent variables with normal responses. I'm still a beginner with statistics and I'm finding the theoretical treatment of the subject in the above book quite daunting to understand. I was wondering if there are other "easier to digest" resources (including books) out there to try and get to grips with the theory. I understand how to use it in R but I still find that there are some fundamental concepts and gaps in my knowledge.


title: Do inputs to a Neural Network need to be in [-1,1]?

text: I want to create an weather forecast using Neural Networks. All the examples I saw used only values [-1,1] as an input. Is it also possible to use bigger values (like air pressure, degree calcius of the last couple of days,...) as inputs and get an number as output? Thank you


title: About the derivation of group Lasso

text: I've been reading the paper of group lasso, "Model selection and estimation in regression with grouped variables". http://www.stat.washington.edu/courses/stat527/s13/readings/yuanlin07.pdf In page 53 of the above paper, I don't know how to obtain `Eqn 2.3` and derive `Eqn. 2.4` from `Eqn. 2.2`: ![enter image description here](http://i.stack.imgur.com/id5CW.jpg) The following is my derivation, from `Eqn. 2.2`: $$ -X_j^TY+X_j^TX_j\beta_j+\sum_{i\neq j}X_j^TX_i\beta_i + \frac{\lambda \beta_j \sqrt{p_j}}{\|\beta_j\|}=0 $$ then, we have $$ \left( X_j^TX_j+\frac{\lambda \sqrt{p_j}}{\|\beta_j\|} \right)\beta_j=X_j^T(Y-\sum_{i\neq j}X_i^T\beta_i)=S_j $$ So, $$ \beta_j={\left( X_j^TX_j+\frac{\lambda\sqrt{p_j}}{\|\beta_j\|} \right)}^{-1}S_j $$ The result above has a little difference from `Eqn. 2.4` in the paper. Could you please indicate what's wrong with my derivation? In addition, in another document http://statweb.stanford.edu/~tibs/ftp/sparse-grlasso.pdf, the result of following is also different. ![enter image description here](http://i.stack.imgur.com/T604B.jpg) **UPDATED:** Yes! The first problem is solved. $$ S_j=X_j^T(Y-X\beta_{-j})=X_j^Tr_j $$ then, $$ \|\beta_j\|=\|S_j\|{\left( 1+\frac{\lambda\sqrt{p_j}}{\|\beta_j\|} \right)}^{-1} $$ Then, we could get $$ \|\beta_j\|=\|S_j\|-\lambda\sqrt{p_j} $$ then re-substitute into the formulation above above, then we get $$ \beta_j={\left( 1-\frac{\lambda\sqrt{p_j}}{\|S_j\|} \right)}S_j $$ But, where is the $+$ in the above formula?


title: Correlate continuous data with dichotomous data, with correction for an ordinal variable?

text: I want to study association b/w continuous data (age) and dichotomous data (complication: yes/no), with "correction" for an ordinal variable (intra- operative nerve laesion). I have a database of 427 patients, of whom 50 have suffered a post operative complication (altered sensation). The goal of this research is _to analyse the influence_ of "age at operation" on the chance of the postoperative complication "altered sensation". "Intra-operative nerve laesion" was scored (0/1/2/3/4/5) and is an imported factor influencing "altered sensation". So this must be included in the analysis. Hopefully somebody could help me with this analysis in SPSS.


title: Question about the relationship between the variance of residuals and variances of the data points

text: I was reading about regression. I understand one of the how the variance of each $x_i$ needs to be normal distributed, but why doesn't that related to the $s_b$? It seems like it should. I was reading the AP College Board paper about Inference. > Two more assumptions are that the errors are normally distributed and that > they have the same standard deviation for all x values. Although it is not > at all obvious, this is not the same as saying that the residuals are > normally distributed and all have the same standard deviation. The > difference is that the “errors” are the deviations between the response > variables and the “actual” underlying (invisible) linear phenomenon, whereas > the “residuals” are the deviations between the response variables and the > least-squares regression line, which is only an estimate of the “actual” > line. I understand (I think) that just because the error of the regression is normal doesn't mean that the variance of $x_i$ is normal because the sample could just randomly follow normality. (don't know if I said that statistically correct), but if all the $x_i$ are normal and the variance are equal, then why doesn't that variance equal the variance that one would get for the $\beta$ I have gotten some great answers here that are very clear and I am hoping that I could get another answer to help clear up my confusion. The kernel of my question is **What exactly is the relationship between the variance of the residuals and the $x_i$ and is there any connection to the _error, variance_ of $\beta$** ![the graph](http://i.stack.imgur.com/i1Qzy.jpg)


title: If the most common leading digit is one, is Benford's Law satisfied?

text: Suppose that we have a list of relevant numbers (e.g. maximum revenues for various years). If the leading digits are $1$ for all the years, can we say that trivially the maximum revenues satisfy Benford's Law?


title: Latent class model with both continuous and categorical indicators in R

text: I have question regarding which R-package to use to create a latent class/mixture model with both categorical and continuous indicator variables. I have not yet found a good example of this using R, even though there are a lot of mixture and latent class analysis packages in R. The software package Mplus have this capability and on the UCLA website there is an example of how to create such a model with Mplus (http://www.ats.ucla.edu/stat/mplus/seminars/introMplus_part2/lca.htm/ the second example). The only way I have been able to replicate this example in R is using the depmixS4 package. Using the following code:               library(depmixS4)     lca<-read.table("http://www.ats.ucla.edu/stat/mplus/seminars/introMplus_part2/lca.dat", sep=",")     names(lca)<-c( "hm", "he", "voc", "nocol" ,"ach9", "ach10", "ach11","ach12")          mod <- mix(list(hm~1,he~1,voc~1,nocol~1,ach9~1,ach10~1 ,ach11~1 ,ach12~1  ), data=lca, nstates=2,                   family=list(multinomial(),multinomial(),multinomial(),multinomial(), gaussian(),gaussian(),gaussian(),gaussian()),                   respstart=runif(32))          fmod<-fit(mod)     summary(fmod)      This code seems to do the job of replicating the example with some minor differences between the output in the depmixS4 variant and the Mplus. **First** , is this code correct for reproducing the Mplus example? **Second** , does anyone have any other packages they would recommend when creating this type of model in R?


title: Showing deviation of data on a fixed scale

text: I'm creating a chart in D3 which shows the Mode from a large (1mil +) dataset; and colors it according to the standard deviation of the dataset. But the issue is that I need to map the range (Red color to Green color) to a fixed domain. Is there a statistical measure which helps to represent information about the deviation of the data on a fixed scale (0-1; -1to1; 0-100 etc)?


title: Is there a KS test for panel data?

text: Wikipedias says > In statistics, the Kolmogorov–Smirnov test (K–S test) is a nonparametric > test for the equality of continuous, one-dimensional probability > distributions that can be used to compare a sample with a reference > probability distribution (one-sample K–S test), or to compare two samples > (two-sample K–S test). The Kolmogorov–Smirnov statistic quantifies a > distance between the empirical distribution function of the sample and the > cumulative distribution function of the reference distribution, or between > the empirical distribution functions of two samples. The null distribution > of this statistic is calculated under the null hypothesis that the samples > are drawn from the same distribution (in the two-sample case) or that the > sample is drawn from the reference distribution (in the one-sample case). In > each case, the distributions considered under the null hypothesis are > continuous distributions but are otherwise unrestricted. Is there a version of KS which could be applied to panel data i.e. where it is known that observations are not independent - same individuals measured over time. PS. I'm looking at comparing the distributions of two groups of subjects (survivals and deaths) although I have multiple measurements for individual subjects over time... Thanks EDIT: Potentialy useful page http://www.itl.nist.gov/div898/handbook/eda/section3/eda35g.htm


title: Is the "hybrid" between Fisher and Neyman-Pearson approaches to statistical testing really an "incoherent mishmash"?

text: There exists a certain school of thought according to which the most widespread approach to statistical testing is a "hybrid" between two approaches: that of Fisher and that of Neyman-Pearson; these two approaches, the claim goes, are "incompatible" and hence the resulting "hybrid" is an "incoherent mishmash". I will provide a bibliography and some quotes below, but for now suffice it to say that there is _a lot_ written about that in the wikipedia article on Statistical hypothesis testing. Here on CV, this point was repeatedly made by @Michael Lew. I have two questions about this line of thought:   1. Why are F and N-P approaches claimed to be incompatible and why is the hybrid claimed to be incoherent? Note that I read at least six anti-hybrid papers (see below), but still fail to understand the problem or the argument. Note also, that I am not suggesting to debate if F or N-P is a better approach; neither am I offering to discuss frequentist vs. Bayesian frameworks. Instead, the question is: accepting that both F and N-P are valid and meaningful approaches, what is so bad about their hybrid?   2. How accepted is this anti-hybrid point of view among professional statisticians? Is it rather a freak position, or maybe a trivial commonplace? Here is how I understand the situation. Fisher's approach is to compute the $p$-value and take it as an evidence against the null hypothesis. The smaller the $p$, the more convincing the evidence. The researcher is supposed to combine this evidence with his background knowledge, decide if it is convincing _enough_ , and proceed accordingly. (Note that Fisher's views changed over the years, but this is what he seems to have eventually converged to.) In contrast, Neyman-Pearson approach is to choose $\alpha$ ahead of time and then to check if $p\le\alpha$; if so, call it significant and reject the null hypothesis (here I omit large part of the N-P story that has no relevance for the current discussion). See also an excellent reply by @gung in When to use Fisher and Neyman-Pearson framework? The hybrid approach is to compute the $p$-value, report it (implicitly assuming that the smaller the better), and also call the results significant if $p\le\alpha$ (usually $\alpha=0.05$) and nonsignificant otherwise. This is supposed to be incoherent. How can it be invalid to do two valid things simultaneously, beats me. As particularly incoherent the anti-hybridists view the widespread practice of reporting $p$-values as $p<0.05$, $p<0.01$, or $p<0.001$ (or even $p\ll0.0001$), where always the strongest inequality is chosen. The argument seems to be that (a) the strength of evidence cannot be properly assessed as exact $p$ is not reported, and (b) people tend to interpret the right-hand number in the inequality as $\alpha$ and view it as type I error rate, and that is wrong. I fail to see a big problem here. First, reporting exact $p$ is certainly a better practice, but nobody really cares if $p$ is e.g. $0.02$ or $0.03$, so rounding it on a log scale is not soooo bad (and going below $\sim 0.0001$ does not make sense anyway, see How should tiny p-values be reported?). Second, if the consensus is to call everything below $0.05$ significant, then error rate will be $\alpha=0.05$ and $p \ne \alpha$, as @gung explains in Interpretation of p-value in hypothesis testing. Even though this is potentially a confusing issue, it does not strike me as more confusing than other issues in statistical testing (outside of the hybrid). Also, every reader can have her own favourite $\alpha$ in mind when reading a hybrid paper, and her own error rate as a consequence. **So what is the big deal?** One of the reasons I want to ask this question is because it literally _hurts_ to see how much of the wikipedia article on Statistical hypothesis testing is devoted to lambasting hybrid. Following Halpin & Stam, it claims that a a certain Lindquist is to blame (there is even a large scan of his textbook with "errors" highlighted in yellow), and of course the wiki article about Lindquist himself starts with the same accusation. But then, maybe I am missing something. * * * ## References   * Gigerenzer, 1993, The superego, the ego, and the id in statistical reasoning \-- _introduced the term "hybrid" and called it "incoherent mishmash"_   * Cohen, 1994, The Earth Is Round ($p<.05$) \-- _a very popular paper with almost 3k citations, mostly about different issues but favourably citing Gigerenzer_   * Goodman, 1999, Toward evidence-based medical statistics. 1: The P value fallacy   * Hubbard & Bayarri, 2003, Confusion over measures of evidence ($p$'s) versus errors ($\alpha$'s) in classical statistical testing \-- _one of the more eloquent papers arguing against "hybrid"_   * Halpin & Stam, 2006, Inductive Inference or Inductive Behavior: Fisher and Neyman-Pearson Approaches to Statistical Testing in Psychological Research (1940-1960) _[free after registration] -- blames Lindquist's 1940 textbook for introducing the "hybrid" approach_   * @Michael Lew, 2006, Bad statistical practice in pharmacology (and other basic biomedical disciplines): you probably don't know P \-- _a nice review and overview_ ## Quotes > **Gigerenzer:** What has become institutionalized as inferential statistics > in psychology is not Fisherian statistics. It is an incoherent mishmash of > some of Fisher's ideas on one hand, and some of the ideas of Neyman and E. > S. Pearson on the other. I refer to this blend as the "hybrid logic" of > statistical inference. > > **Goodman:** The [Neyman-Pearson] hypothesis test approach offered > scientists a Faustian bargain -- a seemingly automatic way to limit the > number of mistaken conclusions in the long run, but only by abandoning the > ability to measure evidence [a la Fisher] and assess truth from a single > experiment. > > **Hubbard & Bayarri:** Classical statistical testing is an anonymous hybrid > of the competing and frequently contradictory approaches [...]. In > particular, there is a widespread failure to appreciate the incompatibility > of Fisher's evidential $p$ value with the Type I error rate, $\alpha$, of > Neyman-Pearson statistical orthodoxy. [...] As a prime example of the > bewilderment arising from [this] mixing [...], consider the widely > unappreciated fact that the former's $p$ value is _incompatible_ with the > Neyman-Pearson hypothesis test in which it has become embedded. [...] For > example, Gibbons and Pratt [...] erroneously stated: "Reporting a P-value, > whether exact or within an interval, in effect permits each individual to > choose his own level of significance as the maximum tolerable probability of > a Type I error." > > **Halpin & Stam:** Lindquist's 1940 text was an original source of the > hybridization of the Fisher and Neyman-Pearson approaches. [...] rather than > adhering to any particular interpretation of statistical testing, > psychologists have remained ambivalent about, and indeed largely unaware of, > the conceptual difficulties implicated by the Fisher and Neyman-Pearson > controversy. > > **Lew:** What we have is a hybrid approach that neither controls error rates > nor allows assessment of the strength of evidence.


title: Minitab and output

text: When you conduct a 1 sample t-test why does minitab only give an upper bound or a lower bound for a confidence interval? So if $H_0: \mu = \mu_0$ and we test $H_a: \mu < \mu_0$ or $H_a: \mu > \mu_0$ why is only an upper bound or lower bound of a confidence interval returned respectively? Whereas testing $H_a: \mu \neq \mu_0$ returns both a lower and upper bound.


title: Chi square test on non-normal distributions

text: I'm a beginner in statistics and would like to understand one of the requirements for the chi-square test to work. It is that the random variable being measured has to follow a normal distribution. My question then is, what other models should I use if the random variables I'm measuring do not follow a normal distribution? This is similar to how the Student-t test should be used instead of the Standard Normal for a small sample, is there an equivalent for chi-square test for a small sample?


title: R: Most efficient way to bootstrap fan charts

text: I want to bootstrap a timeseries prediction to generate fan charts of marginal prediction intervals. What is the most efficient way to do this? Something like this?               library(boot)     data <- load(mydata)     boot_data <- tsboot(data, NULL, ...)     prediction_paths <- data(boot.array(boot_data))     ...      I am thinking that this is likely to be slow: first it generates random numbers to created the unneeded point estimate, then it re-creates the same random numbers from the stored seed to create the matrix of prediction paths. Would it be better to write my own bootstrap algorithm that skips the point estimate and creates directly the needed matric of prediction paths? E.g. along the lines of:               indices <- sample(1:length(data), R*N)     dim(indices) <- c(R,N)     prediction_paths <- data(indices)      Or is tsboot still faster due to a better implementation that I can hope to achieve?


title: Bernoulli and Normal random variable

text: Let $Y$ be a mixture random variable that is distributed $\text{Normal}(\mu_1, \sigma_1^2)$ with probability $p$ and $\text{Normal}(\mu_2, \sigma_2^2)$ with probability $1-p$. Can I estimate $p$ as well as $\mu_1, \mu_2, \sigma_1^2, \sigma_2^2$ from a random sample? How?


title: Regressing a differenced variable on a lagged variable. How can I fix the error in R?

text: I have a time series (std) of 324 observations with no missing values, starting from January 1987 and ending in December 2013. I want to regress via OLS the one in the question. In R, the code:                lm(dstd ~ lstd)      Where 'dstd' is the differenced variable, via:                dstd<-diff(std,differences=1)      And lstd is the lagged variable via:                lstd<-lag(std, k=1)      Does not work, as it gives the following error:                Error in model.frame.default(formula = dstd ~ lstd, drop.unused.levels = TRUE) :       variable lengths differ (found for 'lstd')      Which is, of course, natural as the length of 'dstd' is 323 while the length of 'lstd' is 324. Searching for the error, na.omit has been suggested as the answer for other models. However, no combinations of 'na.omit' or 'na.exclude' would work. This video (at the relevant time) says ' you need to change these', referring to these variables. I have tried excluding while forming another object, as well as subsetting in the regression model, yet none of these seem to work. What code would I need to fix this? Here is a Minimal Reproducible Example: An object 'stdshort' is created with the same characteristics of the parent object, via this code:                > stdshort<-dput(window(std, 1987, 1989, 12))      structure(c(4.5, 4.7, 4.2, 4.4, 3.9, 3.9, 3.7, 3.7, 3.4, 3.6,       3.5, 3.1, 3.5, 3.3, 3.7, 3.7, 3.7, 3.8, 3.6, 3.5, 3.5, 3.3, 3.5,       3.5, 3.4), .Tsp = c(1987, 1989, 12), class = "ts")      The following code replicates the steps previously mentioned: differencing, lagging and linear regression. First, a differenced object is created. The object dtstdshort's structure looks like this:                > dstdshort<-diff(stdshort,differences=1)      > dput(dstdshort)       structure(c(0.2, -0.5, 0.2, -0.5, 0, -0.2, 0, -0.3, 0.2, -0.1,       -0.4, 0.4, -0.2, 0.4, 0, 0, 0.0999999999999996, -0.2, -0.1, 0,       -0.2, 0.2, 0, -0.1), .Tsp = c(1987.08333333333, 1989, 12), class = "ts")      A following object, lstdshort, for 1 period lagging is constructed. This is the structure:                > lstdshort<-lag(stdshort,k=1)      > dput(lstdshort)      structure(c(4.5, 4.7, 4.2, 4.4, 3.9, 3.9, 3.7, 3.7, 3.4, 3.6,       3.5, 3.1, 3.5, 3.3, 3.7, 3.7, 3.7, 3.8, 3.6, 3.5, 3.5, 3.3, 3.5,       3.5, 3.4), .Tsp = c(1986.91666666667, 1988.91666666667, 12), class = "ts")      The regression of those two objects as explained before, gives the following error (with code):                > lm(dstdshort ~ lstdshort)      Error in model.frame.default(formula = dstdshort ~ lstdshort, drop.unused.levels = TRUE) :         variable lengths differ (found for 'lstdshort')


title: Is this a case of semi-supervised classification?

text: I have a rule-based classifier and I know for sure I can classify my data in a number of N classes. And I have also a "small" dictionary where my rules check (during classification process) for every observation iterated. Is this a semi-supervised classification? If the dictionary was built by-hand from a curated process, is this dictionary considered a "training set"?


title: Compare linear regression models (same and different response variable)

text: How can I (1) compare two linear models between years and (2) Can I compare 2 models with different response variables? My data have 4 variables: y_meas, x, year, y_calc. "y_meas" is a lab measured response varibale, and "y_calc" is an estimate of the same variable, using a standard calculation. "x " is a dosage, similar(ish) between two years:               #create dataset     set.seed(100)     dat <- within(data.frame(x = rep(1:10, times=2)),                      {                        year <- rep(1990:1991, each = 10)                        y_meas <- 0.5 * x* (1:20) + rnorm(20)              y_calc <- 0.3 * x* (1:20) + rnorm(20)                        year <- factor(year)                 # convert to a factor                      }                      )      I have two related questions: (1) Is there any difference between slope/intercept of models for 1990 and 1991?               m.1990<-lm(y_meas~x, data=subset(dat, year==1990))     m.1991<-lm(y_meas~x, data=subset(dat, year==1991))     anova(m.1990)     anova(m.1991)      # both models are significant      I can't run `anova(m.1990,m.1991)` because the models are not nested? Do I need to use year as a dummy variable and run ANCOVA? What does this look like (roughly)? (2) Assuming I can combine 1990 and 1991, can I compare the slope/intercept of 'y_meas~x' and 'y_calc~x'? Yes, two different response variables, but are on the same scale.


title: How can I calculate genetic similarities between groups?

text: I know how to calculate genetic similarity between genotypes, but I do not know how to calculate genetic similarity between groups. From cluster analysis, I found seven groups. Now, I need to calculate genetic similarity between these groups. How should I go about doing this?


title: Simple case of MNAR missing data

text: I believe I have a very simple problem with missing data, but I'm a bit lost because all the materials I read seem to be focused on much more complicated cases. I have a random variable $X$ which has a Binomial distribution with parameters $n$ and $p$ (i.e., $X\sim B(n,p)$), where $p$ is unknown ($n$ is known). I have $K$ independent samples of $X$, let them be $X_1,\dotsc,X_K$ but some of them are missing, in the sense that I observe the values $R_iX_i$, where $R_i\in\\{0,1\\}$. I know that $R_i = 0$ when $X_i< C$, for some known constant $C$, otherwise $R_i=1$. This means that my missing data are MNAR (Missing Not At Random), but I know the "mechanism" leading to their being missing, which I believe is what makes my case easier. Edit: it's a case of left-censoring with Binomial data. How can I estimate $p$? I'm also interested in finding a confidence interval for $p$. References and links are much appreciated. Thanks in advance.


title: Expected value of the log-determinant of a Wishart matrix

text: Let $\Lambda \sim \mathcal W_D(\nu, \Psi)$, i.e. distributed according to a $D \times D$ dimensional Wishart distribution with mean $\nu \Psi$ and degrees of freedom $\nu$. I would like an expression for $E(\log |\Lambda|)$ where $|\Lambda|$ is the determinant. I've google'd a bit for the answer to this and have gotten some conflicting information. This paper explicitly states that $$ E(\log|\Lambda|) = D \log 2 + \log |\Psi| + \sum_{i = 1} ^ D \psi\left(\frac{\nu - i + 1} 2\right) $$ where $\psi(\cdot)$ denotes the digamma function $\frac d {dx} \log \Gamma(x)$; the paper does not give a source for this fact as far as I can tell. This is also the formula used on the wikipedia page for the Wishart, which sites Bishop's Pattern Recognition text. On the other hand, google turned up this discussion with a linked paper that states that $$ \nu^D \frac{|\Lambda|}{|\Psi|} \sim \chi^2_\nu \chi^2_{\nu - 1} \cdots\chi^2_{\nu - D + 1}. \qquad (\dagger) $$ They conclude by stating that $$ E(\log | \Lambda|) = D \log 2 - D \log \nu + \log |\Psi| + \sum_{i = 1} ^ D \psi\left(\frac{\nu - i + 1} 2\right) $$ which is derived using the fact that $E(\log \chi^2_\nu) = \log(2) + \psi(\nu /2)$. I checked this calculation starting from $(\dagger)$ and it seems okay, but we have an extra $-D\log\nu$.


title: change in category variable

text: dependent variables (ordinal): credit rating 1970 (cr70) and credit rating 1980 (cr80). Here is what i want to do: regress cr80-cr70 = independent vars. how could this be done and how could you interpret it!? if the dependent variable is continuous it would be simple. but can you make a new var from the difference of two ordinal vars, and have that be the dependent var?


title: Setting a proper prior distribution on sigma for Bayesian Regression

text: I'd like to know how to set a proper prior distribution on sigma when using JAGS in a Bayesian regression. Here is what I have:               model{         for( i in 1 : N ) {             y[i] ~ dnorm(y_hat[i], tau )             y_hat[i] <- inprod(b[1:K], x[i,1:K])         }                 sigma ~ dunif(0,10000)         #tau ~ dgamma(.001,.001)         tau <- pow(sigma,-2)              for(j in 1:K) {           b[j] ~ dnorm(1/K,4*K*K)             }     }      K is the number of predictor variables I have. If I have a y[i] with known range of 0 to 1 (or any known range), then how could I use that knowledge to change the distribution I've chosen on sigma? It seems a bit crazy to choose a prior distribution like I've done when there's no chance in the world that it would be anywhere close to 10,000.


title: Can I use wilcoxon signed rank test for time series data

text: I have a question about wilcoxon signed rank test. I have 2 different speed profiles (speed starts from time 0 to 100 sec) of the same participant that was given 2 different treatment. I would like to ask if I could use wilcoxon signed rank test to compare if both speed profiles are different from each other. I am a little bit confused how to fit the null hypothesis of this test (median difference) to my case. I also think that once I use wilcoxon signed rank test to my time-series data, I may lose the sence of data sequence from my speed profile. Any help from anyone will be appreciated. Thank you, Ayla


title: Interpret Ordinal Data on Likert Scale

text: I am curious to know whether it is current practice in the statistical field to interpret and report data obtained on a Likert Scale (1: Strongly Agree to 5: Strongly Disagree), by finding the median value rather than averaging all the values for a particular item. Please see link here My second question is how do we decide whether we should use 1: for Strongly Disagree rather than Strongly Agree such that: 1: Strongly Disagree 2:Disagree 3:undecided 4:Agree 5:Strongly Agree 1: Strongly Agree 2:Agree 3:undecided 4:Disagree 5:Strongly Disagree


title: Deriving the limit value of an estimtor

text: I kinda need some help with this question. This came from a lecture about Maximum likelihood estimators. ![a busy cat](http://i.imgur.com/PSEGffU.png) I tried to see if I could derive the limit Function for Qn which I think would be 1 if theta = 0 1/2 - |theta-2| for 3/2 <= theta <= 5/2 0 otherwise (i'm sorry for the lack of mathjax. I can't get it to work for some reason) But I'm not sure about how to find the maximum estimate of both theta's in part c.


title: Which test do I use; need to compare tests with one another pre and post intervention

text: I am writing a thesis on the validity of balance training tests. I need to compare the tests and see if they correlate with one another in testing balance in the elderly. There's a control and intervention group, and 9 different balance tests were used to test each subject. Need to see if the tests correlate with one another or not.


title: Uniform convergence of Monte Carlo approximation

text: Usually Monte Carlo method is used to compute integration. For example, let $g(x,\theta)$ be a continuous function about $x$ and $\theta$, $f(x \mid \theta)$ is a continuous pdf with parameter $\theta$. if the integral $E[g(x,\theta)]=\int g(x,\theta)f(x\mid \theta)dx$ has no closed form . Then we can use Monte Carlo approximation $E[g(x,\theta)] \approx \frac{1}{m}\sum_{l=1}^m g(X_i,\theta)$, where $X_1,\dots,X_m$ is a random sample from distribution $f(x \mid \theta)$. By law of large number, $\frac{1}{m}\sum_{l=1}^m g(X_i,\theta)$ converges to $E[g(x,\theta)]$ in probability for each $\theta \in \Theta$. Do we have theorem talking about uniform convergence of Monte Carlo approximation? I mean, under what condition, we have $\sup\limits_{\theta\in\Theta} \left\| \frac1n\sum_{i=1}^n g(X_i,\theta) - \operatorname{E}[g(X,\theta)] \right\| \xrightarrow \ 0$ a.s. or in probability Assume parameter space $\Theta$ is compact


title: Random walks in Markov-Switching states

text: Does anyone know of any literature which discusses the idea of random walk processes in states of a Markov-Switching model? I am struggling to find anything which discusses the theory/fitting behind that (most compare the M-S and RW models with no mention of more specific stochastic process in the M-S model)


title: R, lavaan package, latent change model: Definition of change parameter?

text: I'm currently struggling with the latent change model that I try to realize in R with the lavaan package. R can't estimate all standard errors, which makes me believe that the model is not fully identified. I did not have identification problems with the latent state model, so I guess the critical part may be due to the change parameter.                     Latent_change_prepost <- '                      # first order latent variable definition                           T1_43_1 =~ 1*A43_1_1                                     + 1*A43_1_2                                     + 1*A43_1_3                                     + 1*A43_1_4                                     + 1*A43_1_5                                 T1_43_2 =~ 1*A43_2_1                                     + 1*A43_2_2                                     + 1*A43_2_3                                     + 1*A43_2_4                                     + 1*A43_2_5                                 T1_43_7 =~ 1*A43_7_1                                     + 1*A43_7_2                                     + 1*A43_7_3                                     + 1*A43_7_4                                     + 1*A43_7_5                                 T1_43_3 =~ 1*A43_3_1                                     + 1*A43_3_2                                     + 1*A43_3_3                                     + 1*A43_3_4                                 T1_43_4 =~ 1*A43_4_1                                     + 1*A43_4_2                                     + 1*A43_4_3                                     + 1*A43_4_4                                 T1_43_5 =~ 1*A43_5_1                                     + 1*A43_5_2                                     + 1*A43_5_3                                  # second order latent variable definition                           T1 =~ 1*T1_43_1                               + 1*T1_43_2                               + 1*T1_43_7                                T2 =~ 1*T1_43_3                               + 1*T1_43_4                               + 1*T1_43_5                            # regressions                           change =~ 1*T2 #(no regression, need to introduce "change")                           T2 ~ 1*T1                            # (co)variances                           change ~~ T1                           T2 ~~ 0*T2                            # intercepts                           T1 ~ NA*1                           T2 ~ NA*1                           change ~ NA*1                           T1_43_1 ~ NA*1                           T1_43_2 ~ NA*1                           T1_43_7 ~ NA*1                           T1_43_3 ~ NA*1                           T1_43_5 ~ NA*1                         '            fit2_latent_change <- sem(                                 model = Latent_change_prepost                               , data = diss_prepost_w_interv_kovar                               , missing = "fiml"                               , meanstructure = TRUE                               )      I get estimates of most parameters, but no std. errors. Also if I alternate the model slightly I might get std. errors but they are astronomically high. Any ideas and help would be much appreciated!


title: Covary two dummy variables in SEM?

text: I am running a structural equation model (SEM) in Amos 18, and I want to test the impact of marital status on several latent variables. Marital status is nominal, so I created three dummy variables:   1. Mar_Single: 1 = yes, 0 = no   2. Mar_Married: 1 = yes, 0 = no   3. Mar_Other: 1 = yes, 0 = no I included Mar_Single and Mar_Married in the SEM, so their coefficients will be interpreted against the omitted (reference) group, Mar_Other. The modification indices suggested fit could be improved significantly if I covary Mar_Single and Mar_Married. Should I do this? In a way, this makes sense because they are perfectly correlated: If Mar_Single = 1 then Mar_Married will always = 0. Larry


title: Pull out most important variables from PCA

text: I would like to get the most important variables from a PCA result. I see two clusters in the plot. I now that is possible that there is no only one variable causing this, so maybe I would have to get more than one variable. I'm using "Adegenet" R package. My original data is a matrix where rows are PubMed papers and columns are MeSH keywords. Data has been transformed into SNP-like to adapt the method to the new input data. Please point me the the correct R package if you think I'm doing the incorrect work with this package, I just chose it because I already knew how it works. ![First two dimensions of PCA based on 1359 metagenomcs papers and 3459 MeSH terms](http://i.stack.imgur.com/1LhUv.png)               #R code     library(adegenet)     #Load SNPs     myPath <- "pubmed_result_metagenomics_ALL_parsed.fasta"#core SNPs retrieved with kSNP from 188 H. parasuis strains, removing from the analysis the strains tagged with ‘NK’ phenotype. kSNP k-mer sizes tested were 25, 20 and 15, selecting the run that gave more SNPs, i.e., 15.     core_SNPs_matrix <- fasta2genlight(myPath, chunk=1000, multicore=FALSE)#     core_SNPs_matrix <- as.matrix(core_SNPs_matrix)          # Principal Component Analysis (PCA)     pca1 <- glPca(core_SNPs_matrix) # 10 components saved     pca1          # Draw PCA colorplot     myCol <- colorplot(pca1$scores,pca1$scores, transp=TRUE, cex=4)     abline(h=0,v=0, col="grey")     add.scatter.eig(pca1$eig[1:40],2,1,2, posi="topright", inset=.05, ratio=.3)     title("First two dimensions of PCA \n based on 1359 metagenomcs papers \n and 3459 MeSH terms")     dev.copy2pdf(file = "Figure_12.pdf") #Save as .pdf#


title: T test with one group

text: I think it is a simple question, but I got many different results and now I am not sure which test to use. I have one group and 2 different scores (higher scores means higher either positive or negative associations). I want to know whether I have differences between these variables. I need to consider that these scores have a similar behavior, so I expect differences (one score could be higher than the other). I tried with t-test, but they are used when I have more than 1 group, right? Thanks


title: How to use OLS or NLS on very small sample size?

text: right now I am facing the problem of handling a very small sample size which consists of 8 observations with 24 variables. I would like to do first just simple OLS with one dependent and 23 independent variables and use these coefficients to estimate NLS estimates. These estimates will be used for further studies. The problem obviously is that I got more parameters than observations. I tried to bootstrap my observations and just do simple OLS with a handful of parameter with bootstrap repetitions between 50 and 50,000 to improve my results, however, in fact, the significance of some coefficients even decreased. Is it even possible to derive somewhat significant coefficients here? What else can I do? I am new here, so I dont know yet how to include a table with exemplary data. Sorry. Thanks for any help, TheJoez


title: How to analyze this data (difference among means)?

text: This topic came up today and there are multiple methods that have come forth; I'm looking for the appropriate way to analyze this data. We have five brands of material (cloth), each brand with four different samples taken from a piece of that cloth. Each sample then has between 10-12 observations of durability taken. What is the best way to categorize and analyze this data? One suggestion that came up was grouping all samples together (i.e., combining the four samples' 10 observations to have 40 total observations for brand A). Another was using a general linear model to perform an unbalanced nested ANOVA, with brand as a fixed factor and sample as a random factor. Another idea was calculating the mean and standard error for each of the four different samples using the 10-12 observations, then using a weight to calculate the overall best estimate of the mean for each of the four samples and using that (in other words, each brand would then have four observations, comprised of the weighted average of its observations). Finally, the idea was put forth to sample one of the samples randomly from the four samples and use a one way ANOVA and post-hoc test to examine difference of means. What is the appropriate method (even if none of the above) to deal with data in this format? Thanks for your time. EDIT: Forgot to mention, ideally what we're looking for is differences among brands (if any exist).


title: How do I calculate confidence intervals for a non-normal distribution?

text: I have 383 samples that have heavy bias for some common values, how would I calculate the 95% CI for the mean? I the CI that I calculated seem way off, which I assume is because my data does not look like a curve when I make a histogram. So I think I have to use something like bootstrapping, which I don't understand very well.


title: Matrix completion: How to assign names to the completed columns?

text: I am wondering if this the right place to ask this question. Normally it should be :). I am recently reading some papers about matrix completion such as in here, and here. I didn't go through some testing yet, but I am wondering about the following problem, suppose I have some elements $E_1,E_2,....,E_n$ and suppose that I have information about $m$ elements $E_1,...,E_m$ such as $m < n$. Suppose that at this stage I applied one of the matrix completion algorithms and I got my $n \times n$ matrix, How to assign labels (names) to the completed columns? The problem is that the column order can be arbitrary, so it may be mathematically correct but may have no real meaning if the labels are not assigned correctly.


title: How to detect which one is the better study when they give you conflicting results?

text: You so often come across in the press various studies that conclude directionally opposite results. Those can be related to the testing of a new prescription drug or the merit of a specific nutrient or anything else for that matter. When two such studies arrive at conflicting results how can you tell which one of the two is closest to the truth?


title: Using propensity scores from twang in coxph

text: I've used the `twang` package to calculate propensity score weights in my observational data set. The `twang` vignette uses the `survey` package to apply the weighted PS in analyses, but does not provide an example for Cox proportional hazard estimation. The `survey` package provides a function called `svycoxph()` which which seems to do the trick. However, when i use `coxph()` from the `survival` package an set the `weight=` argument to the extracted weighted PS from `twang`'s `ps()` object i get an indentical HR estimate, albeit a much more narrow 95%CI. **Can I use the extracted propensity score weights from a ps object from the`twang` package as a vector in the `weight=` argument of the `survival` package's `coxph()`?** Should i include an example? **Update:** When i added a `cluster(patient ID)` term in the the `coxph` formula, both the HR and the 95% CI became identical to the `syvcoxph` output. So I guess that did the trick. But can i conclude that this is ok to use?


title: What is wrong with this "naive" shuffling algorithm?

text: This is a follow-up to a Stackoverflow question about shuffling an array randomly. There are established algorithms (such as the Knuth-Fisher-Yates Shuffle) that one should use to shuffle an array, rather than relying on "naive" ad-hoc implementations. I am now interested in proving (or disproving) that my naive algorithm is broken (as in: does not generate all possible permutations with equal probability). Here is the algorithm: > Loop a couple of times (length of array should do), and in every iteration, > get two random array indexes and swap the two elements there. Obviously, this needs more random numbers than KFY (twice as much), but aside from that does it work properly? And what would be the appropriate number of iterations (is "length of array" enough)?


title: Explain model adjustment, in plain English

text: Reading about methods and results of statistical analysis, especially in epidemiology, I very often hear about **_adjustment_** or **_controlling_** of the models. How would you explain, to a non-statistician, the purpose of that? How do you interpret your results after controlling for certain variable? Small walk-through in Stata or R, or a pointer to one online, would a true gem.


title: Prediction interval on untransformed scale

text: I am attempting to use simple linear regression to construct a 95% prediction interval for a continuous response variable (Y) using a continuous input variable (X). When examining my data, I realized that if I log-transform both X and Y, the assumptions of equal variance and normality are essentially met, but if I do not log transform both variables, they are not. I'm able to run all of the R code needed to get the prediction interval for the transformed data, and I've read online and elsewhere how to interpret this interval. However, is there a way to get a 95% prediction interval on the untransformed scale using the 95% prediction interval on the log-transformed scale? In other words, I'm interested in talking about the conditional mean response instead of the median response (which, according to what I've read, is how I must interpret my prediction interval). I could exponentiate the interval endpoints, but I'm fairly sure that's incorrect because E(e^{Y}) is usually not equal to e^{E(Y)}.


title: Regression discontinuity design parametric versus non-parametric different result

text: I am using the parametric approach and non-parametric (local linear regression) approaches of regression discontinuity design (RDD) to compute the treatment effect using Stata. To get the user-written `rd` and the 102nd Congress data, I do this:               net get rd     use votex      The local linear approach:               rd lne d,bw(0.20) mbw(100) ker(rec)     Two variables specified; treatment is      assumed to jump from zero to one at Z=0.            Assignment variable Z is d      Treatment variable X_T unspecified      Outcome variable y is lne          Estimating for bandwidth .2     ------------------------------------------------------------------------------              lne |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------            lwald |  -.1046939   .1147029    -0.91   0.361    -.3295075    .1201197     -----------------------------------------------------------------------------      As far as I understand this is equivalent to following :               gen win_d=win*d     reg lne d win win_d if d>=-0.2 & d<=0.2                Source |       SS       df       MS              Number of obs =     267     -------------+------------------------------           F(  3,   263) =    0.43            Model |  .271662326     3  .090554109           Prob > F      =  0.7339         Residual |  55.7885045   263  .212123591           R-squared     =  0.0048     -------------+------------------------------           Adj R-squared = -0.0065            Total |  56.0601668   266  .210752507           Root MSE      =  .46057          ------------------------------------------------------------------------------              lne |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------                d |   .8450601   .7855123     1.08   0.283    -.7016333    2.391753              win |  -.1046939   .1257913    -0.83   0.406    -.3523801    .1429923            win_d |  -.8707605   1.048807    -0.83   0.407    -2.935887    1.194366            _cons |   21.44195   .0925378   231.71   0.000     21.25974    21.62415     ------------------------------------------------------------------------------      However, when we use the parametric approach (let's say with the polynomial of order one), we use all the observations. But, I am trying to see how parametric approach can be compared with non-parametric approach with the same number of observation as in non-parametric approach. So, I do as follows:               reg lne d win if d>=-0.2 & d<=0.2                Source |       SS       df       MS              Number of obs =     267     -------------+------------------------------           F(  2,   264) =    0.30            Model |  .125446108     2  .062723054           Prob > F      =  0.7440         Residual |  55.9347207   264  .211873942           R-squared     =  0.0022     -------------+------------------------------           Adj R-squared = -0.0053            Total |  56.0601668   266  .210752507           Root MSE      =   .4603          ------------------------------------------------------------------------------              lne |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------                d |   .3566172   .5201877     0.69   0.494    -.6676274    1.380862              win |  -.0964314   .1253232    -0.77   0.442    -.3431916    .1503288            _cons |   21.39136   .0696112   307.30   0.000      21.2543    21.52843     ------------------------------------------------------------------------------      My concern is why the non-parametric approach result (-.1046939) is not the same as parametric approach (-.0964314), although we are using the same observation for both.


title: Regression analysis to estimate affinity or antipathy among celebrity wedding parties

text: Please consider the following plot of a graph: ![Plot 1](http://i.stack.imgur.com/EIB8r.jpg) Each node represents a wedding party while each edge is the relation between the two parties expressed in terms of the number of people who attended both weddings (with an edge of weight 0 in case of no common participant). ![Plot 2](http://i.stack.imgur.com/IaZGo.jpg) Each node $V_i$ has two attributes: **size** $S_i$, the number of people participating, and **time** $T_i$, the year the wedding took place. Each edge $E_{ij}$ has one attribute: **weight** $W_{ij}$, the number of participants that $V_i$ and $V_j$ have in common. I am interested in comparing each relation $E_{ij}$ against the whole dataset to understand which relation is stronger/weaker than expected. In other words, I want to obtain from the data a benchmark value — or the dataset "average" when controlling for the the two nodes attributes **size** ($S_i$ and $S_j$) and **time** ($T_i$ and $T_j$) — for each edge and then comparing the "average" to the real value of $W_{ij}$. * * * > **Put simply, what I want to infer from my data is whether the community > participating to a weeding party $V_i$ shows an affinity or an antipathy > with the community participating to $V_j$.** * * * Because of my limited statistical tool-box, I though of simply proceeding with a regression analysis with **weight** (`shared_ppl` in the R script below) as response or dependent variable. As predictors I used `mult_size` (as $S_i \times S_j$) and `time_diff`(as $\left|T_i - T_j\right|$). As it is reasonable to assume, I want to test whether the **size** of the parties are correlated with **weight** and the difference in **time** are correlated with **weight** , and I want to estimate sign and quantity of both correlations. In my simple example (see below for complete script) the linear regression > $\mathbf{shared\\_ppl} = \beta_0 + \beta_1 \times \mathbf{mult\\_size} + > > \beta_2 \times \mathbf{time\\_diff}$ will return the following coefficients:               Coefficients:                      Estimate    Std. Error t value Pr(>|t|)       (Intercept) -6395.7169331  3558.2955398  -1.797   0.0891 .     mult_size       0.0016097     0.0008846   1.820   0.0855 .     time_diff      -7.1564129     3.0076486  -2.379   0.0286 *      (NOTE: I know that the p-value for `mult_size` is slightly above the standard .05, but still let's considered it for the sake of the example as significantly correlated to `shared_ppl`) Now I can "predict" for each $E_{ij}$ an expected $\hat{W}_{ij}$ and then compare it with the real $W_{ij}$. Specifically what I did was to calculate the difference (labelled as `offset` in the script) from the .95 confidence interval produced by the model. Then   1. a value of 0 would indicate that $W_{ij}$ of the relation $E_{ij}$ is within the .95 confidence interval for $\hat{W}_{ij}$ as predicted by the regression analysis above   2. a value of -14.37 would indicated that $W_{kl}$ of the relation $E_{kl}$ is short of 14.37 people if compared to the .95 confidence interval for $\hat{W}_{kl}$. The `offset` for each relation is plotted below ![Plot 3](http://i.stack.imgur.com/MOawY.jpg) * * * > Then, given my model above, I can infer from the data that > >   1. $V_i$ and $V_j$ do not show neither reciprocal affinity or antipathy, > while >   2. $V_k$ and $V_l$ do show reciprocal antipathy > **Are my analysis and its conclusions correct from a statistical point of view?** * * * _For those who want to replicate the example in R here is the code_               library(igraph)          # Weddings from http://www.huffingtonpost.com/2011/10/02/famous-married-couples-ou_n_980125.html          # Create a full graph (all vertices are connected)     g <- graph.full(7)          # Assign attribute to graph     V(g)$name <- c("Ozzie\nAnd\nSharon Osbourne",                      "Ali Hewson\nAnd\nBono",                      "Jamie Lee Curtis\nAnd\nChristopher Guest",                     "Kevin Bacon\nAnd\nKyra Sedgwick",                     "Rita Wilson\nAnd\nTom Hanks",                     "John Travolta\nAnd\nKelly Preston",                     "David Bowie\nAnd\nIman")         V(g)$year <- c(2000,2001,2004,2005,2008,2009,2011)     V(g)$size <- c(200,300,80,90,100,400,500)         E(g)$weight <- c(100,20,11,10,20,0,21,6,10,70,40,100,8,3,0,2,14,100,100,150)          # Create layout for plot     layout <- layout.circle(g)          # Plot 1     plot(g,          layout=layout,          vertex.size=V(g)$size*0.1)          # Plot 2     plot(g,          layout=layout,          vertex.label="",          vertex.size=V(g)$size*0.1,              edge.label=E(g)$weight)          # Create data frame for analysis: one row for each pair of vertices     v_pairs <- get.edges(g, E(g))     v_pairs <- cbind(v_pairs, E(g)$weight)         v_pairs <- cbind(v_pairs, V(g)$year[v_pairs[,1]] * V(g)$year[v_pairs[,2]])         v_pairs <- cbind(v_pairs, abs(V(g)$year[v_pairs[,1]] - V(g)$year[v_pairs[,2]]))     v_pairs <- data.frame(v_pairs)          # Assign names to df     names(v_pairs) <- c("wedding_i", "wedding_j", "shared_ppl", "mult_size","time_diff")          # Regression analysis     model <- lm(shared_ppl~mult_size+time_diff,data=v_pairs)     options(scipen=999)     summary(model)          # Create data frame from prediction based on model     prediction <- predict(model,newdata=v_pairs,interval='confidence',level=.95)     prediction <- cbind(prediction, E(g)$weight)          # Function to calculate offset from .95 confidence interval of prediction     calc_offset <- function(x) {       if (x[2]<=x[4]&x[4]<=x[3]) {return(0)}       if (x[4]>x[3]) {return(x[4]-x[3])}       if (x[4]<x[2]) {return(-(x[2]-x[4]))}     }          # Calculate offset based on formula above     offset <- round(as.numeric(apply(prediction, 1, calc_offset)),digits=2)          # Add results to graph     E(g)$offset <- offset          # Plot 3     plot(g,          vertex.label=V(g)$size,              layout=layout,              vertex.size=V(g)$size*0.1,          edge.label=E(g)$offset)      FINAL NOTE: This question is part of a line of thought which includes also these questions here, here, here and here.


title: Coverage rates of confidence intervals in reality

text: One proves mathematically that if assumptions of a model are satisfied, then the coverage rate of a $100p\%$ confidence interval is $100p\%$. But then statistics gets applied to the world, where model assumptions may not be satisfied. Are there any studies comparing the coverage rates of confidence intervals applied to the real world with theoretical coverage rates?


title: Amazon interview question—probability of 2nd interview

text: I got this question during an interview with Amazon:   * 50% of all people who receive a first interview receive a second interview   * 95% of your friends that got a second interview felt they had a good first interview   * 75% of your friends that DID NOT get a second interview felt they had a good first interview If you feel that you had a good first interview, what is the probability you will receive a second interview? Can someone please explain how to solve this? I'm having trouble breaking down the word problem into math (the interview is long over now). I understand there may not be an actual numerical solution, but an explanation of how you would walk through this problem would help. edit: Well I did get a second interview. If anyone is curious I had gone with an explanation that was a combination of a bunch of the responses below: not enough info, friends not representative sample, etc and just talked through some probabilities. The question left me puzzled at the end though, thanks for all of the responses.


title: Plain English explanation of Bernoulli mixture models?

text: Not exactly the most accessible explanation can be found here, but I'm looking for something more intuitive, examples of applications and so on. Help is much appreciated.


title: Is it most beneficial to become fluent in R, or only learn select packages that do what you need?

text: I don't know enough about R. Here are two possibilities that I envision regarding the use of R:   1. You are fluent in R. You are able to readily use R for almost any kind of analysis. If you do not have experience with a package, you can learn it very quickly (minutes) and use it, assuming you are familiar with the analysis in general.   2. You are not fluent in R. You learned and use a few specific R packages based on your needs. You would not be able to easily use R for analyses that do not involve the packages you already know.  I know this is a false dichotomy, but it hope it serves to illustrate my questions. My questions are:   1. Are most R users fluent, or not fluent as described above?   2. I appreciate that becoming fluent in R would be most beneficial, but I also recognize that this would involve a significant time investment. For what types of users is such an investment worth it?    3. How many fluent users of R also resort to using other software?  I know these questions are relatively subjective, but I suppose so are other questions about R (e.g., "What R packages do you find most useful in your daily work?"). Edit: I will mention why I am asking. I have never used R and I perform relatively complex analyses regularly (e.g., propensity score matching, structural equation modelling, latent growth analyses, item response theory). I use three or four different statistical packages to perform these analyses. I have been lucky to find these software packages to do what I need, but I feel like eventually I won't be that fortunate. I am constantly debating learning R, despite being so very busy otherwise.


title: z-transformation of Likert-scale questions before running a pca

text: Does it makes sense to standardize (z-transformation) all ordinal variables (I only have Likert-scale questions) before running a PCA? Or is the z-transformation only used when having variables with different scales?


title: Parametric modelling of customer lifetime value using survival (Cox) analysis

text: I have half a year long data of purchases in e-commerce site (100K purchases by 60K customers). In simplest A/B testing framework random customers got a discount on the next purchase after order completion ( `disc` = 0/1 below). I want to estimate to what extent `disc` influenced an interval between orders (`diff`). Following KISS principle, I just drop all 1-time customers and regress ln(`diff`) on `disc`, but don't observe the effect at all. I have two obvious problems:   1. data is very censored, 60% of customers appeared only once   2. selection bias - frequent buyers had more chances to get a discount To address (1), I turn to Cox-model `coxph(Surv(diff, event) ~ disc + cluster(customer_id))` (and observe the effect!), but can't figure out whether it is the best method to handle multiple failure times(purchases) per customer. For (2), I'm thinking of introducing lagged `diff`, but don't know how to do it in a censored case. There were number of relevant discussions (RFM & customer lifetime value modeling in R , Survival Analysis with Multiple Events), but I fail to find a solution for my problems. There is also BTYD package, but it's not parametric. Guess this is very standard question, but can't find a step-by-step (CrossValidated) guide


title: Best open source data visualization software to use with PowerPoint

text: What is the best open source data visualization software? I require the following:   1. Can import data from Microsoft Excel (importing data from Oracle databases would be good too, but this is not mandatory).   2. Charts generated by the software can be exported to Microsoft PowerPoint (copy and paste is fine with me).   3. Open source & easy to use.


title: Finding motif in an on-line (real time) time series

text: I have a data stream (one integer value per second):               480, 481, 479, 482,...      If I show these values in a 2D-diagram, I have something like (for limited values): ![enter image description here](http://i.stack.imgur.com/1Mf37.png) Now I want to discover a similar motif in this data-stream. I have read about Symbolic Aggregate approXimation, but as I understand, I can use it for offline data (if I have 10,000 time series values) and now my question is:   * How can I find a motif in a data stream?   * Actually I want to this in R. Do you know of any package, with which I can implement this job in R?


title: Testing the accuracy of a string matching algorithm?

text: I have a list of records of people that received a treatment at different points in time; some individuals have received the treatment multiple times. I use a string matching algorithm to identify which records refer to the same person. My question is: how can I test how accurate the matching is? I can look at a random sample of matched records and manually check for errors and calculate a p-value based on this, but that ignores the error from not matching records that should have been matched. I'd be grateful for any pointers!


title: How to estimate an mlogit (R package) model with a fixed (offset) variable

text: I am trying to estimate a multinomial logit model using an offset variable with the mlogit package of R Using the syntax here [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/formula.html] it seems that the correct way to do this is to use something like this               v <- mFormula(choice ~ gc + log_gc + offset(log_attr) + 0)      which you would then pass to the model estimation as               model <- mlogit(v, mlogit_data)      Printing the coefficients seems to indicate that the parameter has not been estimated               model$coefficients               gc       log_gc      -0.003023088 -0.477115780      attr(,"fixed")         gc log_gc       However these parameters are exactly the same as the parameters I get if I don't include the variable log_attr at all... ie               v_no_log_attr <- mFormula(choice ~ gc + log_gc + 0)     model_no_log_attr <- mlogit(v_no_log_attr, mlogit_data)     model_no_log_attr$coefficients               gc       log_gc      -0.003023088 -0.477115780      attr(,"fixed")         gc log_gc       FALSE  FALSE       Given that the data has significant variation in this variable, I feel that this is unexpected and that the mlogit package is simply ignoring the offset variable. ![variation_in_data](http://i.stack.imgur.com/thCvR.png) Is this possible to do with this package? Do I have a misunderstanding of how to perform the estimation?


title: Parametric transition matrix in Markov Chains

text: I am trying to model a discrete-time MC with transition probabilities that depend on some function of parameters i.e $p_{ij} = f(X_0,X_1)$. Suppose we take a log-linear model where $p_{ij} = e^{\beta_0+\beta_1X_0+\beta_2X_1}$. Now I have aggregate observations for the state distribution q. I want to estimate $\beta_0,\beta_1,\beta_2$ which minimizes the expression, $\sum_{t=1}^T||q^tP-q^{t+1}||_2$ where $P$ is the matrix with $p_{ij}$ entries I have been trying to find material which deals with such parametric model for transition probabilities but haven't found anything really useful. I'd appreciate some suggestions or links to some resources about how to go about solving this. Thanks. PS: Have also posted this on math.stackexchange. Was confused which is a more appropriate place.


title: Confused about the basic maths in Incremental/Uplift modelling

text: I have recently gotten to study Incremental modelling, which is used in Marketing to study the incremental impact of a certain action on a treatment group against a control group. However, I am very confused about the very basics of this form of modelling: specifically about the formulae used in its calculations. Consider the following example, please:               |CampaignName | Group Type | TotalPeople | Respondents | AverageSpend |     |-------------|------------|-------------|-------------|--------------|     |ToysRus      | Control    | 80          | 10          | 30           |     |ToysRus      | Treatment  | 800         | 130         | 39           |      Supposing that we were to find out:   1. The response rate for the treatment group   2. The incremental response rate   3. The incremental spend per customer   4. Total incremental spend due to targetting I understand that for **1.** , we simply use the formula:               No. of respondents in the treatment group/Total no. of people in treatment group      And this will be the rate of response in treatment group - the percentage of people from the treatment group that responded to our 'treatment'. Likewise we can find the rate of response for the control group. For **2.** , I am confused. I understand that the formula is:               Response rate from treatment group - response rate from the control group      But I personally used the following formula in my practice:               No. of respondents in the treatment group - no. of respondents in control group / Total no. of people in the treatment group      I know that the answer I get due to this is certainly wrong, but I am unable to convince myself as to why the formula is wrong. It seems intuitive to me on first glance. But I think there is basic Math here at play, which very unfortunately I am unable to see presently. I think perhaps it has to do with qualitative and quantitative differences between the groups? That is, one is a control and the other a treatment group, and both are of different sizes, so it makes less sense to subtract them straight away? Like comparing apples and banana? Even if we do get to apply my formula, I don't see how its interpretation would be possible. Yet, I am unable to convince myself! For **3.** , I think the following formula would apply:               Average Spend per customer in the treatment group - average spend per customer in the control group.      For **4.** , I am unsure, but I think we would use:               The incremental response * the total no. of people in the control group      But how to understand this? Are we extrapolating to the control group with our incremental spend, to project how much the 'untreated' group or the control group would spend if we direct our offers/advertisement/etc. at them? **Lastly** , I would be grateful if people answering this question could also suggest further reading for me. Perhaps you feel I need to study the basics of proportion, or percentages, in order to get myself comfortable with such basic calculation in Incremental modelling? Or perhaps you feel that there are certain books in Marketing or Incremental modelling that would help me? Thanks a lot!


title: What's the likelihood function in this case?

text: I have a data set with a binary response variable as non-zero survival cells (y=0) or zero survival cells (y=1), and I would like to build a model with explanatory variable X. I know that the most straightforward method is to use logistic regression. However, in my case, it is more biologically sound to use a Poisson model. Thus, the number of survived cells follows a Poisson distribution $f(k,\lambda)$, and the response (probability of having zero survival cell) is the Poisson distribution at $k=0$. Therefore, my model can be written as $y=f(0,\lambda)=e^{-\lambda}$, where $\lambda=g(\beta,x)$. $\lambda$ is a function of $x$ with parameter $\beta$. I am trying to estimate $\beta$ using MLE. Now I am bit confused about which likelihood function I should use. Shall I assume that $y$ follows a binomial distribution: $y\sim B(1,e^{-\lambda})$, so that I could conduct MLE on the likelihood function from this binomial distribution? What about the Poisson distribution? It is not used in the estimation? Suppose $\beta^{*}$ is the estimated parameter, which gives $\lambda^{*}$. This means that the probability of having zero survival cells follows a $B(1,e^{-\lambda^{*}})$ distribution and the survival cells follows a Poisson $f(k,\lambda^{*})$ distribution?


title: Combining the Counts and P-Values

text: I have experimental and computational data for n genes in data matrix. The rows are the genes and columns are its values. Now, the data is from two sources:   1. Computational Prediction Results in form of p-values from Fisher Exact Test.   2. Count Data results from the experiments across 4 different time points. This count data is basically is the number of signals you get per time point and it could be any positive number. There is no fixed bound.  I am interested in ranking the genes and therefore, I would like to know if there is a way to combine the computational prediction and experimental count result in the form of a statistic or some through some equation. I thought of sorting them in terms of p-values and then in terms of count data but I think it is very naive to do that way.


title: Interpretation when only the squared term is significant

text: I am using a Poisson model. My theory suggests a positive $X_1$ and a negative $X_1^2$. However, the results show an insignificant, positive estimate for $X_1$ and a significant, positive estimate for $X_1^2$. How should I interpret the result when $X_1$ is not significant but $X_1^2$ is significant? Thanks!


title: Entropy in decision tree classification

text: In the decision tree based classification technique. What is the difference among the different approaches like entropy, gini index? When to use entropy and when to use gini index?.


title: Drawing with repetition

text: I thought I was fairly good at combinatorics, but this puzzle is giving me some trouble! I have a bag with $N$ balls. I pick one at random, mark it and put it back inside the bag. I repeat this operation $D$ times. What is the probability that I will end up with exactly $C$ clean balls?


title: Building training set from unlabeled data

text: I want to build classifier using naive bayes. I know that naive bayes is supervised learning that requires labeled data for training set. However, I only have data with no label on it. Is there any method to give label to unlabeled training data or we only need to label it manually? I want to label my data with two classes. My data is all in continuous variable and has three attributes to determine the class. Edit: > I've found useful java program which uses MVE for multivariate outlier > detection in this link : http://www.kimvdlinde.com/professional/mve.html > > I use that program to label the data which is all in continuous variable. > What kind of method that I can use to label data if my data has three > categorical attribute? Is it still possible to label it manually?


title: svm manual implementation to find missing values

text: how svm algorithm used to impute the missing data in a dataset. I need a manual implementation of svm algorithm to find missing values with example


title: When I can model count data as continuous

text: I have time series of several variables of 60 or so rows of count data. I want to do a regression model `y ~ x`. I've chosen to use a Quasipoisson & Negative Binomial GLMs as there's overdispersion etc.               x     Min.   : 24000       1st Qu.: 72000       Median :117095       Mean   :197607       3rd Qu.:291388       Max.   :607492            y     Min.   : 136345     1st Qu.: 405239     Median : 468296     Mean   : 515937     3rd Qu.: 633089     Max.   :1218937      The data itself are very high and so it may be best to model these as count data (this is what I'm trying to investigate - at which point I can model count data as continuous). It seems to be very common practice, what I want to know is the motivation for this? Are there any texts that actually show the problem of modelling high count data with Poisson distribution? Perhaps something that shows the factorial in the distribution makes things difficult.


title: Does this kind of overlap between in-sample data and forecast cause inflated $R^2$?

text: I am using a simple UIP model to forecast exchange rates using interest rates with a twelve month horizon. The equation I use is: $E(t+12) - E(t) = α + β(I*(t) - I(t))$. I apply OLS linear regression on the in-sample period to get $α$ and $β$ and then use these values to calculate the forecast. For example, if I use data points $t=1$ to $t=24$ to get $α$ and $β$ with OLS, then $E(t=25) - E(t=13) = α + β(I*(t=13) - I(t=13))$. This gives me a forecast for $E(t=25) - E(t=13)$. Once I have a whole time series of forecasts (using an automated rolling regression, but nevermind that), I can use another linear regression to compare the predicted values with the actual values and see how well I did after the fact based on the $R^2$ value of this new regression. The problem I see with this approach is that I am using using $E(t=1)$ to $E(t=24)$ in order to forecast $E(t=25) - E(t=13)$: That's an overlap of 11 months between the forecast and the in-sample period! So I'm worried that this will make my forecasting ability look much better than it actually is. But maybe I'm not thinking about this right. Is this actually an issue or does it not affect the result at all? If it really is an issue, I can get rid of the overlap by reducing the in- sample period by 12 months and forecasting $E(t=25) - E(t=13)$ using $E(t=1)$ to $E(t=12)$. However this will weaken my forecasting ability so I'd rather not. Please tell me what I should do! Many thanks in advance to anyone who can help me out.


title: Comparing GLM models using predict

text: Suppose I have two models created by calling `glm()` on the same data but with different formulas and/or families. Now I want to compare which model is better by predicting on an unknown data. Something like this:               mod1 <- glm(formula1, family1, data)     mod2 <- glm(formula2, family2, data)     mu1 <- predict(mod1, newdata, type = "response")     mu2 <- predict(mod2, newdata, type = "response")        1. How can I tell which of the predictions `mu1` or `mu2` is better?   2. Is there some simple command to compute the log likelihood of a prediction?


title: Is there some theory of SVMs with infinitely many data?

text: I am trying to understand what does it means to have a (linear) SVM classifier (with soft margins) given the generative model of the data. And I realize I have not seen any paper on it, nor can I extend the theory by myself. At least two problems:   * the penalty term $C \sum \xi_i$ will become unbounded. There will be infinitely many finite $\xi_i$.   * there will be infinitely many restrictions $y(w.x_i - b) \ge 1-\xi_i$  Any suggestions or references? thanks


title: Estimating probability or frequency with low N?

text: I am trying to estimate the probability of an event using a low number of observations. The naive estimator $\hat{p} =\frac{\text{number of positive observations}}{\text{total number of observations}}$ works well when the total number of observations is big enough, but if you have only a few observations, there is a decent chance that you will erroneously conclude to a 0 or 100% probability. I suppose you could set a prior distribution on the estimated probability (say, uniform), and look for better estimators. I suppose this problem has already been tackled many times, so where should I look?


title: Binomial / multinomial logistic regression or chi-squared

text: I am currently doing my fourth year thesis examining the moral stages of children/adolescents. The DV is the moral stage (categorical variable, stages 1, 2, 3, or 4) and the IVs are age group (I currently have 2, but may potentially have 3) and gender, so my design will be either a 2 x 2 or a 3 x 2. I have been told by my supervisor that a logistic regression is the way to go, but it hasn't been covered in our previous statistics classes (at least not multinomial logistic regression) and he said I can otherwise use 2 chi square tests, one for each gender. I then wouldn't be able to examine any interaction between the age groups and gender. Could anyone tell me what the best method to use would be, and if I was to use logistic regression, would it be multinomial rather than binomial, considering my DV is not dichotomous?


title: Geometrical interpretation of correlation of a variable and the residual

text: I came across this blog post about least angle regression, and at a point he says: >   * Find the variable $x_1$ most correlated with the residual. (Note that > the variable most correlated with the residual is equivalently the one that > makes the least angle with the residual, whence the name.) > >   * Move in the direction of this variable until some other variable $x_2$ > is just as correlated. > > But I couldn't visualize it. By "residual" he is referring to the vectors minus the average of the vectors, am I right? What could be a geometrical interpretation of "variable correlation with the residual"? What would be a, let's say, 2D example of this interpretation?


title: Model Averging for Negative Binomial GLMM in R

text: Ive been trying to find a way to average negative binomial GLMM using MuMin R Package but it seems not to work for negative binomial GLMM, any alternatives? thanks in advance!


title: What model is most appropriate to estimate the probability of something occurring?

text: e.g. estimate the probability that of 10 (ten) small businesses 3 (three) or less will fail in the first year of operation


title: Interpreting estimates of a Gamma distribution

text: I have 399 claim amount of data, i want to fit them using a Gamma distribution. I have done so in Matlab software However i am getting this as my result: > > [phat,pci] = gamfit(data) phat = 1.0e+004 *               0.0001    2.8951      why i am getting the 1.0e+004*? Thanks


title: R detect increasing/decreasing trend of time series

text: I have lots of time series with periods: day, week or month. With `stl()` function or with `loess(x ~ y)` I can see how trends of particular time series look. I need to detect if trend of time series is increasing or decreasing. How can I manage that? I tried to compute linear regression coefficients with `lm(x ~ y)` and play with slope coefficient. (`If |slope|>2 and slope>0 then` increasing trend, `else if |slope|>2 and slope<0` – decreasing). Maybe there is another and more effective method for trend detection? Thank you! For example: I have `timeserie1`, `timeserie2`. I need a simple algorithm that would tell me that `timeserie2` is an increasing algorithm, and in `timeserie1`, the trend isn't increasing or decreasing. What criteria should I use? `timeserie1`:               1774 1706 1288 1276 2350 1821 1712 1654 1680 1451 1275 2140 1747 1749 1770 1797 1485 1299 2330 1822     1627 1847 1797 1452 1328 2363 1998 1864 2088 2084  594  884 1968 1858 1640 1823 1938 1490 1312 2312     1937 1617 1643 1468 1381 1276 2228 1756 1465 1716 1601 1340 1192 2231 1768 1623 1444 1575 1375 1267     2475 1630 1505 1810 1601 1123 1324 2245 1844 1613 1710 1546 1290 1366 2427 1783 1588 1505 1398 1226     1321 2299 1047 1735 1633 1508 1323 1317 2323 1826 1615 1750 1572 1273 1365 2373 2074 1809 1889 1521     1314 1512 2462 1836 1750 1808 1585 1387 1428 2176 1732 1752 1665 1425 1028 1194 2159 1840 1684 1711     1653 1360 1422 2328 1798 1723 1827 1499 1289 1476 2219 1824 1606 1627 1459 1324 1354 2150 1728 1743     1697 1511 1285 1426 2076 1792 1519 1478 1191 1122 1241 2105 1818 1599 1663 1319 1219 1452 2091 1771     1710 2000 1518 1479 1586 1848 2113 1648 1542 1220 1299 1452 2290 1944 1701 1709 1462 1312 1365 2326     1971 1709 1700 1687 1493 1523 2382 1938 1658 1713 1525 1413 1363 2349 1923 1726 1862 1686 1534 1280     2233 1733 1520 1537 1569 1367 1129 2024 1645 1510 1469 1533 1281 1212 2099 1769 1684 1842 1654 1369     1353 2415 1948 1841 1928 1790 1547 1465 2260 1895 1700 1838 1614 1528 1268 2192 1705 1494 1697 1588     1324 1193 2049 1672 1801 1487 1319 1289 1302 2316 1945 1771 2027 2053 1639 1372 2198 1692 1546 1809     1787 1360 1182 2157 1690 1494 1731 1633 1299 1291 2164 1667 1535 1822 1813 1510 1396 2308 2110 2128     2316 2249 1789 1886 2463 2257 2212 2608 2284 2034 1996 2686 2459 2340 2383 2507 2304 2740 1869  654     1068 1720 1904 1666 1877 2100  504 1482 1686 1707 1306 1417 2135 1787 1675 1934 1931 1456 1363 2027     1740 1544 1727 1620 1232 1199      `timeserie2`:                122  155  124   97  155  134  115  122  162  115  102  163  135  120  139  160  126  122  169  154      121  134  143  100  121  182  139  145  135  147   60   58  153  145  130  126  143  129   98  171      145  107  133  115  113   96  175  128  106  117  124  107  114  172  143  111  104  132  110   80      159  131  113  123  123  104  101  179  127  105  133  127  101   97  164  134  124   90  110  102       90  186   79  145  130  115   79  104  191  137  114  131  109   95  119  173  158  137  128  119      109  120  182  140  133  113  121  110  122  159  129  124  119  109  108   95  167  138  125  105      139  118  115  166  140  112  116  139  121  109  164  135  118  121  112  111  102  169  136  151      132  135  130  112  156  134  121  116  114   91   86  141  160  116  118  112   84  114  165  141      109  123  122  110  100  162  145  121  118  115  107  103  162  142  130  139  134  121  118  164      147  125  120  134  107  130  158  141  144  148  124  135  118  212  178  154  167  155  176  143      201  170  144  138  152  136  123  223  189  160  153  190  136  144  276  213  199  211  196  170      179  460  480  499  550  518  493  557  768  685  637  593  507  611  569  741  635  563  577  498      456  446  677  552  515  441  438  462  530  699  629  555  641  625  544  585  705  584  553  622      506  500  533  777  598  541  532  513  434  510  714  631 1087 1249 1102  913  888 1147 1056 1073     1075 1136  927  922 1066 1074  996 1189 1062  999  974 1174 1097 1055 1053 1097 1065 1171  843  441      552  779  883  773  759  890  404  729  703  810  743  743  946  883  813  876  841  742  715  960      862  743  806  732  669  621


title: Augmented Dickey-Fuller test

text: I am now working with Stata and I found out I have a unit root in my regression. How can I correct for this, because I can read everywhere what the test does but if you have a unit root, what to do next?? Sorry, it is one of my first times working with Stata so it is all really confusing..


title: Calculating 95% CI of a proportion of 0%

text: I have a sample of 72 individuals. 0 of these had the trait I'm investigating. How do I calculate the 95% CI when the proportion is 0%?