title: Are t test and one-way ANOVA both Wald tests?

text: t-test for testing whether the mean of a normally distributed sample equals a constant is said to be a Wald test, by estimating the standard deviation of the sample mean by the fisher's information of the normal distribution at the sample mean. But the test statistic in the t test has a student t distribution, while the test staistic in a Wald test asymptotically has a chi- square distribution. I wonder how to explain that? In one-way ANOVA, the test statistic is defined as the ratio between between- class variance and within-class variance. I was wondering if it is also a Wald test? But the test statistic in one-way ANOVA has a F distribution, and the test statistic in a Wald test asymptotically has a chi-square distribution. I wonder how to explain that? Thanks and regards!


title: Regressing coffee beans in US to coffee beans in EU

text: We're trying to model two time series: a random walk (independent variable) vs. the sum of this random walk and a mean-reverting process. For example: coffee bean 100kg prices (EU) vs. coffee bean 100kg prices (US). The difference between these two can be seen as costs in the transportation process, seen as the mean-reverting part. Now we want to regress these two against each other. Through OLS we achieve, as expected, a coefficient of around 1. However, our results show a large amount of autocorrelation in the residuals due to the fact that transportation costs has high and low periods. To remove this autocorrelation we have tried to make them stationary by taking first differences for x and y. After doing this, our autocorrelation is completely gone, but the regression coefficient now deviates from 1 (for example 0.6), which is counterintuitive as the prices are for the same amounts of coffee beans. As such, we lose the trend in our data. How to correct for autocorrelation without losing the trend in our time series? Would using Newey-West standard errors combined with the OLS coefficients give an accurate alternative to modeling the autocorrelation? PS: there are no known tariffs/taxes for importing coffee beans


title: Correlation, regression and causal modeling

text: This is probably a blindingly obvious answer for any seasoned statistician, but I am still confused as to how correlation differs to regression, technically. I understand that one is a measure of association and one a measure of causation, but how can you actually measure causation mathematically, without actually conducting a real life experiment. If that makes sense. Thank you.


title: Mboxcox, interpreting difficult regressions

text: Mboxcox in Stata suggests transforming my variables using a power of 0.1 for the independent variable, and 0.4 for the dependent variable. I have run the model, and it fixes problems associated with the assumptions of OLS. But certainly, it complicates matters in terms of interpretations. Please outline possible interpretations and solutions. Dependent variable is in millions of dollars, and independent variables are number of years. and a number of dummies which do not require transformation.


title: Structural equations: how to specify interaction effects in R lavaan package

text: I am using R lavaan package to estimate a structural equation model. Let's say the model consists of 1 endogenous manifest variable with 1 latent and 2 manifest explanatory variables:               group = {0,1}     attitude1 = latent,scale     age = respondent's age      The desired lavaan model is then (doesn't work):               model <- '     attitude1 =~ att1 + att2 + att3     outcome ~ age*group + attitude1*group'      My goal is, in the lines of what can be done in linear regression, to **establish main and interaction effects between each variable and group. Can this be done?**


title: Poor parameter estimation in frailty models

text: I fitted a one-parameter gamma shared frailty model with Weibull and exponential baseline hazards for recurrent event data (the dataset contained 80 subjects of which 16 patients had events (13 subjects with just one event, 2 subjects with 2 events and one subject with 3 events)) and I found:   1. The Weibull and exponential scale parameter very large (about 4000) with a standard error of 3000. What could be the meaning of this?   2. The estimate for the variance of the frailty ($\theta$) is also large (about 15), with a standard error of about 7.2. Can I interpret this value of theta (15) as high between cluster variability? Or does the high standard error change its interpretation? Note that model convergence was achieved. Can someone help please.


title: Compare means of different variables

text: I've got more than 20 (10 point likert scale) variables with more than 1000 entries each. What I want to do is compare the means of the answers on the questions. A one-way anova seems suitable for this, but you can only categorize by the values of a variabele. I want to categorize by question, the variable itself. Is there any way to do this without having to put all answer beneath each other and having another variable saying which question it is. Just to be clear, I've got: variable: question1 --> 1000 entries ranging from 1 to 10 variable: question2 --> 1000 entries ranging from 1 to 10 ... variable: question20 --> 1000 entries ranging from 1 to 10 I want to compare the means of the different questions with one-way anova but I can't choose to factor by question.


title: What is a good introduction to  statistical hypothesis testing for computer scientists?

text: I was recently exposed to some statistical hypothesis testing methods (e.g. Friedman test) at work, and I would like to increase my knowledge on the topic. Can you suggest a good introduction to statistical significance / statistical hypothesis testing for a computer scientist? I am thinking of a PDF book or similar, but any other kind of help is welcome. Edit: I've already found this website but I was looking preferably for some which is easily printable. Thank you   Tunnuz


title: Interpretation of 3-way ANOVA

text: I have 2 genes (tf1 and tf2), which are affecting a third gene (tg). By affecting a gene I mean the changes in the value of tf1 and tf2 changes the value of tg. what we measure in the whole experiment is this value. We want to see if the values of tf1 and tf2 are dependent (which means they are regulating each other) and if (tf1,tg) and (tf2,tg) are dependent. The response of the tf1 and tf1 and the tg differs in different conditions. So my 3 factors are tf1, tf2, and conditions.   * First I want to make sure that the design of my 3way-anova is correct. One dimension is the conditions.If we see the second and third factors (tf1 and tf2) as a matrix, I have designed it like this.   |tg |tf1 |   |tf2 |func(tf1,tf2)| so in the first column, there is no tf1, and in the secong column we always see its effect. In the first row there is no tf2, but in the second row we see its effect. I thought that it would be good idea to fill the cell in the first row, first column with the tg. and the last cell has an effect of both of tf1 and tf2 (I combined the values of tf1 and tf1 with a function like average).   I appreciate if you tell me what you think about this design.   * If a:tf1, b:tf2, and c:condition, which of the SSs shows if there is a dependency between tf1 and tf2? Is it SSab?  The reason that I am asking is that I don't see much of a difference in the histogram of the SSab for samples in which there is an interaction between a and b with the samples that there is no interaction between a and b.   * Also, what is the interpretation of SSabc for my particular problem?


title: Evaluating recommender systems with (implicit) binary ratings only

text: I'm analyzing a set of news articles and user libraries. User library is the set of news articles shared by one user. Obviously, the rating is 1 (the article is in user's library) and 0, otherwise. I hypothetically present every user with M articles sorted by their predicted rating and evaluate based on which of these articles were actually in the library. So, my only evaluation metric is > recall@M = # of articles the user shared in top M / total # of articles the > user shared Using precision doesn't make sense as zero ratings are uncertain (a user is not interested in the article or is just unaware of it). If I can't afford the user studies (explicit user feedback on recommendations), what are other evaluation metrics can I use? Thank you in advance!


title: Software for fitting generalized least squares model with errors that follow seasonal ARMA model

text: I need to fit a GLS model, with some known regressors, and where the errors follow an **unknown** ${\rm ARIMA}(1,0,1) \times (1,N,1)$ model. It seems like the main tool out there for such models is the `gls` function in the `nlme` package for `R`. In `gls`, one specifies the correct correlation struction using a `corStruct` object, but I cannot find any `corStruct` objects for specifying my (really simple) seasonal model. I am new to R, so I don't think I am up for coding a new `corStruct` for my purposes. Are there any other packages out there for solving this problem? If not, can you point me to some references about how to create custom `corStruct` objects. Thanks for all your help!


title: Improvements in AIC when adding a variable that is not close to signficant

text: In a linear model (regular OLS regression) I am exploring various models of a continuous dependent variable. When adding some variables, the Akaike Information Criterion (AIC) improves (sometimes dramatically) even though the variable is not close to significant (e.g. p = .49 or even higher). One reason I can see this might happen is that the _other_ variables in the model become better fit upon addition of the new variable. Are there other reasons? As an example: Model 1: $F_{10, 363} = 4.31, R^2 = .11, AIC (\text{smaller is better}) = 1497.5, -2LL = 1495.5$ Model 2: $F_{11, 362} = 3.93, R^2 = .11, AIC = 1496.1, -2LL = 1494.1$ and the new variable added has $t = -0.10, F = 0.92$. (This happened with several models)


title: lme and lmer comparison

text: I was wondering if anyone could enlighten me on the current differences between these two functions. I found the following question: How to choose nlme or lme4 R library for mixed effects models?, but that dates from a couple of years ago. That's a lifetime in software circles. My specific questions are:   * Are there (still) any correlation structures in `lme` that `lmer` doesn't handle?   * Is it possible/recommended to use `lmer` for panel data? Apologies if these are somewhat basic. A bit more detail: panel data is where we have multiple measurements on the same individuals, at different points in time. I generally work in a business context, where you might have data for repeat/long-term customers over a number of years. We want to allow for variation over time, but clearly fitting a dummy variable for each month or year is inefficient. However, I'm unclear whether `lmer` is the appropriate tool for this sort of data, or whether I need the autocorrelation structures that `lme` has.


title: Regularization and feature scaling in online learning?

text: Let's say I have a logistic regression classifier. In normal batch learning, I'd have a regularizer term to prevent overfitting and keep my weights small. I'd also normalize and scale my features. In an online learning setting, I'm getting a continuous stream of data. I do a gradient descent update with each example and then discard it. Am I supposed to use feature scaling and regularization term in online learning? If yes, how can I do that? For example, I don't have a set of training data to scale against. I also don't have validation set to tune my regularization parameter. If no, why not? In my online learning, I get a stream of examples continuously. For each new example, I do a prediction. Then in the next time step, I get the actual target and do the gradient descent update.


title: Population vs sample variance

text: This is perhaps a very basic question, but I am a bit confused. I have a time series of day some time varying data, say the amount of sunlight as the day passes. I bucket the data in say 10 bins say according to just time, so later bins occur later in time and calculate the variance of each of the bins. Is there a relationship between the calculated variances of the 10 bins and the overall population variance (variance of the whole data set)?


title: What is the correct way to choose tests for pairwise comparison after ANOVA?

text: I am measuring secretion of insulin upon glucose stimulation of pancreatic islets after various treatments. So I have two independent variables/factors:   1. **Treatment** (three different compounds and a control treatment)   2. **Glucose concentration** ( _low_ and _high_ ) My sample size is only 3 to 4 individuals per experimental group. I analyzed the data by two-way ANOVA (using Prism for Mac), confirming that the effects of the two factors are not independent. Now I am unsure about what test I should choose to do pairwise comparisons of each treatment vs. the control. Prism allows for including **Bonferroni** 's post-hoc tests, after quickly checking the Wikipedia article about ANOVA, I believe **Dunnett** 's test might be appropriate, but I have seen publications using **Duncan** 's test, **Tukey** 's test as well as **Student's t test** for pairwise comparisons in similar setups. Having no suitable textbook at hand, and since I didn't find a good descriptive example from the life sciences/biology field online, I am not sure what test is appropriate. Can someone point me to some material (or recommend a good textbook) that guides the experimenter which (post-hoc/a priori) test to use, and ideally explains why?


title: Is R limited for matrix?

text: I am using R to find the value of lambda by cross-validation for ridge regression but the problem is my code as follows does not work for such a big matrix. Any help would be highly appreciated. Thanks in advance. My own work is               # X<- my own dataset which has n=205 observations by 18856 variables.      I = diag(18856)     cv = 1.0e30   #arbitrary value     for (i in 1:100) {         lambdai = i/100*1 # value of lambda cvi = 0        for (j in 1:n) {           bR = solve( t(X[-j,]) %*% X[-j,] + lambdai*I) %*% (t(X[-j,]) %*% y[-j])           yje = X[j,] %*% bR # fitted            values cvi = cvi +(yje-y[j])^2         }     if (cvi < cv) {      cv = cvi      lambda = lambdai      }        } Sorry for inconvenience.


title: Estimating parameters of categorical distribution from sum-of-outcome data

text: Let $X$ be a categorical random variable with possible outcomes $o_1,...,o_n \subset [l, u]$ (real numbers with a known lower bound $l>0$ and known upper bound $u$) that occur with probability $p(X = o_i) = p_i$ for all $i$ and the probabilities sum to one, $\sum_{i=1}^n p_i = 1$. This random variable cannot be observed. Instead, we observe $Y$ which is a sum of outcomes, $$ Y = \sum_{i=1}^{k} X. $$ For every sample from $Y$, we also known the number of elements in the sum $k$. How can we estimate   * the set of outcomes $o_1,...,o_n$, and    * the corresponding probabilities $p_i, ..., p_n$,  from a sample of $Y$ denoted by $(y_1, k_1), ..., (y_m, k_m)$ where each $k_i$ is the number of elements in the sum? Also, given an estimate, how can we update it as soon as more data becomes available? What do we know and what can we assume:   * the total number of outcomes $n$ is known   * the upper and lower bound of the outcomes $l$ and $u$ is known   * we can place a well-informed prior on each outcome $o_i$   * we have no a-priori knowledge about the probabilities $p_i$


title: Incidence density

text: Suppose there are $6$ people in a population. During $2$ weeks $3$ people get the flu. Cases of the flu last $2$ days. Also people will get the flu only once during this period. What is the incidence density of the flu? Would it be $\frac{3}{84 \text{person days}}$ since each person is observed for $14$ days?


title: Necessary sample size to get at least one object from each of 10 groups

text: I saw this slide with no explanation: http://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CC4QFjAA&url=http%3A%2F%2Fwww.ist.temple.edu%2F~vucetic%2Fcis526fall2004%2Flecture3.ppt&ei=yrSvUp6rG6PIsAS61oLgAw&usg=AFQjCNFVA4ElIWuMdOCF6pqzL_ACl8bEFg&sig2=FCWfx3azhpbxyJOYSdQj5g&bvm=bv.57967247,d.cWc ![enter image description here](http://i.stack.imgur.com/6ouqk.png) With a sample size of 2, there is a .9 probability of getting one object from each of two groups. For a sample size of three, there would be (0.9 * 0.8) chance of one object from each of three separate groups. Things start to get interesting beginning with a sample size of 11. The graph seem to say that a sample size of 27 has a 50% probability of having representation from each of the 10 groups. This looks superficially similar to the Birthday Problem or some sort of reverse probilistic pigeonhole problem if there is such a thing, but I haven't been able to take those ideas anywhere.


title: Estimating carryover constants

text: I am wondering how can I estimate $ \alpha $ in the following carryover recursive estimation: $ \lambda_i = \alpha_\lambda \cdot \lambda_{i-1} + (1-\alpha_\lambda) \cdot y_{i-1} $ I have read Estimating Within- and Between-Method Carryover Evaluated Statistically, but I do not fully comprehend I can use this 'MCSS' algorithm to compute the carryover constant $ \alpha $. I was wondering if there are simpler methods/ preprogrammed solutions to estimate this carryover constant in e.g. EVIEWS or MATLAB.


title: Overfitting naive Bayes

text: I understand naive Bayes is used largely in text classification. However, the number of features tend to outnumber the number of documents. Does this not result in overfitting where the number of parameters outweigh the number of samples. I am trying to classify a set of documents using multinomial naive Bayes but however I have only achieved 80% accuracy using cross validations of 10 folds. Could overfitting be a problem in this? I have about 25,000 features before features subset selection and after feature selection about 5000+ features. However, some of my samples set per class does not exceeds 1000 documents and I have about 7 classes.


title: Proving that the variance observed in time series A is due to the variance observed in time series B

text: I'm testing the hypothesis that the **the variance observed in the number of medical consultations felt from 2009 through 2013 is due to the variance in user fees price** , and not due to something else. I have 3 time series (from Jan. 2009 to Dec. 2013, measured quarterly) corresponding to the following variables (same population):   1. User fees price   2. Ratio of exempt consultations by exempt population   3. Ratio of not exempt consultations by not exempt population In theory, with the exception of the user fee price, the variables that influence 2 and 3 are the same. I've already established that there is a strong correlation between 1 and 3 (Pearson _r_ = -.83; _df_ = 18; _p_ < 0.01). Also, between 1 and 2 I've found a weak, insignificant correlation ( _r_ = -.4). **Question:** Is this enough to statistically prove my theory? If not, how should I approach this? * * * **Notes:** Healthcare demand is influenced by several variables, some of them easily measurable (such as population ageing, number of doctors, mean price of consultation, etc.), but others not so much (such as technology advances in healthcare), which makes it difficult to create a good regression model. In this context, a user fee is a small fixed contribution payed by users for the medical consultation. Some people are exempt from this fee (the poor, children, patients with chronic diseases, pregnant women, ...).


title: PCA for predictors in lmm or glmm

text: Imagine we have the following mixed model               mod <- lmer(y ~ x + z +(1|g), data = dat)      where:               y <- 1:12     x <- 1:12 + rnorm(12)*2     z <- 1:12 + rnorm (12)^2     g <- c(letters[1:4], letters[1:4], letters[c(1,2,5,6)])     dat <- data.frame(y=y,x=x,z=z,g=g)      One problem with this model is that, `cor(x,z)` is very high. So imagine we want to do a PCA with `x` and `z` and then pit `PC1` and `PCA2` to `y`. However this is not a simple `glm` as we have complex `group` structure `g`. Can we do a simple `PCA` on `x` and `z` then include this in `mod` or do we have to do something more complex accounting for the group structure `g`?


title: Standard error computation with interaction term

text: I aimed at studying the effects of an exposure (E) on fetal growth estimated by repeated ultrasound measurements (n = 2 measures per participant in the following example). I used interaction terms of gestational age with exposure, in order to compute the effect of exposure at different age (using `lincom` in Stata). I am surprised because the standard error of the effect estimates computed with lincomest (performed after xtmixed) are similar at 22.5 and 32.6 gestational weeks (both SEs = 8, while Beta = - 1 and - 19 g at 22.5 and 32.6 gestational weeks, respectively). I am wondering if I correctly used lincomest.                   gen expoAG =  expo*AG       /* Gestational Age */     gen expoAGsq =  expo*AGsq  /* squared of Gestational Age */           xi : xtmixed EFW expo AG AGsq expoAG expoAGsq  if go_pen==1 ||n_ident:,  mle     lincomest _b[expo] +  _b[expoAG]*(22.5) + _b[expoAGsq]*(22.5*22.5)  /* Prediction at 22.5 GW */           xi : xtmixed EFW expo AG AGsq expoAG expoAGsq  if go_pen==1 ||n_ident:,  mle     lincomest _b[expo] +  _b[expoAG]*(32.6) + _b[expoAGsq]*(32.6*32.6)  /* Prediction at 32.6 GW */


title: Generating IMA(1,1) series

text: I'd like to generate a series that follows an IMA(1,1) process, where $θ$ is the moving average parameter. I generated the series based on different representations and I got different results, I'm wondering which one is correct? $$ \begin{aligned} d_1 &=μ+e_1 \\ d_t &=d_{t−1}−θ\cdot e_{t−1}+e_t \\ \\ \text{or} \\ \\ d_1 &=e_1 \\ d_t &=μ+d_{t−1}−θ\cdot e_{t−1}+e_t \end{aligned} $$


title: How is the MA part of ARMA solved for?

text: In an AR model the coefficients on the lags can be solved for using least squares. How is the MA part of ARMA solved for? Since the MA part is a sum of white noise terms I imagine that it is not solved using least squares.


title: Citation for Continuous Space Hill Climbing Algorithm pseudocode on Wikipedia?

text: Can anyone provide a reference for the Continuous Space Hill Climbing Algorithm pseudocode in the Wikipedia article on Hill Climbing? The Russell and Norvig text is cited, but they only provide the discrete case. I'm wondering where this particular version of the continuous case is first described. http://en.wikipedia.org/wiki/Hill_climbing


title: Likelihood ratio test for MLE (Markov Chain)

text: As suggested in Calculating log-likelihood for given MLE (Markov Chains) I want to perform a likelihood ratio test for two fitted models (i.e., first and second order markov chains). Simply comparing the resulting log-likelihood values is as suggested in the other thread not enough. I know how to calculate the likelihood ratio, but I am unsure about how to determine the statistical significane. I need the degree of freedoms for both models. What exactly is this in the case of my Markov Chain MLEs. The number of non-zero probabilities in the MLE, or the number of states?


title: Density of robots doing random walk in an infinite random geometric graph

text: Consider an infinite random geometric graph in which the node locations follow a Poisson point process with density $\rho$ and edges are placed between the nodes that are closer than $d$. Therefore, the length of the edges follow the following PDF: $$ f(l)= \begin{cases} \frac{2 l}{d^2} \;\quad l \le d \\ 0 \qquad\; l > d \end{cases} $$ In the above graph, consider the nodes inside the circle of radius $r$ centered at the origin. Assume, at time $t=0$, we place a tiny robot inside each of the mentioned nodes. That is, the density of the robots on the plane is given by: $$ g(l)= \begin{cases} \rho \quad l \le r \\ 0 \quad\; l > d \end{cases} $$ where $l$ is the distance from the origin. The following figure shows an example of the initial placement of the robots. ![example](http://i.stack.imgur.com/0u5n7.png) At each time step, the robots go to one of the neighbours randomly. Now, my question is that: what is the density function of the robots at $t>0$? Is that possible to compute the density function when $t \rightarrow \infty$? Sorry guys, I am by no means a mathematician. Please let me know if anything is unclear.


title: Calibrated p-value

text: Why do people very rarely report the calibrated p-value which is: $$\mathcal{P}(p) = \frac{-ep \log(p)}{1-ep \log(p)}$$ This gives the probability that one falsely rejects the null hypothesis (i.e. Type 1 error). It is a lower bound.


title: Forecast with STL and regressors (using R package forecast)

text: I would like to use the stlf forecast function from the R package forecast (http://cran.r-project.org/web/packages/forecast/index.html) and include regressors in the model. Question 1: This can only be done using "" method="arima" "" - right? Question 2: For the model calibration the parameter "xreg" should work but how can I enter the new data for th forecast? "newxreg" did not work for me. Thank you


title: Probability of an order statistic

text: If I randomize time points $$T=\{266.14646, 107.28526, 108.89631, 593.03129,\\ 118.10284, 425.60470, 39.02817, 291.90210\}$$ into two groups of equal size, then what is the probability that $T_i$ (the $i$-th time point) will be the 3rd order statistic of the second group? Given that $T_i$ follows an exponential with rate $\lambda_i$. I am getting confused to calculate this probability. Especially with what should be the probability that $T_i$ is the 3rd largest value of the 2nd group.


title: Beta regression of proportion data including 1 and 0

text: I am trying to produce a model for which I have a response variable which is a proportion between 0 and 1, this includes quite a few 0s and 1s but also many values in between. I am thinking about attempting a beta regression. The package I have found for R (betareg) only allows values in between 0 and 1 but not including 0 or 1 them selves. I have read elsewhere that theoretically the beta distribution should be able to handle values of 0 or 1 but I do not know how to handle this in R.I have seen some people add 0.001 to the zeros and take 0.001 from the ones, but I am not sure this is a good idea? Alternatively I could logit transform the response variable and use linear regression. In this case I have the same problem with the 0 and 1's which cannot be log transformed. Thanks Jonathan


title: Post-hoc tests on linear mixed model give mixed results.

text: I am quite new to R so apologies if I fail to ask properly. I have done a test comparing bat species richness in five habitats as assessed by three methods. I used a linear mixed model in lme4 and got habitat, method and the interaction between the two as significant, with the random effects explaining little variation. I then ran Tukey's post hoc tests as pairwise comparisons in three ways: Firstly in lsmeans:               lsmeans(LMM.richness, pairwise~Habitat*Method, adjust="tukey")      Then in ‘agricolae’:               tx <- with(diversity, interaction(Method, Habitat))     amod <- aov(Richness ~ tx, data=diversity)     library(agricolae)     interaction <-HSD.test(amod, "tx", group=TRUE)     interaction      Then in ghlt 'multcomp':               summary(glht(LMM.richness, linfct=mcp(Habitat="Tukey")))     summary(glht(LMM.richness, linfct=mcp(Method="Tukey")))          tuk <- glht(amod, linfct = mcp(tx = "Tukey"))     summary(tuk)          # standard display     tuk.cld <- cld(tuk)   # letter-based display     opar <- par(mai=c(1,1,1.5,1))     par(mfrow=c(1,1))     plot(tuk.cld)     par(opar)      I got somewhat different levels of significance from each method, with ghlt giving me the greatest number of significant results and lsmeans the least. When I ran this on a generalized linear model of abundance data from the same study, ghlt was then the most conservative. All the results from all packages make sense based on the graphs of the data. Can anyone tell me if there are underlying reasons why these tests might be more or less conservative, whether in any case I have failed to specify anything correctly or whether any of these tests are not suitable for linear mixed models? It seems from the code to me that the 'agricolae' and 'multcomp' approaches fit their own model rather than the one I have specified, but I found these codes as suggested post-hoc tests on forums.


title: How to compare the strength of two Pearson correlations?

text: I've been asked by a reviewer if the Pearson correlations (r-values) presented in a table can be compared to each other so one can claim that one is "stronger" than another (other than just eyeballing the actual r-values). How would you go about this? I've found this method http://vassarstats.net/rdiff.html but am unsure if this applies.


title: How to analyze interdependent interaction terms of lmer model

text: Assume I test a number of patients repeatedly over time to see how a certain treatment changes their skin conductance in response to a certain colour (`cond`) after 2 months, 4 months, ... etc. I test the skin conductance on the palm, wrist or arm simultaneously to see whether the place matters. I have a longitudinal dataset of the form: `SC` = skin conductance: dependent variable (metric)   `time` = fixed factor with 6 levels (time1, time2, ..., time6) for each measurement in time   `place` = fixed factor with 3 levels (place A, place B, place C)   `cond` = fixed factor with 2 levels (cond1=blue, cond2=red)   `ID` = random factor with subject ID (because not each subject could be tested at each time point) I put the data into the following form:               data$time  <- relevel(data$time, "time1")     data$cond  <- relevel(data$time, "blue")     data$place <- relevel(data$time, "placeA")      and do some model comparison:               model_1  <- lmer(sc ~ time * cond * place + (1 | ID), data)     model_1x <- update(model_1, REML = F)          model_2  <- update(model_1, .~. - time:cond:place)     model_2x <- update(model_2, REML = F)     anova(model_2x, model_1x)     # three-way interaction is not significant (p=0.93) --> leave three-way      #  interaction out and continue with model_2          # remove two-way interactions:     model_3  <- update(model_2, .~. - time:cond)     model_3x <- update(model_3, REML = F)     model_4  <- update(model_2, .~. - time:place)     model_4x <- update(model_4, REML = F)     model_5  <- update(model_2, .~. - cond:place)     model_5x <- update(model_5, REML = F)     anova(model_3x, model_2x) # interaction is significant (p=0.003)**     anova(model_4x, model_2x) # interaction is not significant (p=0.46)     anova(model_5x, model_2x) # interaction is significant (p=0.039) *      In the end I end up with the final model being `sc ~ time + cond + place + time:cond + cond:place + (1|ID), data)` **Here are my questions:**   * How can I look into the interaction terms? If e.g. `place` would have been dropped from the model altogether, looking at the following post hoc test would take the average over the factor `place`, right?              posthoc_test <- glht(model_final, c("condred == 0", "condred + time2:condred == 0",                        ..., "condred + time6:condred == 0"))      Because, however, the factor `place` is still in the model, the above post hoc test is applied only to the baseline level of `place`, i.e. `placeA` What if I want to see how `time` and `cond` interact, regardless of `place`? Is this even possible given that `place` is itself "captured" in the interaction `cond:place`?   * How do I report the significant interactions? Is the p-value I get from the model selection procedure (i.e., p=0.003 for `time:cond`) also the one I can later report for `time:cond`?  P.S.   By now I'm almost convinced that the answer to question 2 is, that the values I get from the model comparison are the ones I can report for `time:cond` in general. In short: it doesn't matter that I got the values from a model comparison instead of getting them from a preformulated ANOVA. Please correct me if I'm wrong!


title: Benchmark Data-sets for Concept Drift where important predictors (independent variables) change with time or stream of observations

text: I'm currently searching the web and literature for streaming classification datasets with concept drift. I've found a number of synthetic datasets where over time the important predictors either change in their "predictive" nature. For example here is a paper with a synthetic data-set with four blocks of step style concept drift. I've also been checking out the references on the wiki article on concept drift. My question is.. what other _REAL_ streaming datasets out there express concept drift in classification problems. In particular I'm looking for datasets where the set of important features completely changes with time. I've already simulated my own dataset which I vaguely discuss here. Any help is appreciated.


title: Difference between series with drift and series with trend

text: A series with drift can be modeled as $y_t = c + \phi y_{t-1} + \epsilon_t$ where $c$ is the drift(constant), and $\phi=1$ A series with trend can be modeled as $y_t = c + \delta t + \phi y_{t-1} + \epsilon_t$ where $c$ is the drift(constant), $\delta t$ is the deteministic time trend and $\phi=1$ Both series are $I(1)$ and I think both exhibit an increasing behavior. If I have a new series that exhibit an increasing behavior, how do I know this series is a series with drift or with trend? Can I do two ADF tests: ADF test 1. Null hypothesis is the series is I(1) with drift ADF test 2. Null hypothesis is the series is I(1) with trend But what if for both tests, the null hypothesis is not rejected?


title: What is "ANOVA"?

text: Wikipedia says: > Analysis of variance (ANOVA) is a collection of statistical models used to > analyze the differences between group means and their associated procedures > (such as "variation" among and between groups). In ANOVA setting, the > observed variance in a particular variable is partitioned into components > attributable to different sources of variation. In one way or two way ANOVA, (which are the types of ANOVA I have seen), the input variable $X$ is categorical, whose value represents the group of the sample. in The Elements of Statistical Learning: Data Mining, Inference, and Prediction. by Trevor Hastie, Robert Tibshirani, and Jerome Friedman, ANOVA seems to be a way of modeling $E(Y|X)$ as sum of functions of various number of components of $X$. It doesn't relate to variance or partition of variance into some form. $X$ is not necessarily categorical either and may be continuous valued. Or am I missing something? Thanks! ![enter image description here](http://i.stack.imgur.com/Nm7HA.png)


title: equation after running a regression analysis

text: I performed a regression analysis about the impact population increase, precipitation, and difference in registered voters in 2012 versus 2008 against the total of ballots cast in 2012 versus 2008. The ratio showed Y= -889.8x - 11837 R2= 0.0751 Does that mean the independent variables greatly affected the dependent variable? Despite an increase in pop and reg voters, noting the rain each county received, NJ still experienced a major decrease in ballots cast after Hurricane Sandy and that is what I am trying to depict here. Thank you, Connie Grams


title: Strange result of post-hoc test

text: I have data for a test on three groups. The measured variable is ratio scaled. The R code is               g1a<-c(7, 3, 40)     g2a<-c(1,1,2)     g3a<-c(0,0,0)      Since the sample is small and normality cannot be guaranteed, I run a Kruskal Wallis test to check for significance:               l<-list(g1a,g2a,g3a)     kruskal.test(l)      The p-value is 0.02336, which is nice. Now I run a post-hoc test, using the Mann-Whitney U:               wilcox.test(g1a,g2a,paired=FALSE,exact=TRUE)     wilcox.test(g2a,g3a,paired=FALSE,exact=TRUE)     wilcox.test(g1a,g3a,paired=FALSE,exact=TRUE)      All the resulting p-values are above 0.05 (0.07652, 0.0636, 0.05935). This is very strange. Shouldn't one of these tests give a much lower p-value? Especially since I'd have to use some sort of correction to account for the multiple comparisons in the post-hoc test. In other words: how can I interpret this result?


title: Structural equation modeling using Stata/SE 12.1 for weighted, non-normal data with missing values

text: I am currently working with a dataset where the data are weighted, contain missing values and violate assumptions of multivariate normality. Specifically, I have built a structural model using Stata - and am wondering if I have attended to these issues properly by using the following techniques.   1. ADF was not an option for me because my sample size was less than 1000.   2. Therefore I used the Quasi Maximum likelihood technique with robust standard errors to account for the non-normality. I did not use the mlmv option because the Stata 12 manual suggests that when observed variables do not meet assumptions of normality, ml, qml or adf are more appropriate.   3. This method (QML) casewise deletes 20% of my data - thus 80% of the responses had a value on all variables and were thus retained. Note: Prior tests determined no dominant patterns of missingness. My questions: A. Have I sufficiently dealt with the problems (non normality, missing data) or are there better methods that I can use in stata or other programs that can better handle weighted data with these issues? B. Is there any way to incorporate observed constructs that are binary?


title: How can I perform post-hoc tests to the ANOVA performed during K means clustering?

text: If I have highly dimensional data and perform k means clustering, I can end up with clusters that differ more on some axes than others. SPSS provides the option to return an ANOVA for each variable, which shows if the variance is significant between clusters. Is it possible to break this down into pairs of clusters? For instance, I might have 8 variables. If I end up with 3 clusters, I can see that there is a significance difference between them with regards to several variables. What I can't see is whether this effect is driven primarily by the difference between cluster A and B, B and C or A and C. Is there a way to do this in SPSS? If not, what tools will allow me to? (I have dabbled in R previously.)


title: Small sample, many variables

text: I am working on some applied research. I have a small non-random sample ($n=90$), and many variables $>50$ IVs (mix of nominal, ordinal, interval) and 6DVs. The data collection instrument itself was obviously broad (but items chosen were due to known relationships or theorised relationships in the past) so as to make the most use of a time-poor population. It was also envisaged that groups of items would represent a common construct (no scale currently exists). What I am left with now, however, is a wealth of data for such a small sample. In an effort to minimise the data, I have conducted principal components analysis on those variables (checking factorability statistics) that were thought to tap common constructs and produced composite variables. These were checked for internal consistency, some of these were reliable $a>.70$, some were not. Composite scores were averaged to apply equal weights (conservative as no prior research done using these previously). I would like to know how the frequency of performance or use of IVs relates to methods of communication (DVs). Ideally, I'd like to know how each IV relates to the DVs. However, multiple testing correction is obviously an issue here. But if I correct, my p value will be so small that I will have no chance of finding a significant effect. I feel as though I am already taking a 'shortcut' by using ordinal measures (likert) as interval (under direction) - another reason for the creation of composite scores. **My question/s are (finally): What should I do now in light of my original goals (correlation/regression)? Should I run them all, not correct and emphasise its exploratory nature OR should I select a small number of variables (probably the composites for higher reliability) and develop a new, smaller model to test?** I obviously want the result to be as rigorous as possible, not misleading in it's use of and reporting of statistics and make some sort of contribution to my area. However, every corner I turn I hit a wall. The majority of literature in my area often uses likert as interval and rarely reports multiple testing correction, but that doesn't make it sit better with me.


title: How to compare two sides of patients with multiple measures per patient?

text: Let's say we have $n$ patients and we did $l_i$ measures of some parameter at the left side of the patients and $k_i$ measures of the same parameter at the right side of them. So we got a table like this:               i    left side                  right side     1    0.1  0.2  0.9  0.5 -0.3    0.5  1.3  0.9 -0.6  1.3  1.3 -0.3     2    -0.4                       -1.2 -0.9 -1.1 -0.9  0.8     3    -0.7                       0.2     4    -0.1 -0.2  0.8             -0.2     5    -0.6 -0.3                  1.1 -1.3      $i$ denotes the number of the patient, i.e. each line refers to another patient. The numbers of measurements of each patient were chosen randomly. We do not have any information about the distribution of the values. The values in the table above are rounded, in reality we have more precise values, so rounding won't be an issue. We want to know whether there is a significant difference between both sides. Normally we could do a paired wilcoxon-test to compare both sides, but in this case I doubt we can do that because of having patients of which we have multiple measures.


title: Density on $\mathbb{R}^{n}$.

text: We consider the model of nonparametric regression: $$Y_{i}=f(X_{i})+\xi_{i}$$ with $f:[0,1]\rightarrow \mathbb{R}$, $\xi_{i}$ are i.i.d with density $p_{\xi} $ with respect to the Lebesgue measure on $\mathbb{R}$ and $X_{i}$ are deterministic. Somebody can explain me why $P_{j}$(the distribution of $Y_{1},\ldots,Y_{n}$, for $f=f_{nj}$) admits $p_{j}(u_{1},\ldots,u_{n})=\Pi^{n}_{i=1}(p_{\xi}(u_{i}-f_{nj}(X_{i})))$, $j=0,1$ as density with respect to the Lebesgue measure on $\mathbb{R}^{n}$.


title: Calculating standard error for a Normal population

text: I'm a programmer, with a decent but not-expert knowledge of stats, and I'm working through these instructions for how to create funnel plots. I understand from a previous question that the standard error in these instructions is being calculated using `p(1-p)` (in step 4), because these are events and therefore follow a Bernoulli distribution. However, if the data was about individuals in the population, rather than events, how would I calculate the standard error? I know it would be a Normal distribution and therefore standard error would be the square root of the standard deviation, but what formula would I calculate that, given the data I have in the example? (i.e. a set of samples with denominator and numerator population). Thanks for your help.


title: Statistics - hypothesis testing

text: (1) (3 points each) John has come up with a unique new remedy for insomnia. He knows that in the absence of his remedy insomniacs sleep an average of 250 minutes a night with a standard deviation of 30 minutes. He believes that his remedy can increase the average minutes slept per night to 300 minutes and that his remedy will not change the standard deviation of the minutes slept. For each participant, he will measure the average number of minutes slept per over an extended period. He wishes to test the null hypothesis that his remedy has no impact on average minutes slept vs. the alternative hypothesis that it increases the average minutes slept at level of significance 0.05. He believes it is biologically impossible that his remedy could decrease the average number of minutes slept. (a) Suppose John will recruit 20 people to be treated with his remedy. What is the type II error rate of John’s test? (b) How many individuals should John recruit for his study in order to be sure to have at least 0.90 power to reject the null hypothesis? (c) Suppose John is incorrect and his remedy can only increase the average minutes slept by insomniacs to 280 minutes. Now how many total subjects does he need in order to have 0.90 power?


title: Is the p-postulate true?

text: The p-postulate is the notion that equal p-values provide equal evidence against the null hypothesis. Royall, R. N. (1986). The Effect of Sample Size on the Meaning of Significance Tests, The American Statistician, 40:4, 313-315   Wagenmakers, E.-J., Lee, M. D., Lodewyckx, T., & Iverson, G. (2008). Bayesian versus frequentist inference. In H. Hoijtink, I. Klugkist, and P. A. Boelen (Eds.), Bayesian Evaluation of Informative Hypotheses, pp. 181-207. Springer: New York.


title: R: Is it possible to estimate the poisson noise?

text: I have a dataset of many discrete counts (RNAseq read counts per base), which contain both real signals and background noise. The noise is random, and should be poisson distributed. What I would like to do is to subtract (or set a threshold for) a baseline that roughly equals the mean of the noise + a CI, based on the data. However, if I just take the mean from the entire data, I include the signal as well. Is there a way to somehow estimate the noise mean, given that one already assumes there is a random noise + specific signal peak distribution in the data? I know there are specific peak calling algorithms, but for my purposes I try to have an estimate of the baseline as well.


title: How can I write up Tukey HSD findings?

text: the result of SPSS is uploaded in image form.![enter image description here](http://i.stack.imgur.com/bFutx.png) This is exactly a part of my HW. And it make me confused for like 1 hour. Use the Tukey HSD procedure in SPSS to determine which levels of the factor differ. State your conclusions briefly. This question has two parts and I just uploaded the result of one part to get some tips or solutions from you guys then finish the left part. Thank you.


title: When would one want to use AdaBoost?

text: As I've heard of the AdaBoost classifier repeatedly mentioned at work, I wanted to get a better feel for how it works and when one might want to use it. I've gone ahead and read a number of papers and tutorials on it which I found on Google, but there are aspects of the classifier which I'm still having trouble understanding:   1. Most tutorials I've seen speak of AdaBoost as finding the best weighted combination of many classifiers. This makes sense to me. What does not make sense are implementations (i.e. MALLET) where AdaBoost seems to only accept one weak learner. How does this make any sense? If there's only one classifier provided to AdaBoost, shouldn't it just return back that same classifier with a weight of 1? How does it produce new classifiers from the first classifier?   2. When would one actually want to use AdaBoost? I've read that it's supposed to be one of the best out-of-the-box classifiers, but when I try boosting a MaxEnt classifier I was getting f-scores of 70%+ with, AdaBoost murders it and gives me f-scores of something like 15% with very high recall and very low precision instead. So now I'm confused. When would I ever want to use AdaBoost? I'm looking for more of an intuitive rather than strictly statistical answer, if possible.


title: Bayesian probabilty

text: I need to know how to find the Bayesian probability of two discrete distributions. For example the distributions are given as follows:               p(y/x) = (start= 100, values = 0.1, 0.4, 0.5);;     p(x) = (start = 150, values = 0.1,0.1,0.1,0.3,0.3,0.1)      and the Bayesian formula is given `p(x/y) = (p(y/x)*p(x))/(summation(p(y/x`)*p(x`)))`. Basically I need to know how to multiply these unequal distribution in python.   Thank you.


title: Preparing survey data for analysis

text: I am preparing a complex survey data for analysis. The survey involved a two- stage stratified cluster sampling. In the first stage, I have 7 strata and 14 sampling units and in the second stage I have 13 strata and 29 sampling units. I have made sure that there are at least 2 sampling units per sampling stratum. However, after declaring the survey design in Stata, the `svydescribe` command yields strata that are duplicates of already existing strata but with single units. The following table depicts the Stata output after submitting the command `svydescribe, final`. ![enter image description here](http://i.stack.imgur.com/epKoz.png) As can be seen from the table, in the second stage of sampling, unit/cluster 7 belongs to stratum 4. However, one observation which belongs to the same stratum and unit/cluster has been considered as belonging to a different stratum and unit/cluster. The same holds with stratum 7 and several other strata. What causes such problems and how are they handled? Update: after serially issuing the commands               egen _svy_tag = tag( PSEUDOPSU PSUSTRATA SSU PSEUDOSSUST FPC1 FPC2), missing     replace _svy_tag = _svy_tag | missing(NWEIGHT)     list PSEUDOPSU PSUSTRATA SSU PSEUDOSSUST FPC1 FPC2 NWEIGHT if _svy_tag        as suggested by StasK I obtained the following output. ![enter image description here](http://i.stack.imgur.com/JuPJA.png) ![enter image description here](http://i.stack.imgur.com/B6xz4.png)


title: Excel's confidence interval function throws #NUM! when standard deviation is 0

text: I have some samples from eight people who all gave the same answer to a question. Now, obviously the sample's mean is the answer all people gave, and the standard dev is 0. Excel throws a `#NUM!` error when I call the function               CONFIDENCE.T(0.05, K33, COUNTA(B33:I33))      where `K33` is the standard dev (0). What would be the correct interpretation of this? Can I even calculate a confidence interval? For clarification: People are asked to give their opinion on a scale of (1, 2, 3, 4, 5), which is ordinal. Nevertheless, one always calculates the arithmetic mean of all judgements (according to ITU-T P.800, see also: Wikipedia), so that's why I also want to get a confidence interval.


title: Quadratic discriminant analysis (QDA) with qualitative predictors in R

text: I need your help with a Statistical Learning homework in R. I have to perform classification over this dataset: mammographic masses predicting Severity (0="not severe",1 = "severe) using these predictors:   * Age (quantitative)   * Margin (qualitative)   * Shape (qualitative) Everything is fine and understandable when I use logistic regression, but I don't know if it's possible to run QDA (or linear discriminant analysis either), since two of the variables are qualitative.


title: Disadvantages of the Kullback-Leibler divergence

text: I'm working on a calibration problem which involves the usage of the Kullback- Leibler divergence as an error between some empirical distribution $p$ and a theoretical distribution $q$. In the model, the $q$ distribution is normal with some fixed parameters. I have two questions:   1. Is the Kullback-Leibler divergence the best f-divergence to consider as error?    2. Does the usage of the Kullback-Leibler divergence entail any kind of issue?


title: Can the EM algorithm be applied to my problem? Input data set is based on a function of parameter

text: I understand EM algorithm is often used for missing data/mixture problem. But can it be used to optimize a particular type of likelihood based on jointly fitting variables and transformations of those variables? I am fitting a logistic regression model: $Y$ is the usual binary dependent variable, $\mathbf{X}'$ is a matrix of independent variables. However, the data $\mathbf{X}' = f_{\alpha}(\mathbf{X})$ are a function of some parameters, $\alpha$. Specifically $X_{.,j}' = \alpha_j^{X_{.,j}}$. We observe only $\mathbf{X}$ and $Y$. This gives the linear model: $\mbox{logit} (Y|X) = \beta_0 + \beta_1 \alpha_1 ^{X_1} + \ldots + \beta_k \alpha_k ^{X_k} $ Is there any other method that jointly estimate this model? The one that I can think of is to use a global search. Steps are:   1. Randomly generate the three alpha parameters within its definition space, fit the model to get coefficient and AUC.    2. Repeat step 1 for 1000 times.    3. Choose the alphas that gives the best AUC. Does this approach make sense? What I am trying to do is to get the optimized parameter. Thank you.


title: Cross validation with first and second order effects

text: Our project involves looking at different chemical properties and developing a model which incorporates first and second order effects to predict properties of new chemicals. The data set consists of 100 obs. and 40 variables. The following approach worked when there were no second order effects involved. How can I include the second order effects and verify the model?               fit1 <- glm(property0 ~ property1 + property2 + property3 + (property4 + property5)^2), data=REORG)     summary(fit1)     step1 <- stepAIC(fit1, direction="both")     summary(step1)     model.1 <- ols(step1, REORG, x=TRUE, y=TRUE)     model.1     val.cross1 <- validate(model.1, method = "crossvalidation", B = 40, bw = FALSE, rule = "aic",                             type = "residual", sls = 0.05, aic = 0, pr = FALSE)     val.cross1


title: Analysing when tunes are alike: testing ordinal and categorical similarity

text: I'd like to compare some Jazz tunes. Jazz tunes are characterized by a succession of chords. The interval between each succeeding chord (number of tones) and the order in which they appear matter. At first, I thought of a Mann-Whitney test that could compare two tunes with a different number of chords using the ranks of their intervals. Another thing matters, though: the quality of the chord. A major chord can follow a minor chord, a dominant one, or a semi-diminished one. So, there is a mix of ordinal and categorical variables of which to take account. What methodology should be used in testing the similarity of two jazz tunes?


title: Practical implication of using constructed sum variables vs. factor scores for regression

text: I have a survey where respondents rate a product on different scales such as design, quality etc, which several questions for each category. I've run a factor analysis and found six meaningful factors from which I've derived factor scores. How would my regression differ if I used these factor scores for a multiple regression rather than sum variables constructed from the question groups, i.e design, quality etc.?


title: Transformation of binary data

text: I have a function that I am trying to optimize using Particle Swarm Optimization. The objective function gets a binary string. These binary strings are candidate solutions of the subject function. I can measure the dissimilarity between two solutions by the Jaccard Coefficient, but I need to move one solution towards another by a factor, meaning I need a vector. Are there any methods to do that? Example:               a = 0011001      b = 0011010      I want to move `a` towards `b` by the factor of 0.4. So I need to find another string that is similarity between `a` is 0.6 and similarity between `b` is 0.4. Actually strings are much longer. I am trying to somehow represent binary data n-dimensional space and then represent a vector from one point to another. I hope I have explained myself clearly.


title: Can we use random forest for classification in combination with distance matrix between classes?

text: With a colleague, we are working on a dataset containing ~5000 continuous variables for 120 individuals belonging to 8 classes. We want to **estimate the relative importance of each variable** to explain the classes. We have used a random forest approach with some success. Now, we could like to go deeper by considering the fact that **the 8 classes we fit are unequally distant from each other**. In fact, in our case **we can _a priori_ generate a distance matrix ( _i.e._ cost matrix) for all possible pairs of classes**. My (very limited) understanding of random forest is that, for regression problems, the error $E$ is computed by the mean square difference between the OOB sample and the prediction for the same sample: $E = n^{-1}\sum\limits_{i=1}^n{{(y_i-\hat{y}_i)}^2}$ Where $y_i$ is the predicted value and $\hat{y}_i$ the real value of an out-of bag-sample $i$. Ultimately, the calculation of the variable importance depends on how the error is computed (right?). In our case, I would like to use a modified loss function, for instance: $E = n^{-1}\sum\limits_{i=1}^n{M_{y_i,\hat{y}_i}}$ Where $M$ is predefined a distance matrix; so $M_{a,b}$ is the distance between class $a$ and $b$. In this way the misclassification error would be more important if $y_i$ and $\hat{y}_i$ represent distant classes and, ultimately, **the variable importance should be more relevant**. My questions are:   1. Does this approach make sense to you, or am I missing something?   2. Can you think of any study that has used something similar.   3. We have so far used the `randomForest` package in R. It does not seem possible to use it in combination with an a priori distance matrix between classes. Do you know if this is already implemented somewhere? **EDIT** I believe this is a very frequent problem in my field, biology, because we deal with classes for which relations can be represented and quantified by trees (dendrograms), often because of there lineage. After some research, it appears that my question is about using a **cost- sensitive** version of random forest. In this respect, it is very similar to this question. I specificity want to use a **cost matrix** rather than a cost vector though. It there any ontological reason why it is not possible or is it simply not implemented?


title: How to use Hartemink's discretization algorithm?

text: From the help documentation of the `discretize` function of the R package `bnlearn`: > Hartemink's algorithm has been designed to deal with sets of homogeneous, > continuous variables Does "homogeneous" here mean that the variables should share a comparable variability ? Then for nonhomogeneous variable is it a good practice to scale them and apply Hartemink's algorithm ?


title: Model selection between exponential and gamma distribution using cross validation

text: I have fitted the exponential and gamma model to my univariate data and obtained the MLE estimates using the R package "fitdistr", now I'm trying to do model selection based on leave-one-out cross validation, are there any readily avaliable R package to do this. Thanks!


title: Effect of delivery mode of questionnaire (online/hardcopy) on results

text: I conducted a survey research where I have used two delivery modes for the questionnaire - online and offline (hardcopy mode). The survey questions measure perceptions of employees on certain organisational behavior variables such as - organisational commitment, organisational justice, perceived organisational support, job satisfaction etc. All items are measured with a likert scale (1-7) The research objective is to test for moderation hypotheses (a total of 15). I observe that different set of variables are emerging as moderators when I run moderation test separately for online and offline(hardcopy) data. How should I proceed? Should I just ignore these differences and focus on overall data or report data according to mode of delivery?


title: Different outputs for linear regression in WEKA and Excel

text: I was double-checking my regression outputs in WEKA and Excel. All coefficients are the same, however I get a considerably different intercept. How come? Also, is it possible to run a t- or p-test or any other significance test in weka (as is included in the Excel regression)? My dependent variable are petrol prices (`Preis`), the independent variables are different weeks, while week 2 (`Woche2`) is my reference category.   * WEKA Result:              Scheme:weka.classifiers.functions.LinearRegression -S 0 -R 1.0E-8     Relation:     Ferien_Magdeburg     Instances:    1600     Attributes:   5                   Preis                   Woche1                   Woche4                   Woche3                   Woche2     Test mode:10-fold cross-validation          === Classifier model (full training set) ===          Linear Regression Model          Preis =                0.0166 * Woche1=1 +           0.0295 * Woche4=0 +           0.0136 * Woche3=0 +           1.4903          Time taken to build model: 0.01 seconds          === Cross-validation ===     === Summary ===          Correlation coefficient                  0.425      Mean absolute error                      0.0265     Root mean squared error                  0.0358     Relative absolute error                 89.6461 %     Root relative squared error             90.4602 %     Total Number of Instances             1600        * Corresponding output in Excel:              Intercept   1.533375          Coefficients          Woche 1 0.0166          Woche 4 -0.0295          Woche 3 -0.013591667          Woche 2 0


title: Validity Index Pseudo F for K-Means Clustering

text: The Validity Index "Pseudo F" is described as: (between-cluster-sum-of-squares / (c-1)) / (within-cluster-sum-of-squares / (n-c)) with c beeing the number of clusters and n beeing the number of ovservations. It measures the seperation between all the clusters and should be high. I want to apply K-Means for a large dataset. (about 3000 observations with about 200 values each). When computing the Pseudo F for say cluster sizes 2 till 20, it decreases continously with the increase in the number of clusters. How should I interpret these values? Possibilities might be.   * The validity index is inappropriate for this type of clustering   * K-Means is not very applicable for this dataset   * The clustering solution with 2 clusters it the best WSS of cluster: ![enter image description here](http://i.stack.imgur.com/B43kV.jpg)


title: Including (and interpreting!) random intercepts and/or slopes in linear mixed models

text: I'm new to linear mixed modeling, and have some theory-driven questions that I'm not sure how to analytically resolve. I am analyzing experimental data with a within-subjects factor (`discount`). My theory hypothesizes that the effect of this within-subjects factor is contingent upon a between-subjects characteristic of respondents (`iipm`). Because my data is in long form, and respondents are making 8 choices over time, I model my data as follows:               library(lme4)     m1 <- glmer(chose ~ iipm*discount + product + (1|id) + (1|time), data=long1,              family="binomial")      All that I'm trying to do here is fit a simple model that accounts for the dependence between observations for a single subject (`1|id`), and the potential effect of making several choices in a row (`1|time`). However, my theory further specifies that this relationship **should not be affected** by the inclusion of other demographic variables in the model (let's say `ideology` and `partisanship`). So, based on some reading I've done (as well as previous answers on this site), I fit the following model:               m2 <- glmer(chose ~ iipm*discount + product + (1 |id) + (1|time) + (1|partisanship) +            (1|ideology) , data=long1, family="binomial")      Because random slopes goes beyond my expertise, I'm just using random intercepts to see what happens when I account for baseline variation amongst individuals attributable to their partisanship and ideology. However, if I **were** to use random slopes to essentially say that the effects of partisanship and ideology vary on an individual basis, even after accounting for baseline variability, I should specify the following model:               m3 <- glmer(chose ~ iipm*discount + product + (1 + partisanship +ideology |id)  +       (1|time) , data=long1, family="binomial")      To test the hypothesis that this baseline variability **_doesn't matter_** , I then run a likelihood ratio test comparing the two models:               library(lmtest)     lrtest(m2,m1) # p=.349     lrtest(m3,m1) # p=.416      If there's no improvement in fit (p>.05), I (very tentatively) interpret this as support for my hypothesis that `demographics` and `ideology` don't matter. Is this a right way to approach the data, or is there a more sophisticated way to test this hypothesis using multilevel modeling? Any expertise and advice is greatly appreciated.


title: Dirichlet processes for supervised learning?

text: It seems when I glance around here at the fashionable learning algorithms, things like neural networks, boosted trees, support vector machines, random forests, and friends are promoted for supervised learning problems. Dirichlet processes and their ilk seem to be mentioned mostly in unsupervised learning problems, such as document or image clustering. I do see them get used for regression problems, or as general purpose priors when one wants to do Bayesian statistics in a nonparametric or semiparametric way (e.g. as a flexible prior on the distribution of random effects in certain models) but my limited experience suggests that this doesn't come as much from the machine learning crowd as it does from more traditional statisticians. I've done a small amount of googling on this and I've found a few definitive uses in machine learning of DPs for supervised learning, but it seems like whenever I need to look up something about (say) hierarchical DPs, whatever paper I find the answer in is using it for unsupervised learning. So, are Dirichlet processes and their cousins most effective as priors for flexible clustering models? Are they not competitive with boosting, SVMs, and neural networks for supervised learning problems? Are they useful only in certain situations for these problems? Or is my general impression incorrect?


title: Modelling when the dependent variable has a "cut-off"

text: Apologies in advance if any of the terminology I use is incorrect. I'd welcome any correction. If what I describe as a "cut-off" goes by a different name, let me know and I can update the question. The situation I'm interested in is this: you have independent variables $\bf{x}$ and a single dependent variable $y$. I'll leave it vague, but assume that it would be relatively straightforward getting a good regression model for these variables. However, the model you're aiming to create is for independent variables $\bf{x}$ and dependent variable $w = min(y,a)$, where $a$ is some fixed value within the range of $y$. Equally, the data you have access to does not include $y$, only $w$. A (somewhat unrealistic) example of this would be if you were trying to model how many years people will collect their pension for. In this case, $\bf{x}$ could be relevant information such as gender, weight, hours of exercise per week, etc. The 'underlying' variable $y$ would be life expectancy. However the variable you'd have access to and be trying to predict in your model would be $w = min(0, y-r)$ where r is the retirement age (assuming for simplicity it's fixed). Is there a good approach for dealing with this in regression modelling?


title: How to find a suitable association of color with data value in a visualization?

text: I'm working on a software project that involves creating a visualizer for flood simulations. As part of this project, I've created a water gradient that shows water depth at particular points. To set what values will represent what colors, I go through the data and get the minimum and maximum values that occur and evenly distribute the colors according to that scale. However, there are often times points in these simulations that have significantly deeper water at them than anywhere else in the simulation. This causes most of the points on the map to have very similar colors and this is not very informative and makes the areas where the water is deeper very hard to see. My goal is to dedicate a larger range of colors to depths that occur more frequently. For example, if depths go from 0 to 12 but most depths are between 1 and 2, I want more color variation to occur within that range than does between say 11 and 12 or 4 and 5. It seems I need to use standard deviation or something involving normal distribution to do this, but I'm a bit fuzzy on how these things work and how I can use them to accomplish my goal. Any help that can be provided will be appreciated. Thank you.


title: Analyse within-subjects data where participants get exposed to some but not all treatments

text: I want to conduct a within-subjects experiment with one independent variable with four levels. However, I ran a few informal pilot tests, and exposing participants to all four conditions seems to negatively affect results on the later conditions. That is, participants get bored of the experiment by the second or third treatment, and this shows in their responses. If possible, I would like to expose participants to only two of the treatments (properly cross-balanced) but have no clue how to analyse the data then as, I think, a normal RM ANOVA would expect participants to have values for all four treatments. * * * NB: This might have been asked before, but the closest question I could find had different underlying assumptions, namely, that different "treatments" (in this case, different questions on a questionnaire) do not have any influence on each other).


title: Comparing incidence rates

text: I want to compare to incidence rate's between two groups (one without a disease and one with). I was planning to calculate the incidence rate ratio (IRR), i.e. incidence rate group B/ incidence rate group A, and then test if this rate equals to 1, and finally calculate 95% CI intervals for the IRR. I found a method for calculation the 95% CI in a book (Rosner's _Fundamentals of Biostatistics_ ): $$\exp\left[\log(\text{IRR}) \pm 1.96\sqrt{(1/a_1)+(1/a_2)}\right]$$ where $a_1$ and $a_2$ are the number of events. But this approximation is only valid for large enough sample sizes and i think the numer of event I have is to small (maybe for the total comparison it's okay.) So I think I should use another method. Im using R and the exactci package and found that I could maybe use `poisson.test()`. But this function has 3 methods for defining the two sided p-values: central, minlike and blaker. So **my questions** are:   1. Is it correct that to compare two incidence rate ratios im using a test for comparing poisson rates?    2. When in use the poisson.test function in R from the exactci package what method is best? The vignette for exactci says: > **central:** is 2 times the minimum of the one-sided p-values bounded above > by 1. The name 'central' is motivated by the associated inversion con dence > intervals which are central intervals, i.e., they guarantee that the true > parameter has less than $\alpha/2$ probability of being less (more) than the > lower (upper) tail of the 100(1-$\alpha$)% confidence interval. This is > called the TST (twice the smaller tail method) by Hirji (2006). > > **minlike:** is the sum of probabilities of outcomes with likelihoods less > than or equal to the observed likelihood. This is called the PB (probability > based) method by Hirji (2006). > > **blaker:** combines the probability of the smaller observed tail with the > smallest probability of the opposite tail that does not exceed that observed > tail probability. The name 'blaker' is motivated by Blaker (2000) which > comprehensively studies the associated method for con dence intervals. This > is called the CT (combined tail) method by Hirji (2006). My data is:               Group A:      Age group 1: 3 cases    in 10459 person yrs.   Incidence rate: 0.29      Age group 2: 7 cases    in 2279 person yrs.    Incidence rate: 3.07     Age group 3: 4 cases    in 1990 person yrs.    Incidence rate: 2.01     Age group 4: 9 cases    in 1618 person yrs.    Incidence rate: 5.56     Age group 5: 11 cases   in 1357 person yrs.    Incidence rate: 8.11     Age group 6: 11 cases   in 1090 person yrs.    Incidence rate: 10.09     Age group 7: 9 cases    in 819 person yrs.     Incidence rate: 10.99       Total:    54 cases in 19612 person yrs.      Incidence rate: 2.75          Group B:      Age group 1: 3 cases    in 3088 person yrs.   Incidence rate: 0.97      Age group 2: 1 cases    in 707 person yrs.    Incidence rate: 1.41     Age group 3: 2 cases    in 630 person yrs.    Incidence rate: 3.17     Age group 4: 6 cases    in 441 person yrs.    Incidence rate: 13.59     Age group 5: 10 cases   in 365 person yrs.    Incidence rate: 27.4     Age group 6: 6 cases   in 249 person yrs.    Incidence rate: 24.06     Age group 7: 0 cases    in 116 person yrs.     Incidence rate: 0       Total:    28 cases in 5597 person yrs.      Incidence rate: 5.0


title: Are there non-uniform random generator directly requiring random seeds?

text: I wonder if my following understandings about random seeds are correct:   1. Uniform random generators require random seeds.   2. Non-uniform random generators depend on random seeds, only through uniform random generators.   3. Or are there (common) non-uniform random generators use random seeds directly, without indirectly via uniform random generators? Thanks!


title: Pseudorandom number generator and pseudo-random number sampling

text: I am curious about some glossaries.   1. Does "pseudorandom number generator (PRNG)" only mean for uniform random number generator? I didn't see one that is not in list of random number generators.   2. Does "Pseudo-random number sampling" only mean for non-uniform random number generator? > Pseudo-random number sampling or non-uniform pseudo-random variate > generation is the numerical practice of generating pseudo-random numbers > that are distributed according to a given probability distribution. Thanks!


title: R square test in matlab

text: I did the code for my $R^{2}$ (R square) test in MATLAB but it is not working accordingly. I want to test the Weibull distribution against my raw data, hence I want to do an $R^{2}$ (R square) test. Below is my code, but the results I obtain are negative. I tried with Gamma and Rayleigh to test, but the code was not even working with other distributions. Could you please verify this code:               sample= N;          h1=histfit(sample,30,'weibull');      xdata1 = get(h1(2), 'XData');      ydata1 = get(h1(2), 'YData');      f = fittype('weibull');      [c2, gof] = fit(xdata1',ydata1',f)      Output > c2 = >   >   >      General model Weibull: >      c2(x) = a*b*x^(b-1)*exp(-a*x^b) >      Coefficients (with 95% confidence bounds): >        a =   6.787e-05  (-0.3653, 0.3654) >        b =       9.961  (-5429, 5449) >   > > gof = >   >   >            sse: 4.9879e+07 >        rsquare: -1.6634 >            dfe: 98 >     adjrsquare: -1.6906 >           rmse: 713.4243 >


title: Binomial probabilities in games when players can lie

text: An oracle can be asked a binary question. If the real answer is _true_ , the oracle says true with probability _tt_ (and false with probability _tf=1-tt_ ). If the real answer is _false_ , the oracle says false with probability _ff_ (and true with probability_ft=1-ff_). So the oracle is defined as P(q) = { tt, ff}. When the oracle is asked _the same question_ _n_ times and answers _true_ _at least m times_ , what is the probability that the real answer is _true_?


title: When should I use FDR or Bonferroni in multiple comparisons?

text: I have the following situations : (1) **Very large number of variables** - say 100 variables and 200 samples (observations). Although level of each variable is just 2 or 3 or 4 max. I am doing linear ANOVA to test the hypothesis. Say, for the first variable, `Ho: l1 = l2 = l3 = l4`, here l1, l2,l3, and l4 are means at levels 1 to 4. The hypothesis is tested in similar way at all 100 variables independently (not a multiple regression rather talking one variable at time). (2) **a single variable with a high number of levels** (lets say 100 levels) Now the null hypothesis is single :               Ho: l1 = l2 = l3 = l4= ..... = l100       In the both situations should I need some short of **multiple comparison** type of corrections such as the Bonferroni correction or the false discovery rate correction? **Why?**


title: Identify significant deviations (gaps) of points in a 2-D distribution

text: Consider an $X_i$, $Y_i$ distribution of points in a 2D plane. How can one assess if a given region in this 2D plane has a significant lack of points when comparing with the rest through a non-parametric test, without making any assumption on the underlying (unknown) distribution? I would be tempted to mirror the distribution against different axes, but I'm sure if I try enough axes I'll end up finding a "significant" difference between the two sides. I'm afraid this just means I'll be affected by multiple hypothesis testing in a way I do not know how to correct for.


title: Standard Errors of Transformed Variables

text: I am carrying out an MLE where some I use a log transformation on the variance parameters which are being optimized. When I calculate the standard errors (se) the se of the transformed variables is massively too high!!! There is a nice explanation of the delta method here: http://www.statlect.com/delta_method.htm I have the jacobian from the optimizer. Given the jacobian and the expression for the delta method, should I take the derivative of the transformed input variable wrt to the input variable? i.e. $\frac{d}{dx} log(x)=\frac{1}{x}$ $(\frac{d}{dx}log(x))^2=\frac{1}{x^2}$ So I should use this transformed version of the jacobian in the se calculations? Baz


title: Post-hoc test for slopes - phia package for R

text: I am trying to do a post-hoc test to make pairwise comparisons between slopes. I found a package called `phia` but is not available for the R version 2.13.2. Does anyone know if there is a similar package for this version of R? Thanks!


title: Extract standard errors of coefficient linear regression R

text: > **Possible Duplicate:**   >  How do I reference a regression model's coefficient's standard errors? If I have a dataset:               data = data.frame(xdata = 1:10,ydata = 6:15)      and I run a linear regression:               fit = lm(ydata~.,data = data)     out = summary(fit)          Call:     lm(formula = ydata ~ ., data = data)          Residuals:            Min         1Q     Median         3Q        Max      -5.661e-16 -1.157e-16  4.273e-17  2.153e-16  4.167e-16           Coefficients:                  Estimate Std. Error   t value Pr(>|t|)         (Intercept) 5.000e+00  2.458e-16 2.035e+16   <2e-16 ***     xdata       1.000e+00  3.961e-17 2.525e+16   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 3.598e-16 on 8 degrees of freedom     Multiple R-squared:     1,  Adjusted R-squared:     1      F-statistic: 6.374e+32 on 1 and 8 DF,  p-value: < 2.2e-16       How do I extract the standard errors of the regression coefficients from either fit or out? I can't seem to figure it out. Thanks!


title: Logistic regression on categorical data

text: I have large dataset (around 2 million records and 300 features) with a lot of missing data. Most of the independent variables are categorical (some of these variables have more than 40 valid values). The outcome is either Y or N. The Y outcome is a rare event: around 98% of outcomes are N. I'm supposed to fit a logistic regression model to these data. I took random sample of them, keeping the same distribution. I am working in R, but I'm new to both R and logistic regression modeling and I have some questions:   1. I applied `factor` to the outcome. Do I need to apply it on every categorical independent variable? I have more than 200 variables, some of them have only 2 valid values while others have 40! Will it affect the size of the data?    2. Is there any advice about attribute selection? Should it be done before fitting the logistic regression model or after, depending on the results?   3. Is it recommended to take biased sample data where the outcome Y is more than the original distribution in the large data?   4. There are fields like `userId`, `groupId`, etc. What type of data do we consider these to be? How to deal with them?   5. What other predictive models are suitable for this kind of data?


title: how does rpart handle missing values in predictors?

text: From the `?rpart` documentation - > na.action : the default action deletes all observations for which y is > missing, but keeps those in which one or more predictors are missing. How does it impute missing values in predictors?


title: Residuals correlated positively with response variable strongly in linear regression

text: I did the multiple linear regression on a dataset of 412 observations, with one response variable (Y) and 25 explanatory variables(X1-X25). Y and most of Xs are not normally distributed. Besides, there are some correlation between several Xs. The plot show that the residuals strongly correlated with Y positively and weakly correlated with fitted Y negatively.(Sorry.As I'm newer in this website, I am n't allowed to post images.) To address these problems, I have tried the principle component regression, the weighted least squared regression and the ridge regression. They all didn't work. I want to know what's wrong with the regression. Why did the residuals correlate with observed Y so obvious?


title: Sampling variance is reduced when removing duplicates -- why?

text: I have carried out a simulation experiment where $n$ out of $N$ items $x_i$ were drawn with replacement with unequal probabilities $z_i$ from a finite population. (If it matters: $z_i = 1/(Nx_i)$ in my case, but perhaps there are more general results.) Then, I estimated the population totals of the indicators $$ y_i := \begin{cases}1: y_i = m \\ 0: \text{otherwise}\end{cases} $$ for each $m$ in two ways:   1. By applying the Hansen-Hurwitz estimator: $\hat{Y}_I = n^{-1} \sum_i (y_i/z_i)$.   2. By first removing the duplicates from the sample and then applying the Horvitz-Thompson estimator for the selection probabilites $\pi_i = 1 - (1 - n/N)^{y_i}$ (which is the probability of selecting the item $i$ at least once, i.e. one minus the probability of not selecting it at all, approximable by $x_i$ Bernoulli trials with success probability $1 - n/N$ each if $x_i \ll n$): $\hat{Y}_{II} = \sum_i(y_i/\pi_i)$. This was repeated 1000 times. I compared the estimates with the true values and observed the ratio between estimate and true value. I can clearly see from my experiments that the second estimate is better in terms of sampling variance. How to support this with theoretical results? From Cochran (1977), Chapter 9, I see the variance of the HH estimator as $$ V(\hat{Y}_I) = n^{-1}\sum_i^N z_i (y_i/z_i - Y)^2 $$ and that of the HT estimator as $$ V(\hat{Y}_{II}) = \sum_i^N \frac{1-\pi_i}{\pi_i}y_i^2 + 2\sum_i^N\sum_{j>i}^N \frac{\pi_{ij}-\pi_i\pi_j}{\pi_i\pi_j}y_iy_j $$ but I have trouble applying this to my case. First, does it mean that for HH the variance depends on $Y$ and for HT it doesn't? Second, I would assume $\pi_{ij} = \pi_i\pi_j$ since the inclusion of items is independent, and then the second term of the HT variance vanishes -- or am I wrong? I appreciate any help.


title: Comparison of cophenetic correlation coefficients on different data sets

text: I applied the same hierarchical clustering (weighted) on two data sets: The first is a 'raw' data set, on which I didn't do anything, and the second on the same data set after I filtered it by removing some items I believe mislead the clustering. I computed the cophenetic correlation coefficients on both methods. But is this relevant? Given that the data sets are linked, but different. Also, I read that 'This coefficient has also been proposed for use as a test for nested clusters'. Is there a method using the cophenetic correlation coefficient to know whether a cluster could be further divided into more clusters?


title: Identifying automotive makes and models in online forums

text: I am working on a project in which I am trying to identify mentions of car makes and models in forum posts. I would like to say what I have tried so far and then ask for any feedback or suggestions of other approaches. So far I have tried two approaches - a named entity recognition method, and a dictionary-based method. This can be thought of as a named entity recognition problem where the named entities are car makes and models. There are many CRF-based tools available for identifying named entities (e.g. Stanford NER tool). The issue is that the pre-trained models that come with these tools were trained to identify NEs like PERSON, LOCATION, ORGANIZATION, etc. but not necessarily car products. Furthermore, these pre-trained models were not trained on forum text which is what I am interested in. I have tried to use my own annotated training data (and gazetteer) but I don't think I have enough. As a result, the precision of the models that I trained on my own data have high precision but very low recall. The other approach that I tried deals with a catalogue of car makes and models. My strategy is to do an exhaustive search through the forum text and match the unigrams and bigrams with catalogue names (using a custom-weighted minimum edit distance measure). I try to reduce the number of searches by removing words that could not be car mentions (i.e. common dictionary words). I am still in the process of refining this model but so far the results show good recall but low precision. I strongly suspect that the low precision is do to the fact that the catalogue that I am using is very large and the vast majority of makes and models in the catalog will never be mentioned. So I need a way to assess the prior probability of each model being mentioned but I haven't figured this out yet. Furthermore, this method is relatively inefficient. What do you think of these approaches? Can you think of any other methods that might be successful? I thought about using Labelled-LDA but I can't quite figure out how it would work for this problem (especially since I don't think that I have enough annotated training data).


title: ISPSS: Cross-sectional time series analysis

text: I am trying to run a multiple regression in SPSS. I am using panel data with 4 independent variables. Of these, one is a dummy variable. Can someone please guide me through the process or give me some hints since I'm a beginner in this. I have 350 companies. These companies have a "score" (dependent variable) for 4 different years. There are 4 different independent variables (1 dummy) that are affecting the dependent variable. When i go to Analyze - Mixed model - Linear, i don't know how to proceed. I don't know what to put as "subjects" and "repeated", and i don't know which variables should be used fixed. ![enter image description here](http://i.stack.imgur.com/qowzp.png)


title: How to avoid 0 determinant when sample covariance matrix has very small values

text: I have $n$, $p$-dimensional vectors and I am constructing the $p \times p$ covariance matrix using the following formula: $$\mathrm{Cov}(j,k) = \frac{1}{n-1} {\sum^n_{i=1}} (x_i(j) - {\mu}(j))(x_i(k) - {\mu(k)}) $$ However, the values in the $p$-dimensional vectors are between 0 and 1, as it is a normalized vector, and $n$ can be large which results in $\mathrm{Cov}(j,k)$ having rather small elements. I am trying to calculate the determinant of the covariance matrix to determine the differential entropy of the underlying multivariate Gaussian. However, because of the small numbers, the determinant is effectively 0 (and in fact rounds to 0 in computation). Are there any methods to avoid this? Perhaps I can avoid the division by $(n-1)$?


title: Bootstrapping t-test for nonparametric data

text: I have a two heavily **skewed** distributions and am trying to use **bootstrapping** to compare their means via the **t-test**. I am concerned about the appropriateness of using the standard error of the original/observed data in the final step when I know that this is not normally distributed. **My process goes like this:**   * Bootstrap - randomly sample with replacement (N=1000)   * Calculate t-statistic for each bootstrap to create a t-distribution: $$ T(b) = \frac{(\overline{X}_{b1}-\overline{X}_{b2})-(\overline{X}_1-\overline{X}_2) }{\sqrt{ \sigma^2_{xb1}/n + \sigma^2_{xb2}/n }} $$   * Estimate t confidence intervals by getting $\alpha/2$ and $1-\alpha/2$ percentiles of t-distribution   * Get confidence intervals via: $$ CI_L = (\overline{X}_1-\overline{X}_2) - T\_{CI_L}.SE_{original} $$ $$ CI_U = (\overline{X}_1-\overline{X}_2) + T\_{CI_U}.SE_{original} $$ where $$ SE = \sqrt{ \sigma^2_{X1}/n + \sigma^2_{X2}/n } $$   * Look where the confidence intervals fall to determine if there is a significant difference in means (i.e. non-zero) I have also looked at the Wilcoxon rank-sum but it is not giving very reasonable results due to the very heavily skewed distribution (e.g. the 75th == 95th percentile). For this reason I would like to explore the bootstrapped t-test further. **So my questions are:**   1. Is this an appropriate methodology?   2. Is it appropriate to use the SE of observed data when I know it is heavily skewed? * * * Possible duplicate: What's the difference between bootstrapping t-test and just using a nonparametric test?


title: Rating system taking account of number of votes

text: I need an idea for a new rating system.. the problem in the ordinary one (just average votes) is that it does not count how many votes... For example consider these 2 cases: 3 people voted 5/5 500 people voted 4/5 The ordinary voting systems just take the average, leading to the first one to be better.. However, I want the second one to get a higher rating, because many people have voted for the second... Any help?


title: Question about foundations of the uniform shrinkage prior

text: I am collecting papers about the uniform shrinkage prior for hierarchical Bayesian model. In "A prior for the variance in hierarchical models" of Michael J. Daniels it is stated at the end of page two that: > An alternative proper prior, shown to give a minimax estimator for the k > means for the above model by Strawderman (1971) under certain conditions, > would be to place a uniform prior on the shrinkage parameter. The mentioned paper by Strawderman (1971) is "Proper Bayes Minimax Estimators of the Multivariate Normal Mean". I do not understand in what sense this article explicates the quote of Daniels (maybe the paper is too technical for me and I missed the point). If somebody has some hints regarding the article or more generally can explain the quoted sentence, I will appreciate it.


title: How to limit impact of a predictor when doing multiple regression?

text: I am doing some multiple regression using R's `lm()` method. I wonder whether there is an easy way to limit the impact of some predictors. I think `lm()` results give too much weight to certain predictors and I would like to somehow cap coefficients of certain predictors. Is it doable?


title: Procedure to sample from an integrated marginal

text: I would like to sample from: $$ p(\theta_2|x)=\int p(\theta_2|\theta_1,x) . p(\theta_1|x) . d\theta_1 $$ Several persons suggest me to use the following procedure to draw $(\theta_2^i)_{i=1:N}$:               for i=1:N         draw $\theta_1^i$ from $p(\theta_1|x)$         draw $\theta_2^i$ from $p(\theta_2|\theta_1^i,x)$      this appears intuitively satifying but I do not see how to show that this procedure is valid. It appears to be linked to Rao-Blackwellisation but it is not clear to me. Thanks for your help. (this question is linked to an older one that I posted Sampling from marginal using integrated conditional)


title: Low probability levels when doing logistic regression

text: I am building a Logistic regression model for a churn problem. When I scored the out of sample data set, I find very low probability levels as the output probability. Conventionally, I would look for .5 as the cut off but this scored population doesnt have many customers above .5 ( say just 1%). Seeing the business cause, we need atleast 5% people to be approached for the impact. I therefore reduced the cut off probability to judge the scored dataset. So now, I am defining .1 probability as the cut off. The model is very good at that level, in that it is perfectly distinguishing my target from non target.   1. Is there any problem with this approach, given that at .1 level, model has very good accuracy.   2. what in general is the cause of low probabilities at scored population level.


title: How do I detect shifts in sales mix?

text: Suppose products A, B and C are sold on markets 1, 2 and 3. I have access to three years worth of monthly historical data. I'd like to detect/foresee possible shifts in sales mix for the three products over the three markets (e.g. I'd like to know if market 1 is buying more of product A than it used to, or if product B is cannibalizing the others on market 3). Does a 5/10/15% increase mean that I'm facing a shift in consumer behaviour? Does it depend on the intrinsic variability of the behavioural pattern? Also how do I choose the optimal time-frame on which to perform such analysis (i.e. how do I know if comparing semester over semester increase makes more sense than comparing quarter over quarter increase)? Would t-test/ANOVA work here? Or maybe a control chart? Or maybe some intervention detection method to check if there was a level shift in my time series? If so, how do I go about it? Any input is greatly appreciated!


title: Categorical variable with a lot of levels and interaction

text: Suppose we have a categorical variable $X_1$ with $100$ levels. Should we **not** test interaction between $X$ and some other variable $X_2$? Because if, for example, $X_1(3)*X_2$ is not significant then the interaction between $X_1$ and $X_2$ is not significant?


title: How to apply coefficient term for factors and interactive terms in a linear equation?

text: Using R, I have fitted a linear model for a single response variable from a mix of continuous and discrete predictors. This is uber-basic, but I'm having trouble grasping how a coefficient for a discrete factor works. **Concept:** Obviously, the coefficient of the continuous variable 'x' is applied in the form `y = coefx(varx) + intercept` but how does that work for a factor z if the factor is non-numeric? `y = coefx(varx) + coefz(factorz???) + intercept` **Specific:** I have fitted a model in R as `lm(log(c) ~ log(d) + h + a + f + h:a)` where `h` and `f` are discrete, non-numeric factors. The coefficients are:               Coefficients:                   Estimate      (Intercept)  -0.679695      log(d)        1.791294      h1            0.870735       h2           -0.447570       h3            0.542033        a             0.037362       f1           -0.588362       f2            0.816825      f3            0.534440     h1:a         -0.085658     h2:a         -0.034970      h3:a         -0.040637      How do I use these to create the predictive equation:               log(c) =  1.791294(log(d)) + 0.037362(a) + h??? + f???? + h:a???? + -0.679695      Or am I doing it wrong? I THINK that that concept is if the subject falls in category `h1` and `f2`, the equation becomes:               log(c) =  1.791294(log(d)) + 0.037362(a) +  0.870735  + 0.816825  + h:a???? + -0.679695      But I'm really not clear on how the `h:a` interactive term gets parsed. Thanks for going easy on me.


title: how to group time series data into stationary & non-stationary?

text: i have time series data (download at:-->http://www.4shared.com/rar/ft3YCN2H/Time_series_data.html) consider potential values as time series data points(4th column) `file<-read.table(filename.txt,skip=1)` `ts_data<-file$V4 #final time seris vetor` please, **give R script** /functions to be used from any package,for follwing   1. to cluster these time seris into 2 groups--stationary & non-stationary(by programming & not by visually seeing ACF/PACF graphs)   2. how to find best model(AR/ARMA/ARIMA/SARIMA/GARCH,etc) to fit over these data for future forecasting ?,i don't want manually for each file,kindly suggest automated script   3. how to transform non-stationary models to get best fit model(i dont know how to get best P(autoregressive order),Q(moving average),D(differencing order))


title: Is it possible to directly read CSV columns as categorical data?

text: I need to analyze with R the data from a medical survey (with 100+ coded columns) that comes in a CSV. I will use rattle for some initial analysis but behind the scenes it's still R. If I _read.csv()_ the file, columns with numerical codes are treated as numerical data. I'm aware I could create categorical columns from them with _factor()_ but doing it for 100+ columns is a pain. I hope there is a better way to tell R to import the columns directly as factors. Or to at least to convert them in place afterwards. Thank you!


title: Entering dummy variables and other covariates into the regression model

text: I am controlling for education-level (two dummy variables) and another covariate in my regression model. I'm conducting a hierarchical regression analysis. I was confused about how I should be entering all the predictors. Do I enter the dummy variables and other covariate in the first block and then the main predictors in the next block? Or do I enter the dummy variables in the first block, the other covariate in the second block and the main predictors in the third block? Thank you!


title: Poisson vs Negative Binomial

text: In the case of count data where my variance is 6.34 and mean is 5.44, which would be a better fit to the data - Poisson or negative binomial? And what will be the assumptions ?


title: Visualization software for clustering

text: I want to cluster ~22000 points. Many clustering algorithms work better with higher quality initial guesses. What tools exist that can give me a good idea of the rough shape of the data? I do want to be able to choose my own distance metric, so a program I can feed a list of pairwise distances to would be just fine. I would like to be able to do something like highlight a region or cluster on the display and get a list of which data points are in that area. Free software preferred, but I do already have SAS and MATLAB.