title: How to draw B-spline plots via the PROC GLIMMIX procedure in SAS?

text: Recently I found that the `PROC GLIMMIX` procedure in SAS added a statement `effect`, which can handle B-spline in models. I tried it, and the program looks like:               proc glimmix data=temp;         effect t = spline(time);         model y = x1 x2 t / solution;         random intercept / subject=id type=un;     run;      It is ok to get outputs, but what I am really interested is the B-spline plot. I went through the entire manual, but did not find any information. Its `ODS GRAPH` function cannot generate B-spline plots automatically. Also, its `ODS OUTPUT` function does not provide any output to draw B-spline plots outside. I am wondering whether any one used to use the `PROC GLIMMIX` procedure to draw a B-spline plot? Any information is pretty appreciated!


title: Wilcoxon rank sum test in R

text: I have results from the same test applied to two independent samples:               x <- c(17, 12, 13, 16, 9, 19, 21, 12, 18, 17)     y <- c(10, 6, 15, 9, 8, 11, 8, 16, 13, 7, 5, 14)      And I want to compute a Wilcoxon rank sum test. When I calculate the statistic $T_{W}$ by hand, I get: $$ T_{W}=\sum\text{rank}(X_{i}) = 156.5 $$ When I let R perform a `wilcox.test(x, y, correct = F)`, I get:               W = 101.5      Why is that? Shouldn't the statistic $W^{+}$ only be returned when I perform a _signed_ rank test with `paired = T`? Or do I misunderstand the rank sum test? **How can I tell R to output $T_{W}$** * * * As part of the test results, not through something like:               dat <- data.frame(v = c(x, y), s = factor(rep(c("x", "y"), c(10, 12))))     dat$r <- rank(dat$v)     T.W <- sum(dat$r[dat$s == "x"])      * * * _I asked a follow up question about the meaning of theDifferent ways to calculate the test statistic for the Wilcoxon rank sum test_


title: Why does Lasso do better than SVM?

text: This is a soft-question: I have been evaluation various regression techniques over a regression dataset that I have. I am surprised by the fact that cross- validated RMSE of Lasso is better than SVM and Random Forest in my case. Can this happen?- I mean, I believed that a non-linear modeling technique like random forest or SVM would do better than a linear model like Lasso. Is that really possible!?


title: How to combine a known and unknown covariance in a linear model

text: I am faced with a problem where I have a situation that can be fit with a linear multivariate regression type model. But, there are two types residual error that will occur from the lack of fit. Right now, my $\mathbf{Y}$ data is uncorrupted by noise source $\mathbf{R}_{2}$, but will be in the "future". Therefore, I need to incorporate the extra uncertainty it will cause into my test statistic. Is this a multivariate random effects model? How do I deal with these two sources of variance? So far, what I have done is as follows. I can write: $\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{R}_{1} + \mathbf{R}_{2}$, where the $\mathbf{R}$'s are a matrices of the residuals, such that of course, $\mathbf{Y} - \mathbf{X}\mathbf{B} = \mathbf{R}_{1} + \mathbf{R}_{2} = 0$ if I have a perfectly fitting model. So, the distributional assumption on each $\mathbf{R}$ is $\mathbf{R} \sim N(\mathbf{0},\mathbf{I} \otimes \mathbf{\Sigma_{R}}) $. Except, $\mathbf{\Sigma_{R_{1}}}$ is unknown while $\mathbf{\Sigma_{R_{2}}}$ is known. I can get an unbiased estimate of $\mathbf{\Sigma_{R_{1}}}$ by computing $\mathbf{\hat{\Sigma}_{R_{1}}} = \frac{1}{n-p} \mathbf{Y} [\mathbf{I} - \mathbf{X}(\mathbf{X}^{'}\mathbf{X})^{-1}\mathbf{X}^{'}]\mathbf{Y}^{'}$. This looks odd, but it is because $\mathbf{Y}$ has "not yet" been corrupted by the noise due to the known $\mathbf{R}_{2}$. However, when I compute test statistics based on this model (eg. null hypothesis: $\mathbf{B}_{1,2} = 0$ and so on), I want this "future" noise to be taken into consideration, where its covariance structure is known. aka, I have at this moment the matrix $\mathbf{\Sigma_{R_{2}}}$ on hand. So, how can I combine an known and unknown covariance matrix, one I know, and one I can easily estimate? Thank you for your time!


title: How to do logistic regression when one of the factors has unequal number of levels across the other factor?

text: I am analyzing an experiment that has the following design:                   Factor1  Factor2           Level1  Level1                  Level2          Level2  Level1.1                  Level1.2                  Level2      In other words, Factor2 has different number of levels across Factor1. I'd like to use logistic regression, but how should I set it up? Is it legitimate to run two regressions ignoring one of the levels in each (e.g. Level1.1 or Level1.2)?


title: What does it indicate when the Spearman correlation is a definite amount less than Pearson?

text: I have a bunch of related datasets. The pearson correlations between pairs of them are typically definitely larger than the spearman correlations. That suggests any correlation is linear, but one might expect that even if the pearson and spearman were the same. What does it mean when there is a definite gap between the pearson and the spearman correlation and the pearson is larger? This seems to be a consistent feature across my datasets.


title: What's the common way to measure the "stability" or "volatility" of a distribution? Variance seems insufficient

text: for example, in a discrete distribution A               P(X) = 0.999  when x = 0     P(X) = 0.001  when x = 1000          E(X) = 1     Var(X) = 1000-1=999      compared to B               P(X) = 0.6  when x = 0     P(X) = 0.2  when x = 1     P(X) = 0.2/2  when x = 2     P(X) = 0.2/4  when x = 4     P(X) = 0.2/4  when x = 8          E(X) = 1 =A     Var(X) = 0.2+0.4+0.8+3.2-1=3.6 <<A      However B seems more volatile(instable) than A. What is the common correct stats to measure this sort of thing? Update: Entropy might be a good candidate, however it doesn't look into the information itself (it only cares about the chance). I feel the time series graph of a fair coin which could score (0 or 1) is less volatile / more stable than a fair coin which could score (0 or 100).


title: Are there any simple techniques to remove/detect fluctuations during training in matlab plot?

text: Image below is a Matlab plot from 3-layer perceptron learning with 2 hidden nodes in both hidden layers. Are there simple technique or simple coding to remove or detect these fluctuations during training? What I've tried only is to reduce step size but I had to compensate that with increasing number of iterations. Please advice ![probabilistic plot](http://i.stack.imgur.com/9Pc0e.png)


title: What is the difference between spatial dependence and spatial (auto)correlation?

text: It seems to me that in literature these terms are used synonymously. Does spatial autocorrelation strictly refer to linear dependence, or does spatial autocorrelation refer to the set of measures of spatial dependence (moran's i,...)? Do you know of any precise definition of both terms? Thank you very much


title: Using offset in binomial model to account for increased numbers of patients

text: Two related questions from me. I have a data frame which contains numbers of patients in one column (range 10 - 17 patients) and 0s and 1s showing whether an incident happened that day. I'm using a binomial model to regress probability of incident on number of patients. However, I would like to adjust for the fact that when there are more patients, there will inevitably be more incidents because the total amount of patient time on the ward is higher on that day. So I'm using an offset binomial model like this (R-code):               glm(Incident~Numbers, offset=Numbers, family=binomial, data=threatdata)      My questions are:   1. Is it okay to have exactly the same variables predicting and in the offset? I want to partial out the tonic increase in incident probability and see if there's anything left, essentially. It makes sense to me but I'm a little cautious in case I'm wrong.   2. Is the offset specified correctly? I know that in poisson models it would read              offset=log(Numbers)      I don't know if there's an equivalent here and I can't seem to find any binomial offsets with Google (major problem being that I keep getting negative binomial which of course is no good).


Let $X$ and $X_{i}$, $i=1,...,n$ be random variables on a probability space $(\Omega , \mathcal F,P)$. Show that the best mean square estimator of $X$ given $(X_{1},...,X_{n})$ is $\hat X =E[X|\sigma(X_{1},...,X_{n})]$


title: Correlation of one binary input series and continous response series

text: I have one input series of type 0/1 (it is an intervention time series) and some metric response series. I need to have a look at cross correlations. Normally I would use Box-Jenkins-technique. But now I have to prewhiten an input series which is binary. Can I handle it like a metric one? Does it make sense to look at the autocorrelations and fit an ARMA-model in the "normal" way? Or would it be possible to fit ARMA-models for the metric response variables and then fit the same model to the input variable and then look at the cross correlations? Nearly all of my response series are white noise - but as I understand I have to fit an ARMA-model to the input if this is not white noise anyway. Is this correct or is it sufficiently for the response to be white noise? Would be very happy if anyone can help me! Thank you very much.


title: When is the median more affected by sampling error than the mean?

text: I'm writing a paper on making probability estimates, and it's been asserted to me that I should take the median of the estimates given by my participants, rather than the mean. I've been told I should do this because the mean is more affected by sampling error than the median. Why is this? Is this something that is always true, or which only holds under certain circumstances?


title: Pearson correlation between discrete variable that's mostly 0 and a standard normal variable

text: Suppose I want to estimate the correlation between $X\sim N(0,1)$ and $Y$, where $Y \in \\{-1,0,1\\}$ and is equal to zero for 99 per cent of the sample. Sample size is 10 million. What are the properties of the Pearson correlation in this instance (bias, usefulness, etc)? Would an alternative correlation estimator be better (e.g. a rank correlation)?


title: Which analysis to use for ratio of events after/before diagnose, between positive diagnoses and negative?

text: I am trying to figure out what kind of analysis is needed to answer this kind of question: > **Is having a positive diagnosis in "D" changing the number of events "E" > after the diagnosis?** I want to find the ratio between number of events after and before diagnosis ($R_{ab}$) and find out if the ones who had "positive" diagnosis had $R_{ab}$ greater than the ones diagnosed "negative". An example could be that $D \equiv \text{cancer}$ and $E \equiv\text{ number of visits to hospital in 2 years}$. Then I would classify people as having cancer (positive) or not (negative). For a positive case $i$ we can look for negatives in similar conditions (age, date diagnosis, region.). Say we have a positive $C_i$ and the "similar" negatives $C_{Ni}$. We would measure number of events in two years after diagnosis ($N_{ai}$), and number of events in two years before diagnosis ($N_{bi}$). We are interested in $R_{abi} = N_{ai}/N_{bi}$ and compare with $R_{ab}$ for the negative diagnoses: $$ RR_i = \frac{R_{abi} (\text{positive diagnosis})}{ R_{abi} (\text{similar negative diagnosis}) } $$ There are many possible $C_{Ni}$ for each $C_{i}$, so we could choose randomly a few or take all of them. One approach would be to do a _case control study_ , making _exposure = diagnosis_ and $\text{"disease"} = R_{ab} > 1$ . Then we could treat the data as a case control study. But I am not sure if this is a correct approach, because for example, time of diagnosis is influencing both the case-control definition and the exposure, etc. Any suggestion or reference to similar analysis?


title: Calculating confidence intervals when binary input variable equals zero

text: I am building a multivariate regression model and trying to use the model to predict a High, Medium, and Low estimated outcome for each individual in a group of people. I am using the 95% confidence interval to get the High/Medium/Low estimates. This means constructing a High, Medium, and Low estimate for each variable in the individual's profile, and then summing the estimates together in the appropriate category to get a High, Medium, and Low overall outcome estimate. To be explicit, here's how I see High/Medium/Low mapping to the 95% confidence interval, where B is the coefficient and SE is the standard error: High = (B*X)+(1.96*SE)   Medium = (B*X)   Low = (B*X)-(1.96*SE) Again, the High estimated outcome is the sum of all the High predictions for the input variables, the Medium outcome is the sum of all the Medium predictions, and so forth. **The "Problem":** Some of my input variables are binary indicators (with values of 0 or 1). I'm not sure how to handle the binary indicators with values of 0 when summing the input variables to generate the High and Low predictions. **My Question:** When a binary indicator variable equals 0, can I leave their High and Low estimates as simply 0, or do I have to use the High/Low formulas listed above [e.g., (B*0)+(1.96*SE)]? Put another way, is the prediction for a 0-value indicator variable always 0, regardless of High/Medium/Low? **Why It Matters:** If I assume that the predictions for the 0-value binary indicators always equal 0, I get a much more accurate High/Medium/Low range of estimated outcomes. If I adhere to the formulas listed above regardless of the binary indicator's value, I get a much broader range of predictions. (This is probably because the first case basically assumes that SE=0 for 0-value binary indicators). This problem has me puzzled, so I would be greatly appreciative of any insights. Happy to illustrate by example as well.


title: R-squared to compare forecasting techniques

text: Is it appropriate when forecasting to use $R^2$ as the measure of how well exponential smoothing fits a data set for the purpose of time-series forecasting? I understand that it is appropriate for regression analysis; could I compare $R^2$ values between a linear regression, logarithmic regression and exponential smoothing as a method of evaluating the best technique for forecasting?


title: Upper 5% value of distribution of the mean

text: The density function of a random variable x is $f(x)=ke^{-2x^{2}+10x}$. Find the upper 5% point of the distribution of the means of the random sample of size 25 from the above population. I need hints to pursue this. My thoughts:   1. Assume normal population. Find the value of k by integrating $\displaystyle \int_{-\infty}^{\infty} ke^{-2x^{2}+10x}dx=1$   2. Find $\mu$ by integrating $\displaystyle \int_{-\infty}^{\infty} xke^{-2x^{2}+10x}dx$   3. Find $\sigma^2$ by integrating $\displaystyle \frac{1}{24}\int_{-\infty}^{\infty} (x-\mu)^2ke^{-2x^{2}+10x}dx$   4. The upper 5% point of the distribution of the means of the random sample will correspond to $\displaystyle z_{0.95}=\frac{(x-\mu)}{\sigma}$.  Should $ z_{0.95}$ be $1.96$ (two-sided) or $1.645$(one-sided)? Finally, replace values of $z_{0.95}$, $\mu$ and $\sigma$ to get the value of x, which is the required answer. Are these steps correct?


title: Model selection on variance parameters (and about REML)

text: Various descriptions on model selection on random effects of Linear Mixed Models instruct to use REML. I know difference between REML and ML at some level, but I don't understand why REML should be used because ML is biased. For example, is it wrong to conduct a LRT on a variance parameter of a normal distribution model using ML (see the code below)? I don't understand why it is more important to be unbiased than to be ML, in model selection. I think the ultimate answer has to be "because model selection works better with REML than with ML" but I would like to know a little more than that. I did not read the derivations of LRT and AIC (I am not good enough to understand them thoroughly), but if REML is explicitly used in the derivations, just knowing that will be actually sufficient (e.g., Chi-square approximation of the test statistic in LRT does not work if a MLE is biased).               n <- 100     a <- 10     b <- 1     alpha <- 5     beta <- 1     x <- runif(n,0,10)     y <- rnorm(n,a+b*x,alpha+beta*x)          loglik1 <- function(p,x,y){        a <- p[1]        b <- p[2]        alpha <- p[3]       -sum(dnorm(y,a+b*x,alpha,log=T))     }          loglik2 <- function(p,x,y){        a <- p[1]        b <- p[2]        alpha <- p[3]        beta <- p[4]       -sum(dnorm(y,a+b*x,alpha+beta*x,log=T))     }          m1 <- optim(c(a,b,alpha),loglik1,x=x,y=y)$value     m2 <- optim(c(a,b,alpha,beta),loglik2,x=x,y=y)$value     D <- 2*(m1-m2)     1-pchisq(D,df=1) # p-value


title: I don't understand what is going on with my t-test and Mann-Whitney results

text: I have 2 samples, A and B. I subjected both samples to extraction and got the amount of yield (the dependent variable). I just want to see if there is a significant difference between the yield for A and B. Each sample had 3 replicates. Sample A:               0.003,0.0029, 0.003      Sample B:               4.832, 4.828, 4.662      Since it's just 2 samples, I went for a two-sample t-test. I tested for normality, which indicated the data was not normal. I did a log transform, but it was still not normal. Anyway I subject the log yield to the t-test and got equal variance and p value 0.00 And I tried non-parametric test (Mann-Whitney) and got p-value 0.046. I don't get why is the p-value of Mann-Whitney so close to 0.05. Can it be that the means of sample A and B do not differ significantly when the values seem so far apart? What am I doing wrong?


What does "model-based" mean for multivariate statistics? E.g. mvabund– an R package for model-based analysis of multivariate abundance data. How is "non-model-based multivariate statistics" then? Thanks!


title: Multivariate statistics vs machine learning?

text: Are multivariate statistics and machine learning solving the same problems? I saw that their books are about the same topics, so I have the impression that they are solving the same problems and probably using the same methods. What are their relations and differences? Can machine learning solve the same problems in univariate statistics? Thanks!


title: What is an interpenetration design?

text: I have come about a quote in Groves & Lyberg (2010, p. 859) saying that "(...) both interpenetration and reinterview designs are required to estimate these terms" (referring to estimation of response deviations, which are the variance component in the total variance of a mean that is not due to sampling variance). I tried to find literature on what is meant by an interpenetration design. Search on CrossValidated gives zero findings and referring to a classic on experimental design as, for example, 'The Theory of the Design of Experiments (by Cox & Reid, 2000) does not give any reference either. What is an interpenetration design? What are its applications? Does it also go under different names? Reference: Groves, R. M., & Lyberg, L. (2010). Total Survey Error: Past, Present, and Future. _Public Opinion Quarterly, 74_ (5), 849–879. doi:10.1093/poq/nfq065


title: Averaging weights learned during backpropogation

text: Is it possible to train a neural network N times, using backpropogation and then average the weights learned to produce a more accurate classifier? My tests are indicating no, but I'm unsure if it is an implementation error? Thank you


title: Is there a quick way to convert z-scores into percentile scores?

text: Does anyone know of a function or R package that can help me turn z scores into percentile scores? The final goal is to classify or rank a group of respondents into four categories based on the heights of their z scores (20% lowest scores, 30%, 30%, 20% highest scores). Any ideas? Many thanks in advance!


title: Wald test in regression (OLS and GLMs): t- vs. z-distribution

text: I understand that the Wald test for regression coefficients is based on the following property that holds asymptotically (e.g. Wasserman (2006): All of Statistics, pages 153, 214-215): $$ \frac{(\hat{\beta}-\beta_{0})}{\widehat{\operatorname{se}}(\hat{\beta})}\sim \mathcal{N}(0,1) $$ Where $\hat{\beta}$ denotes the estimated regression coefficient, $\widehat{\operatorname{se}}(\hat{\beta})$ denotes the standard error of the regression coefficient and $\beta_{0}$ is the value of interest ($\beta_{0}$ is usually 0 to test whether the coefficient is significantly different from 0). So the size $\alpha$ Wald test is: reject $H_{0}$ when $|W|> z_{\alpha/2}$ where $$ W=\frac{\hat{\beta}}{\widehat{\operatorname{se}}(\hat{\beta})}. $$ But when you perform a linear regression with `lm` in R, a $t$-value instead of a $z$-value is used to test if the regression coefficients differ significantly from 0 (with `summary.lm`). Moreover, the output of `glm` in R sometimes gives $z$- and sometimes $t$-values as test statistics. Apparently, $z$-values are used when the dispersion parameter is assumed to be known and $t$-values are used when the dispersion parameter is esimated (see this link). Could someone explain, why a $t$-distribution is sometimes used for a Wald test even though the ratio of the coefficient and its standard error is assumed to be distributed as standard normal? **Edit after the question was answered** This post also provides useful information to the question.


title: Principal component analysis (PCA) on long-tailed data

text: (1) When doing PCA, do you assume the variables to be bell-shaped? Say if I have a bunch of variables, some are bell-shaped but some have characteristic long (right) tails (highly skewed and asymmetric with the median close to the minimum and and a few extremely large values). Should I do log transformation on those long-tailed variables or leave them like that? What are the implications to the results of the PCA with and without transformation? (2) The second question is, when I do principal component regression (PCR), the response variable have a long tail & one of the PC scores also have long tail. Should I log transform either or both of them? NB: Tails in this question are all one-sided.


title: Correlation Clustering

text: Correlation Clustering : Given a signed graph where the edge label indicates whether two nodes are similar (+) or different (−), the task is to cluster the vertices so that similar objects are grouped together. Papers report that task is NP-Complete and suggest Approximation Algorithms. Can one recommend me a fast-implementation of a naive algorithm to do the task?


title: Counterfactuals for Variables with Negative Values

text: Lets imagine I have estimated the following simple linear regression model: $y_{i} = 10 + 0.5x_{i} + \varepsilon_{i} $, and want to work out the counter- factual, or what would $ y_{i}$ be in the absence of $ x_{i}$. The most simple example of this is a simple treatment control case. So $x_{i} = T_{i}$ and $T \in \\{0,1\\}$. In this case there are two possible counter- factual outcomes: 1: Get the treatment, so $ T=1$ and the counter-factual is $y_{i}-0.5T_{i}=y_{i}-0.5$, 2: Don't, so $ T=0$ and the counter-factual is simply $y_{i}-0.5T_{i}=y_{i}$. This is easy enough to understand. However, I am having a lot of trouble with calculating such counter-factuals in the case where $ x_{i}$ takes negative values. In my case, I have a continuous regressor $ x_{i}$, that has values roughly ranging from $-8$ to $+ 1$. To see my issue, lets assume the same regression model as in the above, and try to work out the counter-factual values of $ y_{i}$ for two potential values of $ x_{i} \in \\{-2, -1\\}$. For $x_{1}=-2$, the lower value, the predicted impact of this on the outcome is: $0.5(-2)=-1$. For $x_{2}=-1$, the greater value, the predicted impact of this on the outcome is: $0.5(-1)=-0.5$. Obviously, the lower value of $x$ has less of an impact on the outcome as $-1<-0.5$. However, this is the point where my brain goes dead. How do you translate this into a counter-factual outcome? From before, we had a formula: $y_{i}-0.5x_{i}$. Although, in this case using this formula predicts the opposite. $y_{1}-0.5(-1)=y_{1}+0.5$ $y_{2}-0.5(-2)=y_{2}+1$ So in the absence of $x_{i}$ both values of $y$ would be, greater even though the estimated coefficient is positive. I am really struggling with this. Is there some kind of transformation I can do that allows the counter-factual values of the higher values of $x_{i}$ to shrink while the lower values of $x_{i} $ shrink less (or not at all) as in the $x_{i}=T_{i}\in\\{0,1\\}$ case? To further elaborate, I want to be able to predict a lot of counter-factual values for a continuous treatment effect, but I want to be able to use these counter-factual values to work out how much of the variation in $y$ is caused by $x$. Let us imagine we are talking about the outcome being economic growth. We want to find out how much good political institutions ($x$) explain of economic growth ($y$). What I would like to do is see what these growth rates were like if institutions did not matter. For example, imagine we have estimated the regression model given in the above for 2 countries: A and B. In A, there are good institutions so $x_{A}=4$, and B is bad so $x_{B}=1$. Now let us say that economic growth in A is $y_{A}=4$ and $y_{B}=2$. So in this case, the counter-factuals (what growth would be if $x$ didn't matter) would be $4-0.5(4)=2$ and $2-0.5(1)=1.5$ for A and B respectively. In the counter-factual scenario the sum of growth rates is 3.5, as opposed to the observed 6. This means that in the absence of institutions we would have had ($\frac{3.5}{6} \approx 58\%$) of the the growth. So institutions explain $42\%$ of the growth difference between these two countries. With negative values of $x$ I am not sure how to go about performing such a calculation.


title: Stream classification of time series

text: I have a set of time series $\mathcal{Y}$, and a test time series $T$ for which I need to find the closest matching time series $Y_i \in \mathcal{Y}$. This has to be done online, i.e., $T$ is a stream of measurements and the detection has to happen as each measurement of $T$ is taken. The time series in $\mathcal{Y}$ are of different lengths, and since $T$ is a stream, the length of $T$ is unknown a priori. Moreover, $T$ can match ${Y_i}^{1+}$, i.e., once a best match is found, it cannot hop to other series, but it can repeat itself. I am wondering what distance metric or set of algorithms best suit this kind of application? I'd also appreciate if you can refer me to any references or papers on that topic.


title: Naive Bayes feature probabilities (Do I double count words?)

text: I'm prototyping my own Naive Bayes bag o' words model, and I had a question about calculating the feature probabilities. Let's say I've got two classes, I'll just use spam and not-spam since that's what everyone uses. And let's take the word "viagra" as an example. I have 10 emails in my training set, 5 spam and 5 non-spam. "viagra" appears in all 5 spam documents. In one of the training documents it appears 3 times (this is what my question is about), so that's 7 appearances in spam total. In the non- spam training set, it appears 1 time. If I want to estimate p(viagra | spam) is it simply: **p(viagra | spam) = 5 spam documents contain viagra / 5 spam documents total = 1** In other words, **does the fact that one document mentioned viagra 3 times instead of once really not matter?** * * * Edit: Here's a blog post where the author uses the approach I just laid out: http://ebiquity.umbc.edu/blogger/2010/12/07/naive-bayes-classifier- in-50-lines/ And here's a blog post where the author says: p(viagra | spam) = 7 viagra spam mentions / 8 total mentions http://www.nils-haldenwang.de/computer- science/machine-learning/how-to-apply-naive-bayes-classifiers-to-document- classification-problems And then one of the answers below says it should be: p(viagra | spam) = 7 viagra mentions in spam / total term count in spam Can anyone link to a source that gives an opinion on this?


title: Line with discontinuous sigmoid function in half-normal residuals plot

text: In R, I've plotted the half-normal residuals for a few different models, i.e.               halfnorm(residuals(model_object)))      I notice that one plot takes a distinct S-shape (i.e. it's Sigmoidal) and has an area of discontinuity (a gap) in the center. I'm curious about what such a result means. My understanding of the half-normal residual plot is that it compares the ordered residuals from the data to the expected values of ordered observations from a normal distribution. Whereas full-normal plots use the signed residuals, the half-normal plots use the absolute values of the residuals. Outliers appear at the top right of the plot as distinct points and departures from a straight line means that the model is not satisfactory. Does this mean that the line with a sigmoid function is simply indicative of a poorly specified, perhaps over-dispersed model? Or is there more to it? Does the line have to be a straight line to indicate that the model is satisfactory (regardless of whether there are outliers, since in this case there are none)?


title: Test data results does not match with cross validation results

text: I'm confused with my data I'm currently playing with. I have a data set which holds 58 attributes in 10000 instances. Attributes are 56 float values typically within 0 to 1. Then there is nominal attribute which tells "data class". Last one is result of data, which is nominal attribute with 3 different values. When playing with Weka with this data set and if i do cross-validation with 1000 folds and try to classify that data with random forest, I get pretty good results, somewhat ~85% of instances are classified correctly. Same applies with Multilayer perceptron (however I haven't done it with that many folds since one fold takes ~5 minutes of time). Problem comes when I test this data with own test set, for example last 10 instances from model, and make a new model with 9990 instances (actually I have plenty of more data, so model is still 10000 instances wide, and no test data is involved there). Now my result is highly worse. ~40% of instances are classified now correctly. How can this be? What I don't understand here and how can I improve this results with my own test data?


title: Pre to Post Repeated measures for categorical dependent variable and more than one independent variables both continues and catgorical

text: I am doing Prenatal (before delivery) and Postnatal (after delivery) depression assessment (2 time points).   Dependent variable is Depressed/Non depressed mother on a measure of depression.   When as there a number of predictor variables both continues and categorical in nature.   Clearly speaking I am using 2 questionnaires at 2 time points: i) dependent variable is dichotomous, ii) dependent variable has continuous predictors. Can i use Generalized estimating equation? What exactly constitutes repeated measures logistic regression?


title: Need to use a percentage as a independent var in regression

text: This is my first question! My firm takes boxes of documents, preps them (removing staples, taping torn documents), and then scans the documents. I want to calculate what effect staples and torn documents have on the time it takes to prep a box of documents. We measured a series of boxes for: the # of total documents, the percentage of staples (# of staples / # of paper), the percentage of torn documents (# of torn docs / # of total docs), and the time it took to prep the box. The issue I'm having is- how do you treat the percentages? (Staples%, TornDoc%) I tried calculating the Staples% * #ofDocuments to get the number of staples, but that results in a high correlation w/ #ofDocuments (multicollinearity). However, if you use it as a percentage, since each box has a completely different number of documents, it seems to me that you first have to convert all % to a common / nominal variable. Please advise!


title: Relation between minimum contrast estimate and minimum distance estimate?

text: What relation are between minimum contrast estimate and minimum distance estimate? If I understand correctly, these two are different methods? or are they equivalent? Thanks and regards! * * * Minimum distance estimate from Wikipedia > Let $\displaystyle X_1,\ldots,X_n$ be an independent and identically > distributed (iid) random sample from a population with distribution > $F(x;\theta)\colon \theta\in\Theta$ and $\Theta\subseteq\mathbb{R}^k (k\geq > 1)$. > > Let $\displaystyle F_n(x)$ be the empirical distribution function based on > the sample. > > Let $\hat{\theta}$ be an estimator for $\displaystyle \theta$. Then > $F(x;\hat{\theta})$ is an estimator for $\displaystyle F(x;\theta)$. > > Let $d[\cdot,\cdot]$ be a functional returning some measure of "distance" > between the two arguments. The functional $\displaystyle d$ is also called > the criterion function. > > If there exists a $\hat{\theta}\in\Theta$ such that > $d[F(x;\hat{\theta}),F_n(x)]=\inf\\{d[F(x;\theta),F_n(x)]; > \theta\in\Theta\\},$ then $\hat{\theta}$ is called the minimum distance > estimate of $\displaystyle \theta$. * * * Minimum contrast estimate from Bickel and Doksum's Mathematical Statistics Vol1 ![enter image description here](http://i.stack.imgur.com/GDin4.png)


title: R-code question: model selection based on individual significance in regression?

text: I'm trying to generate an R function that keeps relevant variables based on their absolute t-value (or p, whichever is easier in code). Basically what I want is to run one regression (1), retain all variables that are significant (based on the t-value, or p). Then run another regression (2), retain all variables that are significant (based on t-value, or p). Then take all retained variables from the regressions (1 & 2) above and combine them in a final regression (3). Currently I am struggling with the step of retaining variables based on their t-value, I can extract the t-values of the variables but don't know how to code that they should be "retained" and kept for the later regression. What I have currently (example edited for better understanding):               summaryI1 <- summary(lm(y~x1+x2+x3+x4)) #x1-x4 are the vars included in regression 1     (coefI1 <- coef(summaryI1))     tI1 <- coefI1[, "t value"] #extracts the t-values for all the x's          cutoff <- which(tI2>cut, arr.in=TRUE) #checks that t-stats of vars are above the cutoff level      Suppose x1 and x2 are significant and I would like to keep them "aside" for the next regression. How would I code that r should create a new matrix of x1 and x2 based on x1 and x2 having t-values higher than my cutoff? I'm not that experienced in r programming yet but learning. Can anyone help?


title: Equivalent of a boxplot for normally distributed data?

text: If I plot some data in function of a categorical variable in R, I get the standard boxplot. However, the boxplot displays non-parametric statistics (quantiles) that don't seem appropriate for normally distributed data. What is the equivalent of the boxplot for normally distributed data (ideally, something based on confidence intervals?) ? How do I get these plots with R?


title: Lognormal distribution from world bank quintiles PPP data

text: I am not a stats person. The world bank has data giving PPP (personal purchasing parity, or something like that) for quintiles (actually first 10%, second 10%, 2nd, 3rd, 4th 20%, 9th 10% and 10th 10%) of a country's population. I've been reading a good deal about the Gini index and also modeling income distribution curves for populations, and I would like to turn this quintiles data into a lognormal distribution representing the same properties (e.g. that if China has 1,300 USD PPP for the 1st 10%, that after modeling it, the weighted average PPP (integration under the curve?) for the poorest 10% of the lognormal distribution would come out to 1,300). Any thoughts on how to do this, tactically? I am cannot think my way through this - but I am a reasonable programmer of simple scripts and would use python's scipy and numpy to fit curves. Given some help on how to proceed. * * * @mpiktas, you are right, they do not give maximum income. The quintiles/decile information must be averages. I do not think you can straightforwardly fit this data as if it is raw data, to a distribution, if that were the case I would not have asked the question! @Michael, given that the quintile/decile data points are average values for the entire quintile/decile bracket, is it possible still to come up with a best fit? Is this what you meant by least squares fitting?


title: How to do factorial analysis for a non-normal and heteroscedastic data?

text: I designed a field experiment with 4 independent factors but data is not normal and heteroscedastic. Friedman test (agricolae package) from R only fits for rbd. Can anybody suggest how to analyze my data please? Many thanks.


title: Fitting a linear model with non gaussian noise

text: I am trying to evaluate the elasticity of prices of some goods. I am concerned about the gaussianity of the noise in the prices. With non gaussianity I am referring to the non existence of the firt/second moment fo the distribution of the error. Is there a way to fit linear model without the assumption of gaussianity of the error?


title: Is there a way to compute confidence intervals for regression estimates of simulated data?

text: I've got a dataset where someone counted birds in the breeding season over 10 years. For each year (x site), we want to see how reduced sampling might affect our ability to detect a trend. So to that end, I have simulated various datasets from the original where we cut down sampling to once every 45 days (and 60, 90,120). I fit the same negative binomial model to the original and each simulated dataset (5oo datasets for each of the different sampling intervals). So now I have 500 regressions (and coefficients) for each dataset x location. Is there some way to throw a confidence band around these? Is it something very trivial, like computing CI around a mean?


title: Interpretation of intercept in linear regression of two-way factorial data without interaction term

text: Suppose I have a two-way factorial design and perform a linear model analysis such as (using _R_ ): `utils::data(npk, package = "MASS")`               lm(yield ~ block*P, npk)       Here, the intercept is the mean of yield of the first level of factors _block_ (with levels 1-6) and _P_ (with levels 0 & 1), i.e. the first field ( _block_ = 1 & _p_ = 0) in a factorial table showing all combinations of factor levels. `mean(npk$yield[npk$P=="0" & npk$block=="1"])` The corresponding t-test appears to be the test against the null that this mean is different from zero. Now suppose I don't model an interaction term:                  lm(yield ~ block+P, npk)      Here, I don't understand where the intercept comes from. I thought it could be:               (mean(npk$yield[npk$block=="1"]) + mean(npk$yield[npk$P=="0"]))/2      i.e. in a factorial table, the mean of the first row and the first column, but this does not appear to be the case. How is the intercept calculated in a linear model of a two way factorial design without an interaction term? How is interpreted and what null is tested with the t-test?


title: Empirical relationship between mean, median and mode

text: For a unimodal distribution that is moderately skewed, we have the following empirical relationship between the mean, median and mode: $$ \text{(Mean - Mode)}\sim \text{3(Mean - Median)} $$ How was this relationship derived? Did Karl Pearson plot thousands of these relationships before forming this conclusion, or is there a logical line of reasoning behind this relationship?


title: How does knnimpute of the preprocess function work?

text: I am new to R and I use a script I do not completely understand. It preprocesses a dataset for data mining. At one point, the data (stored in `fil`) should be centered and scaled and the NA values are replaced by values calculated with `knnImpute`.               preProcValues <- preProcess(fil, method = c("center", "scale","knnImpute"))     filpred <- predict(preProcValues, fil)      I get the following Error:               Error in nn2(old[, cols, drop = FALSE], new[, cols, drop = FALSE], k = k) :     NA/NaN/Inf in foreign function call (arg 1)      It seems to me that this error has to do with the way how `knnImpute` works. Is it possible that `knnImpute` uses surrounding columns to impute NA values? Because in the column where this stops to work there are a lot of NA values in one row of several columns. Can anybody explain if this might be the cause of the error? Is anybody able to explain how `knnImpute` works or what else might be the cause?


title: How to perform multiple regression analysis with a small sample size? (n=20)

text: I have a sample of 20 years, I want to perform multiple regression analysis to estimate economic growth. I have 9 independent variables, but because my sample is small i get insignificant results. I do not have data to maken it a sample of 100 years or anything near that. What can I do? Can I incorporate the sample size into my problem? How do I enter a sample size in SPSS?


title: How can the IID assumption be checked in a given dataset?

text: **1-** How can I check if a set of data can be assumed as IID data? I'm not so familiar with statistics, but I guess I should look at the first lag of autocorrelation for independent distribution. Have no idea about identical distribution condition! **2-** It seems that I was not clear enough! I'm trying to detect outliers in a series of records (turbulent flow velocity in a river). I transform data into wavelet space and then I shrink the wavelets over a certain threshold. Since standard deviation is the worst option as an scale estimator, I looked for a new estimator. Rousseeuw and Croux developed new robust estimators for measuring dispersion in iid random variables, Sn and Qn. I don't know offhand if the high breakdown properties they enjoy carry over to the time-series case or not. From the answer given by kwak, I can infer that wavelets do NOT follow independent distribution property. Since after shrinkage, location of non-zero elements indicates the spike location in the original time series. Am I true? (shuffling the indices results in losing location of spikes) If so, other scale estimators like median absolute deviation (MAD) are not valid in case of time series as we calculate the median. How about identical distribution assumption requirements? **3-** OK, let me ask my question in simple manner: I want to use robust scale estimators Sn and Qn for shrinking a series of wavelets. the wavelets are obtained from decomposing observations of a turbulent flow field velocity vectors collected at 1 Hz sampling rate. if the data can be assumed as iid e.g. Qn has breakpoint of 50% and efficiency of 82% (Gaussian distribution). My question is **the high breakdown properties they enjoy carry over to the time-series case or not.** Or how can i approve that the wavelets follow iid characteristics.


title: The Sleeping Beauty Paradox

text: ### The situation Some researchers would like to put you to sleep. Depending on the secret toss of a fair coin, they will briefly awaken you either once (Heads) or twice (Tails). After each waking, they will put you back to sleep with a drug that makes you forget that awakening. When you are awakened, to what degree should _you_ believe that the outcome of the coin toss was Heads? _(OK, maybe you don’t want to be the subject of this experiment! Suppose instead that Sleeping Beauty (SB) agrees to it (with the full approval of the Magic Kingdom’s Institutional Review Board, of course). She’s about to go to sleep for one hundred years, so what are one or two more days, anyway?)_ ![Maxfield Parrish illustration](http://i.stack.imgur.com/zLmrR.png) _[Detail of aMaxfield Parrish illustration.]_ ### Are you a Halfer or a Thirder? **The Halfer position.** Simple! The coin is fair--and SB knows it--so she should believe there's a one-half chance of heads. **The Thirder position.** Were this experiment to be repeated many times, then the coin will be heads only one third of the time SB is awakened. Her probability for heads will be one third. ### Thirders have a problem Most, but not all, people who have written about this are thirders. But:   * On Sunday evening, just before SB falls asleep, she must believe the chance of heads is one-half: that’s what it means to be a fair coin.   * Whenever SB awakens, _she has learned absolutely nothing she did not know Sunday night._ What rational argument can she give, then, for stating that her belief in heads is now one-third and not one-half? ### Some attempted explanations   * SB would necessarily lose money if she were to bet on heads with any odds other than 1/3. (Vineberg, _inter alios_ )   * One-half really is correct: just use the Everettian “many-worlds” interpretation of Quantum Mechanics! (Lewis).   * SB updates her belief based on self-perception of her “temporal location” in the world. (Elga, _i.a._ )   * SB is confused: “[It] seems more plausible to say that her epistemic state upon waking up should not include a definite degree of belief in heads. … The real issue is how one deals with known, unavoidable, cognitive malfunction.” [Arntzenius] * * * ### The question Accounting for what has already been written on this subject (see the references as well as a previous post), how can this paradox be resolved in a statistically rigorous way? Is this even possible? * * * ### References Arntzenius, Frank (2002). _Reflections on Sleeping Beauty_ Analysis 62.1 pp 53-62. Bradley, DJ (2010). _Confirmation in a Branching World: The Everett Interpretation and Sleeping Beauty_. Brit. J. Phil. Sci. 0 (2010), 1–21. Elga, Adam (2000). Self-locating belief and the Sleeping Beauty Problem. Analysis 60 pp 143-7. Franceschi, Paul (2005). _Sleeping Beauty and the Problem of World Reduction_. Preprint. Groisman, Berry (2007). _The end of Sleeping Beauty’s nightmare_. Preprint. Lewis, D (2001). _Sleeping Beauty: reply to Elga_. Analysis 61.3 pp 171-6. Papineau, David and Victor Dura-Vila (2008). _A Thirder and an Everettian: a reply to Lewis’s ‘Quantum Sleeping Beauty’_. Pust, Joel (2008). _Horgan on Sleeping Beauty_. Synthese 160 pp 97-101. Vineberg, Susan (undated, perhaps 2003). _Beauty’s Cautionary Tale_.


title: Can AIC compare across different types of model?

text: I'm using AIC (Akaike's Information Criterion) to compare non-linear models in R. Is it valid to compare the AICs of different types of model? Specifically, I'm comparing a model fitted by glm versus a model with a random effect term fitted by glmer (lme4). If not, is there a way such a comparison can be done? Or is the idea completely invalid?


title: Estimating parameters of a dynamic linear model

text: I want to implement (in R) the following very simple Dynamic Linear Model for which I have 2 **unknown time varying parameters** (the variance of the observation error $\epsilon^1_t$ and the variance of the state error $\epsilon^2_t$). $ \begin{matrix} Y_t & = & \theta_t + \epsilon^1_t\\\ \theta_{t+1} & = & \theta_{t}+\epsilon^2_t \end{matrix} $ I want to estimate these parameters for each time, **without any look ahead bias**. From what I understand, I can use either a MCMC (on a rolling window to avoid the look ahead bias), or a particule filter (or Sequencial Monte Carlo - SMC). **Which method would you use** , and   **What are the pros and cons of these two methods?** Bonus question: In these methods, how do you select the speed of change of the parameters? I guess we have to input an information here, because there is a bargain between using a lot of data to estimate the parameters and using less data to react more quickly to a change in the parameter?


title: Sampling Distribution

text: Suppose $X \sim N(\mu, \sigma^2)$. After $n$ independent draws $\bar{X} \sim N(\mu, \frac{\sigma^{2}}{n})$. I am assuming that each draw contains $n$ numbers? Because if each draw only had one number then we would only have one mean.


title: How to calculate 95% CI for OR for a different reference category without running the SAS logistics again?

text: My question is about calculation of confidence interval (CI) for odds ratio (OR) from a SAS output of a logistic regression model for a different reference category without running the SAS program again. I’ve run SAS logistics modeling dichotomized income level (low coded as 1 versus high coded as 0) with 4 dummy coded race variables as covariates- White, Chinese, Black and Filipinos with White as the reference category. While looking at the OR and CI for Filipinos’ income level against White from a SAS output that only furnished intercept, beta, SE, OR and 95% Wald CI for OR, I’ve become interested to know about the OR and 95% CI for Filipinos against Chinese (instead of White). I know how to calculate the new Beta and OR by calculating the predicted values for Chinese and Filipinos and then taking the difference of the log pie for Chinese and Filipinos to get the new beta and OR, but I’m completely at a loss about calculation of the new CI. Can anyone HELP please?


title: Sample size for t-test when the outcome is percentage point

text: first time asking question here so I hope I will not mess up. I am in the process of calculating the sample size for a study that utilizes an independent _t_ -test design: comparing some outcome before and after a certain policy, subjects in the pre and post periods are not linked. Simply put, the study is about infant weaning food distributed in developing countries. The outcome is percent weight of cooking ingredients that is oil, which is fortified with nutrients. The current recommendation is 9% weight of the baby's meal should be the fortified oil (about 1 part of oil to 10 part of powered grains). The organization is trying to promote a new amount of 16% oil, which was shown to be associated with better anthropometric outcomes. I could use the usual _t_ -test sample size formula, but I'm concerned that the SD of the mean 9% point may not be meaningful because the outcome is bound between 0% to 100%, the distribution of the percent point at the extremes (like 9%) is likely to be non-normal, and additionally not all caretakers follow the cooking guideline so the SD is probably huge. So my question is: do you know of any formula is applicable to this situation of comparing mean when the **individual** (not group) level data are proportion? Many thanks for reading. I appreciate any advice or reference.


title: Stochastic coordinate descent for $\ell_1$ regularization

text: I recently came across the following paper: "Stochastic Methods for $\ell_1$ Regularized Loss Minimization" by Shai Shalev-Shwartz and Ambuj Tewari, ICML 2009. In the paper, the authors propose a modification of the coordinate descent algorithm for the LASSO in which the coordinates (the $\beta$s) are updated in a random order. This modification appears to have better run-time performance than deterministic coordinate descent. Can you offer any intuition for why such a modification would make the algorithm faster in practice?


title: classification algorithms that return confidences?

text: Given a machine learning model built on top of scikit-learn, how can I classify new instances but then choose only those with the highest confidence? How do we define confidence in machine learning and how to generate it (if not generated automatically by scikit-learn)? What should I change in this approach if I had more that 2 potential classes? This is what I have done so far:               # load libraries     from sklearn import neighbors     # initialize NearestNeighbor classifier     knn = neighbors.KNeighborsClassifier(n_neighbors=3)     # train model     knn.fit([[1],[2],[3],[4],[5],[6]], [0,0,0,1,1,1])     # predict ::: get class probabilities     print(knn.predict_proba(1.5))     print(knn.predict_proba(37))     print(knn.predict_proba(3.5))


title: Is it acceptable to not include high-order interactions (3-way and above) in the model when they are not by themselves of interest?

text: Is it acceptable to not include high-order interactions (3-way and above) in the model when they are not of interest and not part of the hypothesis that is being tested? NB. I am not talking about model reduction when a full model was fitted and higher interactions later excluded. I am asking about not including them at the outset, because their effects are not of interest and often uninterpretable. For example, in a situation where the response is modelled as a function of one experimental factor and several covariables and the data is balanced only with respect to the experimental factor. I heard somebody argue that it is acceptable, and if so, is there a reference to a credible source? As far as I understand, this is somewhat similar to ANCOVA, where the covariate is included as a main effect term, but not part of an interation.


title: Pizza statistics for the masses

text: A short entry on the NY Times website provides the Facts and Figures of pizza consumption in the United States. I have a casual interest in how statistics are used (or abused) to provide information to general audiences, and a couple questions have arisen based on the statistics presented:   * If 1 out of 8 Americans will eat pizza today, does that mean that the average American will eat pizza once every 8 days? There's an assumption here that _every_ american eats pizza, which is not the case; however, that raises the question of how to make a valid assumption of how many Americans eat pizza.   * It is reported that 25% of a child's caloric intake is pizza. I'll define a child as a 9 year old who is moderately active and therefore needs a 2000 calorie daily intake. If we trust Google's estimate that the number of calories in a pizza slice is 285, then does that suggest that a child consumes 12 slices of pizza per week on average? (2000 * 7 * 0.25 / 285) I suspect my interpretation of the statistics is flawed; it does not seem to me that a child could be part of the 1-out-of-8 Americans eating pizza today while also eating approximately 1.7 slices per day to achieve the 25% caloric intake number.


title: Bootstrap to evaluate variance of AUC ROC

text: I have a toy dataset and want to eval the AUC - ROC with bootstrap                   set.seed(117)         test <- c(rnorm(70, 40, 10), rnorm(30, 55, 10))         dx <- c(rep(0, 65), rep(1,35) )         df <- data.frame(cbind(test, dx))          r <- 1000         boot.f <- function(d, i){           data <- d[i,]           p <- unlist(subset(data, data$dx==1)["test"])           n <- unlist(subset(data, data$dx==0)["test"])           mean(sample(p, r, replace=T) > sample(n, r, replace=T))         }          roc.boot<-boot(df, boot.f, r)          roc.boot          boot.ci(roc.boot, type="bca")      It seems that bootstrapping needs two samplings (in `mean` and in `boot`). Is there a more efficient way? I'd like to obtain a `boot` object, so prefer to use the function.


title: R lme4 1.1-7: REML=FALSE giving error "extra argument(s) ‘REML’ disregarded "

text: I'm trying to fit a model with the function `glmer` (`lmer4` 1.1-7 package) in R using REML but I just get an error saying `extra argument(s) ‘REML’ disregarded` (see below). How else can I fit it with REML to optimize the random effects structure with this package?               a1<-glmer(Infected~(1|River),family=binomial,data=covertfrag,REML=FALSE)     Warning message:     extra argument(s) ‘REML’ disregarded


title: Rejection sampling - picking a g(x) distribution

text: Let's say I have to sample from a pdf $\pi(x) = 3x^3+\frac{3}{4}x^2, 0 \le x \le 1$. I know we have to pick some pdf $g(x)$ such that $cg(x) \ge \pi(x)$ for all x. The only thing that came to mind was using $g(x) = 1$ with $c = 3.75$. So $cg(x) = 3.75$. No better distributions come to mind that are easy to sample from. Is this a good option? Also two conceptual things - what makes a distribution easy to sample from? Why do we want to pick c as small as possible. Why not use $c = 50$? Is this for algorithm efficiency? i.e. more of the randomly drawn samples from $cg(x)$ will fall under the curve of $\pi(x)$?


title: Why the Breusch-Pagan rejects H0 on apparently non-heteroskedastic data?

text: Breusch-Pagan rejects the H0 on this residuals:               > length(model$residuals)         [1] 515959         > summary(model$residuals)         Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      -205.000   -4.420   -0.451    0.000    4.130  196.000      > quantile(model$residuals, seq(0, 1, 1/10))              0%         10%         20%         30%         40%         50%         60%         70%         80%         90%        100%      -205.228344  -10.137559   -5.705923   -3.400891   -1.814776   -0.450546    1.036849    2.914375    5.648389   10.733350  196.011350      > bptest(model)              studentized Breusch-Pagan test          data:  model     BP = 1385.57, df = 14, p-value < 0.00000000000000022      Looking at the residuals variance with @whuber visualization function I don't see any heteroskedasticity issue: ![ResidualsvsPredicted](http://i.stack.imgur.com/wvYlh.png) So yes, there are some extreme observations at the tail of distribution, but is this enough to reject the H0? **Manual Breusch-Pagan test**               summary(lm(I(model$residuals^2) ~ data$X1 + data$X2 + data$data$X3 + data$X4))          Call:     lm(formula = I(model$residuals^2) ~ data$X1 + data$X2 +              data$dataX3 + data$X4)          Residuals:             Min     1Q Median     3Q    Max        -356   -100    -79    -23  41965           Coefficients:                         Estimate Std. Error t value             Pr(>|t|)         (Intercept)           176.52       4.82   36.66 < 0.0000000000000002 ***     data$X1                -5.21       0.28  -18.60 < 0.0000000000000002 ***     data$X2               -75.36       4.75  -15.88 < 0.0000000000000002 ***     data$X3                15.70       0.85   18.48 < 0.0000000000000002 ***     data$X4fact2           10.84       3.15    3.44              0.00059 ***     data$X4fact3            1.14       2.15    0.53              0.59735         data$X4fact4           96.31       6.98   13.80 < 0.0000000000000002 ***     data$X4fact5          191.39      12.73   15.04 < 0.0000000000000002 ***     data$X4fact6            3.44       4.03    0.85              0.39280         data$X4fact7            8.30       2.79    2.97              0.00298 **      data$X4fact8           19.38       2.66    7.29     0.00000000000031 ***     data$X4fact9            1.44       2.37    0.61              0.54468         data$X4fact10          46.35      11.19    4.14     0.00003435124602 ***     data$X4fact11          -6.23       2.71   -2.30              0.02168 *       data$X4fact12          18.07       3.46    5.22     0.00000017455269 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 502 on 515944 degrees of freedom     Multiple R-squared:  0.00269,   Adjusted R-squared:  0.00266      F-statistic: 99.2 on 14 and 515944 DF,  p-value: <0.0000000000000002


title: Do you use a chi-squared test or a t-test for equality of variances?

text: How can I find if I should t-test or chi-squared test if I am given a problem like the following? Consider testing $H_0: \sigma^2_X = \sigma^2_Y$ against $H_1: \sigma^2_X ≠ \sigma^2_Y$ from two independent samples from normal populations with unknown means $\mu_X$ and $\mu_Y$ and standard deviations $\sigma_X$ and $\sigma_Y$. The $X$'s are 11.4, 9.7, 11.4, 13.3, 7.4, 8.5, 13.4, 17.4, 12.7. The $Y$'s are 3.2, 2.7, 5.5, -0.9, -1.8. Find the value of the test statistic. P.S.: I know how to do the `chisq.test` and `t.test` when I just one hypothesis ($H_0$)! How should I write R script to do the above problem when I have more than one hypothesis? What are some good external R related script to this question that I can cover for seeing similar example?               > X = c( 11.4, 9.7, 11.4, 13.3, 7.4, 8.5, 13.4, 17.4, 12.7)     > Y = c(3.2, 2.7, 5.5, -0.9, -1.8)     > ?t.test     > t.test(X, Y)              Welch Two Sample t-test          data:  X and Y     t = 5.9114, df = 8.306, p-value = 0.0003089     alternative hypothesis: true difference in means is not equal to 0     95 percent confidence interval:       6.092637 13.805141     sample estimates:     mean of x mean of y       11.68889   1.74000           > chisq.test(X, Y)     Error in chisq.test(X, Y) : 'x' and 'y' must have the same length


title: Multivariate time series simulation in RapidMiner?

text: I'm actually "Getting started" with RapidMiner (RM). I'm an R expert but totally newbie to RM. The problem involved is analysing and forecasting the dispersion of multivariate timeseries (in finance). The R codes we have are very difficult for economists. I found RM as a real alternative because of its visual "self-documenting" and simplcity in ETL tasks. To be honest, I fitted some models to my data as an experiment , it was also very expressive. What are the experiences of others? Is it possible to **SIMULATE** a Series model in RM (or in other visual and simple tool)? The Looping operators were not enaugh to do this because (for example in a VAR1 setup) this kind of model simulation requires peturbating the error terms with some noise that affects the following predictions. To be extreme it would be neccesary to modify the VAR1 parameters (to modfy the base cases, e.g. the point prediction's trajectory). It’s possible in R, but not too intiutive, so It would be great implement a this kind of model in RM. (The problem in short: /multivariate historical data/ -> /few factors/ -> /multivariate timeseries model fitted to these factors with indogeneous and exogeneous variables and lags/ -> /simulated factors/ -> /simulated data/ -> /transform and visualize dispersion of data/) I know that RM is not well suited for this task by design, it is designed for others. Any case: I'm sure that we could use it much more than an (on site) ETL tool. Thanks, regards


title: Time series modeling in R on a weekly basis over multiple years featuring different number of weeks in each year

text: My problem is that I have weekly data from 2008, 2009 and 2010 but the number of weeks in each of these years are not the same since in we (in Denmark) code the last week in 2009 as 53. This means that I have 52 weeks in 2008 and 2010 but 53 weeks in 2009. Is there some way to make R understand this and treat the time series accordingly. Typically I want it to integrate well with ts, stl, etc. Any help or pointers on this subject would be greatly appreciated.


title: EM algorithm - the values of likelihood and complete likelihood for different starting values

text: I use finite mixture models. For two sets $A$ and $B$ of starting values (the initial guess) I get the following values of $Q(\bf{\hat{\Theta}})$ function and $\log L(\bf{\hat{\Theta}})$: Set $A$: $Q(\bf{\hat{\Theta}_A})=110$ and $\log L(\bf{\hat{\Theta}_A})=7$ Set $b$: $Q({\bf{\hat\Theta}_B})=78$ and $\log L(\bf{\hat{\Theta}_B})=20$ What makes me confused is that for these local probability maxima I should choose $\bf{\hat{\Theta}_A}$ when referring to $Q$ function and conversely choose $\bf{\hat{\Theta}_B}$ if I decided to rely on $\log L$. McLachlan on page 78 (The EM Algorithm and Extensions - McLachlan, Krishnan 2ed. (2008)) allows me to derive that if $\bf{\hat{\Theta}_B}$ doesn’t globally maximizes likelihood function, there exists $\bf{\hat{\Theta}_A}$ that must satisfy: $Q(\bf{\hat{\Theta}_B}) < Q(\bf{\hat{\Theta}_A})$ (which is true for my case) implying that $\log L(\bf{\hat{\Theta}_B}) < \log L(\bf{\hat{\Theta}_A})$. The last inequality doesn't hold for the above values. It means that based on theoretical grounds that something is wrong with my implemented algorithm. Do you think that above derivations are appropriate?


title: Sufficiency in Lehmann Scheffe

text: We are wondering what sufficiency in the Lehmann Scheffe Theorem is needed for. Our reasoning was:   1. If an unbiased estimator is uncorrelated with all unbiased estimators of 0, it is UMVUE   2. If the estimator is from a complete family, it is uncorrelated with all unbiased estimators of 0   3. Therefore an unbiased estimator from a complete family is UMVUE   4. UMVUE is unique So what exactly is sufficiency needed for? Why is completeness not enough? Would be nice if anyone could point out the flaws in our reasoning and explain!


title: Phillips–Perron unit root test instead of ADF test?

text: I have read that PP unit root is often used in economy. Is it sensible to do both tests (PP & ADF) or is PP test enough?


title: Examples of wrapping open source machine learning software in PMML?

text: Does anyone know of examples of wrapping an existing piece of supervised learning software to output models in PMML format? Of particular interest are learners that just take in labeled vectors of numbers as training data and put out models that are pretty much just coefficient vectors (liblinear, SVMlight, BXRtrain, BOW, etc.). That is, they don't have any smarts about data types, ranges of legal values of features, etc.: something else is assumed to deal with that and present the learner with appropriate numeric vectors. For such software, all the interesting data dictionary stuff would need to be supplied alongside the input data if it's going to show up in PMML model that's output. There's nothing conceptually difficult about this: what I'm curious to see is if there's any conventions, design patterns, etc. that have grown up in the PMML community for doing this. The PMML website list of software that either consumes or produces models in PMML, but this is mostly commercial closed source software. The programs with open source versions listed there (Rapidminer and WEKA that I can spot) are rather complex data mining suites. What I'd like to see an example of is a a minimalist wrapping of a simple one-trick pony kind of learner.


title: Confidence Interval Application to a time series

text: So I'm not sure if I 100% grasp confidence intervals. Say I have a huge data set of a bond prices from 1996 to present in MINUTES. Suppose I separate each data by day. If I were to use a Dickey Fuller Test (or any other t-test), on a particular day (thus a sample is drawn based on all points in one day) and the value I got rejected the null hypothesis, would it be conceptually right for me to conclude that because based on my sample I was able to reject the null, I am able to make a conclusion about the whole dataset (population)? In other words, would picking one day and not randomly picking datapoints throughout matter?


title: Adjusting for covariates in ROC curve analysis

text: This question is about estimating cut-off scores on a multi-dimensional screening questionnaire to predict a binary endpoint, in the presence of correlated scales. I was asked about the interest of controlling for associated subscores when devising cut-off scores on each dimension of a measurement scale (personality traits) which might be used for alcoholism screening. That is, in this particular case, the person was not interested in adjusting on external covariates (predictors) -- which leads to (partial) area under covariate- adjusted ROC curve, e.g. (1-2) -- but essentially on other scores from the same questionnaire because they correlate one to each other (e.g. "impulsivity" with "sensation seeking"). It amounts to build an GLM which includes on the left-side the score of interest (for which we seek a cut-off) and another score computed from the same questionnaire, while on the right- hand side the outcome may be drinking status. To clarify (per @robin request), suppose we have $j=4$ scores, say $x_j$ (e.g., anxiety, impulsivity, neuroticism, sensation seeking), and we want to find a cut-off value $t_j$ (i.e. "positive case" if $x_j>t_j$, "negative case" otherwise) for each of them. We usually adjust for other risk factors like gender or age when devising such cut-off (using ROC curve analysis). Now, what about adjusting impulsivity (IMP) on gender, age, and sensation seeking (SS) since SS is known to correlate with IMP? In other words, we would have a cut- off value for IMP where effect of age, gender and anxiety level are removed. Apart from saying that a cut-off must remain as simple as possible, my response was > About covariates, I would recommend estimating the AUCs with and without > adjustment, just to see if the predictive performance increase. Here, your > covariates are merely other subscores defined from the same measurement > instrument and I never faced such a situation (usually, I adjust on known > risk factors, like Age or Gender). [...] Also, since you are interested in > prognostic issues (i.e. screening efficacy of the questionnaire), you may > also be interested in estimating the positive predictive value (PPV, > probability of patients with positive test results who are correctly > classified) provided you are able to classify subjects as "positive" or > "negative" depending on their subscores on your questionnaire. Note, > however, that it is necessary to know the prevalence of this disorder to > correctly interpret the PPV in turn... Do you have a more thorough understanding of this particular situation, with link to relevant papers when possible? **References**   1. Janes, H and Pepe, MS (2008). Adjusting for Covariates in Studies of Diagnostic, Screening, or Prognostic Markers: An Old Concept in a New Setting. _American Journal of Epidemiology_ , 168(1): 89-97.   2. Janes, H and Pepe, MS (2008). Accommodating Covariates in ROC Analysis. _UW Biostatistics Working Paper Series_ , Paper 322.


title: Nested ANOVA in SPSS / Science Education

text: I am analyzing student performance on a pre/post test designed to assess gains for a week long high school biology curriculum on natural selection. After comparing pre/post performance using a variety of different metrics (test score, Rasch measure, chi-square per item) between the paired samples, I conducted an ANOVA to gauge whether there were differences between teachers. Two of the (seven) teachers are significantly different from the others (their students have drastically lower mean Rash measures and did worse overall), so I'd like to account for these differences in the analysis. Is a nested ANOVA model appropriate here? I know nested models aren't in the SPSS dialogs, but I'm not much of a programmer and thus do not know what syntax to produce. Moreover, each teacher has a different number of students. Would a random sample (ensuring the same amount from each teacher) be useful for paired comparisons of means (like pre vs. post), or does it not matter since it's a comparison of means between groups anyway? Apologies for any obvious questions. Education is not my original field (I'm a bio-anthropologist), and it's been a bit since I took multivariate stats in grad school.


title: How to get rid of bias in data?

text: I have been trying to classify a set of data into one of four classes. The data has already been generated and I have set aside 10,000 for training and 2,000 for testing. I have also generated the labels for each of the data. Let's call the classes - 0,1,2 and 3. Now when I observe the classification, I notice that there are a lot of 0s in the training data and hence in most cases, the classifier is just learning to predict 0 no matter what the features are. (I am using random forests for classification) Generating the data again to ensure uniformity, takes a lot of time and I prefer to avoid that. Is there anyway I can still use the data that I have? Thanks in advance


title: Predicting the missing data out of three values in each of the two vectors

text: I have 2 vectors of rural and urban populations of the same country. (years from 1975 to 2020) with only three values (1980, 1990 and 2001 years) in each. And I need to predict the missing data. My first approach was to interpolate the data in between the known years(using na.approx in R) and then use make an ARIMA model, however I am not sure if its a good idea, since the graphics look kind of awkward. Also the values are "waving", for example, urban population - 4000 in 1980, 3800 in 1990 and 3950 in 2001. Any ideas on how to make a better prediction would be really helpful.


title: Statistical Power of ROC/AUC Test with non-IID Samples :: To how many IID Samples are my non-IID Samples Equivalent?

text: I've been assigned to solve the following problem as part of a serious, biological research project. I think I have a tentative solution, but I'm wondering whether the approach I've picked is the best. Your help/input would be appreciated. We monitored the vital signs (V) of ten subjects as they travelled up and down a mountain. From each subject, we took approximately 100 samples, one every minute --- fifty when they were at the top, and fifty at the bottom of the mountain. A sample contained a measurement of five vital signs of the subject (respiration rate, heart rate, etc.). Along with each sample of vital signs, we also recorded whether the patient was on the top (T) or on the bottom (B) of the mountain. For each sample, we computed Stress (S), where S is a function the vital signs, so S = f(V). (The computational details of f are not relevant.) 1000 samples were drawn from 10 subjects, whereby we had 1000 values for Stress, each of which was coupled to a piece of binary information (B or T). The purpose of collecting the binary information is to test our (previously written) _binary classifier_ : based on a Stress-value, the classifier decides whether or not the subject is on the top or bottom of the mountain. By comparing this to the actual binary value (top / bottom), we can populate a confusion matrix. We then iteratively vary the discrimination threshold of the binary classifier, to finally construct an ROC curve. This gives us the Area Under the Curve (AUC), which we use as to assess the quality of our classifier. That is the end result we are interested in. The important question has been apparent from the start: what's the statistical power of our end result (the AUC), given that we're really only sampling from 10 subjects? More concretely, our samples are obviously not IID. We're asking ourselves: approximately how many IID samples are our 1000 samples equivalent to? We have been thinking about the power of our test with 1000 non-IID samples as equivalent to a test with $n$ IID samples and the same power. Is there any statistical methodology that will let us find $n$? This is the general direction of the question I'm trying to answer. My tentative solution involves the use of General Estimating Equations. A similar approach seems to have been taken by the authors of the CLEPSYDRA study (available here). My reasoning is that the covariance structure of the samples is essentially unknown, which is a problem that GEEs would let me overcome. But I'm not sure if there won't be further complications, i.e. if this is really the right tool for the task, or if there isn't a simpler approach. Is there a simpler approach? Are there further issues that using a GEE won't solve?


title: What would be an effective strategy to correlate congressional attendance and precipitation?

text: I have data sets of the attendance of congresspeople, and of the weather in Washington, DC. I would like to examine whether and to what degree precipitation is correlated with the likelihood of individual congresspeople being present or absent. For example, the 111th congress was in session 298 days of the 388 days between 2009-01-06 and 2010-01-29. Precipitation was recorded during 187 of those 388 days. (I also have the amount of precipitation, the temperature, the pressure, lots of weather measurements) Congresswoman A was absent 52 times; of those 52, there was precipitation 39 times. Et cetera for each congressperson. The null hypothesis would be that there is no relationship between precipitation and attendance. I'd like to create some sort of visualization; perhaps a histogram of number of congresspeople on the y axis and... er... something Bayesian? Chances they will be absent, given the fact it's raining?


title: Get fitted values estimated in ARIMA in Matlab

text: I am using Matlab to forecast time series data using ARIMA algorithm. I am able to get forecasted values, but unable to get the fitted values. This is what I mean. In R:                model <- auto.arima(data)");      forecasted <- forecast.Arima(model, h=36)      forecasted$fitted      `forecasted$fitted` will give the data which has been fitted to the original data by the ARIMA algorithm. But I cannot see anything like that in Matlab. This is what I have done so far:               model = arima(1,1,0);     [est, ~, logL, info] = estimate(model, data);     [output, YMSE] = forecast(est, 10);      output gives me the 10 values Arima has forecasted, but not the original data values it has fitted. I would like to know those. Any help is appreciated.


title: How to cope with exploratory data analysis and data dredging in small-sample studies?

text: Exploratory data analysis (EDA) often leads to explore other "tracks" that do not necessarily belong to the initial set of hypotheses. I face such a situation in the case of studies with a limited sample size and a lot of data gathered through different questionnaires (socio-demographics data, neuropsychological or medical scales -- e.g., mental or physical functioning, depression/anxiety level, symptoms checklist). It happens that EDA helps to highlight some unexpected relationships ("unexpected" meaning that they were not included in the initial analysis plan) that translates into additional questions/hypothesis. As is the case for overfitting, data dredging or snooping does lead to results that do not generalize. However, when a lot of data is available, it is quite difficult (for the researcher or physician) to postulate a limited set of hypotheses. I would like to know if there are well-acknowledged methods, recommendations, or rules of thumb that may help to delineate EDA in the case of small-sample studies.


title: Linear mixed effect models with two independent variables

text: I am estimating a random intercept and a random slope model using the following R code. My dependent and independent variable are both continuous.               randominterceptfixedslope<-lmer(y ~ x + (1|state),data=data,method="ML") # model with fixed slope but random intercept     randominterceptrandomslope<-lmer(y ~ x + (1+x|state),data=data,method="ML") # model with random slope and random intercept          anova(randominterceptfixedslope,randominterceptrandomslope)      Anova tells me that my randominterceptrandomslope model is a better fit on the data. So far good, please correct me if I am wrong. My question is: If I have another independent variable $x_1$, can I put two independent variables in the above model i.e. can I have a randominterceptfixedslope and randominterceptrandomslope model with two independent variables. If yes, how do I do that? As in what the code should look like? Thanks for your response. I got a second query. lets say my full model is this:                    randominterceptrandomslope<-lmer(y ~ x1 + x2 + x3 + x4 (1+x1+x2+x3+x4|state),data=data,method="ML")      If some of my independent variables are correlated, what is the procedure of reducing the collinearity issue in a linear mixed effect model? I could spot collinerity using VIF and retain the most significant independent variables but I can do this for each factor level (levels of state) individually. But won't it result in retaining some independent variables in one factor level while deleting the same in other factor level? I guess the main question is how to spot collinearity in a mixed effect model and what to do with it when you have 5 or 6 independent variables?


title: How to interpret a hazard ratio from a continuous variable -- unit of difference?

text: I am reading an article which shows Hazard Ratios for continuous variables, but I'm not sure how to interpret the given values. My current understanding of hazard ratios is that the number represents the relative likelihood of [event] given some condition. E.g: if the hazard ratio for death from lung cancer given smoking (a binary event) is 2, then smokers were twice as likely to die in the monitored time period than non-smokers. Looking on wikipedia, the interpretation for continuous variables is that the hazard ratio applies to a unit of difference. This makes sense to me for ordinal variables (e.g number of cigarettes smoked a day), but I don't know how to apply this concept to continuous variables (e.g. grams of nicotine smoked a day?)


title: What is the likelihood of an observed standard deviation given a known normal distribution?

text: Given a sample X of size n from a normal distribution $N(\mu,\sigma)$ one can estimate $\sigma$ by $\hat{\sigma}$ from X. Then we know that: $p(\sigma|\hat{\sigma},n)\propto\hat{\sigma}\sqrt{\frac{n-1}{\chi^2(\nu=n-1)}}$ Now I am looking for the likelihood function: $p(\hat{\sigma}|\sigma,n)$. Is there an analytic solution for this or is it necessary to derive the likelihood function from case-by-case simulations? Can someone help me with this relation? Or better, point to a reference where this relation has been used before?


title: Detect equation from data

text: I know R is a very powerful tool for data mining and predictive model building. However, I'm finding it difficult to extract an equation from any sort of modeling that I do. For example, using something like:               model <- lm(y~var1+var2+var3)      For data mining, I can use something like:               analysis <- rpart(y~var1+var2+var3)      to understand my data. However, it would be nice if R had the capability of displaying the model it uses, whether it is linear or something much more complicated. There is a program called Eureka from Nutonian that is able to take data and try and fit an equation to it. Is R able to do this?


title: Statistical test for spatial difference of incidents

text: I have an area divided into grids. I'm looking for a method to determine statistically, if the occurence of criminal incidents for two time periods is significant equal. I have two datasets with equal periods and grid locations. The target values are only assigned to a specific time period. TARGET = 0 -> No incident, TARGET = 1 -> Incident   **Timeperiod 08AM-16PM**   Day, XCOR, YCOR, TARGET   1, 1, 1, 0   1, 1, 2, 0   2, 1, 1, 0   2, 1, 2, 0 **Timeperiod 16PM-24PM**   Day, XCOR, YCOR, TARGET   1, 1, 1, 1   1, 1, 2, 0   2, 1, 1, 0   2, 1, 2, 1   I'm looking for a statistical method to validate/test if the occurence of an incident for both time periods is spatial different? But, how can I do so? Any tips, links or advice much appreciated


title: Multiplicative regression coefficient

text: The expected model is: $$ Y = X1*X2*X3 $$ Taking logs we should get a coefficient of $1$ on all three when we regress. That is: $$ \ln (Y) = b + b1 \ln(X1) + b2 \ln(X2) + b3 \ln(X3), \\\ \text{gives: }b1 = b2 = b3 = 1. $$ What could be the reasons when we collect the data for all the variables and run the regression but get coefficients not equal to $1$? The only information more than what is already mentioned is that the measure of Y comes from a separate dataset than the other variables. There is no measurement error and this is observational data. I guess this would be more of a theoretical question.


title: PDF of the normal distribution and probability values

text: I have one conceptual problem about a standard normal variate (SNV), $X$. If $X \sim N(0, 1^2)$, then the probability density function (pdf) is given as $f(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{1}{2} x^2);\quad -\infty<x<\infty$ _**My Problem :_** I understand $f(x)$ as defined above gives the probability that SNV takes the value $X = x$. Thus, $P(X=0)=f(0)$, and $f(0)= 0.398942$ which I understand is very high and it should be high only as mostly all values assumed by $X$ will cluster around $0$. Similarly, $f(-1) = f(1) = 0.241971$ and $f(-3) = f(3) = 0.004432$ However, as I understand it, since $X$ is a continuous variable, $P(X = 0)$ should be $1 / \infty = 0$. I am very confused now about my understanding about my concept of pdf. Does the bell shaped curve gives me plot of values assumed by random variable $X$ against the corresponding probability that $(X =x)$. I am not able to understand if my thinking process is right or wrong.


Can the mean deviation about mean exceed the standard deviation for the Pareto distribution? I just went through some books and found they are claiming that it cannot. How can I prove that? What is the underlying logic that confirms this conclusion?


title: Is it correct to use a paired t-test for this question?

text: this is an exam practice question for my 1st year stats class and I just want to see if I have approached this problem correctly and to clarify my understanding. Here is the complete question: A professor is deciding whether to schedule an end of semester Stats exam in the morning (9am) or in the afternoon (2pm). He's not sure if attention levels are lower or higher in the afternoon than in the morning. To test this he selects a group of 20 students and exposed them to a psych task. In this task lights randomly ash over a period of 20 minutes and the subject is required to press a button within a certain time of a light ashing. Failure to react to a flashing light is called a `lapse'. The more lapses the poorer the attention of the student. The students were given the task at 9am (AM) and again at 2pm (PM) and the number of lapses for each student is recorded in R as follows:               AM = c(1,5,4,2,6,6, 7, 11,5,3,4,5,8,3,4,5,7,7,4,5)         PM = c(4,5,5,8,9,10,12,15,4,8,6,6,9,5,9,8,7,8,7,6)      (a) State an appropriate null and alternative hypothesis, stating any parameters used. (b) State an appropriate test statistic to test this hypothesis and the null distribution of this test statistic, stating and checking any assumptions required. (c) Calculate the test statistic chosen in (b), the corresponding p-value and form an appropriate conclusion. (d) Calculate an appropriate confidence interval and provide a formal interpretation. * * * **Here is my answer to this question:** (a) My null hypothesis is that there is no difference in the average number of lapses for morning and afternoon. I.e. $H_{0}: U_{d} = 0 $ The alternative hypothesis is that the difference is greater than 0. i.e. $H_{1}: U_{d} > 0 $ The parameter of interest is the population mean difference (actually I don't completely understand what the question is asking here - when it says state any parameters used?) (b) My assumption is that the differences are independent and identically distributed normal random variables. i.e.: $D_{i} \sim \mathcal{N}(U_{d}, \sigma^2),\mathrm{iid} $ Upon checking the distribution in R - it does appear normal. The test statistic used is: $ t = \frac{\bar{D}-U_{d}}{S_{d}\sqrt(n)} $ Where $n$ is the number of pairs, $\bar{D}$ is the mean of the sample differences, $S_{d}$ the standard deviation of the differences, and $U_{d}$ is the mean of the population differences. Also assumed is that the null distribution of the test statistic is normal. (c) Calculating the data in R, I have $\bar{D} = -2.45$, $S_{d} = 1.959$, $n=20$. Therefore the observed test statistic is: $ t_{obs} = \frac{-2.45-0}{1.959\sqrt(20)} \approx -5.5917 $ Using a 5% level of significance, the critical value is $t^* = $ `qt(0.95,19)` $\approx 1.729133$ Since $ -5.5917 < -1.729133$, we reject the null hypothesis and conclude that there is evidence to suggest that there are less lapses (better attention) in students in the AM. (Am I correct to take the negative for the critical value here?) The corresponding p-value is $\approx 1.079061\times 10^{-5}$ (`pt(-5.5917,19)` in R) and since this is much smaller than $0.05$ we reject the null hypothesis. (d) Finally, looking for a 95% confidence interval - I observe that the alternative hypothesis is one sided $(H_{1}: U_{d} > 0)$, so a one sided confidence interval is used. The lower bound of the one sided CI is: $\bar{D}-t^{*}\frac{S_{d}}{\sqrt{n}} = -2.45 - 1.729133\frac{1.959457}{\sqrt{20}} \approx -3.207616002$ So the 95% CI is $(-3.207616002, \infty)$ since $ -5.5917$ falls outside this interval - so I conclude that there is a significant difference in student's attention levels between AM and PM (specifically, that it is better in the AM). Have I made any errors here in my logic here and are the calculations reasonable? Many thanks! Just want to make sure that I am on the right track!


title: Train / Validate / Test sets in Caret

text: I want to use caret to compare two different classification algorithms. For example SVM and Elastic net. I want to put aside some samples for test set and then use the rest of the samples for training the model, which involves tuning some parameters (like alpha and lambda for elastic net) for which I use cross validation as well. But once I fixed the model, I want to calculate it's performance on the test set (that was not used for parameter tuning). I want to do 10 fold cross validation to calculate performance of each method (SVM vs Enet). Note that this cross validation should be in an outer loop than the cross validation used for parameter tuning for each method. I was wondering if I can use Caret to do that. In the pseudo code given by Caret (http://caret.r-forge.r-project.org/training.html) the first loop (line 2) is for each parameter set, and the next loop (line 3) is for each resampling iteration. Is it possible to add another outer loop such as (Line 0: for each resampling iteration) to this schema ? * * * EDIT: @topepo: yes we have a limited data, so we want to have multiple test sets. Use 90% data for training (including parameter tuning) and then calculate the performance on the 10% remaining data. then repeat this 10 times with different partitionings of the data.


title: Can you use a moving average as an instrumental variable?

text: I have panel data and am interested in changes in total expenditures. I would like to consider an instrumental variable approach to deal with an endogenous regressor – the short run elasticity of income, lagged by 3 years. There are two ways I can think to do this:   1. I could calculate the elasticity using all years and then run a cross-sectional analysis.    2. I could calculate the elasticity by using a moving average of the last five years and then comparing that to expenditures in time _t_ \+ 3. However, this would make the instrument correlated over time. I don't think this violates the exclusion restriction in my case, but I am wondering if this violates other assumptions, such as randomness.


title: Probability questions for statistics

text: A set of final examination grades in a course is normally distributed with a mean of 73 and a standard deviation of 8.   1. What is the probability of getting a grade below 91 on the exam?    2. What is the probability that a student scored between 65 and 89?    3. If the professor grades on a curve (gives A’s to the top 10% of the class, regardless of the score), are you better off with a grade of 81 on this exam or a grade of 68 on a different exam, where the mean is 62 and the standard deviation is 3? Explain why.


title: Calculating the standard error for a ratio estimator in a stratified sample

text: I need to use Lohr's (2d Ed) method of calculating a standard error for an estimate generated from a ratio estimator based on a stratified sample. In working out her example on page 145-146, I find that the fraction in her example is calculated opposite to how she shows it in her MSE equations on page 144. Specifically, on p144, the fraction has the estimated $x$ as numerator, and the population $x$ in the denominator. but in her example she reverses them. Neither method produces her example result. I'm trying to figure out which version is correct. I would appreciate some help. (The Cochran method uses a sampling fraction in place of the fraction of pop and sample $x$'s).


title: The Monty Hall Problem - where does our intuition fail us?

text: From Wikipedia : > Suppose you're on a game show, and you're given the choice of three doors: > Behind one door is a car; behind the others, goats. You pick a door, say No. > 1, and the host, who knows what's behind the doors, opens another door, say > No. 3, which has a goat. He then says to you, "Do you want to pick door No. > 2?" Is it to your advantage to switch your choice? The answer is, of course, yes - but it's incredibly un-inituitive. What misunderstanding do most people have about probability that leads to us scratching our heads -- or better put; what general rule can we take away from this puzzle to better train our intuition in the future?


title: How can I calculate Conditional Probability of the following case?

text: Recently I attended a workshop on probability, I was asked the following question but I did not find a way to solve it. A discussion on this question can help to learn something new. We were told that, We are guests in a game show and close to win a great fortune. The quiz master asks us to choose one of three (closed) doors. She explains that behind one of them awaits you a million Euros. Once you fixed your choice the quiz mastress opens one of the other doors and shows you that this was only a goat. She gives you a final chance: you may either retain your door or switch to the remaining closed one. > (i) Say door 3 is opened. How can we Calculate the conditional probability > that our door is the winning one given that the door 3 is a fail, and its > complement. > > (ii) How to Calculate the unconditional probability that our door is the > winning one, and its complement. Thank you in advance.


title: What causes extremely small p values in Granger causality test?

text: I am generating some Granger causality tests for various time-series. I have noticed something that appears odd. I have timeseries $A$ and $B$. Sample size is about 1500 elements in each. $A$ and $B$ are highly correlated about 70%. When I do a 1 ,2,or 3 lag Granger test I get p values like 0.0000000000034. That seems way too low - I would be happy if it was around 0.01. Is there something else I need to take into account?


title: How to compare two groups on a measure of social skills that includes 5 subscales where each subscale is number correct out of 12?

text: I have two groups of participants. Each group contains 30 children (15 girls and 15 boys). The first group has been sampled from a population of children being raised under some kind of constitutional care, while the second group represents children living in their families. I want to measure their social skills and compare the results of the two groups. This is how the measurement goes: each child is presented a series of 12 pictures. From what it says about a single picture, I get 5 true/false values - each value represents a single social skill; true means the child proved the particular social skill when confronted with that picture and false means it didn't. In this way, I get 12 $\times$ 5 = 60 true/false values for each child. So to visualize my data, I have something like this $$ \begin{array} & & \mbox{Ch}_1 & \mbox{Ch}_2 & ... & \mbox{Ch}_{30} & \mbox{Ch}_{31} & ... & \mbox{Ch}_{60} \\\ \mbox{Skill}_1 & 12 & 11 & ... & 10 & 11 & ... & 11 \\\ \mbox{Skill}_2 & 9 & 10 & ... & 8 & 10 & ... & 12 \\\ \mbox{Skill}_3 & 10 & 10 & ... & 9 & 11 & ... & 10 \\\ \mbox{Skill}_4 & 7 & 7 & ... & 8 & 7 & ... & 7 \\\ \mbox{Skill}_5 & 5 & 6 & ... & 5 & 5 & ... & 6 \\\ \end{array} $$ Where an entry at $[i,j]$ denotes a number of $true$ values in the result list for the $i$-th social skill and the $j$-th child. The children 1 to 30 are from the first group and the children 31 to 60 are from the second group. Now, I would like, for each of the 5 social skills, to evaluate whether these two groups of children scored equally or if there was a statistically significant difference. From what I gathered from various sources (e.g. this), it seems like the two most likely usable tests for me would be the Fisher's test or the Mann-Whitney test. So my question is basically: **What kind of test would be appropriate for the given data and objective?**


title: How to calculate similarity in gene expression for each gene in two conditions and rank them?

text: I have expression values (log2) for 200 genes in two conditions treated and untreated and for each condition I have 20 replicates. The dataset looks like this:               Gene      UT1             UT2             T1              T2       DDR1      8.111795978   7.7606511867    7.9362235824    7.5974674936     RFC2      10.2418824097 9.7752152714    10.0085488406   9.5723427524     HSPA6     6.5850239731  6.7916563534    6.6883401632    7.3659252344     PAX8      9.2965160827  9.2031177653    9.249816924     8.667772504     GUCA1A    5.4828021059  5.3797749957    5.4312885508    5.1297319374      I have shown only two replicates for each sample in the sample data. **How can I calculate similarity in gene expression for each gene in two conditions and rank them from highest to lowest (in terms of similarity)?** I am looking for a solution in R or python. The `cor` function in R does not give me what i want.


title: Formative model identification problem using SEM

text: I have problem in identifying a formative structural equation model using AMOS.Does anyone have any suggestion for dealing with the identification problem in formative measurement models?


title: Sufficient statistic for function of parameter

text: Suppose we are estimating $\tau(\theta)=\theta e^{-\theta}$ from $X_1,...,X_n \sim \mathrm{G}(\theta,r)$ (G is the gamma distribution) then it is easily shown that $T=\sum_{i=1}^n \ln(X_i)$ is sufficient or $\exp(T)$ also. I'm wondering, is a continuous function of a sufficient statistic also sufficient and if $T$ is sufficient for $\theta$ is it also sufficient for $\tau(\theta)$?