title: Lmer model design

text: My data is a series of repeated measures in time (14 measures). I am trying to model the variable `HbA1c` which is a blood test performed at each visit to measure the global blood glucose level. Other predictors are recorded at each visit are BMI, age, insulin dose, and insulin type, etc.... By looking at the data in facetted (per Patient ID) scatter plots ($Y$: parameters, i.e. BMI, age, insulin dose and insulin type; $X$: Visit NR), I see different intercepts and slopes for the included parameters in the model. So is this the model I am looking for?               HbA1cfit <- lmer(HbA1c ~ VisitNR + BMI + Age + Insulindose + ... +          (1|PatientId) + (BMI|PatientId) + (Age|PatientId) + (...|PatientId),         data=type2diabetes ,REML=F)      The plot of the fitted values vs. the residuals shows a normal distribution with no clear trend. The `anova(HbA1cfit1, HbA1cfit2 [or more, each time model + extra parameter])` is significantly better every time. The scatter plots `(y=HbA1cRAW, x=visitNR)` and `(y=fitted(HbA1cfit), x=visitNR)` look very similar to me. One more question: how should I order the parameters in the hierarchy of significance in their determining power towards `HbA1c`? Let me rephrase the whole problem: Data set: PatientID Sex VisitNr Age Insulindose Insulintype C-peptide PO-drugs HbA1c   1 f 1 35 50 2 1.5 1 65   1 f 2 36 55 2 1.6 1 66   ...   1 f 14 42 60 3 0.2 2 70 2 m 1 60 50 4 2.5 2 80   ...   2 m 14 67 40 4 1.3 3 75 ... 485 m 1 50 20 3 2.5 2 50   ...   485 m 14 57 30 3 2.5 3 55 So data for 485 patients, most of them for 14 visits, if not completed: marked af missing data. Variables : PatientId : 1 to 485 Sex : categorical : f or m Visit nr : 1 to 14 Age : continuous integer, visits are typical 6 months appart , so age goes up with about 7 years from visit 1 to visit 14 Insulinedose: continuous integer Insulintype: factorial ordened : 1 = 1xlongacting , 2 = 2xmix , 3 = 3xmix, 4 = basal-bolus C-peptide: continuous double   PO-drugs : factorial ordened : 1 = none, 2 = 1 drug, 3 = 2 drugs, 4= 3 drugs HbA1c : continuous integer Questions to solve: which factors determine the outcome : HbA1c ? Time ? (as visit Nr ?) , Age, Others ??? In what order : most to least to not significant ? Is HbA1c going up or down in time significantly ? Since there are missing values all over the database and since I am dealing clearly with repeated measures, I tought it might be adressed by lmer in R. Everything seems to be nested in PtientId , only F or m stays the same in all visits, the other values can all change and do clearly not group the data. My effort to model this (if appropriate at all) is on top of this post. Please help , the more I read about lmer the more I get confused. Jan


title: Seasonality of treatment and Average Treatment Effect

text: I have panel data of sales for many stores in two comparable cities. One of the cities holds a special event once a month ( _the treatment_ ) which is expected to boost sales across the board on that day. I would like to estimate the average treatment effect on the treated, i.e. identify whether or not the special event does what it should. While the two cities are quite comparable (and I can find a number of controls for the stores), sampling is clearly not random. Furthermore, there is autocorrelation in sales between dates at each individual store, and it seems reasonable to think the monthly seasonality of the treatment has some impact. My question is twofold. 1/ Do you think running a difference-in-differences analysis in this “experimental” setup makes sense at all?   2/ The examples of diff-in-diffs I have seen are all pre-treatment/post- treatment. Do you have any idea of how I could account for the periodicity of the treatment? It seems wrong to consider each month separately because of the autocorrelation. Thank you very much for your help!


title: Using controls as basis for quantile/tertile categorization

text: I have a number of continuous but really skewed variables in the statistical analyses I am currently performing. Even though I am strongly in favour of keeping them just as they are and relaxing linearity assumptions by using splines, several faculty suggestions favour categorization for various reasons, including ease of interpretation, tradition, comparison with former studies in the field and some other unsupported claims. The most recent suggestion has been to categorize the variables in the case population based on the tertile/quantile distribution among the control population. The argument is that the control population is the "normal" standard against which the case population should be viewed. It appears to be fairly common in epidemiological studies. However, no one ever justifies the logic. See for example, strategy outlined in this SAS instructional site. a) Isn't the approach of categorizing variables in the case population using the distribution of the variables in the control population the same as using the outcome (case/control status) to determine categories and hence, strongly biases effect estimates? b) If I must categorize, is it better to categorize on the basis of the whole population (cases and controls) taken together or base the categorization on distribution in the controls alone? Any thoughts, references, insights on the pros and cons of this suggested control-based categorization method would be really appreciated! Thanks!


title: Clustering Survey Responses Based on Dichotomous Responses

text: I have a set of about ten questions that I would like to use to create groupings from. The responses are all dichotomous (responses are in the form of 1 or 2 where 1 and 2 represent differences in preference discovered through qualitative research). The questions, were provided in the form:               Which describes you best:     1. I prefer apples     2. I prefer pears      So far, I've looked at latent class analysis (challenging to interpret and not consistently reproducible), linear discriminant analysis, kmeans clustering (not consistently reproducible), and multiple correspondence analysis (MCA provided the most interpret-able results, but I'm unclear IF or how one could classify respondents using the results). What would be the most reasonable clustering method to use? If you want to play around with my data feel free to do so (n=799):               dat <- read.csv("http://www.bertelsen.ca/media/stackoverflow/cluster.csv")


title: Nonlinear regression: heteroskedasticity and correlated residuals

text: I'm performing regression analysis of some data. I believe I have to use a non-linear model with the form $y = at^b + c$ where $t$ is time. A log transform won't linearize the data here because of the constant term $c$. When I perform the regression, I get a model that fits the data very well. My R-squared value is 0.996 with an RMSE value 1/100th the scale length of the problem. However, my residuals are growing in time and they exhibit a sinusoidal pattern. I believe this means they are correlated with non-constant variance. How do I deal with this? All of the theory I've seen developed deals with serial correlation of errors or heteroskedasticity in linear regression only. My best guess is a weighted autoregressive model, but I'm not sure how to pick weights since I've tried the common suggestions and they aren't helping the problem.


title: Does more variables mean tighter confidence intervals?

text: Assume that the true (but unknown) relationship in a population between $Y$ and $X1, X2, X3, X4$ is $$Y=\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4.$$ Further assume that I have a sample from the population that I wish to model using linear regression and I fit the following two regression models: \begin{align} Model \;1: \quad Y & =\beta_0 + \beta_1 X_1 + \beta_2 X_2 \\ Model \;2: \quad Y & =\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 + \beta_4 X_4 \end{align} After obtaining my estimates $\hat{\beta}_0, \ldots, \hat{\beta}_4$ I then compute the predicted value of $Y$ for some covariate values $\tilde{x}_1, \ldots,\tilde{x}_4$. Additionally, I compute the confidence interval of my prediction, e.g. for model 1, $$ \tilde{y} = \hat{\beta}_0 + \hat{\beta}_1 \tilde{x}_1+\hat{\beta}_2 \tilde{x}_2 $$ with confidence interval in the usual way, i.e., $$\tilde{y} \pm t_{n-3} \hat{\sigma} \sqrt{\tilde{x}^T(X^TX)^{-1}\tilde{x}}$$ where $\tilde{x}=(\tilde{x}_1, \tilde{x}_2)^T$. Now I do the same thing with model 2, which is the more correct model and compare the interval width of the confidence interval for the prediction. Do I get a tighter confidence interval for my prediction? I tried it in R for a simple simulated data set and the answer was yes. In this example I varied the value of $\tilde{x}_1$ and fixed the values of the other covariates.               library(survival)     set.seed(13)     n <- 200          x1 <- rnorm(n)     x2 <- rbinom(n, 1, 0.5)     x3 <- rnorm(n)     x4 <- rbinom(n, 1, 0.5)     b0 <- 2     b1 <- 1.2     b2 <- 1.8     b3 <- -2.2     b4 <- 0.8     mu <- b0+b1*x1+b2*x2+b3*x3+b4*x4     eps <- rnorm(n,0,1.5)     y <- mu+eps          fit2 <- lm(y~x1+x2)     summary(fit2)     fit4 <- lm(y~x1+x2+x3+x4)     summary(fit4)          newdata <- data.frame(x1=seq(-2,2,0.1), x2=1, x3=0.5, x4=1)          pred2 <- predict(fit2, newdata, se.fit=TRUE,interval="confidence")     pred4 <- predict(fit4, newdata, se.fit=TRUE,interval="confidence")          plot(newdata$x1, pred2$fit[,3]-pred2$fit[,2],type="l", ylim=c(0,2), lwd=2, main="Interval     width (upper-lower bound)", ylab="Interval Width", xlab="x1")         lines(newdata$x1, pred4$fit[,3]-pred4$fit[,2], lty=2, lwd=2)     legend(-2,0.6, lty=c(1,2), lwd=2, legend=c("2 covariates", "4 covariates"))      ![enter image description here](http://i.stack.imgur.com/IAIgA.png) My question: Is this always the case? Or more general:   1. If the variables have an effect on my model, do I ALWAYS get a tighter confidence interval if I include more variables in my model?   2. On the other hand, if the variables do not have an effect, do I get wider confidence intervals for the prediction?   3. How about the standard errors of the parameter estiamtes?   4. Is this also true when confounding is present?   5. How about other regression models such as logistic regression or mixed models? Here I only considered confidence intervals for the prediction, but my question also applies to prediction intervals.


title: Testing for significant binary attributes

text: I have binary vectors, each with a binary outcome: for example:               |A|B|C| outcome      1 0 0 -> 0      0 0 1 -> 0      1 0 1 -> 0      0 1 0 -> 1      0 1 1 -> 1      How can I test which of my attributes (A,B and C in this example) have a significant impact on my outcome? Test the summed quadratic difference (in this binary case the same as the absolute difference) between the outcome and each variable over all vectors against 0? If yes, how is this test statistic distributed? Many attributes will be 1 in less than 1% of all cases including the outcome, but I have ca 1 million of those vectors. How do I deal with this while testing for significance? Do I even need to with this sample size? To test for covariance I would just add constructed variables like A*B, A*C, B*C and A*B*C. Am I right with this approach?


title: What is a suitable probability distribution for monetary amounts?

text: In my Bayesian model, I have a random variable representing the amount of money received due to a user's click on an affiliate link on a website. There are several such links with different payouts; in my model I want to capture the uncertainty due to the unknown behaviour of the user. The amounts are typically in the range of 0.30 EUR to 0.80 EUR. I am wondering: what is a suitable probability distribution for monetary amounts?


title: Conditional independence isn't guaranteed when specifying the marginal distributions?

text: It was mentioned very briefly in a lecture related to graphical models that two random variables $X_3$ and $X_4$ are both dependent on $X_2$. But even when conditioned on $X_2$, the two variables $X_3$ and $X_4$ could be dependent through some other means. I wasn't sure what the person meant by this, below is an example provided: Random variables: $X_1$, $X_2$, $X_3$, $X_4$ with the following density functions. $X_1 \sim \mathrm{Bernoulli}(1/2)$   $X_2 \sim \mathrm{Bernoulli}(1/2)$   $X_3 \mid (X_1,X_2) \sim \mathcal N(X_1+X_2,\sigma^2)$   $X_4 \mid X_2 \sim \mathcal N(aX_2+b,1)$ The above specification doesn't necessarily mean that: $p(x_3,x_4|x_2)=p(x_3|x_2)p(x_4|x_2)$ Meaning the conditional independence does not necessarily hold between $X_3$ and $X_4$ given $X_2$. The graphical model shows that the conditional independence will hold, but the specification I provided above does not guarantee it.


title: symmetric r.v. raised to an odd power

text: My prof claims that raising a symmetric r.v., like N(0,1), to an odd power gives a distribution with expectation 0. What's the best way to see this?


title: How can I interpret the results of LSA?

text: I implemented LSA on MATLAB. I have a $D\times N$ term-document matrix, where $D$: # of words, $N$: # of docs. I did low-rank approximation using SVD, and got $$X_k = U_k \cdot S_k \cdot V_k' (D=1000, N=600, K=4)$$ Now I want to classify the documents into 4 classes, and I know I have to use the col-vectors of $V_k'$. Each column of it has 4 values. I think each value indicates how much the document is related to the topic (in latent spaces). Am I right? But when I see the column's value, it has both positive and negative values. How can I interpret it?


title: Discrete functions: Confidence interval coverage?

text: How to calculate discrete interval coverage? What I know how to do: If I had a continuous model, I could define a 95% confidence interval for each of my predicted values, and then see how often the actual values were within the confidence interval. I might find that only 88% of the time did my 95% confidence interval cover the actual values. What I don't know how to do: How do I do this for a discrete model, such as poisson or gamma-poisson? What I have for this model is as follows, taking a single observation (out of over 100,000 I plan to generate:) Observation #: (arbitrary) Predicted value: 1.5 Predicted probability of 0: .223 Predicted probability of 1: .335 Predicted probability of 2: .251 Predicted probability of 3: .126 Predicted probability of 4: .048 Predicted probability of 5: .014 [and 5 or more is .019] ...(etc) Predicted probability of 100 (or to some otherwise unrealistic figure): .000 Actual value (an integer such as "4") Note that while I've given poisson values above, in the actual model a predicted value of 1.5 may have different predicted probabilities of 0,1,...100 across observations. I'm confused by the discreteness of the values. A "5" is obviously outside the 95% interval, since there's only .019 at 5 and above, which is less than .025. But there will be a lot of 4's -- individually they are within, but how do I jointly evaluate the number of 4's more appropriately? Why do I care? The models I'm looking at have been criticized for being accurate at the aggregate level but giving poor individual predictions. I want to see how much worse the poor individual predictions are than the inherently wide confidence intervals predicted by the model. I'm expecting the empirical coverage to be worse (e.g. I might find 88% of the values lie within the 95% confidence interval), but I hope only a bit worse.


title: What does $d$ mean in this notation of the "usual noninformative prior of $\mu_i$ and $\sigma_i$?"

text: Samiuddin, (1976) states: ![enter image description here](http://i.stack.imgur.com/jvBNh.png) or, typset with $\LaTeX$ as originally posted > We start with the usual noninformative prior distribution of $\mu_i$ and > $\sigma_i (i = 1,2,\ldots, k)$ > > $$\pi(\mu_1, \mu_2, \ldots, \mu_k; \sigma_1, \sigma_2, \ldots, \sigma_k) > \propto d\mu_1d\mu_2\ldots d\mu_k \frac{d\sigma_1d\sigma_2\ldots d\sigma_k} > {\sigma_1\sigma_2\ldots \sigma_k}$$ What does this notation mean? * * * Samiuddin, M. 1976. Bayesian Test of Homogeneity of Variance. Journal of the American Statistical Assoc. Vol. 71, No. 354


title: Kruskal-Wallis test data considerations

text: I'm about to apply Kruskal-Wallis test (non-parametric ANOVA) on three groups of unequal length. I was taught/advised to apply Krukal-Wallis only if:   * dependent variable is at least at ordinal level of measurement   * group's $ n > 5 $ (otherwise H statistic is not $ \chi ^2 $ distributed, so exact p-value cannot be calculated, etc.)   * Wikipedia page says that it requires _an identically-shaped and scaled distribution for each group_ Apart from that, what are general considerations when applying Kruskal-Wallis test? Should I check homoscedascity with Fligner-Killeen's test, and what should I do in case of small $ n $, or unequal group sizes? There are `exactRankTests::wilcox.exact` and `PASWR::wilcoxE.test` functions in R, but I can't seem to find an analogous one for Kruskal-Wallis... Advice, anyone?


title: Is it possible to perform k-NN classification using 5 fold cross validation?

text: I am going to validate my feature selection algorithm using k-nn classifer.


title: GBM package: Why there is a missing node?

text: Why there is a missingNode as 3 as there are no missing values? I have the data in the following form:               Y Score ot     1    10 3     0     7 2     0     9 2     0    22 1     0     5 2      Then               rawData <- read.delim("file.txt", header=TRUE, sep=" ");     data <- rawData[,c("Score","ot")];     y <- data[,c("Y")];     train.data <- as.data.frame(data);     train.data$y <- y;          fit.adaboost <- gbm(y ~ ., data=train.data, dist="adaboost", n.tree = 200,         shrinkage = 1, train.fraction = (2/12));     pretty.gbm.tree(fit.adaboost,i.tree =2);      When I do pretty I see                   SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight Prediction       0        0     6.5000000        1         2           3       3007.974  49087 -0.2998575       1       -1     0.8330843       -1        -1          -1          0.000   3061  0.8330843       2       -1    -0.3752048       -1        -1          -1          0.000  46026 -0.3752048       3       -1    -0.2998575       -1        -1          -1          0.000  49087 -0.2998575


title: Why am I getting different results for my logistic regression when performed by different software?

text: My data is simple, my independant variable is continous from 0-1000 and the response is either a 1 or a 0. I'm performing a logistic regression to determine the 50% inflection point. When I put this data into Statgraphics 5.1 I get a different respose from my collegaue who is using Statgraphics XVI (my organization is very behind on its software purchases). Is there a reason why these two software versions should give a different response when fed the exact same variables? Quick excel plot of the two models, and the green dots are the data ![enter image description here](http://i.stack.imgur.com/bOecz.png)


title: Can Neyman-Pearson lemma apply to the case when simple null and alternative don't belong to the same family of distributions?

text: 1. Can the Neyman-Pearson lemma apply to the case when a simple null and a simple alternative don't belong to the same family of distributions? From its proof, I don't see why it can't. For example, when the simple null is a normal distribution and the simple alternative is a exponential distribution.   2. Is the likelihood ratio test a good way to test a composite null against a composite alternative when both belong to different families of distributions? Thanks and regards!


title: hyperspherical nature of K means (and other squared error clustering method)

text: This famous paper points out that any squared error based clustering method (K means being the most common example) tends to generate hyperspherical clusters. However, it does not point to a theoretical justification. I would like to know: 1. If there is a theoretical justification for this observation. 2. If the K-Means algorithm can be changed to generate rectangular clusters.


title: What is model adequacy checking in ANOVA?

text: Does model adequacy checking mean "Checking the normality assumption that $\epsilon_{ij}\sim \mathcal{N}(0,\sigma^2)$"? More specifically, does it mean checking that the residuals, $e_{ij}=y_{ij}−\hat{y}_{ij}$, are distributed in this way? And is an independence check performed by showing that the mean square due to treatment and mean square due to error independently follow a chi-squared distribution according to Cochran's theorem?


title: Variance-covariance matrix of the parameter estimates wrongly calculated?

text: I fitted an hyperbolic distribution to my data with the `hyperbFit(mydata,hessian=TRUE)` command (package HyperbolicDist). The hessian looks like:               > hyperbfitmymodel$hessian                 hyperbPi      lZeta     lDelta            mu     hyperbPi   536.61654  -23.82800   25.62345   26153.16561     lZeta      -23.82800  250.74196 -261.20570     -35.58481     lDelta      25.62345 -261.20570  272.77771     182.75927     mu       26153.16561  -35.58481  182.75927 2028904.75586      Now I want to calculate the variance-covariance matrix of the parameter estimates, according to this page 2: > The asymptotic covariance matrix of $\hat{\theta}$ is given by the inverse > of the negative of the Hessian matrix evaluated at $\hat{\theta}$. I therefore calculate:               solve(-hyperbfitalv$hessian)      which gives                        hyperbPi         lZeta        lDelta            mu     hyperbPi -5.113433e-03 -0.0091511819 -0.0083271877  6.650321e-05     lZeta    -9.151182e-03 -1.6617499980 -1.5905496996  2.320893e-04     lDelta   -8.327188e-03 -1.5905496996 -1.5261031428  2.169113e-04     mu        6.650321e-05  0.0002320893  0.0002169113 -1.365591e-06      This looks clearly wrong to me, because there are negative values for the variance, but a variance cannot be negative? The covariance yes, but not the variance? EDIT: The complete output of `hyperbFit(mydata,hessian=TRUE)`:               Data:     mydata     Parameter estimates:            pi       zeta      delta         mu        0.090747   0.204827   0.002035  -0.002494       Likelihood:         756.911      Method:             Nelder-Mead      Convergence code:   0      Iterations:         365       2nd EDIT: If I use `solve(hyperbfitalv$hessian)` I get                            hyperbPi         lZeta        lDelta            mu     hyperbPi  5.113433e-03  0.0091511819  0.0083271877 -6.650321e-05     lZeta     9.151182e-03  1.6617499980  1.5905496996 -2.320893e-04     lDelta    8.327188e-03  1.5905496996  1.5261031428 -2.169113e-04     mu       -6.650321e-05 -0.0002320893 -0.0002169113  1.365591e-06      3rd EDIT: The output of `summary(hyperbfitalv)`:               Data:      mydata     Parameter estimates:            pi          zeta         delta         mu             0.090747     0.204827     0.002035    -0.002494        ( 0.071508)  ( 0.264040)  ( 0.002514)  ( 0.001169)     Likelihood:         756.911      Method:             Nelder-Mead      Convergence code:   0      Iterations:         365       4th EDIT: Ok, this is the hessian of pi, log(zeta), log(delta), and mu but how can I get the hessian of pi, zeta, delta and mu?


title: Statistical softwares providing multiple seasonality forecasting

text: Does anybody know which software provide forecasting with more than two seasonalities? I know R includes tools for double seasonal holt winter forecasting but I am looking for a software which can forecast triple seasonal data.


title: Statistically comparing results from two methods

text: I am not exactly sure if CrossValidated is the place to ask this question, but decided to give it a shot. This may be thought as an abstract question. I am trying to compare the feasibility of two different methods, let's call them:               Method A     Method B      Both methods take the same parameters and yield a real number greater than 0 as the result:               Method A(parameters) -> a real number greater than 0     Method B(parameters) -> a real number greater than 0      If the result is greater than 1.0, this is a "pass". If the result is less than 1.0, this is a "fail". The closer the number is to 0, the more severe is the situation. Assuming both methods are reliable, I am trying to figure out which method would give me better results (let's say more "passes"). First, I thought normalizing the results from one method with the other, I could create some kind of index (let's call it `i`, which would let me achieve the goal mentioned above:               i = Result from Method A / Result from Method B > 1.0 -> Method A gives a better result.      However, using this approach I loose the severity aspect of the problem due to the normalization, meaning: **Example 1:**               Result from method A = 5.65     Result from method B = 4.00         i = Result from method A / Result from Method B = 1.41 (Method A gives a better result)      **Example 2:**               Result from method A = 1.10     Result from method B = 0.78         i = Result from method A / Result from Method B = 1.41 (Method A gives a better result)      As seen above both indices yield the same value, 1,41. However, using method A or B does not affect the severity of the result in Example 1 (both > 1.0). Using any method would be fine. In the meanwhile, using method A in Example 2 is crucial, since it yields a result > 1.0 (pass), whereas the result from method B is < 1.0 (fail). I am wondering what would be the best way to compare these two methods? A somehow weighted index (this is just intuition, I am not sure how to construct something like this) This question may belong to the field of decision making, however, I am not very knowledgeable about that field (as I am not in statistics, other than basics). Any pointers would be welcome! Thanks! PS: I am not sure if I used correct tags for this question. Please let me know if there is a suitable tag, then I can edit and correct them.


title: What is the meaning of "finite sample error control"?

text: I encountered this phrase while reading a paper which goes like this -- "These methods lack finite sample error control due to instability". Although it might not be important, the paper deals with a feature selection strategy for classification.


title: a representation of ARMA(1,1) process

text: let $X_{t}=\frac{4}{5}X_{t-1}+u_{t}$ with $u_t$ a white noise with $\textbf{E}(u^2_{t})=\frac{9}{50}$. assume that $X_{t}$ is affected by ameasurement error,and let $Y_{t}=X_{t}+v_{t}$ where $v_{t}$ is a whine noise uncorrelated with $u_{t}$ and with $\textbf{E}(v^2_{t})=\frac{1}{2}$. Show that $Y_{t}$ can be written as an ARMA(1,1) process, and derive a representation of $Y_{t}$ of the type $Y_{t}=-\sum_{i\geq 1}\alpha_{i}Y_{t-i}+\epsilon_{t}$ with $\epsilon_{t}$ a white noise such that $\textbf{E}(\epsilon_{t}Y_{t-i})=0$ $\forall i>0$, $\forall t$ Some help would be appreciated


title: "effects not estimable" error in R aov

text: I have data on the effect of treatment on leg shaking. I reused several individuals so I want to know   1. Is there a significant effect of treatment on shaking   2. Do the reused individuals present some effect or significant variation Nevertheless when I run               model1<-aov(Legshaking~Idmale*Idfemale*tratament)      or               model1<-aov(Legshaking~Idmale+Idfemale+tratament)      I get               13744 out of 13824 effects not estimable     Estimated effects may be unbalanced      What could be wrong here? Any advice about how to handle this data in R will be welcome! I have attached an image picture so you can have an idea of the data set * * * Screenshot of data: ![enter image description here](http://i.stack.imgur.com/5O1aZ.png)


title: How to test differences in multiple pairwise distances? With mixed model?

text: I want to compare two set of variables:   * X1 = All the pairwise distances between individuals from group A and group ref R   * X2 = ALL the pairwise distances between individuals from group B and group ref R is it ok to test if the distribution X1 and X2 are different with a mixed model?               lme(X ~ Y, random = ~ 1 | ID/ID_ref )          X<-c(X1, X2)      Y is a factor describing the type of comparison with two levels A-R/B-R               ID<- c(IDa, IDb) = a list with the identity of the sample involved in each comparison in group A and group B          ID_ref<-a list with the identity of the sample involved in the comparison in group of reference      For example dataset will look like               X  Y  ID  ID_ref          5 A-R A1  R1          7 A-R A1  R2          2 A-R A2  R1          10 A-R A2  R2          1 B-R B1  R1          6 B-R B1  R2          0 B-R B2  R1          7 B-R B2  R2


title: Is it legitimate to present results from univariate and multiple regressions?

text: I recently saw a paper that presented the results of a multiple regression, and then proceeded to also present the results from univariate regressions for the independent variables in which they were most interested. The multiple regression already showed 4 significant independent variables (interactions were not considered), is it then legitimate to present, and interpret, results from univariate regressions for these variables? I would particularly appreciate references to books or published articles that discuss this issue. Thanks!


title: Can I copy output from R to Word or Excel?

text: I would like to put an output from R to my appendix in my thesis; however when I pasted it to Word, the columns were moved. Do you have any idea how to adjust it to be exactly the same as in R? so that it is easier to read


title: Bimodal curve fitting

text: I have been at it for days and simply can't solve this. I would like to fit a curve with two peaks (i.e. bimodal) to my data set. I have both x and y values in my data (i.e. not a probability distribution data set). The curve i fit may or may not be made up of two normal distributions (the distributions might instead have different slopes left and right of the peak) - not sure. y-values are only ever between 0 and 1 (these are probabilities) x-values are not evenly distributed the data look something like this: x: [-0.250 -0.100 -0.050 -0.038 -0.025 0.00 0.025 0.038 0.050 0.100 0.250] y= [ 0.590 0.700 0.890 0.780 0.675 0.68 0.650 0.805 0.825 0.630 0.540] Any ideas would be very much appreciated! I have played around with nlinfit, glmfit, fit (with 'gauss2' specified) and gmdistribution.fit but haven't solved this. In the end, i would like to have the x-values for each peak, and the stds of each, and ideally also the x-value of the valley/trough (between the peaks). I'd prefer to fit with mle, but least squares would be fine too. Vanessa


title: Significance of difference in means

text: In “A Comparative Study of Nursing Home Residents …”, Aigner et al. write (p. 20) > Acute visits were significantly higher for the nurse practiticoner/physician > team at 3.0 visits per year (± 2.4) versus 1.2 visits per year (± 1.5) for > the physician-only group (P <0.0001) Furthermore, (p. 19) > The statistical analysis was performed using the chi-squared or Fisher exact > test for comparisons of percent and Student t test for comparison of means I always thought that if confidence intervals overlap by that much, differences can’t be significant. Am I wrong? Aigner, M. J., Drew, S., & Phipps, J. (2004). A comparative study of nursing home resident outcomes between care provided by nurse practitioners/physicians versus physicians only. Journal of the American Medical Directors Association, 5(1), 16-23.


title: Can variance be equal to mean?

text: I understand that variance can be zero if all data points are equal. Can anybody state an example, where sample variance is equal to sample mean and what - if any - is its significance to distributions like the Poisson.


title: Finding the MLE of parameter $\mu$

text: Suppose $X_1,\ldots,X_n$ are a random sample of a continuous and strictly increasing distribution $F(x)$ with mean $\mu$. If $$ Y_i = \begin {cases} 2 & \text{if}\ X_i>\mu \\ 1 & \text{if} \ X_i\leq\mu \end {cases} $$ Determine the MLE of parameter $\mu$


title: Why isn't the Anova() function in the Car package returning an F statistic?

text: I ran a simple psychology experiment that included 4 conditions, each containing 8 blocks of training. There were different participants in each condition. Hence, condition and block and subject and block are crossed, but subject is nested in condition. I'm trying to do a basic repeated measures anova to test for effects of block and condition. I have an unbalanced design, and a mixed effects model. My approach was to use the lme4 package in conjunction with the car package. I am running R version 2.14 on mac os x lion. Here is what I've done:               library(nlme)     library(car)          rm( list = ls() )          data <- read.table("anova_data", header = TRUE)          condition <- factor(data$condition)     block <- factor(data$block)     subject <- factor(data$subject)     accuracy <- data$accuracy          fm1 <- lmer( accuracy ~ block*condition + (1|subject %in% condition) )          Anova(fm1)      My problem is that this returns a summary table without F values, like such:               Analysis of Deviance Table (Type II tests)          Response: accuracy                           Chisq Df Pr(>Chisq)          block           17.169  7    0.01634 *            condition       68.294  3  9.897e-15 ***          block:condition 26.481 21    0.18869      Any help is greatly appreciated.


title: How much does the value of a car depreciate due to time (only), and mileage (only)?

text: I have a database of second hand cars. It contains (among other things) the asking price, the mileage, the fuel consumption and the year it was built. I would like to know how much the value depreciates due to time (only), and due to mileage (only). So for car model X, I'd like to know: - how much value the car loses with each year - how much value the car loses with each mile My dated statistical knowledge from college tells me that I can do an OLS regression with price as dependent variable and year as independent variable, controlling for mileage. And a regression analysis with price as dependent variable and mileage as independent variable, controlling for year. Would this work? How much observations would I need per car model? Year seems to have an exponential effect on value, and is it an ordinal variable that creates problems in such a model? Are there any other solutions? Any help would be greatly appreciated.


title: Regularized logistic regression and support vector machine

text: L2 regularized logistic regression differs with L2 regularized support vector machine with their loss function. Are there more deep differences for these two models? I tried several data sets, and found L2 regularized logistic regression are always better than L2 regulared support vector machine. Are there any references/researches for discussion their differences in terms of classification performance?


title: ANOVA with Repeated Measurements based on randomly assigned repetitions

text: first of all the basic idea of my experiment: changes/differences of reaction time (RT) after 3 different types of sport (a: no jogging=SP1; b: 10min of jogging=SP2 and c: 30min of jogging=SP3); 36 subjects (SUB) are randomly assigned to one of the 6 possible combinations of “sport” (abc, acb, bac, bca, cab and cba, see variable FCT). So the first subject has "a" on day 1, "b" on day 2 and "c" on day 3; the second one might have "a" on day 1, "c" on day 2 and "b" on day 3 [etc...] Therefore I have the following table containing my raw data:               SUB    SP1    SP2    SP3    RT1    RT2    RT3    FCT     s01     a      b      c     0.1    0.4    0.8    abc     s02     a      c      b     0.6    0.3    0.8    acb     s03     b      a      c     1.0    1.0    0.7    bac     s04     b      c      a     0.3    0.7    0.1    bca     s05     c      a      b     0.3    0.5    0.4    cab     s06     c      b      a     0.0    1.0    0.3    cba     [.]    [.]    [.]    [.]    [.]    [.]    [.]    [.]     s36     c      b      a     0.1    0.9    0.0    cba      In general I think I have to use an ANOVA based on a repeated measurement design - but I have no clue on how to deal with the factor “sport” since it also changes with the repeated measures per subject; Is it OK to use the variable FCT defined by 6 possible combinations of “sport” as one factor with 6 categories with the repeated measurements-ANOVA? Thank you for comments!


title: How to define confusion matrix of classification rules

text: How to define confusion matrix of the database and the classification rules are found below. And calculation precision and recall. MinSup = 3% and MinConf = 30%               @attribute temperature real      @attribute humidity real     @attribute windy {TRUE, FALSE}     @attribute play {yes, no}            No.  outlook temperature humidity windy play         1   sunny       85.0    85.0    FALSE   no         2   sunny       80.0    90.0    TRUE    no         3   overcast    83.0    86.0    FALSE   yes         4   rainy       70.0    96.0    FALSE   yes         5   rainy       68.0    80.0    FALSE   yes         6   rainy       65.0    70.0    TRUE    no         7   overcast    64.0    65.0    TRUE    yes         8   sunny       72.0    95.0    FALSE   no         9   sunny       69.0    70.0    FALSE   yes         10  rainy       75.0    80.0    FALSE   yes         11  sunny       75.0    70.0    TRUE    yes         12  overcast    72.0    90.0    TRUE    yes         13  overcast    81.0    75.0    FALSE   yes         14  rainy       71.0    91.0    TRUE    no          Classify rules found:     1: (outlook,overcast) -> (play,yes)     2: (outlook,rainy), (windy,FALSE) -> (play,yes)     3: (outlook,rainy), (windy,TRUE) -> (play,no)     4: (windy,FALSE) -> (play,yes)     5: (outlook,sunny), (windy,FALSE) -> (play,no)     6: (outlook,sunny) -> (play,no)     7: (windy,TRUE) -> (play,yes)


title: Where to find datasets related to US elections at the individual level?

text: I am currently looking for datasets that I can use to predict how an individual voter will vote in a US election. Ideally, the dataset that I am looking for includes socioeconomic features such as age, sex, income, county, and include an outcome variable for the party/candidate that this person cast their vote for. Given that voting is highly confidential, it has been hard to really find such data - though I'm wondering if anyone on CrossValidated knows of anything out there or knows of a different dataset that I could use to predict how a particular individual will vote.


title: Finding statistically significant favorites among objects

text: I'm designing an experiment where subjects are asked to choose their favorite object among 3 different objects. Each subject gets the same 3 objects to choose from. How can I determine if one of the objects is favorited over the other objects in this experiment? The data looks like this:               Subject  Object 1  Object 2   Object 3     1        0         0          1     2        1         0          0     3        1         0          0     4        0         1          0     5        1         0          0     ...      A 1 denotes that a particular object was a favorite. As an example, Subject 4's favorite object was Object 2. I also have a small sample size, ~20-35 subjects (class room size), so any techniques for dealing with small samples would be desirable.


title: Estimation based Distance observations

text: Let $X_1,\dots,X_n$ are i.i.d. copies of $X$, however I don't know the value of them. I only know the distance for every pair $(i,j)$. Now, assume that we are given $Y$ and $Z$ which are also i.i.d. copies of $X$, and we will be interested to estimate the distance between $Y$ and $Z$. Is there a literature relating to my problem above? Can anyone help me?


title: Scaling a series of numbers

text: Suppose I have 30 numbers that vary between 0 and 1.0 and which sum to 1.0. The mean is obviously 0.033. A client wants these scaled to lie between 0 and 1.0 but to have a mean of 0.5. By the way, I'll probably have to do this for any set of numbers (positive, negative, outside the range 0 - 1) so that they lie in the range 0 - 1 and have a mean 0.5. Any suggestions?


title: Clustering groups based on frequency and last date of choice

text: I have data coming in the following format: UserId, GroupId, Frequency (how many times the user chose the group), Max timestamp (the last time the user chose the group). Based on this dataset we need to figure out which groups look alike. I was thinking of using collaborative filtering but I have 3 concerns:   1. We don't have any knowledge of "Dislikes"; we only know which groups users "Like".   2. We don't really have ratings. In other words, users don't rate the group. Either they choose OR they don't.   3. Frequency doesn't really imply interest level. Should I use value '1' to imply interest & '0' to imply 'don’t know the interest level'?  Should I use collaborative filtering or is there a more appropriate machine learning algorithm? Thanks for your help.


title: Is Ye Shiwen's 400m IM performance in Olympics "anomalous" statistically?

text: The 16-year-old Chinese swimmer Ye Shiwen swam more than 5 seconds faster than her personal best in women's 400m individual medley in London 2012 Olympics, winning her a gold. Her performance received high exposure to media and mixed comments, even some controversy. For example, Nature published a news report titled "Why great Olympic feats raise suspicions -- 'Performance profiling' could help to dispel doubts" and immediately raised great attention and a great number of reader comments. As an effort for myself to better understand the issue and an exercise, I attempted a 'performance profiling' as what the Nature news report's subtitle says. What I did is to analyze the short-term performance increase (or decrease) of top swimmers in 400m IM worldwide during 2011 and 2012, and try to see how Ye Shiwen compares to other good performance improvers. I found there are more than 20 swimmers who are arguably better performance improvers than Ye. What's also worth noting is that Ye is at a particularly interesting spot that her performance is between the majority male and majority female swimmers. Probably not exactly a surprise given she's world's No. 1 in the event. See my summary: http://bit.ly/swimdata. But I wonder if some of the critical eyes here can "cross validate" this approach. What kind of factors should be examined in this kind of "performance profiling" for sports, and specifically for the Ye case?


title: From confidence interval to standard deviation - what am I missing?

text: Given: Crude incidence rate per 100,000: 164 population: 211,101 95% CI: 147-182 What is the standard deviation of the incidence rate : My calc: Standard error:(182-147) / (2 x 1.96) = 8.93 Standard deviation = 8.93 x sqrt(211,101) = 4103 This can't be right, and I think it must be to do with the 100,000 Is this correct: Standard deviation = 8.93 x sqrt(211,101/100000) = 12.98 ? Many thanks !


title: What is the most suitable statistical test for a series of data?

text: Sorry, this is a homework help question: > In a bacterial growth experiment, ten μ (exponential growth rate constant) > data points have been collected: one each of nine different concentrations > of an antibiotic, plus one control (no antibiotic). What is the most > appropriate statistical treatment for such a data set? The background of the question is also that each data point has been collected by running a single experiment of measuring population against time, and calculating μ from the exponential phase of the growth curve (probably using linear regression of log(population) against time to get slope μ). I'm thinking that running a chi squared test of each data point against the control would be sufficient to show significance, but I'm not really sure, and would greatly appreciate the input of the learned statisticians here!


title: Parametric or non-parametric tests? ... or transformation?

text: I am using R to examine the relationship between two variables in a small data set ($n=16$). My problem is that I'm not really sure how to handle the analysis (read: I'm in deeper waters than I've traveled before). Do I use Spearman's Rho to calculate the rank coefficient? Or do I assume normality and use standard parametric tests (Pearson's)? Or do I transform the data, and then do the parametric tests? Or should I do something else? My initial feeling is that I should use non-parametric tests on this data set for various reasons:   1. I feel strongly that there is an upper and lower bound to the possible values observed for either the dependent or independent variable.    2. I do not feel that the relationship between the variables is linear. I've provided the QQ Plot of the dependent variable, and the scatterplot of the dep & indep vars below, and neither seems to adhere to the normality assumptions: The thing is, I'm getting hung up on several things:   1. the Pearson's coefficient seems really high 0.94366 which makes me wonder if I'm sacrificing something by using Spearman's instead (for the sake of full disclosure, the Spearman's coefficient of 0.86765 is also significant at $\alpha = .001$). **Now Answered**   2. If Spearman's is the answer, what is the next step in building a predictive model? **Now Answered**   3. Since the sample is relatively small, should I use some kind of resampling to calculate the correlation instead of just using the 16 values?   4. What other major things might I be overlooking? * * * Plots QQ Plot of Indep Var: ![enter image description here](http://i.stack.imgur.com/qc5wV.png) http://dl.dropbox.com/u/27272488/qqPlotY.png Scatterplot of Dep & Indep Vars: ![enter image description here](http://i.stack.imgur.com/2pGTb.png) http://dl.dropbox.com/u/27272488/XvsY.png * * * Data               > df$year         [1] 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009         > df$x     [1] 1.39 1.41 1.48 1.46 1.24 1.34 1.70 1.61 1.47 1.69 1.94 2.30 2.51 2.64 3.01 2.14     > df$y     [1] 2320 2227 2161 2116 2294 2483 2897 3197 3270 3714 4028 4576 4837 5174 5312 4462


title: When is it incorrect to draw stacked area charts?

text: Stacked area charts work great when charting raw values, averages, counts, and percentage of total. The sum of the values across the different groups will add up to that of the entire set. But what about in the case when the value being charted is a Median, Geo Mean, or any Percentile? In these cases the sum of the values (say Media) across the categories will not add up to the Median value of the entire set. Example: We have the data on the consumption of Meat, Vegetables, Fruits, and Dairy for a country. For each category (and the entire set) we have the median pounds consumed per person, each year, for the last 10 years. If we were to draw them in a stacked area chart, the sum of the values for a given year would not equal the median pounds consumed per person across all the categories. Is it fair to assume that the stacked charts should be avoided in these cases as they will cause confusion?


title: "t-Chart" value

text: I wonder what is wrong with my own answer for the problem which I have shared below, could you please help me some about that: " Suppose that we have 12 observation units $ ( n=12 ) $ in a statistical investigation. What is the "t-table" value $ ( {\downarrow}v^{ \backslash ^{\large{p} \, \rightarrow} } \; ) $ for $ \, t_{0,01} $ ? " And, I thought, it would be better to add also a link for a sample of "t-chart" I have tried to find the solution $ (v=n-2=10 \; \land \; t_{0,01}= 3,169 ) $ but, it did not work for me so far because, the answer key for this problem shows $ ' 3,106 ' $ instead of $ ' 3,169 ' $ I am just guessing that there might be something wrong with the answer key, so what do you think in order to get a correct solution for this problem? Please, show me how to do that. Thank you very much.


title: Proportional hazards assumption meaning

text: Suppose we have a Cox-PH model and death is the outcome of interest. Suppose we have the predictors $X_1, X_2$ and $X_3$. Then does the PH assumption mean that for any person (i.e. arbitrary values of $X_1,X_2$ and $X_3$) the hazard of death is $c \cdot h_{0}(t)$ where $c$ is a constant? So the hazard function is the same for everyone?


title: Isn't a test for cointegration the same as testing for heteroskedasticity in the residual error terms?

text: Testing for cointegration tests to see if the residuals from a regression between the two variables is a stationary process. So wouldn't a test on those residuals for heteroskedasticity be the same thing? To detect if the variance of the errors is stable and the error terms are not a random walk? Wouldn't a White test or Breusch-Pagan test for heteroskedasticity be the same as an Augmented Dickey-Fuller test? Both would detect if the errors are getting bigger or have a constant variance, right? An increasing variance would imply a random walk and unit root, correct?


title: "War stories" a la Skiena for clustering analysis?

text: what I mean by war stories is "tales from experience with real problems" where researchers "describe real problems and the process that led to the algorithm that solved the problems". Here are some examples (link). Basically I'm looking for some examples of real world problems that were solved by applying clustering analysis, and how those problems were solved (e.g. why clustering analysis was thought to be a good fit, how the algorithm was chosen & tuned etc.)


title: Alternative name for "incidence rate ratio"

text: I'm using a negative binomial GLM to predict individual level scores on a psychometric test that ranges from 0-27 and presenting the coefficient "incidence rate ratios". It seems a little strange to use the term "incidence rate ratio" when actually it is a psychometric score ratio. Should I still use the term "incidence rate ratio" in my manuscript for easy recognition/interpretation? If not, is there another, more generic term that is used or should I make up my own domain specific term (e.g. PGSI ratio where PGSI is the psychometric index I am modelling).


title: Pointers on how to perform a Monte Carlo Analysis

text: Im looking for some pointers on how to perform a Monte Carlo simulation. Say im conducting a survey of how late flights are for a certain airline. I amass 1000 records, and plot a cumulative sum of these values. I want to determine how likely it is to have several flights grouped together that a decrease in the cumulative sum would exceed X%. Would i 1) calculate the mean and standard deviation of this sample, use these to generate multiple samples, and record the probability of such a decrease? 2) take the original sample, shuffle it, and then select values from it (without replacement), and record the probability of such a decrease?


title: What does correlation between two beta coefficients mean?

text: I am now trying to use PRODCLIN to do the mediation analysis to obtain a confidence interval for binary single mediation analysis. http://ripl.faculty.asu.edu/downloadable-programs/prodclin/ All of the variables I used are binary variables, including IV, mediator and DV. I plugged in the numbers of betas and SEs from the two logit regression models: estimated a and b. [a is the coefficient of IV when regressing the mediator(M) on IV(X). M= constant + aX + residual] [b is the coefficient of the mediator, when r egressing the DV(Y) on both IV(X) and the mediator(M). Y=constant+ c'X + bM + residual] The program asked for the correlation rho between a and b. I am not sure what it means by correlation between two coefficients, especially when they were estimated based on logit regression models. I don't know how to obtain the correlation between these two coefficients. I assume that assuming zero is the most conservative way?


title: Bayesian approach and least-squares approach to multivariate regression with structural design

text: Assume for example a trivariate Gaussian model: $$ {\boldsymbol Y}_1, \ldots, {\boldsymbol Y}_n \sim_{\text{iid}} {\cal N}_3\left({\boldsymbol \mu}, \Sigma\right) \quad (*) $$ with ${\boldsymbol \mu} = (\mu_1,\mu_2,\mu_3)'$. The Bayesian conjugate theory of this model is well known. This model is the most simple case of a multivariate linear regression model. And more generally, there is a well known Bayesian conjugate theory of multivariate linear regression, which is the extension to the case when the multivariate mean ${\boldsymbol \mu}= {\boldsymbol \mu}(x_i)$ is allowed to depend on the covariates $x_i$ of individual $i \in \{1, \ldots, n \}$, with linear constraints about the multivariate means ${\boldsymbol \mu}(x_i)$. See for instance Rossi & al's book accompanied by the crantastic `bayesm` package for `R`. We know in addition that the Jeffreys prior is a limit form of the conjugate prior distributions. Let us come back to the simple multivariate Gaussian model $(*)$. Instead of generalizing this model to the case of linearly dependent multivariate means ${\boldsymbol \mu}(x_i)$ depending on individuals $i=1,\ldots,n$, we can consider a more restrictive model by assuming linear constraints about the components $\mu_1$, $\mu_2$, $\mu_3$ of the multivariate mean ${\boldsymbol \mu}$. ### Example Consider some concentrations $x_{i,t}$ of $4$ blood samples $i=1,2,3,4$ measured at $3$ timepoints $t=t_1,t_2,t_3$. Assume that the samples are independent and that the series of the three measurements $(x_{i,t_1},x_{i,t_2},x_{i,t_3}) \sim {\cal N}_3\left({\boldsymbol \mu}, \Sigma\right)$ for each sample $i$ with a mean ${\boldsymbol \mu} = (\mu_1,\mu_2,\mu_3)'$ whose three components are linearly related to the timepoints: $\mu_j = \alpha + \beta t_j$. This example falls into the context of _Multivariate linear regression with a within-design structure_. See for instance O'Brien & Kaiser 1985 and Fox & Weisberg's appendix So my example is a simple example of this situation because there are only some predictors (the timepoints) for the components of the mean, but there are no predictors for individuals. This example could be written as follows: $$ (**) \left\{\begin{matrix} {\boldsymbol Y}_1, \ldots, {\boldsymbol Y}_n \sim_{\text{iid}} {\cal N}_3\left({\boldsymbol \mu}, \Sigma\right) \\ {\boldsymbol \mu} = X {\boldsymbol \beta} \end{matrix}\right. $$ with ${\boldsymbol \beta}=(\alpha, \beta)'$ and $X=\begin{pmatrix} 1 & t_1 \\ 1 & t_2 \\ 1 & t_3 \end{pmatrix}$ is the matrix of covariates for the components of the multivariate mean ${\boldsymbol \mu}$. The second line of $(**)$ could be termed as the _within design_ , or the _repeated measures design_ , or the _structural design_ (I would appreciate if a specialist had some comments about this vocabulary). I think such a model can be fitted as a generalized least-squares model, as follows in `R` :               gls(response ~ "between covariates" , data=dat,        correlation=corSymm(form=  ~ "within covariates" | individual ))      (after stacking the data in long format). **My first question** is Bayesian: what about the Bayesian analysis of model $(**)$ and more generally the Bayesian analysis of the multivariate linear regression model with a structural design ? Is there a conjugate family ? What about the Jeffreys prior ? Is there an appropriate R package to perform this Bayesian analysis ? **My second question** is not Bayesian: I have recently discovered some possibilities of John Fox's great `car` package to analyse such models with ordinary least squares theory (the `Anova()` function with the `idesign` argument --- see Fox & Weisberg's appendix). Perhaps I'm wrong, but I am under the impression that this package only allows to get the MANOVA table (sum of squares analysis) with an orthogonal matrix $X$, and I'd like to get (exact) confidence intervals about the within-design parameters for an arbitrary matrix $X$, as well as prediction intervals. Is there a way to do so with `R` using ordinary least squares ?


title: Rank data from 1 to N (N may not be equal to number of elements in the dataset)

text: I calculate `percentile ranks` based on the method at http://www.psychstat.missouristate.edu/introbook/sbk14m.htm I need to assign ranks to a dataset where ranks are from `1-to-N`. Since N may or may not be equal to the size of data, ordinal ranking is not the solution. I wanted to know if the following formula would be right. If not, please suggest. Thanks                        Fb + 0.5 Fw     Rank =  ------------- x N   (notations are similar to those in the URL)               SizeOfData


title: How do I determine if a survival model with missing data is appropriate?

text: Oversimplifying a bit, I have about a million records that record the entry time and exit time of people in a system spanning about ten years. Every record has an entry time, but not every record has an exit time. The mean time in the system is ~1 year. The missing exit times happen for two reasons:   1. The person has not left the system at the time the data was captured.   2. The person's exit time was not recorded. This happens to say 50% of the records The questions of interest are:   1. Are people spending less time in the system, and how much less time.   2. Are more exit times being recorded, and how many. We can model this by saying that the probability that an exit gets recorded varies linearly with time, and that the time in the system has a Weibull whose parameters vary linearly with time. We can then make a maximum likelihood estimate of the various parameters and eyeball the results and deem them plausible. We chose the Weibull distribution because it seems to be used in measuring lifetimes and is fun to say as opposed to fitting the data better than say a gamma distribution. Where should I look to get a clue as to how to do this correctly? We are somewhat mathematically savvy, but not extremely statistically savvy.


title: Cholesky ordering of variables in PVAR in Stata

text: I'm currently working on a panel of 4 countries and 29 years data. I want to arrange my variables _GDP_ , _government consumption_ , _government investment_ , and _revenues_ according to Cholesky ordering. Can anyone tell me that how I do this in Stata? I start my panel estimation in Stata by using the `xtset` command.


title: Analysis of multiple time series

text: i'm not sure how to google for this as i am not very familiar with time series analysis. i have 500 websites, and i am measuring the number of visitors to each website each day. at some point, i turn on SEO (search engine optimization) for each of the websites. this happens on different days for different websites. the distribution of visitors by account for any given day has a long tail to the right. SEO may not have an immediate effect; it might take a few days/weeks to really start to see some results. i want to measure something like the "average" lift in the number of visitors, but an average is probably not going to do the trick because of the mix of websites. (a daily/weekly/monthly trend curve would be really cool, but the average problem will exist there, too) i can probably average the number of visitors per day for any given account, but i can't do it across accounts. do i simply need to segment the websites into "number of visitor" groups? what other kinds of approaches should i read about?


title: The Lagrange multipliers of SVM

text: Actually the solve the SVM is to solve the following Lagrangian Equation:![enter image description here](http://i.stack.imgur.com/mF7lk.jpg) If we don't use kernel function, $\langle x^{(i)},x^{(j)}\rangle$ is just the vector vector inner product. The $a_ia_j\langle x^{(j)},x^{(j)}\rangle $ is the same to $\langle a_ix^{(i)}, a_jx^{(j)}\rangle $ in the formula.   But if we use the `kernel function`. **That's totally different, right?** **Q** : So Can I get the same train result between $\langle a_ix^{(i)}, a_jx^{(j)}\rangle$ and $a_ia_j\langle x^{(j)},x^{(j)}\rangle $?


title: What is the expected MINIMUM value drawn from a uniform distribution between 0 and 1 after n trials?

text: Assume you draw a uniformly distributed random number between 0 and 1 _n_ times. How would one go about calculating the expected _minimum_ number drawn after n trials? In addition, how would one go about calculating a confidence interval to state that the minimum number drawn is in the interval [a,b] with m% confidence?


title: How to calculate the probability that the maximum of two independent standard normal distributions is greater than zero?

text: The R command `hist(pmax(rnorm(10),rnorm(10)))` gives a histogram of random variable Y, where Y is the maximum observation of two independent Gaussians.   * How can I mathematically determine the probability that Y is greater than zero?   * How can I use a simulation to determine the probability that Y is greater than zero?


title: Nonparametric test for comparing trends in two time series

text: I am interested in doing a nonparametric hypothesis test in Stata. I have 2 series of data, each time series and that overlap in time period. 10 time periods each. I am interested in a test that checks whether the series are trending similarly. I am not interested in the levels. For example, suppose I have a series that happened to obey the equation: y = 2x + 100 and one that obeyed y = 2x + 1. They are clearly very different in level but identical in how they trend. I would want the test to fail to reject any difference in how they are trending in that case. If, however, the series obeyed y = 100, that is trending very differently from y = 2x + 100 despite the common point. Even if it is y = 150 intersects it. I would want the test to tell me that these two series are clearly trending differently. Is there a nonparametric test for this? I was thinking of something like K-Smirnov, but I cannot seem to figure out how to apply it in Stata to the case at hand. The alternative I've been doing is a linear regression of y on x for each series and then testing whether the coefficients are statistically different from one another just with the test command (storing the coefficients and then testing). Thank you very much for reading and for any help!


title: How to associate one predictor and 10 responses?

text: I have to test one predictor associate with 10 responses, the predictor is nominal, and response are continuous, so what kind of test can be used to test the association between them?


title: Function with variable having gaussian distribution

text: If I have a variable $X$ whose Gaussian distribution is known and let $f$ be a known function, is there a way to compute $f(X)$ _i.e._ the resulting Gaussian distribution from this? Is the result actually a Gaussian distribution?


title: Type of probability used by insurance company (and influence of gender)

text: I am currently working on analyzing a case and I need some help with the statistics. When a certain woman tried to find an insurance for her car, she was told that the premium she had to pay would be more than average, based on the fact that she was a woman. According to the insurance company, on average, women were a higher risk. She won the court case: the court decided that her individual risk could not be based on averages of other women. Now my questions are: (a) What kind of probability is used by the insurance company? (b) Can you comment on the kind of probability that is suggested by the verdict of the court? I am really interested in the difference of the type of probability used by the insurance company and the court.


title: Which distribution to use with MCMC and empirical data?

text: In class, we've been learning a myriad of really interesting techniques to sample from a given distribution, filter online data, particle filters, etc. My issue is that when I take some real-world data and plot it, the distribution is clearly not Gaussian. So, I need to estimate some distribution. Or, in the case of an online filter (particle, etc.) I need to estimate some form of transition kernel. How do people normally do this? What would be considered "best practices" for developing some distribution to fit empirical data? What are some reliable "goodness of fit" tests?


title: Can I derive the PDF ( f(C) ) from a spatial function defined as C(x)?

text: Say I have a random variable (such as concentration) that is defined spatially by some function $C(x)$. I would like to derive a PDF $(f(C))$ from the concentration function $C(x)$ over some domain of x, so that if I pick a random value of $x$, $f(C)$ describes the probability of $C$. Here's an example (from http://www.cs.ucr.edu/~ciardo/teaching/CS177/section7.1.pdf). The line segment joining two points on a circle has length $Y$ defined by $Y = 2r \sin{ (Θ/2)}$ if $Θ$ is Uniform$(0,2π)$, the pdf of Θ is $f (θ) = 1/2π$. The expected (mean) value of Y can be defined as: $E[Y] = \int^{2π}_0 Y(θ) f(θ) dθ$ How could i define the probability density function of $Y$, $f(Y)$? (without assuming that f(Y) has a normal or some other distribution).


title: Dimensionality reduction with locality-sensitive hashing

text: I have a conceptual understanding of locality-sensitive hashing (LSH) in relation to near-neighbour search. However, the articles that I have read so far seem to gloss over the details of how LSH can be used for probabilistic dimensionality reduction (e.g. Wikipedia). **Questions**   1. How does LSH reduce the dimensionality of data?   2. Can I specify the dimensionality of the subspace or is this predetermined?


title: pdf of Beta distribution with parameters (a+1, b)

text: Let $A$ be a Beta distributed random variable, $\mathcal B(a,b)$, with probability density function: $$ f(x)= \frac{Γ(a+b)}{Γ(a)Γ(b)}*x^{(a-1)}*(1-x)^{(b-1)}, $$ where $0<x<1$. Suppose that $B$ is a random variable, $\mathcal B(a+1,b)$. Write down its probability density function $g(x)$.


title: ets() in R returns a flat forecast for intermittent demand

text: In my attempt to forecast sales demand by month utilizing the last 3 years of history to predict balance of the year, `ets() from forecast() package` yields an answer for some time series but not others. I'm following Dr Hyndman's presentation on best fist ets() function. I understand while the flat forecast (average of the intermittent/volatile inputs) might be the best based on the intermittent history. My question is more of whether the process below has any validity or if it's tottay wonky. For the ones where ets() provides a flat forecast I tried the procedure outlined below to fix it, please let me know if I'm totally lost in my way of thinking:   1. ETS() flat forecast when there are too many 0's (or volatile)   2. I increase every observation of the time series by 100   3. ETS() forecast a curve   4. Can I just now reduce everything by 100 to get my forecast?


title: hierarchical Bayesian models vs. empirical Bayes

text: Would you consider the HBM vs EB to be two alternatives in which the hyperparameters are "in the game" of being sampled/estimated/etc.? There is clearly a connection between these two. Would you consider HBM more "fully Bayesian" than EB? Is there some place where I can see what are the differences between being "fully Bayesian" and the other alternatives? Thanks.


title: How to I compare the growth differences between two time based datasets?

text: (See note**) I am attempting to put together an infograph and am almost done collecting all the data I need. In addition to the base raw data collected, I'd like to add something I'm calling parallel stability and differential volatility. Here's what I mean:   * First, take a look at these graphs: ![enter image description here](http://i.stack.imgur.com/hJTN4.jpg)   * While interesting, I'd like to add data for growth/recession realized between the two parties for a trailing 14 days.   * Also like to be able to add the "parallel stability and differential volatility" measurement, which basically says if both party experience the same growth/recession during the same time period they have "parallel stability". On the other hand, if one party grows, and the other experiences a recess in growth, that is differential volatility. I'm looking to compare YouTube-to-Youtube, Searches-to-Searches, News-to-News data between the two parties, not across the datasets. If it's of use I can attempt to post the data, but it's basically one data point per day for the past 365 days per dataset (YouTube, Searches, News) per person with Y-values between 1-100; meaning you could just as easily create a dummy dataset based on how you describe how to find the solution. ** Let me clearly state that my understanding of math, stats, etc - is very limited; so assume that I know almost nothing and use plain-English where possible.


title: Generating random nos based on 'k' moments

text: How do I generate random nos based on say _k_ moments? (no other constraints on support) When _k = 2_ , we generate random nos. from a normal distribution defined by the 2 moments. Can we generalize this idea to an arbitrary _k_. A related question: Is normal distribution the only (known) distribution which is defined by just 1st 2 moments? I wouldn't mind suggestions for books/links, instead of direct answers to these questions. **(Editing the question based on a couple of answers I got)** I'll try to provide the context in which I thought of this question. I have a uni variate sample and no more information about anything. Now, I want to try and simulate nos. from the distribution from which this sample was drawn. I could just find its 1st 2 moments, assume this sample to be normally distributed, and be done. But if I want do away with 'Normal' assumption, I can try to find it's first _k_ moments and somehow draw from the distribution defined by these _k_ moments. I am interested in its theoretical underpinnings as well as practical implementation.


title: Repeated measures ANOVA for patient counselling data

text: I am trying to check the significance of counselling in the technique for inhaler usage (rotahaler only) at different points of time in same subjects. The subjects' baseline technique will be assessed using a checklist scoring : zero for incorrect/ missed step and 1 for correct step. The total score is 8 (i.e., the checklist has 8 steps). Then subject will be counselled and an immediate post intervention score is taken followed by another assessment after 2 weeks. I want to compare the means at 3 points in time. Can I use a repeated measures ANOVA here? If not, which test is suitable?


title: Markov Switching Forecast. How can I derive this?

text: Consider the autoregressive model, $\left[ \begin{array}{l} y^{\ast}_t\\ x_t^{\ast} \end{array} \right] = \left[ \begin{array}{l} a_{11}\\ a_{21} \end{array} \begin{array}{l} a_{12}\\ a_{22} \end{array} \right] \left[ \begin{array}{l} y^{\ast}_{t - 1}\\ x^{\ast}_{t - 1} \end{array} \right] + \text{} \left[ \begin{array}{l} \varepsilon_t\\ \upsilon_t \end{array} \right],$ where $\{ \varepsilon_t \}$ and $\{ v_t \}$ are white-noise processes with zero mean, $y_t^{\ast}$ and $x_t^{\ast}$ are given by $\begin{array}{lll} y^{\ast}_t & = & y_t - \alpha_1 - \alpha_2 S_t,\\ x^{\ast}_t & = & x_t - \alpha_3 - \alpha_4 S_t, \end{array}$ and $\{ S_t \}$ follows a two-state Markov process with transition probabilities $\begin{array}{lll} p & = & P ( S_t = 1 | S_{t - 1} = 1),\\ q & = & P ( S_t = 0 | S_{t - 1} = 0 ) . \end{array}$ Derive the expected value of $y_{t + n}$ conditional on information available at time $t$ about the current and past values of $( y_t, x_t)$ and the current value of $S_t$, i.e., $E ( y_{t + n} |_{} y_t, y_{t - 1}, \ldots ., y_1, x_t, x_{t - 1}, \ldots ., x_1, S_t)$.


title: How to integrate cdf(ppf(x)-A) for standard normal ppf and cdf

text: I'm looking for a way of integrating the following formula where _ppf()_ is the percentile point function for the standard normal distribution, _cdf()_ is its inverse, and _A_ is a constant: \begin{equation} \int_{0}^{1} cdf(ppf(x)-A)dx \end{equation} I can do it with a monte carlo technique but I'm hoping there's a faster way! Thanks in advance if anyone can help...


title: Is there autocorrelation in volatile series?

text: I read in many books or notes online that (1)volatile series do not differ significantly from white noise and that (2)their squared values will exhibit correlation. Although I agree with (2), I can't seem to agree with (1). I used the tsdisplay function from forecast package to plot the acf of SP500(MASS) and bmw(evir) data and found that there are a few significant autocorrelation. This isn't so obvious when plotting the acf using the built in acf function as the correlation at lag 0 mask the correlation at subsequent lags. I also generated a white noise series and although the acf seem the same there were not really any significant autocorrelation. Can anyone explain this to me? Thanks


title: How does logistic regression use the binomial distribution?

text: I'm trying to understand how logistic regression uses the binomial distribution. Let's say I'm studying nest success in birds. The probability of a nest being successful is 0.6. Using the binomial distribution, I can calculate the probability of r successes given n trials (numbers of nests studied). But how is the binomial distribution used in a modelling context? Let's say I want to know how mean daily temperature affects nest success and I use logistic regression to explore this question. Within the context I've described, how does the logistic regression used the binomial distribution? I'm looking for an intuitive answer, therefore an answer without equations! I think equations are only useful once understanding has been achieved at an intuitive level.


title: Invalid parent values, Adaptation incomplete, and print problem in Jags

text: I am trying to fit a GLMM for binary responses and I want to model the covariance matrix of random effects, say bigE, with some covariates. To do this, we decomposed the bigE into two different matrices, bigT and bigD, as seen in the following code. When I run the following code, it gives me the following error. Error in node (inverse(bigT[1:5,1:5,100,2])) Invalid parent values I couldn't figure out my errors, although I spent days on it. Furthermore, when I delete the bigT and gamma, and assume that bigE=bigD, this time it gives me that "adapation is incomplete". At this point, I need your help to see the problem with bigT and MCMC input arguments, such as n.adapt and so on. My last question is also, is there any "cat" kind of function in Jags. Since, if there is, it will be easy to see what's going on every line by printing them. I would be very much thankful if you could give me some suggestions for this problem. mymodel<-function(){ for ( k in 1:K) { for ( i in 1:n) { for ( t in 1:nt) { Y[t,i,k] ~ dbern(cp[t,i,k]) cp[t,i,k]<-min(max(cp1[t,i,k],0.0000001),0.9999999) cp1[t,i,k]<- 1/(1+exp(-eta[t,i,k])) eta[t,i,k]<-inprod(X[t,,i,k],beta[])+ b[t,i,k] }#t b[1:nt,i,k]~dmnorm(mu[1:nt], InvbigE[,,i,k] ) InvbigE[1:nt,1:nt,i,k]<-inverse(bigE[,,i,k]) bigE[1:nt,1:nt,i,k]<-inverse(bigT[1:nt,1:nt,i,k])% _%bigD[1:nt,1:nt,i,k]%_ %inverse(t(bigT[1:nt,1:nt,i,k])) } # i } # k for ( t in 1:nt) {mu[t]<-0 } for ( k in 1:K) { for ( i in 1:n) { for ( t in 1:nt) { for ( tt in 1:nt) { bigT[t,tt,i,k]<-equals(t,tt)*1 + equals(t,tt+1)*inprod(-Xg[t,,i,k],gamma[]) bigD[t,tt,i,k]<-equals(t,tt)*exp(inprod( Xl[t,,i,k],lambda[])) } # tt } # t } # i } # k ## Priors for ( p in 1:np) {beta[p]~dnorm(0,0.00001) } for ( g in 1:ng) {gamma[g]~dnorm(0,0.00001) } for ( l in 1:nl) {lambda[l]~dnorm(0,0.00001) } } #mymodel


title: Multiple comparisons with ANOVA including one between and one within-subject effect

text: I'm trying to look at multiple comparisons (across levels of my between subjects factor) for a model with one between-subject factor and one within- subject factor. I was trying to use TukeyHSD on an aovlist object, and then found out that I couldn't do that. After some browsing, I found how to do what I wanted using `lme()` and `glht()`:               options(contrasts=c("contr.sum", "contr.poly"))     lme_model=lme(dv ~ between*within, data frame, random=~1|ID,                   correlation=corCompSymm(form=~1|ID))     anova(lme_model, type="marginal")     summary(glht(lme_model, linfct=mcp(between="Tukey")))      It seems to work (the `anova()` results seem to jive with what I was getting from `aov()`), but I guess the problem is that I don't really understand what `lme()` is, or what it's doing. Does anyone have any documentation that I can look at?


title: Sum of very low probability

text: I have a score of some feature, $F_1$ and $F_2$ where this score is the logarithm of probability. This score is very low and very sparse, for example i have: $F_1$ with score $-800$ $F_2$ with score $-3000$ I need to get the probability with the exponential function and make the sum of probability like this: $e^{-800} + e^{-3000}$ But if I make the exponential of this number i get 0 for the underflow problem. With log-sum-exp trick i can do: $\log(e^{-800}+e^{-3000})$ $= \log(e^{-800+800}+e^{-3000+800})-800$ $= \log(e^{0}+e^{-2200})-800$ But now i have $e^{-2200}$ that is too small. What's wrong? how I can make my $e^{-800} + e^{-3000}$ sum ?


title: How to specify/restrict the sign of coefficients in a GLM or similar model in R

text: **The situation:** I'm struggling with a predictive analysis of food sales prices using a generalized linear model. My dataset contains different kinds of food (cheeses, vegetables, meats, spices etc.) and hence I am splitting the dataset completely by these kinds when doing the analysis, because they are very different by nature. **The current model:** The dataset/model contains both factors such as "country of production" and numeric variables such as "transport distance" which is all used in the gamma based GLM i R. **The problem:** Now in general my model fits pretty well, however sometimes in rare cases some of the metric variables gets the opposite sign (+/-) than you would expect it to have, because the model somehow catches other influences. **An example:** An example would be spices. All spices have a relative long "transport distance" and a relative long shelf life and hence a pretty small impact on the sales price compared to e.g. meat. So in this case the model might by accident end up giving the "transport distance" variable a small but negative value - which is of cause wrong because it would mean that the longer the distance the food was transported the lower the price would be. **My question:** What kind of model should I use in R if I wan't something similar to a GLM model but I want to be able to specify restrictions on some of the variables/coefficients? E.g. if I want to say that an increased "transport distance" should ALWAYS have a positive impact on the sales price? **Ideas:** I have heard something about both "Bayesian GLM" models or using a so called "prior distribution" but I have no idea which one, if any, would be the best to use..? **UPDATE** The answer below by @ACD is not, exactly what I'm looking for. I don't need an explanation of WHY this occurs, I need a solution to restricting the coefficient signs :-)


title: Why does correlation come out the same on raw data and z-scored (standardized) data?

text: I just discovered by mistake that raw data and the same z scored data produces the same correlation. Why is this? Can someone walk me through the logic?               > x = c(2,4,6,7,8,3,6,8,9)     > y = c(3,3,5,7,6,8,9,10,5)     > cor(x,y)     [1] 0.3780247     > zx = (x-mean(x))/sd(x)     > zy = (y-mean(y))/sd(y)     > cor(zx,zy)     [1] 0.3780247


title: Linear model fit seems off in R

text: I'm posting this in this forum because I think the question is related to stats. If it's not, I'm fine with it being moved over to Stackoverflow. I have data in R and am trying to fit a linear model to it. Here's what the data looks like (sorry that it's not reproducible, it's just too much data to type out). The colored dots are based on density of points. The black line is what a linear model is returning as the best fit line (via `lm(y~x)` which spits out a slope of `0.67` and an intercept of `0.002`) and the dashed line is what I would think the best fit line should be (slope of `5`, intercept of `-3`). Why is R's `lm` method giving me a line which looks like it doesn't fit at all. Is it true that that line is better fitting than the dashed line I propose? ![enter image description here](http://i.stack.imgur.com/uRUu9.png)


title: How to compute whether the mean of a non-random sub-sample is statistically different from the mean of the sample?

text: I have a variable called "obs" and from this variable I generated a new variable called "obs_sub" by excluding all observations for which a dummy variable is equal to one. Now what I want to know is if I remove these observations whether the mean of the sub-sample is equal to the mean of the sample. This is my null hypothesis. I could run a regression of $E[obs|d] = \alpha + \beta \cdot d$, but the regression with the dummy $(d)$ will only tell me whether the coefficient on the dummy is statistically significant from zero. It will not tell me whether the data are unlikely given $E[obs|d=1] = E[obs]$ or is the null hypothesis $\beta=0$ the same as the null hypothesis as $E[obs|d=1] = E[obs]$? Cheers, Martin


title: Mixed Models : Effect of Adding Individual Level Predictor on Between Group Residual Variance

text: This is a portion of a tutorial from the Centre for Multi-Level Modeling regarding random intercept models. Can anyone explain why this scenario is so? It seems like by the logic given and that the u_j (group effects) are moved toward zero by adding the new predictor, the variance of the u_j would DECREASE not increase. ![enter image description here](http://i.stack.imgur.com/0iocN.jpg) ![enter image description here](http://i.stack.imgur.com/s0qXz.jpg)


title: Graphical Multiple Linear Regression in Stata

text: Consider this standard linear regression model: $Y = \beta_0+\beta_1X_1+\cdots+\beta_pX_p+\epsilon$ I've loaded such a dataset into Stata 12.0, so I have some variables $y,x_1,x_2,\dots,x_p$. How do I produce a plot, like I would with `. scatter y x` for a simple linear regression model?


title: How to select the best classification scheme in survival analysis (SurvivalROC, R2, Concordance, AIC)?

text: I'm studying several schemes on classifying patients about their survival time. Let me illustrate the problem with supposing I have just two schemes. Let's suppose that Scheme 1 put the patients in 5 groups, and Scheme 2 put the patients in 5 groups too (but the groups have different compositions, for example patient 1 in Scheme 1 may be in the group 3 and in Scheme 2 in the group 4). So, it is my problem but with several schemes. KM plots are useful to analysis, but it is not practical, and some plots are very similar in visual terms. The logrank test is not the case I think, because I think it is just useful to verify if I really have any difference in the survival experience between the curves and no more. In really, what I want is to verify which scheme give me a better prediction in some sense of accuracy, in really a measure for that. One way I think was to just use the survivalROC package and apply the survivalROC function for times 60, 120 and 180 for example, and whatch which scheme give the best AUC. Another way I thought is to use these schemes with cox proportional hazard models, is it really an alternative? Using the coxph function I can have R2 and concordance in its return, and I can compute the AIC with the extractAIC function too and choose the best scheme using the less AIC. In the AIC case I have some doubts beacuse I thought that AIC was good a time ago, but when was I reading recently about it I just found examples with nested models for COX PH, is it applicable to non-nested models too? So, in summary I thought about 4 options: 1. survivalROC and use AUC 2. coxph and use R2 3. coxph and use concordance 4. coxph and use extractAIC I want to automate my analysis using one numerical measure to aumotatic selection of the best scheme (remembering the non-nested models nature). Is there any problem in this approach I am thinking? Which could be the best measure for that? Thanks, Felipe Prata.


title: Logistic regression and discrepant sample sizes between 0/1 groups

text: I currently working in a multivariate logistic model but I have a problem regarding the sample size of my observations: -The "success" (1) event group has a sample size of 249 distinct observations - The "non success" (0) event group has a sample size of 48,957, and it's a significant part of the population, and many times larger than the "success" group. So when I fit a multivariate logistic regression model, two things happen: - Because of the big sample size, The p-values of the model coefficients become always significant at alpha 1% - The fitted predicted probability variation between the groups becomes very tiny, even if the independent variables are good predictors. So I was suggested making a bootstrap of the model. So here is my doubt : -Should I do it the more usual way, taking smaller arbitrary size random samples of the complete sample (both groups) with replacement ( or without replacement in this case?) OR -Should I keep the the small "success" group constant and them add an equal number of different "non success" randomly picked cases in each sample. Is this "cheating"? OR -Should I just stay with the original whole sample model, and bear with it's little predictive power (less than 1% fitted predicted difference between the groups ), even when using significant predictors that differentiate the groups. In both cases this is not exactly the usual bootstrap as I wish to make sub- samples of a larger group, the original big sample, and perhaps what I want is not bootstrap at all. My idea is to minimize the discrepancy between the two groups sizes, the inflated p-values, and the small probability fitted difference between groups due to tho fact that the "success" event is rare in my sample. Is this theoretically correct? The final objective is to obtain an "avarage" model with the mean coefficients and use it to calculate the probability of the "non success" cases actualy being "successes" Or does anybody has other idea for this question? Feel free to post any suggestion in R language. I know that there is a similar question but I want to know also if my ideas about this are wrong.


title: Spearman's Correlation using Stata

text: I am trying to use Stata for the first time to calculate Spearman's Rank Correlation for differences in weight and hemoglobin between annual visits for a small dataset. Do I enter the differences with the negative sign directly or do I square the differences first?


title: PCA and component scores based on a mix of continuous, binary and categorical variables

text: My question is strongly related to this one: PCA and component scores based on a mix of continuous and binary variables. I will basically use the same code, but add a new nominal feature (x6) to the data set. I want to apply a PCA on a dataset consisting of continuous, binary and categorical variables.               # Generate synthetic dataset     set.seed(12345)     n <- 100     x1 <- rnorm(n)     x2 <- runif(n, -2, 2)     x3 <- x1 + x2 + rnorm(n)     x4 <- rbinom(n, 1, 0.5)     x5 <- rbinom(n, 1, 0.6)     x6 <- c(rep('A', 25), rep('B', 25), rep('C', 25), rep('D', 25))     data <- data.frame(x1, x2, x3, x4, x5, x6)          # Correlation matrix with appropriate coefficients     # Pearson product-moment: 2 continuous variables     # Point-biserial: 1 continuous and 1 binary variable     # Phi: 2 binary variables     # For testing purposes use hetcor function     library(polycor)     C <- as.matrix(hetcor(data=data))          # Run PCA     pca <- princomp(covmat=C)     L <- loadings(pca)      Now in order to calculate the component scores, it was suggested to multiply the data set with the loadings L, which works fine for numerical and binary variables, but not on categorical data. The following computation causes the categorical feature to be a vector of NA´s.               scores <- data * L      How can I obtain the scores for this feature? Do I have to split it up into dummy variables to make this work?


title: Trends in noisy data and the applicability of the 95% confidence interval

text: I am performing simulations while measuring a quantity A which depends on the parameter B. I make N independent measurements of A for given values of B. I can then calculate the mean to get an estimation of what the real value of A is. My question is about the error. I can calculate the 95% confidence interval as 1.96xstddev(A)/sqrt(N) for each value of B. What I often find however is that this is not really representative of the error. For example, if I have the following data:               B   A(B)  95%CI     1   2     10     2   4     10     3   6     10     4   8     10      Without knowledge of the error, one would conclude that there is a linear relationship between B and A. And in my experience, as I make more measurements, this relation doesn't change, despite the noise in the measurement. But if I show this result to someone else, they say "given those error bars (+-10 on each point), you can't say anything about the trend". Since confidence intervals assume a normal distribution, are my measurements somehow more-sharply peaked somehow? Or am I using the 95%CI wrongly? Is there a better way to show the confidence in my data and/or the trend? Any advice is greatly appreciated.


title: Application of Poisson distribution to image processing

text: I'm trying to write a program to detect water bubbles in heated oil. I've applied a canny edge detection filter to the image and the results look like the following: No bubbles: http://imageshack.us/photo/my-images/20/nobubbles.png/ Bubbles: http://imageshack.us/photo/my-images/35/swirlsg.png/ I'm trying to find a way to determine when the bubbles are present. Would it be appropriate to check if the white pixels are Poisson distributed?


title: R as an alternative to SAS for large data

text: I know that R is not particularly helpful for analysing large datasets given that R loads all the data in memory whereas something like SAS does sequential analysis. That said, there are packages like bigmemory that allows users to perform large data analysis (statistical analysis) more efficiently in R. I wanted to know, apart from all the theoretical information, has anyone used / is using R for analysing large datasets in an enterprise environment and what are the typical issues that could arise. By large datasets I am referring to datasets that are ~ 200 GB in size. Also, any thoughts on real-life examples of migrating from SAS to R in such use cases would be helpful.


title: Building a time series model using more than independent variables

text: I am working on a project, and I am totally new to statistics. I have sales data for last two years at week level, along with other variables like temperature, holiday (TRUE/FALSE), where holiday are nominal variables. I have to do forecasting for the next 52 weeks. I have the following questions:   1. Can I use time series regression model where sales would be dependent, and temperature and   holiday would be independent variables?   2. How to decide which independent variable would have more impact on the sales?    3. Can we do forecasting using nominal variables? Will dummy coding work?   4. Can we do it in R/SPSS? I would appreciate any kind of help. Thanks in advance.


title: Meaning of partial derivative in gradient descent

text: Sorry for a newbie question. I haven't quite understood the correct meaning of what it means to take the partial derivative of a cost function with respect to the parameters, say theta. Suppose these parameters include 2 arrays of dimensions axb and cxd. When we do stochastic gradient descent, then what will be the partial derivative of J wrt theta (let's call it x)? Mathematically this is a symbolic expression. So when I program SGD, what should be the term x. What happens if J is a scalar, and what happens if J is a matrix? If in any case it is a matrix, what will be its dimensions? Thanks in advance! Hope I can clarify this issue soon


title: Hierarchical clustering of repeated measures in a few locations

text: I have a water quality data ( **value** ) measured 10 times (every month - **data** ) on three **depths** (shallow, medium, deep) in five **location** (A, B, C, D, E): ![enter image description here](http://i.stack.imgur.com/9ANC6.png) I want to find in which location values are similar, etc. I don't care about changes in depth or time, but overall characteristics. I've thinking about using hierarchical clustering, but what is the best way to prepare this data for that method? For example - using only location's mean and sd, some properties could be lost.


title: How do you choose the order p and q for ARMA(p,q) process for modeling a time series?

text: How do you choose the order p and q for a ARMA(p,q) process for modeling a time series? Can that be told from ACF or PACF alone, just like for a AR(p) or MA(q) process?


title: Why is Gamma(0,0) equivalent to the Jeffreys prior

text: I'm trying to use some code that includes Gamma priors for Poisson (rate) and Exponential (rate) distributions. I want to make the priors noninformative. I read that using a Gamma(0,0) is equivalent to the Jeffreys prior, which is noninformative. Is this true? Can anyone suggest any references for this? I thought the shape and scale parameters of a Gamma distribution had to be positive. Thanks.


title: Estimating means of correlated distributions with long tails

text: Suppose I have a relatively large number of samples (~1k) drawn from a series (~40) of increasingly long-tailed distributions (going from approximately normal to approximately log-normal). I want to estimate the mean and its uncertainty for these distributions, which I do using jackknife resampling because the samples are highly correlated between adjacent distributions in the series. However, as the distributions become increasingly long-tailed, the estimated uncertainty of the mean becomes so large that the data are useless. I might think to apply robust estimators, but my review of the literature seems to indicate that these methods assume that the "outliers" are erroneous data, drawn from some other distribution than the one of interest. For this case, they are samples from the true distribution that appear to be outlying merely because the region in which they lie is sparsely sampled.   * Are the usual robust estimators valid for this case? If so, how should they be applied if an outlier in one distribution does not correspond to outliers in other distributions?   * If not, are there other suitable methods? I experimented with power transformations, but could not work out how they could be used while retaining the correlations, since the transformation parameters would differ between distributions.


title: Standardized beta for logistic regressionin R

text: For my survey data analysis, I ran an Ordinal Logistic regression using the 'polr' function. The summary of the regression is as follows: ![enter image description here](http://i.stack.imgur.com/csKGq.png) My question is:   1. Do I need to standardize my beta values?   2. If so, is lm.beta the right approach (as per my understanding, it only works for linear models)? And if not, could you please provide a method to do so. Thanks everyone!


title: Gaussian Process Regression with additional Basis Functions

text: I'm working through Rasmussen's Gaussian Processes book, and I have a question about the possibility of optimizing additional basis function hyperparameters (in section 2.7 http://www.gaussianprocess.org/gpml/chapters/RW2.pdf). The text explains that the hyperparameters can be varied to maximize the log marginal likelihood, which is given in eq. 2.44 and 2.45. The derivative for the general case is given in eq. 5.9. I've looked in the gpml code that goes with the book and relevant literature, but I haven't been able to find the derivative wrt the hyperparameters including the basis function hyperparameters. In other words, for anyone familiar with Gaussian Process Regression, what is the derivative of eq. 2.44 wrt the hyperprior? (again, GPR Book, Ch. 2).


title: Naive forecasting

text: I used to know that Naive forecast is equal ( ft is the same like the previous year) Someone told me that there is another equation used in the sales which is equal ( current year - previous year)/ previous year Is it true?


title: mixed effects model output

text: lets say we have this: > model2 <- lmer(milk.amount~(1|cow), data=milk, REML=FALSE) > > model1 <- lmer(milk.amount~(1|cow), data=milk) > >> summary(model2) >> >> Linear mixed model fit by maximum likelihood ['lmerMod'] >> >> Formula: milk.amount ~ (1 | cow) >> >> Data: milk >> >> AIC BIC logLik deviance df.resid >> >> 186.5 191.6 -90.2 180.5 37 >> >> Scaled residuals: >> >> Min 1Q Median 3Q Max >> >> -2.0244 -0.4104 0.1795 0.6621 1.3879 >> >> Random effects: >> >> Groups Name Variance Std.Dev. >> >> cow (Intercept) 6.755 2.599 >> >> Residual 2.999 1.732 >> >> Number of obs: 40, groups: cow, 10 >> >> Fixed effects: >> >> Estimate Std. Error t value >> >> (Intercept) 27.0150 0.8663 31.18 then: > > summary(model1) >> >> Linear mixed model fit by REML ['lmerMod'] >> >> Formula: milk.amount ~ (1 | cow) >> >> Data: milk >> >> REML criterion at convergence: 178.9 >> >> Scaled residuals: >> >> Min 1Q Median 3Q Max >> >> -1.9981 -0.4136 0.1775 0.6561 1.4021 >> >> Random effects: >> >> Groups Name Variance Std.Dev. >> >> cow (Intercept) 7.589 2.755 >> >> Residual 3.000 1.732 >> >> Number of obs: 40, groups: cow, 10 >> >> Fixed effects: >> >> Estimate Std. Error t value >> >> (Intercept) 27.0150 0.9132 29.58 question is why model1 (with REML) doesn't show AIC, BIC, logLik, deviance coefficients? is it possible due to some kind of software dependecy? thanks in advance for reply