title: What is the qualitative difference between a Michaelis-Menten model and a log-linear model?

text: In a book chapter I am reading1 the authors compare the fit of the relationship between biodiversity and ecosystem function using both a log- linear model and a Michaelis-Menten model. The authors refer to the Michaelis- Menten model as "saturating" and the log-linear model as non-saturating. > ...authors of the publications used in [Meta-analysis 1] did not distinguish > between non-linear curves that do (Michaelis-Menten) and do not saturate > (log-linear). In [Meta-analysis 2] this distinction could be made and showed > that the average R2 value was 0.690 and 0.682 for the saturating Michaelis- > Menten and log-linear relationship, respectively. My initial investigation into these models (via Wikipedia etc...) appears to show very similar behavior in that they both approach an asymptote with increasing X-values. Furthermore, I am not sure what the authors are referring to with respect to the saturating and non-saturating descriptions because both models appear to saturate. **My question is, how should I interpret the relationship differently if it fits one of these models versus the other?** NOTE: I have looked up the models and while I can somewhat understand each, I lack the training to relate them mathematically. However, I am looking for a non-mathematical comparison of what fitting one or the other would indicate about the relationship they describe. * * *   1. Schmidt, B. et al. 2008. "Consequences of species loss for ecosystem functioning: meta-analyses of data from biodiversity experiments" in Naeem et al. eds. 2008. Biodiversity, Ecosystem Functioning & Human Wellbeing: an ecological and economic perspective. Oxford University Press.


title: How to do Linear Discriminant Analysis (LDA) with only theoretical distributions

text: Keeping things in short. I want to distinguish two classes on a feature space. All what I know are:   1. the theoretical distributions of the class 0 data on each feature dimension.   2. the expected distance between the mean of the class 1 data and the mean of the class 0 data on each dimension.  but I do not have any labeled data. I wonder whether LDA can help in this case and how to do it. Previously, I tried hypothesis testing, but it is very difficult to decide the significance level on each dimension: when I plotted the distribution of accepted data, their sample distributions are far away from the theoretical, often due to the fact that the variance of each class is of a similar scale of the difference between the means of two classes. PS: all features are correlated.


title: Random variables

text: I've just come across a few questions on the net which I have no idea how to answer, can anyone point me towards the right direction?   1. Random variables $X$ and $Y$ are such that $X$ has mean 1 and variance 4, $Y$ has mean 2 and variance 9, and $\text{Corr}(X, Y ) = 1/3$ What is the variance of $3X − 2Y + 1$?   2. In question 1 above, what is the covariance between $X + 2Y$ and $X − Y$ ?   3. In question 1 above, if $Z$ is another random variable satisfying $E[3X −2Y +Z] = 0$, what does the mean of the random variable $Z$ equal?


title: Significance of an impulse response function

text: I've read several paper that all compare different cumulative IRF of the same VAR equation for statistically significant difference. The IRF they use are simply the sum of the coefficients of the VMA representation resulting from Cholesky decomposition. I wonder how one can construct tests for IRF. The authors mention that they use Newey-West standard errors. Could it be that they don't just take the sum of the coefficients but run two different simulations, one with and one without a shock? I still don't know how to apply NW to it and how to derive significance from that though. Here you can find one of the paper (p.16). Thanks for your hints!


title: Given a frequency distribution, calculate probability of A given B

text: My test recorded the number of times a variable reached a threshold value when a predetermined event occurred. Then, the test recorded the number of times the variable surpassed the threshold value by 1 level when the same event occurred. The test progressed in this manner with the highest observation at level 22. The resulting frequency distribution follows (sample size was 9,970 observations):   * 0:3468 <= the threshold level   * 1:2163 <= level 1   * 2:1408   * 3:1005   * 4:641   * 5:487   * 6:266   * 7:196   * 8:135   * 9:71   * 10:49   * 11:30   * 12:16   * 13:6   * 14:10   * 15:8   * 16:3   * 17:4   * 18:0   * 19:1   * 20:1   * 21:1   * 22:1 <= level 22 Note that the frequency of Level 1 observations do not include the frequency of Level 0 (the threshold) observations. Finally, in order for the variable to reach subsequently higher levels, it must pass through lower levels. During runtime, given that the variable has reached the threshold level(0), what's the probability that the variable will continue to level 1? to level 10?


title: Linear regression: b1 has the minimum variance among all unbiased linear estimators of beta1

text: There is a proof provided in _Applied Linear Regression Models (1983)_ by **Kutner et al.** (Page 64), which is quite clear and easy to understand, except one point, namely, it assumes that $\sum k_i d_i = 0$, from the "restrictions on $k_i$ and $c_i$", without expounding on what restrictions these are. Could anyone please explain? Thanks! ![enter image description here](http://i.stack.imgur.com/ffJVS.png)


title: How to combine results of logistic regression and random forest?

text: I am new to machine learning. I applied logistic regression and random forest on a same dataset. So I get variable importance (absolute coefficient for logistic regression and variable importance for random forest). I am thinking to combine the two to get a final variable importance. Can anyone share his/her experience? I've checked bagging, boosting, ensemble modeling, but they are not what I need. They are more of combining information for the same model across replicates. What I am looking for is to combine result of multiple models.


title: How do I calculate prediction intervals for random forest predictions?

text: In regression problems random forests yield a prediction for each case by averaging the results of each single tree in the forest. I would like to draw a 95% prediction interval around each prediction. There is a nice option `interval="prediction"` for the prediction method of linear models, but there's no such thing for random forests (at least in the R package randomForest). Can I presume that the casewise predictions of the single trees are normally distributed and apply the formula http://mathurl.com/bbvuvx9 or should I determine it by bootstrapping? If so, how can this be done?


title: Using ANOVA to analyse a distribution

text: I found a method in an paper, which I do not quite understand. In this experiment the influence of a within subject factor was analysed. In order to analyse the influence of the factor not only on general parameters, such as mean or standard deviation, but on the whole distribution a method was used that I cannot find in the literature. For this the data was sorted into quartiles and the the average for each quartile was calculated. Then using these averages, a standard ANOVA was calculated using the within-subject factor as well as the quartile number was performed. The only article I could find, which includes something similar is this one. However, as I am not really familiar with fmrt and similar techniques, I do not understand most of this article and cannot be sure if it actually describes the method I am looking for. Is this method common for analysing data distributions? If so, where can I find some references explaining this method and how to interpret the results?


title: how do I interpret the following hypothesis test?

text: Let's say I have two hypotheses for a coin with probability $p$ for heads: $H_0$ - the null hypothesis - the coin is fair $p = 0.5$. $H_1$ - coin is unfair $p \neq 0.5$. Say the test is $|X-n/2| > r$ for some $r$ where $n$ is the number of tosses. It is possible to choose $r$ for a desired significance level quite easily (theoretically, all the details are there for computation). But what if I want to design a test with a certain type II error? It seems odd to me that one cannot compute $p(\textit{accept } H_0 | H_1)$ because $H_1$ doesn't specify what is the value of $p$. When I looked around, I noticed that for this simple example, it is requested to calculate type II error for a certain value of $p \neq 0.5$ (such as $p = 0.7$). Is there a meaning to type II error when $p$ is unknown? Also, how is that related to the idea of "frequentist repeated experiments?" Thanks.


title: Standard deviation of a particular dimension in a multivariate Gaussian distribution

text: I have a set (cluster) of vectors in dimension d. From this I have calculated the sample mean and covariance matrix ( I make the assumption that they are from a multivariate Gaussian). My question is, given a new vector (in dimension d) I am trying to decide if it belongs to this cluster by checking if the distance from the mean is less than 2 standard deviations. In the one dimensional case I would simply check if x-x_bar > 2*sigma. How does this extend to the multivariate case? Thanks


title: PYMC Confusion: are observed nodes fixed or stochastic?

text: I've been trying to gain a better understanding of factor potentials in PYMC. In reading this article by Cam Davidson-Pilon on Yhat, I got confused about how observed nodes are understood by PYMC. Are they fixed in PYMC when `observed=True`? Or are they random like a stochastic node? The factor potential is specified as:               @mc.potential     def censorfactor(obs=obs):          if np.any((obs + birth < 10)[lifetime_.mask] ):              return -100000         else:             return 0      The variable birth and lifetime_ are fixed numpy arrays. Hence, it seems, the only changeable variable is `obs`. But I thought that was fixed in PYMC because it is defined by:               obs = mc.Weibull( 'obs', alpha, beta, value = lifetime_, observed = True )      If obs is fixed then the potential function will return a log probability of -100000 every time. I suspect I'm missing a fundamental point of MCMC here. Appreciate if someone can help me understand. **UPDATE** I did some additional research and it appears that observed nodes can indeed be (partly) stochastic. If you pass a masked array as observed data to a stochastic node, each of the masked values will be individually represented as a stochastic node. Chris Fonnesbeck's tutorial from this year's Scipy conference covers this topic at the end of this IPython notebook. Here he describes how the posterior predictive distribution is used to estimate each of the missing data points: $$p(\tilde{y}|y) = \int p(\tilde{y}|\theta) f(\theta|y) d\theta$$ The ability to impute missing data conditional on the data **and** the model parameters (simultaneously) a powerful technique. Plus we get traces, distributions, HDI etc for each missing value.


title: What is the probability of randomly selecting n random numbers in the range 1-m in sorted order?

text: For example say m is 10 and n is 3. I am selecting 3 numbers from 1-10 and want to know the probability with which they will be in sorted order. 1, 1, 4 would be fine. 1, 6, 3 would not be fine. For m=10 and n=2, the probability of selecting the 2 numbers in sorted order is 55/100. How could I calculate the probability for any arbitrary n and m?


title: R - Linear Regression - Control for a variable

text: I have a computer science background & I am trying to teach myself data science by solving the problems available on the internet I have a smallish data set which has 3 variables - race, gender and annual income. There are about 10,000 sample observations. I am trying to predict income from race & gender. I have divided the data into 2 parts - one for each gender & now I am trying to create 2 regression models. Is this possible in R? Can some one provide example syntax.


title: What type of predictive analysis should I use?

text: I have a pure math background with knowledge of basic statistics (random variables, inference, etc.) but am new to predictive modeling. Here is my situation: I have a bunch of independent variables that I am trying to use to predict a single dependent variable. The independent variables are categorical (ordinal), and the dependent variable is a metric. I'm not sure what areas I should be looking into/learning about. Multiple correspondence analysis, multiple regression, and factor analysis all seem to pop up in my searches. Which of these kinds of analysis would best serve my situation? I'd welcome any alternatives, but there's no need to cover all available options. Thanks!


title: Big Data state of art

text: What are the current tools used to analysis big data other than Hadoop and Hadoop related technologies like HIVE, PIG and etc?


title: Random Forest Regression - Coping with extreme values

text: I'm not sure if I used the concept "extreme values" right. Anyhow, I'm trying to produce a model that estimates maximum tree heights / $\text{km}^2$. I have a database of around 24000 points ($\text{km}^2$), each has the max tree height value and 33 predictors. After playing around with random forest I manage to achieve a correlation of 0.67 between the real height and the estimated height on the test sample (20%). A MSE of around 1.6 meters. But Maximum errors of up to 33 meters. What I can see is that patches with very tall trees or very short trees (50 meters - 1 meters) are out of the scope of the model. Thinking in linear regression it is analogous to losing prediction power as you move away from the center of gravity of the observations. Right? How can I cope with this if at all? p.s. this was implemented in R


title: Is there an equivalent to Kruskal Wallis one-way test for a two-way model?

text: If the model does not satisfy ANOVA assumptions (normality in particular), if one-way, Kruskal-Wallis non-parametric test is recommended. But, what if you have multiple factors?


title: Analysis of variance not statistically significant... but is there still a pattern to the data?

text: **Problem** I have a computer algorithm, and am attempting to measure the affect of a parameter (NS) on the output of the algorithm (SC). My hypothesis is that if the level of NS is lowered, there should be a corresponding decrease in SC as well. I'm doing all of my statistical analysis in the R language (see below). I've performed an analysis of variance on my data, and though the results aren't statistically significant, it looks like there is a pattern to the data (that lowering NS does decrease SC). I tried increasing my sample sizes, which does result in more pairwise comparisons becoming statistically significant. However when I do this, the results of Levene's Test show that my variances are no longer homoscedastic, and thus the analysis of variance isn't valid... If you look at a graph of the mean values of the output (sc) you can see what I'm talking about when I say "it looks like there is a pattern". I've generated different data sets several times and each time, the means of SC are scattered similarly to graph below (despite the pairwise comparisons being statistically insignificant). ![Mean SC values](http://i.stack.imgur.com/mx51N.png) **What I am Doing** Here is my R code and some output:               t <- read.table("output.dat")      names(t) <- c("sc", "ns")     leveneTest(t$sc, group=t$ns, center=median)          Levene's Test for Homogeneity of Variance (center = median)               Df F value Pr(>F)     group     24  1.0447 0.4018       124975                    Warning message:     In leveneTest.default(t$sc, group = t$ne, center = median) :     t$ne coerced to factor.          t.aov <- aov(t$sc ~ as.factor(t$ns))     summary(t.aov)                              Df  Sum Sq Mean Sq F value Pr(>F)         as.factor(t$ne)     24    1448   60.32   5.488 <2e-16 ***     Residuals       124975 1373548   10.99                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           TukeyHSD(t.aov)  #This prints out a HUGE table which I'm not going to include                      #The point is all but a few of the comparisons aren't statistically                               #significant at the p = 0.05 level      plot(TukeyHSD(t.aov))      ![Plot of confidence intervals](http://i.stack.imgur.com/WxWqI.png) **So...** at this point I'd be ready to just say that adjusting the level of NS doesn't affect SC (for all but the 1-24, 2-24 etc, comparisons, its hard to see in the graph, but its there) and thus lowering the level of NS doesn't result in a minimization of SC. However I just can't reconcile the graphs of the means of SC with the statistical implications from the analysis of variance... is my intuition leading my astray here, should I simply reject my hypothesis? Is there a way I can still increase my sample size to get a significant result, even though the Levene Test says my data is no longer homoscedastic? Should I be using a different statistical tool instead of analysis of variance to decide these things? Any suggestions, advice, or criticism is appreciated. **P.S** I'm not a statistician, so if I'm doing something glaringly stupid, please let me know.


title: A doubt on the definition of p-value

text: The p-value is the probability, under the assumption of the null hypothesis $H_0$, of obtaining a result equal to or more extreme than what was observed at given data. This means, if I were to repeat my experiment an infinite number of times, the proportion of times that my observed data would give results equal or more extreme, under $H_0$, would be equal to the p-value. My doubt is the meaning of the «under $H_0$» or «under the assumption of the null hypothesis $H_0$» parts. What do they mean? Only if $H_0$ is true, will the p-value magnitude have some importance in assessing whether the results obtained at this given data study provide evidence against or for $H_0$? If I were to get a low p-value, while $H_0$ was true, what could I conclude? \footnote{To me this seems possible to happen. If we simulate a linear DGP, with the usual assumptions to make sure OLS is BLUE, and if the error terms had great variance, then at an experiment I could some estimates completely different from the DGP, making the p-value very low, no?} What would mean to calculate $P(T(\mathbf{X})\geq T_0)$ when $\theta\in \Theta_1$? Any help would be appreciated.


title: Is it OK to use the CLT to create a normal distribution where there is none?

text: I have some data that looks like this: ![enter image description here](http://i.stack.imgur.com/Cu0GO.png) Procrastinator has come up with one good suggestion for how to test hypotheses under this distribution, but it relies on some guesswork to fit constants. I would therefore feel more comfortable if I had multiple methods, and could check that they all agreed. One thing I can do is find the mean of a bunch of samples, and by CLT this will approximate the normal distribution. In fact, this happens at reasonably small sample sizes. Here is my data with 1000 samples (with replacement) of 100 points each: ![histogram](http://i.stack.imgur.com/GCxya.png) It looks fairly normal: ![qq](http://i.stack.imgur.com/5MIBV.png) So would it be acceptable for me to compare two groups by first transforming them into a set of samples, and then using e.g. a T-test? Or am I introducing an unknown bias here? (In this modified form, comparing their means does seem to be a useful metric, so I think a t test is good from that perspective.)


title: Sample size determination for block design with repeated measurement in SAS

text: I have a design where birds of 1, 3, 5 week of age are placed in 2 chambers(independent) with 4 treatments( light technology) for 4 days. for example, x no of birds of 1 weeks of age are kept for 4 days and then new birds of 3 week age are placed in both chambers and so on...readings are taken for their position, feed consumption/day/bird etc. I need to calculate the sample size required. I am using Proc GLMpower in SAS, using twowayanova to calculate the sample size but I am not sure how to account for repeated measurement for 4 days. Also, Can you suggest any other way to do this type power analysis. I tried in G*Power as well but not quite sure how it works. Thank you.


title: How to characterize symmetric discrete distribution?

text: I have a very basic question on when a discrete distribution might be called a _symmetric_ distribution. Let say I have a r.v. $X$ that can take two possible values $(x1, x2)$ with $x1 \neq x2$ and corresponding probabilities $(0.4, 0.6)$. Then can I say that $X$ is symmetric? Thanks,


title: Pre-whitening of AR(2) time series

text: I'm working on a time-series problem. Consider the model below: $$Y = \beta_0 + \beta_1*t + \beta_2 * \sin(wt) + \beta_3 * \cos(wt) + \epsilon$$ For now, assume that $\epsilon$ is AR(2) with KNOWN parameters $\phi_1$ and $\phi_2$. I want to estimate $\beta$ coefficients above and I want to use pre- whitenning approach. So, what I do is: $\epsilon_t = w_t + \phi_1 \epsilon_{t-1} + \phi_2 \epsilon_{t-2} \rightarrow w_t = \epsilon_t - \phi_1 \epsilon_{t-1} - \phi_2 \epsilon_{t-2}$ So, I define: $Bt = t-1$ and $B^2t = t-2$ and I write: $\Phi(B) = 1 - \phi1B - \phi 2B^2$ and I multiply this operator function to the equation above: $$\Phi(B)Y = \Phi(B)\beta_0 + \Phi(B)\beta_1*t + \Phi(B)\beta_2 * \sin(wt) + \Phi(B)\beta_3 * \cos(wt) + \Phi(B)\epsilon$$ and I get: $$\tilde{Y} = \tilde{\beta_0} + \tilde{\beta_1}*t + \tilde{\beta_2} * \tilde{\sin(wt)} + \tilde{\beta_3} * \cos(wt) + w_t$$ where $w_t$ is white noise. I can then simply fit this recent model since noises are white noise. My question what is the relation between estimated $\tilde{\beta}$ and my original $\beta$'s?


title: How to choose the link function when performing a logistic regression?

text: I am running a logistic model. In SAS Entreprise Miner, I noticed there's a link function that has three possible options: `logit`, `probit` and `cll` (complementary log-log). Can you please shed light on the following questions:   1. Can we use any of these link function to carry out a logistic regression?   2. Are there situations where one would be better than others?   3. Is it intuitively possible to get some insight about which kind of function could be useful in which situation? (By just looking at the formula, complementary log-log function might be good for normalization of data when data does not depart (too much) from a normal distribution.) Any additional pointers would be greatly appreciated.


title: How do I ensure PROC ARIMA is performing the correct parameterization of input variables?

text: I'm trying to forecast using ARIMAX with two exogenous (input) variables. I'm using PROC ARIMA, but I can't figure out from the SAS documentation whether my code is producing the parameterization I want. I want to extend an ARI(12,1) model so that it also includes the last 12 terms of each of the two exogenous variables in my forecast. So, using `VariableX` with the two exogenous variables `VariableY` and `VariableZ`, my best attempt at the code is:               proc arima;       identify var=VariableY(1) nlag=24;       estimate p=12;       identify var=VariableZ(1) nlag=24;       estimate p=12;       identify var=VariableX(1) nlag=24 crosscorr=( VariableY(1) VariableZ(1) );       estimate p=12 input=( VariableY VariableZ );       forecast id=MonthNumber interval=month alpha=.05 lead=24;       run;     quit;      The documentation leads me to believe the first four lines of the procedure are required for setting up the forecast at the end. But when I run the procedure, the output appears to show a forecast using only the last term of each of the two exogenous variables. ## In summary, I'd like to be sure where each of the following are controlled:   * The $p$ of $AR(p)$, and similarly for each of the exogenous variables   * The $d$ of $I(d)$, and similarly for each of the exogenous variables   * The $q$ of $MA(q)$, and similarly for each of the exogenous variables


title: Effect of difference operator on sampling

text: Suppose I have a dataset $D$ and let it be split into 2 smaller datasets $D_1$ and $D_2$ such that $D = D_1 + D_2 $. Thus we can also say that $D_1=D - D_2$. Will we get a simple random sample of $D_1$ if we took a simple random sample of $D$ and removed from it those data points obtained by taking a random sample of $D_2$?


title: Questions about variable selection for classification, and different classification techniques

text: I have a question concerning feature selection and classification. I will be working with R. I should start by saying that I am not very familiar with data mining techniques, aside from a brief glimpse provided by an undergraduate course on multivariate analysis, so forgive me if I am lacking in details regarding my question. I will try my best to describe my problem. First, a little about my project: I am working on an image cytometry project, and the dataset is composed of over 100 quantitative features of histological images of cellular nuclei. All of the variables are continuous variables describing features of the nucleus, such as size, DNA amount, etc. There is currently a manual process and an automatic process for obtaining these cellular images. The manual process is (very) slow, but is done by a technician and yields only images that are usable for further analysis. The automatic process is very fast, but introduces too many unusable images - only about 5% of the images are suitable for further analysis, and there are thousands of nuclear images per sample. As it turns out, cleaning the data obtained from the automatic process is actually more time consuming than the manual process. My goal is to train a classification method, using R, to discriminate between good objects and bad objects from the data obtained from the automatic process. I have an already classified training set that was obtained from the automatic process. It consists of 150,000 rows, of which ~5% are good objects and ~95% are bad objects. My first question deals with feature selection. There are over 100 continuous explanatory features, and I would like possibly get rid of noise variables to (hopefully) help with the classification. What methods are there for dimensionality reduction with the goal of improving classification? I understand that the need for variable reduction may vary depending on the classification technique used. Which leads to my second question. I have been reading about different classification techniques, but I feel that I cannot adequately determine the most suitable method for my problem. My main concerns are having a low misclassification rate of good objects over bad objects, and the fact that the prior probability of the good objects is much lower than the prior probability of the bad objects. Having a bad object classified as good is less of a hassle than recovering a good object from the pool of bad objects, but it would be nice if not too many bad objects were classified as being good. I have read this post and I am currently considering Random Forests as per chl's answer. I would like to explore other methods as well, and I would like to collect the suggestions of the good people here at CV. I also welcome any readings on the subject of classification that may be helpful, and suggestions for R packages to use. Please ask for more details if my post is lacking in details.


title: Proportion of Intervals from t-interval method

text: Here is an excerpt of a book question: > The t-interval method assumes that the underlying distribution is > approximately normal. If we were to generate 800 samples of sample size 15 > from this exponential distribution and the calculate 95% t-intervals, what > is the proportion of intervals that include the true mean 900? Is this what > you would expect? Explain. Here is the R code for 800 samples.               sims = 800     n = 15     tsamples <- replicate(sims, rexp(n, rate = lambda))     mu = 1/lambda          results <- as.numeric(sims)     t.int <- matrix(FALSE, sims, 2)     for (i in 1:sims) {         t.int[i, ] <- t.test(tsamples[, i], conf.level = 0.95)$conf.int         results[i] <- t.int[i, 1] < mu & t.int[i, 2] > mu     }     sum(results)/sims      Output changes every time but a few are 0.9175, .91, .905, .925 My question is, I don't truly understand what it means by asking for "proportion of intervals that include the true mean 900". How can I figure out what this is asking and what the answer is? Is it asking how many times the R result is within a certain amount? Am I supposed to run the code multiple times with sims = 900 and find the lowest and highest number for a range?


title: Variable magnitudes and principal component analysis

text: I am completing some PCA on various sets of data, some with 30 time-series variables, some with about 100 (meaning matrices with between 30 and 100 columns). I know that the covariance matrix is calculated and then the eigenvectors and eigenvalues etc., so we are talking about the relationships between variables. My understanding: we are aiming to reduce the number of dimensions needed to explain the whole data set as accurately as possible, or expressed in other words, we want to know what the driving factors of a model are. I have read explanations trying to visualise PCA as a body of data. E.g if we have 3 variables, we could say the data represents a 3-d shape. we could have a sphere (all dimensions/variables possess equal weight/importance), a sheet of paper (one dimension is pretty much redundant) or even a piece of string (two dimensions are not very explanatory). **I want to know, do the actual magnitudes of the original data play a role at all in the end result in PCA?** Is there a way anyone can explain this with similar examples to those I have given? I am asking because my data (at some point along the time series) underwent a few changes, which meant that only a few of the variables where changed in magnitude - the were scaled down. So the value on one day was 100, it was cut in half to 50, and then all historical data points were retrospectively reduced accordingly. I want to know if these specific changes affect the output of PCA. Thanks in advance!


title: Mean squared error definition

text: I'm currently working through (part of) a textbook on non-parametric regression techniques. Regarding the choice of smoothing parameter the book starts out explaining the MSE which is defined as: $MSE(\hat f(x)) = E[(\hat f(x)-f(x))^{2}]$ Would this really be the MSE? The book also mentions the mean squared error of prediction (separately) but doesn't give a definition other than to explain that MSEP is for new values $(x^{*},y^{*})$. I'm a bit confused as I've seen the term MSEP applied to the definition of MSE given by the book. Thanks!


title: Mixed-model ANOVA vs multilevel modelling

text: Is mixed model ANOVA the same as multilevel modelling? If so, how do they differ? I am trying to compare inter- and intraindividual differences and not sure which one is the better approach.


title: Interpretation of error term in Moving Average (ARIMA)

text: I have an elementary question regarding the error term in MA (ARIMA)-- From where does this error term come from? From what I understood from the question raised earlier in the following link: Moving-average model error terms , the error term is the not the residual of regression yt (t is the subscript) over its previous value y(t-1), y(t-2),... (Please validate my understanding by answering this above sentence in YES/NO) but then how does the model get the estimates of the coefficients of the error terms from the data. Here both the error term and the coefficients are unknown, so how does it get the coefficients? I am definitely missing something. Please help.


title: What is the difference between scores in Princomp vs. factanal?

text: In R the `princomp()`and the `factanal()` are somewhat similar. At least their output looks pretty similar. I learned that this is not surprising since the print function of `princomp` comes from `factanal`. I understand that SS loadings do not make much sense for `princomp` as it is bounded to `1` anyway. Moreover, as Joris stated on nabble, the proportion of variance is only printed because of the common print function, but does not contain valuable information when princomp is used. What I do not understand is rather not an R question but more a multivariate stats question what is the conceptual difference between these PCA and Factor Analysis functions as they are used in R? This question relates particularly to the scores (let's assume "regression" scores for FA) respectively the difference between scores in both concepts? What should I rather use when I want to use to resulting scores in a regression model (for example in order to circumvent multicollinearity)? I also understand that PCA has a fixed number of components while FA has fewer factors than variables. richiemorrisroe's answer in the thread suggested by Rob Hyndman might go into that direction.


title: How to accurately track the 75% quantile in a non-stationary timeseries?

text: I have a non-stationary timeseries with a mean (µ) and standard deviation (SD) which both vary across time. The distribution of the timeseries is skewed, so the left and right sides of the distribution must be treated separately. I am wondering what is the best way to accurately track the 25% and 75% quantile? **What I have tried** I have tried calculating the 25% quantile within a rolling window, then feeding these values into the forecast() package to smooth them out. I'm wondering if there is a better way of doing this, because one of the constraints is that there is not much data: I only have about 50 samples within each window, any more and the SD and mean have changed so much that they are out of date.


title: How can I convert my dataset into a new dataset which is adjusted for confounding covariates?

text: I have used software before to do linear regression and factor in/out the confounding variables, but what I would like to do is generate a new data set which is adjusted for the confounding variables. How would I go about this? Further explanation: I have 1 dependant variable, multiple interesting independant variables and multiple confounding independant variables. I would like to   1. Adjust the dependant variable for the confounding variables.   2. Do PCA analysis on the interesting independant variables.   3. Colour my data points on the PCA plot by the adjusted dependant variable. I can do steps 2 and 3 but don't know how to do step 1. I can use R and have used MiniTab and SPSS once or twice, so any specific instructions relating to those would be amazing but even to be pointed in the right direction would be great.


title: Interpreting LASSO tables in SAS

text: I have been working on LASSO in SAS lately, and I'm still trying to figure out how to work with the options, but my main question for which I have not been able to find an answer on the internet so far is: how do you interpret the SAS output? For example (using SAS website examples): > ![Lasso summary](http://i.stack.imgur.com/wIvFg.png) I know that in this one, you can see what has been entered and removed, but I'm not sure about `CP`. > ![Lasso 2](http://i.stack.imgur.com/sPs5y.png) Here I know it tells me which variables have been retained, but I'm not sure if I need to report anything from this or what it means. > ![enter image description here](http://i.stack.imgur.com/upRsp.png) Do these numbers represent the effect each variable has on the outcome? What numbers tell me if this is a good model or if the variables chosen with LASSO are actually good predictors? Thanks.


title: Determining trends in text

text: I have a bunch of documents (66 quarterly reports on grievances and complaints on health care services provided and growing) and a list of words that I'd like to follow over time. What is the easiest way to do this? I have played with R's text mining library and have got rather frustrated. I have also tried to use RapidMiner and it chokes on three documents (ran out of memory). I would greatly appreciate any suggestions, ideas etc...


title: Profile likelihood confidence intervals

text: The theory behind **profile likelihood** (PL) **confidence intervals** (CIs) is clear to me. (See here, for example). SAS is surprisingly quick in calculating the PL CIs for all the covariates in a given model. Therefore I was wondering if there is another way of calculating PL CIs that does not require the time-consuming loop of maximizing the likelihood (over the nuisance parameters) given different fixed values of the parameter of interest. Example: I run a logistic regression in SAS on ~43000 subjects with 10 covariates (+ the intercept). It took 0.6 seconds to fit this model without asking for PL CIs, and only 3 seconds when I required also the PL CIs for all the 10 covariates (option `clodds=pl` in `proc logisitc`). In Stata, the same model was fit in 0.35 seconds (without PL CIs) and in 10 seconds when asking the PL CIs for just one of the 10 covariates in the model (I used the `pllf` command, which iteratively fixes values of the parameter of interest and maximizes the likelihood over the remaining ones).


title: Optimising control routines after creating a predictive model

text: I have a dataset with 1500 patients with data on recurrence of a disease. Follow-up time varies between 1 and 15 years. Approx 10% have recurrence. What I´d like to do is create a predictive model for recurrence into n similar groups that share the same temporal recurrence risk. This so that I could advise on an optimal control scheme. I can create the predictive model quite allright, but I´m uncertain as to what technique I can use to optimise the interval for controls. I guess a brute force simulation approach could be allright, but I´m sure there are other ways to do this. Advice or pointers to literature would be very welcome. EDIT: By control I mean how often they should be seen and evaluated for recurrence in the outpatient clinic. @rolando: 1) from what is already know on recurrence of this disease- the patients can be broken down into high, medium and low recurrence risks. This aso makes the most sense with respect to advising on a control scheme for the recurrence. Some patients recur within one year (2-3%), some 3-5 years and some after 10 years. Clearly it would be costsaving to have the latter patients be screened for recurrence at 10 years instead of yearly as would be indicated in the high risk group. However, exactly how I can utilise a prediction model and create the different groups so that they are homogenous with respect to time of recurrence is one of my problems. 2) n is the number of different groups with the same recurrence risk-see 1 3)By optimising the control routines (see 1)- A patient with low recurrence risk should not be seen as often as a patient with high recurrence risk. How can the optimal time-interval for controls in each group be calculated? //M


title: Why is the first postulate of the Poisson process that $\lambda dt$ is the probability of exactly one event in $[t,t+dt]$?

text: From the pdf of the Poisson distribution I would expect $\Pr(x=1)$ to be $$\lambda dt \cdot \exp(-\lambda dt)$$ I can see that as $dt$ gets very small, $\exp(-\lambda dt)$ becomes close to $1$, and so suggests $\lambda dt$, but I don't see why in the limit as $dt \to 0$ of $\lambda dt \cdot \exp(-\lambda dt)$, doesn't that make the $\lambda dt$ term zero also?


title: Maximizing number of positive targets in first decile in logistic regression

text: When We estimate logistic regression using MLE we try to minimize $-2ln(likelihood)$, which is equivalent to minimization of sum of squared deviance residual. Is it possible to fit logistic regression in such a way that number of positive targets in first decile is maximized? Decile is calculated using predicted score. The order in other nine deciles should not be relevant at all. **Update 2013-03-28 15:35**   Maybe I will describe this problem from a bit more business perspective. Lets say we have a huge database of potencial customers. Some of them we have already called and we know whether they bought our product or not (this is out target). There is limit on number of phone calls that we can make and we want to maximize the sales. The limit is equal to 10% of customer number. We can build logistic regression model and call 10% of people with highest scores, but intuitively I feel that something better can be done. MLE will also maximize the likelihood of those who probably won't buy. If we focus more on those with highest score then we might obtain better result. From top of my head I can think of following method:   1) Build logistic regression model on whole population.   2) Build logistic regression model using only 20% of people with highest scores in first model. I'm not sure whether this solution makes sense, but I just wanted to show my way of thinking.


title: Explanatory power of a variable

text: I have simple linear regression model. What I want to calculate is how "important" each of my input variables are i.e. to make a statement something like this: "60% of predictive power in this model comes from variable var1, wheres var2 and var3 have30% and 10%" respectively" What I need to do to calculate these percentages?


title: Multiple t-tests

text: Suppose we have some baseline data $X_1$ with $n$ observations and some post- baseline data $X_2$ with $n$ observations. We want to compare whether the post-baseline data differs significantly in mean from the baseline data. Would using multiple t-tests mean that we would need to make $n^2$ comparisons? That is, we compare each data point from post-baseline to every other baseline data point? Would this have a high Type 1 error? Should we be looking at the family-wise error rate?


title: Confusion related to marginal precision

text: I was reading this book where it was mentioned that if my gaussian distribution has unit marginal precision, the covariance matrix equals the correlation matrix. I didn't quite get it. Any clarifications will be much appreciated


title: Adding variables to the model one by one, or at the same time

text: What are the benefits of adding the variables into a model one by one, as compared to adding them all at the same time? As I see in most research, the first model that is tested is consisted of all the control variables, then in the second model the first predictor is added, then the second one, and so on, until all predictors are in the model. Are there any benefits in using this approach, or would it be ok to have just do two models - one with just the controls, and one with controls + all other predictors. (p.s.) - am interested in this with regards to logistic regression


title: What is the effect of transforming the DV on log likelihood?

text: I would like to know which (non-linear) transformation of my DV yields a better model (not necessarily with the aim of producing a DV that better fits the assumptions of the test). My impression is that the Log Likelihood is dependent on the scale of the DV such that I can not compare models based on the Log Likelihood (or AIC/BIC for that matter). I suspect that I could do an inverse transformation on each of the IVs relative to the transformation I was considering on the DV, but a) I'm not entirely sure and b) it seems like there must be a better way. Would the proposed method work? Is there a better way?


title: How to make a one-sided test for proportion in R

text: I've to make some tests for proportions in R because of some A/B duties. I know that given s1, sc the number of success of two samples and n1, nc the total number of experiments, I can do this for H0: s1/n1 = sc/nc               prop.test(c(s1, sc), c(n1, nc), conf.level=.95, correct=FALSE, alternative='two.sided')      to get the confidence interval of the difference p1 - pc (when p=s/n) and the p-value. What about if I need to test H0: p1 - pc <= 0? I tried this but I'm not sure it's correct:                prop.test(c(s1, sc), c(n1, nc), conf.level=.95, correct=FALSE, alternative='greater')      Is it ok to check this way? How can I calculate the test statistic under pooled variance assumption?


title: K-means as limit of Soft K-means algorithm

text: I encountered the following exercise > Show that as the stiffness $\beta$ goes to $\infty$, the soft $K$-means > algorithm becomes identical to the original hard $K$-means algorithm, except > for the way in which means with no assigned points behave. Describe what > those means do instead of sitting still. The notation is as follows: The point $\mathbb{x}^{(n)}$ gives a weight to the $k$-th cluster update through the responsibility $$r_k^{(n)}=\frac{\exp{(-\beta d(m^{(k)},\mathbb{x}^{(n)}))}}{\sum_i \exp{(-\beta d(m^{(i)},\mathbb{x}^{(n)}))}}$$ where $m^{(k)}$ is the $k$-th cluster point and $d(x,y) =\frac{1}{2}\sum_i (x_i-y_i)^2$. Then the update of $m^{(k)}$ is done with $$m^{(k)}=\frac{\sum_i r_k^{(i)}\mathbb{x}^{(i)}}{\sum_i r_k^{(i)}}$$ As stiffness $\beta$ goes to $\infty$ I see how soft $K$-means becomes hard $K$-means. The second part of the problem is asking about means which would have no assigned points in the context of hard $K$-means and thus not move. With soft $K$-means however they would move but does anyone know geometrically what they do?


title: Nested Standard Deviations

text: This is a finance question but most of it is statistics related. A few years ago the Chicago Board Options Exchange created the volatility index VIX, which is really just the 30-day standard deviation of returns on the S&P500 index. Then as people start trading options on the VIX index, they created another index that measures the 30 day standard deviation of of VIX index. By doing this, the concept of "nested volatility" becomes a real-life object with lots money attached behind it. I wonder if there are any academic treatment on the subject. The type of questions I seek answers to are something like: _We know that variances follow chi-square distributions, but what happens when you take the variance n times?_ To take a look at these things, I simulated a driftless Brownian motion time series $$x_{t+1}=x_t+\mathcal N(0,1)$$ for 5,000 days, then I calculated its 10 day volatility, and the 10 day volatility of the volatility, and so on. The process is repeated 8 times. A few interesting observations were made.   1. The average value of the time series decrease by about 65–73% each time you take the volatility. The scale of this decrease also decreases over time.   2. The correlation the between time series and its volatility ranges from 0.4 to 0.6, and this correlation increases each time you take the volatility. This is very strange as CBOE claims that historically VIX has little correlation with its volatility index VVIX. It would be nice to find some theoretical explanations for these observations. I really hope there is something similar to power series expansions where we can break down an object into an infinite number of ever-decreasing smaller objects.


title: Modelling non-independent y through a naive model by using shifted values of y as input

text: Excuse the unclear title, but I'm trying to ask about the validity of using a shifted dependent variable as a predictor of itself as a crude way of increasing the accuracy of predicting outcomes which are known to be non- independent from each other. If that's not clear I mean something like: $$ [ y_1, y_2, \ldots, y_{n-1}] = \mathbf{X\beta}+ b [ y_2, y_3, \ldots, y_n] $$ I see that that there's an issue here with leakage between my dependent variable and training set but does anyone know of any term / references for doing this? I understand that a modelling technique incorporating non-independence from the outset is a more correct way of going about this, but assuming I am to use something like a Random Forest, is this a statistically valid approach? * * * **edit** : I think my poor notation was confusing, what I mean is that in the case of predicting $y_{n}$, I would use (along with my other predictors) the value of $y_{n+1}$ and/or $y_{n-1}$. Alternatively, given a sequence: `1,2,3,4,5`, when trying to predict the third value in this sequence, I'm looking to treat the outcomes as independent but include `2` and `4` as predictors, and am wondering what the pitfalls of this approach are, if there's a name for it or if it is something that outright "shouldn't be done" and I should just use something akin to an HMM.


title: Repeated measures design and longitudinal study

text: What is the relation and difference between repeated _measures design and longitudinal study? > **The repeated measures design** (also known as a within-subjects design) > uses the same subjects with every condition of the research, including the > control. For instance, repeated measures are collected in a longitudinal > study in which change over time is assessed. Other studies compare the same > measure under two or more different conditions. For instance, to test the > effects of caffeine on cognitive function, a subject's math ability might be > tested once after they consume caffeine and another time when they consume a > placebo. > > **A longitudinal study** is a correlational research study that involves > repeated observations of the same variables over long periods of time They look like the same from Wikipedia. Thanks and regards!


title: Estimating the constant in multiple regression

text: Estimating an OLS model with a single dummy variable in Stata, the constant is exactly equal to the mean of the dependent variable when the dummy is equal to 0. When estimating an OLS model with two dummy variables in Stata, the constant is no longer exactly equal to the mean of the dependent variable when the dummies are set equal to 0. Does any one know how Stata predicts the constant term in univariate and multiple regression?


title: Is there any R package can smooth the coefficients as a function of coefficients sequence?

text: I want to do some regression analysis that constrains the coefficients to vary smoothly as a function of their sequence. It is similar to the "Phoneme Recognition" example in the part 5 "Basis Expansions and Regularization" of _The Elements of Statistical Learning_ (Hastie, Tibshirani & Friedman, 2008). in the book it says: > The smooth red curve was obtained through a very simple use of natural cubic > splines. We can represent the coefficient function as an expansion of > splines $\beta(f)=\sum_{m=1}^M h_m(f)\theta_m$. In practice this means that > $\beta=\mathbf{H}\theta$ where, $\mathbf{H}$ is a $p × M$ basis matrix of > natural cubic splines, defined on the set of frequencies. Here we used $M = > 12$ basis functions, with knots uniformly placed over the integers 1, 2, . . > . , 256 representing the frequencies. Since $x^T\beta=x^T\mathbf{H}\theta$, > we can simply replace the input features $x$ by their filtered versions $x^* > = \mathbf{H}^T x$, and fit $\theta$ by linear logistic regression on the > $x^*$. The red curve is thus $\hat\beta(f) = h(f)^T\hat \theta$. But I am not sure about how to create the basis matrix $\mathbf{H}$. When using ns() function in R, how to set the knots? Any hint will be appreciated.


title: Cronbach's Alpha and Skewed Data

text: I am attempting to assess the internal reliability of a questionnaire, but the results are highly negatively skewed. Range = 6-20   Median = 20   IQR* = 1 Log normalisation doesn't work. Is there an alternative to Cronbach's or can I use ranking in any way (more info on how to do this would be great!)? *IQR = Interquartile range


title: Supervised Latent Dirichlet allocation (LDA) on a sentence level

text: I'm writing a bachelor thesis on formality classification of sentences in English via projecting the formal pronouns from foreign (e.g. French tu/vous) texts to English (you). It's not really relevant, so don't think about it. Therefore I want to use LDA (Latent Dirichlet allocation) on a larger segment than just words, i.e. sentences. Currently I just use labeled LDA where one document is a single labeled sentence. Therefore I completely ignore context, which is really important for this task. Currently I have two ideas for a supervised sentence LDA:   * Based on SLDA: ![enter image description here](http://i.stack.imgur.com/nn82I.png) Based on Chen (2010) An Adaptation of Topic Modeling to Sentences extended with a label set: ![enter image description here](http://i.stack.imgur.com/VO50F.png) Do you guys know of any work concerned with LDA on larger segments, e.g. sentences or paragraphs?


title: Interpretation of ordinal logistic regression

text: I ran this ordinal logistic regression in R:               mtcars_ordinal <- polr(as.factor(carb) ~ mpg, mtcars)      I got this summary of the model:               summary(mtcars_ordinal)          Re-fitting to get Hessian          Call:     polr(formula = as.factor(carb) ~ mpg, data = mtcars)          Coefficients:           Value Std. Error t value     mpg -0.2335    0.06855  -3.406          Intercepts:         Value   Std. Error t value     1|2 -6.4706  1.6443    -3.9352     2|3 -4.4158  1.3634    -3.2388     3|4 -3.8508  1.3087    -2.9425     4|6 -1.2829  1.3254    -0.9679     6|8 -0.5544  1.5018    -0.3692          Residual Deviance: 81.36633      AIC: 93.36633       I can get the log odds of the coefficient for `mpg` like this:               exp(coef(mtcars_ordinal))      mpg      0.7917679       And the the log odds of the thresholds like:               exp(mtcars_ordinal$zeta)                 1|2         2|3         3|4         4|6         6|8      0.001548286 0.012084834 0.021262900 0.277242397 0.574406353       Could someone tell me if my interpretation of this model is correct: > As `mpg` increases by one unit, the odds of moving from category 1 of `carb` > into any of the other 5 categories, decreases by -0.23. If the log odds > crosses the threshold of 0.0015, then the predicted value for a car will be > category 2 of `carb`. If the log odds crosses the threshold of 0.0121, then > the predicted value for a car will be category 3 of `carb`, and so on.


title: Confidence or prediction limits for significant difference between forecast and observation?

text: My task is to to forecast business data without advertising and compare the real observations (with advertising effect) to the forecasted values in order to measure the effect of advertising. I am confused as to whether I should use the confidence limits of the predicted values or the prediction limits to calculate significant difference? I think I should perform a hypothesis test. My null hypothesis would be: there is no effect, the real observations will lie within the interval.


title: How do I extend Inverse distance weighting to take into account uncertanities?

text: From a Multibeam echosounder (MBES) I have $n \times m$ measurements: {$(x_i,y_i,z_i,d_i)$}$_i$ where $(x_i,y_i,z_i)$ is the location of the $i^{th}$ reflection from the seafloor and di is the distance of the $i^{th}$ reflection from the MBES. I want to grid these irregularly points to a $XY$ regular grid: $\ x = x0 + i*dx, 0<= i < nx, y = y0 + j*dy, 0<= j < ny $ So far I am using Inverse distance weighting (IDW). This makes sense: measured points closer to a node in the $xy$ regular grid have more say in deciding the $z$ value in this node than the measured points far away. However the measurements are uncertain, and the larger the distance; the larger the uncertanity. I have a simple model of the uncertanity as a standard deviation in $x$ or $y$ and in $z$: $\ 1.96*\sigma_{xy} = \alpha + \beta*d$ $\ 1.96*\sigma_{z} = \sqrt{a^2+(b*d)^2} $ where $\ \alpha = 2, \beta = 0.02, a = 0.25, b = 0.0075 $ I want somehow to take into account the uncertanities when gridding so that the more certain measurements gets a greater say than the less certain ones. **Any ideas on how to do this?**


title: Confusion related to difference of kriging and gaussian processes

text: I am having a hard time understanding what is the difference between kriging and gaussian processes. I mean wiki says they are the same but their formulas for prediction are so different. I am a bit confused why they are called similar. Clarifications?


title: How to interpret insignificant categorical variables for logistic regression

text: I am trying to interpret categorical variables with more than two classes. Some are significant whilst other classes are not. What can I infer from the insignificant ones? Does this mean the insignificant ones and the reference category equally influence the dependent variable? For example:               ETHNICITY (Reference Category - Indian)     Other Asian: Sig = .273     exp(b) = 1.123      African:     Sig = .000     exp(b) =  .148


title: What is wrong with this chi squared test of independence?

text: I'm trying to learn how to run the chi squared test of independence. I've created two fictional variables, Variable 1 and Variable 2, and I want to know if they are correlated. I'm working with the following values:                       Variable 1  Variable 2  Row Total     Group 1     100         200         300     Group 2     150         300         450     Group 3     200         400         600     Group 4     250         500         750     Group 5     300         600         900     Col. Tot.   400         800         3000      Which gives me the following table:               Observed    Expected    (O – E) (O – E)^2   (O – E)^2/E     100         40          60       3600        90     150         60          90       8100        135     200         80          120      14400       180     250         100         150      22500       225     300         120         180      32400       270     200         80          120      14400       180     300         120         180      32400       270     400         160         240      57600       360     500         200         300      90000       450     600         240         360      129600      540      Which gives me a X^2 of 2700. So, according to this test, Variable 1 and Variable 2 are not correlated. I just don't see how this could be possibly be true. Am I running the test wrong, or is my interpretation of how the test should work incorrect?


title: Is there any tool that can do Vector ARIMA modeling in time series

text: Vector ARIMA model is used in multiple time series analysis. I am just wondering if there is any software or tool can be used to build the model. Some tools,like R, can only be used to predict the value after the Vector ARIMA is already built by the user. So is there any tool can create the vector ARIMA model? Thanks in advance!


title: Three binomial random variables: test null that no two have equal mean

text: I have a dummy variable for health (0 or 1), and for each observation a treatment (one of three possible treatments for each observation). I want to test the null hypothesis that the proportion (mean) healthy within each treatment group is not significantly different than any other treatment group's proportion. \begin{align} H_0\\!: P_j = P_{j'}\quad &\text{for some } j \ne j' \\\ H_A\\!\\!: P_j\ne P_{j'}\quad &\text{for all } j \ne j' \end{align} That is, my alternative hypothesis is that every proportion (mean) is different than every other proportion. My data looks something like this:               treatment, health     1,         1     1,         0     2,         0     2,         0     3,         1     3,         1     ...


title: What are the definitions of semi-conjugate and conditional conjugate priors?

text: What are the definitions of **semi-conjugate priors** and of **conditional conjugate priors**? I found them in Gelman's _Bayesian Data Analysis_ , but I couldn't find their definitions.


title: How to detect covariate shift?

text: According to this site A Simple Machine Learning Method to Detect Covariate Shift, covariate shift happens when the distribution of test and train data are not the same or not similar. The site also provided a way to check if the data has covariate shift. However, I'm not sure of the validity of the method that the site proposed nor do I know how to apply the phi coefficient. If the method happens to be invalid, how then would I determine if the data has covariate shift? **ADDITIONAL:** This site provides a way on what to do if data has a covariate shift.


title: Combining covariances

text: Imagine we have one sample (n=40) from a population on 2 variables, for which we estimate the covariance (say, 0.1). If we would have another sample (n=10) from the same population, and we estimate the covariance to be 0.2, what would we expect to be the covariance in the population? Would it be: (4/5 * 0.1) + (1/5 * 0.2) = 0.12? That is, could we weigh the covariances according to the sample size? If yes, how to combine the data into having a covariance of 0.12, since copying one data below the other won't do.


title: How do you verify a rating system?

text: Let's say I have a database of a lot of (thousands) chess games. In chess, every player have a rating, called ELO that gives an estimation how strong is this player. If we know the the ELO ranking of both players, we can give an estimation of the probabilities of the possible outcomes of the match. In theory, a player whose rating is 100 points greater than their opponent's is expected to score 64%. This is theory. The player's true skill can be different. E.g. if he/she is really tired, or out of practice, his/her true skill will be lower than his/her ELO . So for every game I have an estimation (if we don't count the draw, it's pratcically a simple binomial experiment) of the probability of the different outcomes, and the outcome itself. Is there a statistic method, with which I can verify how accurate is this estimation? I need something that can tell _how accurate_ is the method of estimation not just disprove a null hypotesis that it doesn't have anything to do with player's true skills.


title: Iterative proportional fitting in R

text: **The mission** I am trying to find a way to do Iterative Proportional Fitting in R. The logic of the procedure is like this: one has a table with e.g. sample distribution of some variables. Let us say it is this one:               sample1 <- structure(c(6L, 14L, 46L, 16L, 6L, 21L, 62L, 169L, 327L, 174L,      44L, 72L, 43L, 100L, 186L, 72L, 23L, 42L), .Dim = c(6L, 3L), .Dimnames = list(         c("Primary", "Lowersec", "Highersec", "Highershort", "Higherlong",          "University"), c("B", "F", "W")))      Another table is from some other source, say another sample:               sample2 <- structure(c(171796L, 168191L, 240671L, 69168L, 60079L, 168169L,      954045L, 1040981L, 1872732L, 726410L, 207366L, 425786L, 596239L,      604826L, 991640L, 323215L, 134066L, 221696L), .Dim = c(6L, 3L     ), .Dimnames = list(c("Primary", "Lowerse", "Highersec", "Highershort",      "Higherlong", "University"), c("B", "F", "W")))      Now, we want to preserve the relationships between variables found in `sample1`, but we want to apply these relationships to the marginal distribution we find in `sample2`. Iterative Proportional Fitting does this as is described here (I could not possibly offer a better explanation). I have tried doing it in LEM, with the following result:                                   B         F        W     Primary     124204.64  960173.6 637701.7     Lowerse     119749.12 1081459.0 612789.9     Highersec   336934.21 1792001.6 976107.2     Highershort  90512.27  736464.1 291816.6     Higherlong   43486.91  238593.0 119431.0     University  163186.85  418628.6 233835.5      I am not 100% sure that this result is, but chances are (say 99%) that it is. The odds ratios from the first table are preserved in the resulting table, while the marginal distributions (row and column sums) are identical to the second input table. **The problem** Strangely, this quite useful algorithm is not readily available in R, at least not in a user-friendly form. One function that is likely to be relevant is `cat::ipf()`. However, I cannot figure out how to use the `margins=` argument. I am certainly not alone in this question. The help example uses a 3-dimensional table, which makes things even more confusing. In addition, there exist some user-written functions, one is here and the other is to be found here. The first one gives an erroneous result, unfortunately. The second one is also very non-transparent, requiring specifically pre-formatted CSV files as input, instead of R matrix objects. **The question**   1. Can anyone please explain how to actually use the `cat::ipf()` function?    2. Are there any alternative functions to achieve the IPF adjustment task, using matrices as input?    3. (SOLVED) Can this function be fixed to deliver a proper result?  Thank you. **ADDENDUM: I was able to get a proper output from the function in (3). After some consideration in turned out that the function does not accept matrix as input for marginal distribution, but simply a list of these marginal distributions. So practically the question is solved. However, a proper answer for 1 and 2 would be of use to the larger community, since IPF is quite essential in loglinear models.** Also reproducing a working IPF function here seems like a good idea, for those who search in the future:               ipf <- function(Margins_, seedAry, maxiter=100, closure=0.001) {         #Check to see if the sum of each margin is equal         MarginSums. <- unlist(lapply(Margins_, sum))         if(any(MarginSums. != MarginSums.[1])) warning("sum of each margin                                                        not equal")              #Replace margin values of zero with 0.001         Margins_ <- lapply(Margins_, function(x) {             if(any(x == 0)) warning("zeros in marginsMtx replaced with                                     0.001")              x[x == 0] <- 0.001             x         })              #Check to see if number of dimensions in seed array equals the number of         #margins specified in the marginsMtx         numMargins <- length(dim(seedAry))         if(length(Margins_) != numMargins) {             stop("number of margins in marginsMtx not equal to number of                  margins in seedAry")         }              #Set initial values         resultAry <- seedAry         iter <- 0         marginChecks <- rep(1, numMargins)         margins <- seq(1, numMargins)              #Iteratively proportion margins until closure or iteration criteria are met         while((any(marginChecks > closure)) & (iter < maxiter)) {             for(margin in margins) {                 marginTotal <- apply(resultAry, margin, sum)                 marginCoeff <- Margins_[[margin]]/marginTotal                 marginCoeff[is.infinite(marginCoeff)] <- 0                 resultAry <- sweep(resultAry, margin, marginCoeff, "*")                 marginChecks[margin] <- sum(abs(1 - marginCoeff))             }                 iter <- iter + 1         }              #If IPF stopped due to number of iterations then output info         if(iter == maxiter) cat("IPF stopped due to number of iterations\n")              #Return balanced array         resultAry         }      Note that this function will not accept a matrix for marginal distribution. The following will help:               m1 <- rowSums(sample2)     m2 <- colSums(sample2)     m <- list(m1,m2)      And then supply `m` as the first argument.


title: Hypothesis Testing of endpoint of Uniform distribution

text: Let $c$ be a positive constant. Suppose I have $X$~$Unif(0,\theta)$ and I wish to test the null hypothesis $H_0: \theta \leq c$ against the alternative hypothesis $H_1:\theta > c$. Suppose I did a sample $X_1,...,X_n$ but I obsrerve them with perturbation - that is, I actually have on hand $Y_1:=X_1+\epsilon_1,...,Y_n:=X_n+\epsilon_n$ where $\epsilon_i$ are i.i.d. ~ $N(0,\sigma^2)$. How should I conduct a hypothesis testing? **Idea:** I was thinking of using Likelihood ratio test, since that allows me to do hypothesis testing on a composite null hypothesis. I can take into account the errors in my likelihood function by convolving the uniform with the normal. However, when I want to define critical values, I have difficulty because I have no idea how to approximate the test statistic to any distribution and hence derive rejection regions. **Idea 2:** Rao's score test, since here I know how to approximate the test statistic distribution for large sample size, but the only form of Rao's test I know involves point null hypothesis (i.e. testing $\theta = c$) - sadly, my bound $c$ is conservative upper bound, so there is a large chance that $\theta$ is smaller than $c$ and I dont want this case to be a rejection case. Other ideas would also be appreciated.


title: How to calculate the probability of it raining at least once in a 30 day period?

text: **Practical situation** : I’ve got 120 days of data collected during rainy season. On an average it rained on 52.5 of those days.   1. What is the probability of it raining at least once in 30 days?   2. If it rains at least once in 30 days, then what is the probability that it continues to rain for 4 days in a stretch in those 30 days? My Answer for (1): I considered a Poisson process with $\lambda = 52.5/120 = 0.4375$, calculated the non-occurrence of an event (in this case "no rain") in 30 days and subtracted that value from 1 to get .9999 as the probability of it raining in 30 days. Please advise if I'm on the right track and also how do we go about part (2) of the question.


title: Hypothesis testing: small timeseries changes

text: I have pairs of timeseries that are estimating the same quantity over the years. It is some survival data: number of "dead" subjects during the year over the number of total subjects at the beginning of the year. The underlying data is regularly revised, so I need to assess the effect of the changes. Below an example is given. ![the two series](http://i.stack.imgur.com/SmpCV.png) I want to specifically test if the mean (which is the parameter of interest) is the same for the timeseries pairs (the analysis will be repeated for a lot of pairs, please do not focus too much on the particular image although it is fairly representative). I have several questions:   1. Should I use a paired test or not?   2. If it is paired, how do I deal with the ties, i.e. zero differences (I will always have a lot of ties)?   3. I believe that `bootstrapping` can be usefull in this case. I have bootstrapped the sample mean difference. This gives a distribution centered around the observed mean difference. How can I compute a p-value out of this distribution? Would it be reasonable to bootstrap the t-test instead?


Best statistical method for analysing a single clubs football (soccer) results? I have every result for over 100 years in date order. Is there a recognised or recommended method for analysing results, making some sort of coefficient to help in predicting results or at least showing +/- progress based on each result. apologies for the noob


title: Is there a convenient form for this large covariance matrix?

text: Consider the following bivariate vector autoregression: $$X_t=\mu +X_{t-1}A+\varepsilon_t,\ \varepsilon_t \overset{iid}{\sim} MVN(0, V),\ X_t=(X_{1,t},X_{2,t})',$$ where the assumptions on the coefficient matrix $A$ are such that the process $\\{X_t\\}$ is stationary, i.e. there are no unit roots. _The goal is to find a nice (from a computational perspective) expression for the inverse covariance matrix for the vector $y=(X_{1,1},X_{1,2},\dots , X_{1,T}, X_{2,1} \dots X_{2,T})'$_. In the univariate case, i.e. when $\\{X_t\\}$ is an AR(1) process such as: $$ X_t=\mu + aX_{t-1}+\varepsilon_t,\ \varepsilon_t \sim N(0, v) $$ then Chen and Deo (2009) gives the following convenient expression: $\Omega^{-1}=B'B,$ where $$ B= \ \left( \begin{array}{ccc} \sqrt{1-a^2} & 0 & 0 & \dots & 0\\\ -a & 1 & 0 & \dots & 0\\\ 0 & -a & 1 &\dots & 0\\\ \vdots & \vdots & \vdots & \ddots & \vdots \\\ 0 & 0 & 0 & -a & 1 \end{array} \right)\ $$ and $v\Omega$ is the $T\times T$ covariance matrix for $(X_1,X_2,\dots,X_T)'$. If it helps, and I think it does, we may assume that $A$ is diagonal so that the bivariate autoregression is really two separate AR(1)-processes with possibly correlated innovations. The problem is research relate and I need to to speed up numerical computations of bilinear forms such as $w'\Sigma^{-1} u$, for some vectors $w$ and $u$, where $\Sigma$ is the covariance matrix of $y$. I have a feeling this may be a standard problem solved in some old paper -- I was indeed expecting to find something in a textbook even -- but I can't find any good references. All help, even if just hints and not full answers, is greatly appreciated. **Update** In the case with diagonal $A$, all the elements of $\Sigma$ are easily found. We have $$\mathrm{Cov}(X_{i,t},X_{j,t+h})=A_{j,j}^h \frac{\mathrm{Cov}(\varepsilon_{i,t},\varepsilon_{j,t})}{1-A_{i,i}A_{j,j}}=A_{j,j}^h \frac{V_{i,j}}{1-A_{i,i}A_{j,j}}.$$ Note that the fraction only takes 3 different values, one for $i=j=1$, one for $i=j=2$ and one for $i\neq j$. These values corresponds to different blocks of $\Sigma$. Each of these blocks' inverse can be found directly using Chen and Deo's expression (I think). Another observation is that the power of the $A_{j,j}$ term is just incremented by $\pm 1$ as we move up/down or left/right within these blocks of the matrix.


title: Count data on proportion or different kind of type?

text: I really do not know which kind of variable my response data is. I cannot find any answer to my question in the world wide web. I have the results from a survey of 900 farmers. Each farmer told how on how many hectares he used one pesticide. So the result can be maybe 0.11 = They used a pesticide on 11% of their fields or 1= Use on 100% of fields. I guess its count data on proportion, but I am not sure, because i did not find any similar example. The distribution of the data is shown below (there are also many farmers who used on 0% of their fields a pesticide: ![](http://s7.directupload.net/images/140808/ywnx9muy.jpg) I need to know it, because i try to do an analysis of my data with a Generalized linear model in R (with poisson,binomial or gamma errors). It is necessary to know if it is count data on proportion or something else. What would you guess is it continuous data or count data on proportion or something else?


title: Inference from a sometimes-random time-series

text: Let's say we have two cointegrated time-series, $Y_{1}$ and $Y_{2}$, and I want to assess the causal impact of $Y_{1}$ on $Y_{2}$. There is good reason to think that both variables are influenced by a variety of outside factors. $Y_{1}$ contains some observations equal to 0 (in my case, about 12 from a total 100) that are generated completely independently from the rest of the data (which are never equal to 0). That is, sometimes $Y_1$ is turned off in a quasi-experimental way. What would be the best way of narrowing down the range of possible impacts of $Y_1$ on $Y_2$? **Context** I am helping a non-profit assess the impact of their advertising spend, $Y_1$ on $Y_2$, their website hits. Both series ramp up exponentially towards the campaign's deadline, so the residuals of $\ln(Y_{1, t}) = \beta_{0} + \beta_{1}\ln(Y_{2, t}) + \epsilon_t$ are stationary. Both series have significant autocorrelation also. At several stages during the campaign the NFP renegotiated contracts with their media buyer, resulting in exogenous advertising black-out periods. I'd be very comfortable drawing inference in a cross-section; my experience with economic time-series is that 'causal' parameters are often poorly identified, using bodgy instruments or sign restrictions. But I don't know much about doing quasi-experimental inference in time-series. Any guidance would be greatly appreciated.


title: Can I use ANOVA if a sample is not normally distributed

text: My data are not normally distributed. I have to compare three groups for their achievement in skill learning. $n=4$ in each group. I have used three methods ($M1$, $M2$ and $M3$) to teach one skill to developmentally delayed children. Can I still use ANOVA? Do I have to use repeated measures of ANOVA or a time series? Please throw some light


title: Predicting the Weather

text: Given a tree trunk with concentric circles, can we predict the weather for each year? Each concentric circle accounts for a year that the tree has been on the Earth. The innermost circle is the oldest year and the outermost circle is the current year. I have a data set that tells how much rainfall there was for every year since rainfall could be measured. Given this information can we predict rainfall for the innermost years? I thought maybe the size of the concentric circles were dependent upon the rainfall for that year. This assumes a direct relationship. I was asked to think of this more broadly and in particular to think about this as a Markov Model.   1. How could I go about doing this prediction? Is a Markov model a good place to start, or are there simpler ways to approach the problem?   2. Even if there are other ways, I want to explore the Markov model as well. Where is a good place to start learning about Markov models? Any help and comments would be greatly appreciated!


title: Mean and variance of a general multivariate skew normal distribution

text: I have a problem about a general multivariate skew normal distribution. There is a $p\times 1$ vector, $\mathbf{y}=(\mathbf{y}_1',\mathbf{y}_2',\ldots,\mathbf{y}_n')',p>n$, which has the density as $$f_p(\mathbf{y})\propto\phi_p(\mathbf{y};\boldsymbol{\mu},\Sigma)\prod_{i=1}^n\Phi(c_1+c_2\mathbf{1}'\mathbf{y}_i)=\phi_p(\mathbf{y};\boldsymbol{\mu},\Sigma)\Phi_n(c_1+c_2A\mathbf{y}),$$ where $\phi$ and $\Phi$ denote the pdf and cdf of the multivariate or univariate normal distribution. $c_1,c_2$ are constant scalars. $A$ is a $n\times p$ matrix, $A=\left(\begin{array}{cccc} 1,\ldots,1 & & & 0\\\ & 1,\ldots,1\\\ & & \ddots\\\ 0 & & & 1,\ldots,1 \end{array}\right).$ The goal is to calculate the mean and variance of this distribution. Gupta et al. (2004) have derived them for a specific distribution with $p=n$ and they seem tedious. I would like to know if there exists a simple solution for my problem. Any clue is welcome! Thanks!


title: Why calculating partition function is important in graphical models?

text: Consider a graphical model as specified at [1] Eq (1) and (2). Why calculating the partition function is important in graphical models? Why do we need it? I was under the impression that people calculate the partition in order to find the MAP, but Ermon assumes having a MAP solver, in order to find the partition solver. Doesn't it contradict my initial understanding of the reason for calculating the partition function. Or maybe this: he has assumes having the MAP solver for small subproblems (in the buckets), but needs to have the partition function in order to do the MAP on the whole problem. Is that a correct understanding? [1] http://www.cs.cornell.edu/~ermonste/papers/LPCount.pdf


title: determine independent and dependent variables without apriori knowledge of the nature of data

text: i have a spreadsheet with seven parameters and their values for 9000 samples.I need to determine which are dependent parameters and which are independent. Please let me know which algorithm might be helpful. i've thought of PCA and regression but i did not understand how i can use them to determine the independent and dependent parameters


title: Implication of convergence in rth mean

text: I'm trying to show that if $X_n \rightarrow X $ in $r$-th mean, then $E|X_n|^r \rightarrow E|X|^r$. (Edit: should have said with $r \ge 1$) My question is whether the following steps are sufficient or not: $[E|X_n - X|^r]^{\frac{1}{r}} \rightarrow 0$ as $n \rightarrow \infty$ Setting $r=1$ $E|X_n - X| \rightarrow 0$ as $n \rightarrow \infty$ And so: $E|X_n| \rightarrow E|X|$ as $n \rightarrow \infty$ Hence $E|X_n|^r \rightarrow E|X|^r$ as $n \rightarrow \infty$ Thoughts appreciated.


title: HLM-related inquiry - comparing the magnitude of coefficients

text: I've conducted two separate HLM analyses, each with two IVs predicting to the same DV in the same sample. The problem that I've encountered is that I would like to compare the **magnitude** of the resulting coefficients across each analysis to determine if they are equivalent. I'm using dyadic data, which is why I used HLM. Here is a more concrete example: Analysis 1: DV-1 IV-1 IV-2 Analysis 2: DV-1 IV-3 IV-4 To be concrete, I want to show that a computer-based coding of a test (e.g., which produces IV-1 and IV-2) produces similar coefficients as a human-based coding of the identical test (e.g., which produces IV-3 and IV-4). That is, I would like to be able to compare the coefficients obtained for IV-1 with IV-3 (which are the same variable, but one was produced using a computer and the other by a human), and IV-2 with IV-4 (again, these two should be the same variable, one produced by a computer and the other by a human), to examine whether they are equivalent in magnitude. Is this possible? Is there a statistical test that exists to quantify the difference? To clarify, I'm not interested in mediation per se. Rather, I'm just interested in whether the resulting coefficients are equivalent in magnitude (regardless of whether they are accounting for different variation in the DV).


title: What is the procedure to compare two different period time series

text: I am currently working on the task that I would like to compare two different period time series like Sales in 2012 vs Sales in 2013. Kindly suggest me any statistical procedure.


title: How to detect outliers in skewed data set?

text: I am working on my school datamining project. Within preprocessing stage I need to remove outliers from my data set which is positively skewed (see description). I have an idea to remove all values which are larger than mean + 3 x standard deviation, but I am not sure this is a suitable technique for my case because the data set is not normally distributed. What technique should I use?                 var     n    mean      sd  median trimmed     mad  min     max   range skew kurtosis   se     1   1 41019 1668.99 1107.08 1453.68 1524.22 1026.05 10.9 5920.74 5909.84 1.18 1.33 5.47


title: TelephoneSurvey, when to call?

text: I conduce a prestudy to understand when is the best time to call bakers for a telephone survey. Here is my preliminatory study design:   1. I have a table of all 15000 baker firms in a country    2. I select randomly 20 for each hour of day, all but from 19:00 to 23:00 (7pm- 11pm) where I do not expect them to work. 20h * 20 = 400 calls   3. I call them asking them politely if they would answer my survey.    4. For each call I record start and end time The longer the call the better.    5. Also, I record if the baker accepted the survey, and if he answered the phone at all How would you optimize this design so that I require fewer calls? I think about using bayesian theorem to establish gutt feeling priors. Like, at 03:00 they have to work so the probability is low or at 13:00 the major work is over so they will respond more frequently. Then I intend to update the bayesian priors with measuremed response rates Any ideas how to improve my study design? Particularly that I require fewer samples to come to a good judgment when to call?


title: Problems with the sample (unkown population)

text: I realized a large research with Brazilian physics high-school teachers and I wanted that my sample could be representative of the Brazilian physics high- school teachers population. However, nobody in Brazil knows how many teachers there are, in which cities they are or who they are to provide contact. So, during one year, I contacted teachers, schools, universities, governmental secretaries, labor unions, and event organizers and I reached more than 7000 names of possible teachers. To build my sample, knowing that a lot of them wouldn’t not answer the survey, I invited 33% of this 7000 names to participate of the research. I chose randomly this 33%, but I can’t contact many of them and others didn’t fit with my population (they are no more teachers, they teach only in Universities and so on). In the end of one year trying to contact this 33%, I got about 300 valid responses. But now, with all data collected, I’m not sure if my sample could be considered a simple random sampling (SRS) and how I can modify and fix my results to extend them to all population if my sample is not a SRS (or another previously known sampling method). Someone could help me? Thank you, Alysson Ramos Artuso


title: Creating a meaningful statistic measure for physical geography

text: I'm writing a paper on the collapse of ancient agricultural terraces in the Mediterranean, and have come across a problem when trying to create a meaningful measure. I have sampled with a GPS terrace failure points, and mapped them with GIS. I have assessed the damage in the field in a simplistic ordinal scale called "status", where "1" is partial damage, and "2" complete collapse to the bedrock. I have an hypothesis that the potential damage to terrace stone walls increases/decreases according to two factors:   * Increasing: The number of other terrace failure points in the immediate higher catchment area of the sample point. (Hydrologically speaking, a terrace failure point creates a channel for surface runoff to pour on the lower terraces, increasing saturation, peak infiltration etc., which brings about a higher chance of failure.)   * Decreasing: The number of other adjacent terrace failure points on the **same** terrace wall. (The reasoning this time being other likely preferred runoff channels that decrease the chance of higher damage in the sample point.) I'm trying to find a way to use these opposing "vectors" together, but I hit a snag. The obvious way is to divide the former (Neighbours) with the latter (Brothers) (Neighbours/Brothers), thus increasing one or decreasing the other will give me a higher ratio that represents the potential damage. The way in which I calculate the value of each is summing all the "status" values of the points which match my spatial query. However, this results in some points that have 0 value for Neighbours and/or 0 value for Brothers, so a simple division doesn't work well in all the cases. So I'm asking if there are ways to make this more meaningful. I have no pre-existing assessment on the possible outcomes of such a measure, and I don't know if I should be expecting a positive or negative (or any) correlation.


title: Words that estimate numerical proportions of a group

text: There seem to be more than a few papers that surround the use of certain words to define the probability of an event without actually using the numerical probability itself. For example, if you are certain that something is going to happen, the probability of the event is 100%. If you are almost certain, the probability of the event is 93% +/- 6%. http://en.wikipedia.org/wiki/Words_of_estimative_probability In election terminology, we have words like majority, plurality, minority, some fraction - majority/minority. These words are easy to convert to numeric proportions (in most cases). With those examples in mind are there any papers or generally accepted conventions for words like "many, most, few, some" that act as a guideline for numerical proportions associated with those words?


title: Fitting a distribution to sample data

text: I have completely modified this post for the sake of clarity. To Whuber and David, thank you guys for your posts. I want to perform a MC Simulation to see what will, more likely, be the tendency in the distribution of vessels by TEU Capacity (if there is going to be more vessels in the 1000-2000 TEU range, or in the 300 - 4000 TEU range, etc.). But to do the MC Simulation I need to generate random numbers from a specific distribution and I do not want to assume normality or lognormality or any other distribution if I do not do the goodness of fit check first. The thing is, what how would I go around the distribution fitting taking into account the number of vessels assigned to each TEU Range Class?               CLASS                   Count of Vessel          FREQ     13,000.01 - 14,000 TEUS       22             25%     8,000.01 - 9,000 TEUS         19             21%     4,000.01 - 5,000 TEUS         17             19%     6,000.01 - 7,000 TEUS          9             10%     5,000.01 - 6,000 TEUS          6              7%     9,000.01 - 10,000 TEUS         5              6%     1,000.01 - 2,000 TEUS          4              4%     3,000.01 - 4,000 TEUS          4              4%     16,000.01 - 16,500 TEUS        2              2%     2,000.01 - 3,000 TEUS          1              1%          Grand Total 89


title: Methods for cross-sectional twin data (GEE, mixed effects modeling, GLMM)?

text: I have a twin dataset of 248 pairs (166 MZ, 82 DZ) that was gathered at one time point. I'm interested in looking at the relationship between a physiological biomarker and self-reported insomnia levels. The biomarker is continuous, but I want to break it into categorical cut-offs so it is easier to interpret (this is how it is commonly treated in the literature). Therefore, I have a categorical predictor and continuous outcome. I'd like to first examine this relationship, and then examine if this relationship is also affected by adding physical activity levels (self-reported, continuous) as another predictor. Next, I'll examine if season of measurement (categorical) moderates any significant relationships. I have a number of confounds that I want to adjust for (e.g., gender, BMI), as well. I plan to model these relationships first cross-sectionally (accounting for correlated nature of data), and then separately for MZ and DZ twins to assess for genetic and familial confounding. I'd like to analyze the data in SPSS, but have access to Stata and SAS if necessary. Both my biomarker and insomnia variables are non-normal (positively skewed), but do not fit the non-normal models that are typically analyzed using generalized linear mixed models. Therefore, I'd like to use GEE or linear mixed effects models with transformed variables. However, is there an approach that would be appropriate/better that I'm missing? Also, does anyone have an opinion on whether GEE or mixed effects modeling is more appropriate for these data? Thank you very much for any thoughts and suggestions.


title: Warning message with loglin/dmod function, package gRim

text: after executing following code: > library(gRim)   >  TA=table(ALL)   >  upp.model<\- dmod(~.^3, margin=colnames(ALL),data=TA) I get this warning message:               Warning message:      In loglin(object$datainfo$data, object$glist, fit = TRUE, print = print,  :      algorithm did not converge      Information about the data:                str(TA)        'table' int [1:3, 1:2, 1:2, 1:2, 1:4, 1:2, 1:2, 1:2, 1:3, 1:4, 1:2, 1:3, 1:2, 1:2, 1:2,   1:2, 1:3] 0 0 0 0 0 0 0 0 0 0 ...     - attr(*, "dimnames")=List of 17     ..$ a: chr [1:3] "(vors\xe4tzliche leichte) K\xf6rperverletzung" "Sonst."  "Gef\xe4hrliche K\xf6rperverletzung gem. \xa7 224 StGB auf Stra\xdfen, Wegen oder  Pl\xe4tzen"     ..$ b: chr [1:2] "ja" "nein"     ..$ c: chr [1:2] "ja" "nein"     ..$ d: chr [1:2] "ja" "nein"     ..$ e: chr [1:4] "1-18" "19-25" "26-59" "Sonst."     ..$ f: chr [1:2] "M" "W"     ..$ g: chr [1:2] "Sonst." "keine Vorbeziehung"     ..$ h: chr [1:2] "Sonst." "Keine Beziehung"     ..$ i: chr [1:3] "Sonst." "Sonstige Opfer" "Sonstige Gesch\xe4digtenspezifik"     ..$ j: chr [1:4] "10-18" "19-25" "26-36" "37-88"     ..$ k: chr [1:2] "M" "Sonst."     ..$ l: chr [1:3] "Sonst." "ohne" "Sonstiger legaler Aufenthalt"     ..$ m: chr [1:2] "Sonst." "nein"     ..$ n: chr [1:2] "ja" "Sonst."     ..$ o: chr [1:2] "ja" "nein"     ..$ p: chr [1:2] "Sonst." "nein"     ..$ q: chr [1:3] "Sonst." "im Bundesland" "in der Tatortgemeinde"      I am really not sure, how to handle this problem. Any suggestions? Thank you


title: eBayes() lmFit()

text: I am using the golub data set from R with the labels.               labelgb <- factor(c(rep("ALL",27),rep("AML",11)))          names(golub) <- labelgb      I need to use the empirical Bayes method in the limma package and calculate p-values for the golub dataset using the `eBayes()` function.               fit = lmFit(golub)     e = eBayes(fit)      In the `lmFit` function what does the design argument do and if its crucial for the golub dataset to supply the design argument?


title: Strange residuals interpretation

text: I have run a model and with the data I have, I wasn't expecting to produce the best model ever but my residuals are really strange. The outcome variable is number of days going to a website in a month (so discrete 0-30) and I have applied a standard regression. model <\- step(lm(day_count_jun ~ vars, data=ml), direction="forward"))) #forward-backward I have several independent variables and I have chosen (with function step in R) 5 of them. It is the same result if I include only 5 variables and all the variables I have. The adjusted R2 is 0.25 Here is the residual plot I have: ![enter image description here](http://i.stack.imgur.com/NPepm.png) Should I understand that one of my coefficients should have a negative sign? Besides, my residuals are not normal. ![enter image description here](http://i.stack.imgur.com/L5k9l.png) What can I do to solve this if I don't have more data?


title: Is normality testing 'essentially useless'?

text: A former colleague once argued to me as follows: > "we usually apply normality tests to the results of processes that, under > the null, generate random variables that are only _asymptotically_ or > _nearly_ normal (with the 'asymptotically' part dependent on some quantity > which we cannot make large); In the era of cheap memory, big data, and fast > processors, normality tests should _always_ reject the null of normal > distribution for large (though not insanely large) samples. And so, > perversely, normality tests should only be used for small samples, when they > presumably have lower power and less control over type I rate." Is this a valid argument? Is this a well-known argument? Are there well known tests for a 'fuzzier' null hypothesis than normality?


title: p-value depends of the size of the sample?

text: I'm trying to find the best distribution that fit my data (sample size of 7210).   Some of the distributions that have a low Anderson–Darling value look very similar to my data. But all the p-values show insignificance. Is that because of the sample's size? Anyway if I do a plot for Loglogistic, for example, it looks very similar to my data distribution.   BTW: I'm using MiniTab to do this. Thanks.               Distribution                   AD       P  LRT P     Normal                     23,480  <0,005     Box-Cox Transformation     13,896  <0,005     Lognormal                  15,057  <0,005     3-Parameter Lognormal      13,664       *  0,000     Exponential              2438,006  <0,003     2-Parameter Exponential  1882,128  <0,010  0,000     Weibull                   103,812  <0,010     3-Parameter Weibull        54,502  <0,005  0,000     Smallest Extreme Value    242,163  <0,010     Largest Extreme Value      64,242  <0,010     Gamma                      14,018  <0,005     3-Parameter Gamma          14,076       *  0,524     Logistic                   17,623  <0,005     Loglogistic                13,213  <0,005     3-Parameter Loglogistic    12,822       *  0,004


title: R: Shapiro-Wilk test with relation of W and P-Value?

text: I'm currently learning about the Shapiro-Wilk test for normality and using R for conducting simple experiments. One thing which is confusing me is that I found some results would produce a high W value whilst the p-value is very low. If I understand this correctly, the W value may be between 0 to 1. A small W leads to rejection of normality whilst high W leads to the alternative hypothesis. How then does this relate to the P-Value? In other words, is p-value somehow derived from the W result? Leading on from this, here's an example which is confusing me:               set.seed(1000)     xgamma <- rgamma(5000,1,2)     shapiro.test(xgamma)          Shapiro-Wilk normality test          data:  xgamma     W = 0.811, p-value < 2.2e-16      I'm generating 5000 random gamma values which result in a W of 0.811 (unexpected) and p-value of 22e-16. I thought that the W value should be low since the p-value is tiny in comparison. How should I interpret the W value in comparison with P-Value? Thanks for any comments in advance!