title: LDA with a categorical predictor variable

text: I have a dataset with a categorical treatment (1 if treated, 0 if control) and many columns (34) of response variables. Each column represents a species and its response (some measured abundance) to the presence or absence of the treatment. Each row is a separate site; exactly half of them were treated and the other half were controls. The experiment was replicated over two years. I am interested in finding out which of these species respond to the treatment, so a LDA was suggested, but I am confused because my treatment is categorical and all of the examples I can find of LDA use continuous predictor variables. So the data look like: Treatment(0/1) Site(ID#) Year(1/2) Species 1 Species 2 ... Species 34 Can I use a LDA for this question? If not, what sort of ordination is capable of accounting for categorical predictor variables?


title: Analysis of predictors of change (outcome measured at two time points)

text: I have data where a Normally-distributed outcome was collected at two time points. Covariates were also collected at both time points. Some of the covariates (e.g. sex) didn't change between the two time points. I am interested in modelling the data to see whether changes in predictor variables were associated with changes in the outcome variable, adjusting for other covariates. I am using Stata.


title: Forecasting rain from weather conditions

text: Please bear with me - I have a very limited knowledge in statistics. Here is my case (simplified):   A weather pattern is derived from 1 year of data. We are monitoring:     * sky conditions (sunny/overcast)      * wind direction (E/W/N/S)    Goal: We are trying to predict probability of "rain within 1hr".   Observations:     * if sky="overcast" AND wind="East": 18 out of 20 times it rains within 1hr. - p=0.9   * if sky="overcast" AND wind="South": 15 out of 20 times it rains within 1hr. - p=0.75 My problem is that there is only 20 occurrences on each pattern that gives me very low confidence (approx. +-22% at 95% confidence lvl) - apparently the predominant wind in the area is N/W :(   **Question:** Can I re-formulate my statement/pattern like this:   if sky="overcast" AND (wind="East" OR wind="South") 33 out of 40 times it rains within 1hr. e.g. p=0.83?   If this is the case now that I have 40 occurrences would that give me better confidence interval, say approx. 15% at 95% cl? Thank you in advance


title: Regression Tree when target is a ratio

text: I am learning a regression tree for data of the form $(x_i,y_i)$: $x_i = (1, 0, 1, ...., 1 , 1)$ a multiple input vector and $y$ is a ratio of the number of observations divided by the number of trials for that input. so $y_i = (obs/trials)$. My problem is that most of the data is sparse, meaning that most of the inputs will have 0 observations and some trials and then its corresponding $y_i$ will be zero (let's say no input has trials =0). My problem is that when I learn a tree, in each leaf the tree calculates the average of $y_i$ and that is not what I need since all the $y_i=0$ are being discarded in the average. What I would like to do is to sum up, the number of trials in each leaf, and the number of observations and make the division for each leaf instead of averaging the ratios (or even each one by itself (observation, trials)). I am using sci-kit for the tree model. Is there any intelligent way to do what i want without implementing it from scratch. I mean suing tree libraries as sci-kit?


title: Significant difference

text: I have two answer pools for the same survey. I want to compare these answers that came from 2 different groups. Or I want to somehow see if these two groups share the same ideas about asked questions.#


title: Background subtraction for signal and error analysis

text: I use a CCD to see the split of a energy level due to Zeeman effect. I have a 1 dimensional CCD of 7926 pixel of 7μm each one. My CCD analyze a region 2 dimensional, and then it steps forward 200 times. So, i have a matrix like this !["matrix"](http://i57.tinypic.com/xn5f7s.png). I select the background between 40 and 60 pixel and i calculate for each row the mean value. Then, i project these 7926 mean values over the Y axis. So, for each bins of the row of the matrix(i.e. 3040) , i subtract the value of the content of 3040 pixel of the background This is my background: !["background"](http://i.stack.imgur.com/pZUaX.png) (in x there are the pixels of my CCD and y is the intensity). I have to subtract the background from the signal, subtracting bin for bin. At each bins of the background it is associated an error $\sigma = \sqrt{N}$ where $N$ is the bin content. Each bins of the signal has an error $\sqrt{N}$, but when i subtract the backgroun how can i calculate the exact error of the pure signal? If $Z'_i$ is the bin content after the subtraction, $Z_i$ the bin content before and $Y_i$ the bin content of the background, then $Z'_i = Z_i - Y_i$ and the error (throw the propagation) $\sigma_{Y_i} = \sqrt{\sigma_{Z'_i}^2 + \sigma_{Y_i}^2 +....}$ In fact, I think, there should be a term of covariance. How can i compute it?


title: Algorithms for prediction of consumption, based on previous data

text: I have data for a few stores which sell apples. For each store I have an averaged value of how many kilos of apples users have bought per day for this month. So it looks like this:               Store 1:     -------     January - average of 10 kilos of apples per day     February - average of 18.7 kilos of apples per day     March - average of 24.5 kilos of apples per day     April - ...          Store 2:     --------     January - average of 3.23 kilos of apples per day     February - average of 2.9 kilos of apples per day     March - average of 7.89 kilos of apples per day     ...      For each store I want to make a prediction - how many kilos of apples per day will be sold in the next few months (4, 5 months in the future). The problem is that sometimes the data is like this: ![enter image description here](http://i.stack.imgur.com/PPnDP.jpg) (it is increasing and maybe it's good to use the method of least squares fitting to fit the data values with line and find the function) but other times the data is like this: ![enter image description here](http://i.stack.imgur.com/lo5j1.jpg) and the least squares won't work :( What I would like to ask is:   * to make the prediction, is it always true that I should find a line or a curve to fit the data best or that's not always the case, and the prediction will be very wrong if I use this approach   * for the second type of graphics, what algorithms can I use :(   * how to decide which algorithm to use to fit the data best and if I should use a line for fitting, or a parabola or a third degree polynomial and so on The ideal solution I want to write is to have one algorithm which takes the data and decides which algorithm (from a set) will give the best prediction and then give the data to the chosen algorithm for predictions. (For example give data like (2, 3, 4, 5, 6, 7) and the first algorithm will say "use the method of least squares" and then I will give the data (2, 3, 4, 5, 6, 7...) to the method of least squares to get the next 4-5 values from it.) I know it's a big question, but any kind of information will be very useful! Thank you very much in advance!


title: Statistics to determine the no. of questions/participants in a Quiz

text: I would like to construct a quiz which has N questions (e.g. N=300) and M participants (e.g. M=1000). Each participant would be asked a subset of N questions (e.g. sub=20). In the end, for each question we would get X amount of answers from different participants. How to calculate the optimal ratio of participants/questions, so that the results are meaningful? Is there a particular formula in statistics that takes into acount these parameters? Sorry that the question is not formulated more formally, statistics is not my strong suite. Thanks for your help, any pointers to a web page which explains this or even a search phrase would be appreciated! :) Best, Tamara


title: Oversampling with categorical variables

text: I would like to perform a combination of oversampling and undersampling in order to balance my dataset with roughly 4000 customers divided into two groups, where one of the groups have a proportion of roughly 15%. I've looked into SMOTE (http://www.inside-r.org/packages/cran/DMwR/docs/SMOTE) and ROSE (http://cran.r-project.org/web/packages/ROSE/ROSE.pdf), but both of these create new synthetic samples using existing observations and e.g. kNN. However, as many of the attributes associated with the customers are categorical I don't think this is the right way to go. For instance, a lot of my variables such as Region_A and Region_B are mutually exclusive, but using kNN the new observations may be placed in both Region_A and Region_B. Do you agree that this is an issue? In that case - how do one perform oversampling in R by simply duplicating existing observations? Or is this the wrong way to do it?


title: Conditional distribution in a Bayes net

text: I have this Bayesian network where variable C is dependent on variables A and B. I want to know how come $ P(B|C) = \frac{P(B)}{P(B)+(1-P(B))P(A)} $


title: Probability to find a modern human with two teeth in different developmental stages

text: Each teeth growths from the crown to the root. There are different stages previously described to divide this process, as Crown initiation, one half of the crown, crown complete, root one quarter, root and a half, etc. Each of these stages of each different tooth (incisors, canines, molars, etc) has its own typical age of attainment (years) defined by the mean and the sd. For example: Incisor 1/4. This stage, defined by the acquisition of one quarter of the root, is attained in modern humans at a mean age of 5.28 years (SD=0.91) Molar root complete. This stage, defined by the complete root lenght formation of the first molar, is attained in modern humans at a mean age of 8.45 years (SD=1.36). You can see the normal distribution in the next figure:![Incisor, black; Molar, red](http://i.stack.imgur.com/Nkyq8.png) **Question 1:** Which is the probability to find a modern human with these stages in both teeth? The answer must be in the overlapping area below both normal distributions, but I don't remember the procedure to solve it. **Question 2:** Imagine we have one individual with an age of 5.5 years. Which is the probability that this individual at this age would have the incisor at that stage (1/4)? And which is the probability this individual would have first molar at stage of Root complete?


title: Clinical trial design

text: I have a clinical trial phase 2 3 stage that I need a design for: 3 groups of 12 blood transfusion done at different Hemoglobin levels. It is not a dose escalation study because the transfusion is the same amount of blood for all groups. I need a study design and a statistical test suggestion.


title: Is this mixed ANOVA correctly specified in R?

text: I am trying to test some results I have for significance. It has been recommended that I use R and I am completely new to this. ### Set-up:   * Groups: two groups of 8 subjects (16 total)   * Two conditions: alert and passive   * Measurements: responses for three different stimuli (A, B, and C) measured in each condition ### Experiment: Testing the order of conditions   * Group one: Alert A, B and C followed by Passive A, B, and C    * Group two: Passive A, B, and C followed by Alert A, B, and C Stimuli A, B and C are randomly interleaved in the experiment, does this matter in my ANOVA? I am interested in making a between and within group comparison of responses to A, B, and C ### Here is what I am doing: My data is arranged in the following way               Group Subject Condition Stimulus  Response     One      S1    Alert          A    _Value_     One      S1    Alert          B    _Value_     One      S1    Alert          C    _Value_     One      S1    Passive        A    _Value_     One      S1    Passive        B    _Value_     One      S1    Passive        C    _Value_     One      S2    Alert          A    _Value_     ...          Two      S9    Alert          A     _Value_     Two      S9    Alert          B     _Value_     Two      S9    Alert          C     _Value_     ...      This is the code I used:               My_anova = aov(Response ~ Condition * Stimulus * Group +                 Error(Subject/Condition * Stimulus ), data=My_Data)      ### Question   * Is this model correctly specified in R?


title: What is the impact of windowing function on time series

text: Greeting I would like to know what is the impact of windowing functions like Hanning,... on a time series. Is it possible to finde anomalies using windowing functions? **EDIT** I have a time series, which contains integer values (470,471,472,472,473,471,...) I try to find anomalies in this vector of value. As I have understood, with Wavelet transform , FFT, and moving average . As I know any discontinuity between the 1st and last time sample (i.e., a jump in the value ) will contribute to the frequency spectrum in FFT, but we should use windowing to detect these anomalies. It means, I should calculate FFT for my time series and then I should run a window function like hanning or etc. on the result to find anomalies, Is it right?


title: Understanding expectation of data points estimators

text: I have a small doubt in the way estimators are defined. We have a sequence of data points $X^{(1)}, X^{(2)}, ... , X^{(n)}$. We then define an estimator as a function of these data points. When calculating the bias of an estimator we find the expectation of this function of "random variables". But are $X^{(1)}, X^{(2)} ... $ all random variables? These are just observed values. I know the definition says they are samples coming from some unknown distribution (at least the parameter of the distribution is not known). How can we possibly compute $E[X^{(i)}]$ when we do not know what the distribution is? Also I tend to see each $X^{(i)}$ as an observed data value of the random variable and not a random variable. What is it that I am missing? Many proofs that I have seen (like the one that proves $\bar{X}$ is unbiased) end up using the iid assumption to get rid of the $E[X^{(i)}]$. None of them talks about calculating expectation explicitly. Any more examples to show what is meant by calculating the expectation of values of random variables and saying that these are iid could be helpful in understanding the concepts. Update: One way I am beginning to think of it is the following: We have a set of n data points corresponding to a random variable. Now instead of looking at these as n distinct samples of one random variable, define n iid random variables having the same distribution as the original random variables and consider these to be the n values of these n random variables. Am I thinking in the right direction?


title: How to perform ANCOVA on mixed experimental design where covariate is one level of a condition using R?

text: This is a second question related to previous post: Is this mixed ANOVA correctly specified in R? To reiterate, the setup is as follows:   * Groups: two groups of 8 subjects (16 total)   * Two conditions: alert and passive   * Measurements: responses for three different stimuli (A, B, and C) measured in each condition Experiment: Testing the order of conditions   * Group one: Alert A, B and C followed by Passive A, B, and C    * Group two: Passive A, B, and C followed by Alert A, B, and C ### New Question:   * If I wanted to use the responses measured in the passive condition as a co-variate to the alert responses, how would I perform an ANCOVA analysis in R? I am not sure how to handle this considering the repeated factors and between subject measures.


title: Overmatching bias and confounding variables

text: As I understand it, matching is one way to identify causality in observational studies. By matching observations that are "similar" and comparing ones that did or did not receive treatment, you can consider this as a quasi-experiment of sorts. What is overmatching? What kind of bias does it introduce? I have mostly seen matching from an economics perspective, but have recently seen some papers on epidemiology suggesting that "overmatching" can result in bias. I find it hard to understand the terminology of the papers and would greatly appreciate if someone could help explain some of the main concepts. Below is an article that references the idea: > Over-matching can cause bias. BMJ. 2002 August 10; 325(7359)


title: What is a good index of the degree of violation of normality and what descriptive labels could be attached to that index?

text: ### Context: In a previous question, @Robbie asked in a study with around 600 cases why tests of normality suggested significant non-normality yet the plots suggested normal distributions. Several people made the point that significance tests of normality are not very useful. With small samples, such tests don't have much power to detect mild violations of normality and with large samples, they will detect violations of normality that are sufficiently small not to be of concern. It seems to me that this problem is similar to the debate around significance testing and effect sizes. If you only focus on significance tests, when you have big samples, you can detect small effects that are irrelevant for practical purposes, and with small samples you don't have sufficient power. In a few instances I've even seen textbooks advise people that you can have "too large" a sample, because small effects will be statistically significant. In the context, of significance testing and effect sizes, one simple resolution is to focus on estimating the size of effect of interest, rather than being obsessed with the binary decision rule of whether there is or is not an effect. Confidence intervals on effect sizes is one such approach, or you could adopt some form of Bayesian approach. Furthermore, various research domains build up ideas about what a given effect size means in a practical sense, for better or worse, applying heuristic labels such as "small", "medium", and "large effect". This also leads to the intelligent recommendation of maximising sample size in order to maximise accuracy in estimating a given parameter of interest. This makes me wonder why a similar approach based on confidence intervals of effect sizes is not more widely espoused in relation to assumption testing, and normality testing in particular. ### Question:   * What is the best single index of the degree to which the data violates normality?   * Or is it just better to talk about multiple indices of normality violation (e.g., skewness, kurtosis, outlier prevalence)?   * How can confidence intervals be calculated (or perhaps a Bayesian approach) for the index?   * What kind of verbal labels could you assign to points on that index to indicate the degree of violation of normality (e.g., mild, moderate, strong, extreme, etc.)? The aim of such labels could be to assist analysts with less experience in training their intuition of when violations of normality are problematic.


title: Mice function in R freezing

text: I am attempting to impute my data for missing values using the `mice` function (from the mice package) in R. However everytime I run the function it freezes or lags. I have tried running it over night, and it still does not finish. I am working with 17000 observations across 32 variables. I can't even stop the function after it starts, I have to force quit the program. Here is my code:               imputed.data = mice(data, m = 1, diagnostics = F)


title: A/B testing in Python or R

text: A/B testing: http://20bits.com/articles/statistical-analysis-and-ab-testing/ http://elem.com/~btilly/effective-ab-testing/ I'm not too familiar with A/B testing, but I was wondering if there were any specific packages/libraries in R or Python that can be used to perform A/B testing.


title: Normalization to non-degenerate distribution

text: I am reading de Haan's Extreme Value Theory (2006). In the discussion of distribution of sample maximum, he said "in order to obtain a non-degenerate limit distribution, a normalization is necessary". Then he gave the following example. "Suppose that there exists a sequence of constants $(a_n)>0$ and $(b_n)$ such that \begin{equation} \frac{\max \{X_1, \cdots, X_n\} - b_n}{a_n} (1) \end{equation} has a non-degenerate limit distribution as $n \to \infty$, i.e., $$\lim_n F^n(a_nx + b_n)=G(x), (2)$$ for every continuity point $x$ of $G$, where $G$ is a non-degenerate distribution function." And he also commented that this is a linear normalization. I have three questions here.   1. What does it mean to normalize to a non-degenerate distribution function, please? In my past studies, normalization means to find the constant $c=\frac{1}{\sqrt{2\pi}}$ such that $\int_\mathbb R c e^{-\frac{x^2}{2}} = 1$. It appears that normalization means different things in de Haan's book.   2. What do the two sequences $(a_n)$ and $(b_n)$ mean here, please? Or what role do they play, please? Why is $(1)$ equivalent to $(2)$, please?   3. What are common non-linear normalization, please? Thank you!


title: Compute the probability of a lower value in the elements of a timeseries

text: The original question was confusing, so I've edited it to describe the core of the problem. I have a large time-series (300K items) of integers between 0 and 9. Any two consecutive integers are distinct:               0, 1, 2, 3, 0, 1, 0, 1, 0, 1, 2, 3, 4, 5, 1, 0, 1, 2, 3, 0, 1, 0, 1, 0, 1, 2      I want to compute 10 probabilities: for each possible value (0 through 9), the probability that the next element in the time-series is smaller than this value. For the time-series above, the probabilities would be:   * for 0: 0 for the next value to be smaller   * ...   * for 3: 0.66 for the next value to be smaller   * ...   * for 5: 1 for the next value after to be smaller   * ... Then, I would like to see how that probability changes through the time-series (let's say over a window of 10000 elements) - a rolling probability. As a bonus, I would also like to compute the same probabilities, but for pairs of values: the probability of the next item being smaller than 2 if the current value is 2 or 3. Can you also describe this problem in formal language so that I can do more research on my own?


title: calculate the rate of change

text: I posted this on maths, but seems it would be better here :S http://math.stackexchange.com/questions/49941/calculate-the-rate-of-change Basically I am trying to calculate the change frequency for a set of data. Each bit of data has the date-time it was created. I would like to say for a specific set of data the change frequency is hourly, daily, weekly, monthly or yearly. So far I have tried getting the list of dates and get the min/max which is easy to calculate an average from which can be converted into a human readable label such as hourly, daily etc How would i take into account the age of the last new bit of data. eg: say there were 50 dates all roughly an hour one after the other. This is hourly. but if the last one was 2 weeks ago, its not quite hourly. In this example I am not sure myself what the frequency would be of the list (hourly, daily, weekly, monthly or yearly) so I'm looking for a bit of direction. Maybe someone here has done this before and has a good model or knows a bit more than me :) Thanks


title: Easy resources for funnel plot

text: I'm just facing an issue of finding an easy tutorial or guide to the `funnel plot`. Are there simple resources you recommend for this? Thanks.


title: Structural Equation Model design through AMOS, interpretation help

text: I have calculate the estimates through Maximum Likelihood procedure with model of six latent variables(presented in big ellipse) with respective observed variables (presented in rectangle) in AMOS 20.0. The relationship between latent and observed variable represent factor loadings which should be greater than 0.5 in my case. Now issue is perceived playfulness two factor loadings are problematic. One is negative (PP_2: -0.14) and other one is too low (PP_1: 0.20). What is the reason for these abnormal factor loadings. If i remove these two factor loadings problem still remains my PP_4 gets negative. I don't have justification grounds in terms of reference paper or supportive evidence? From the diagram below, a clear illustration of my problem is available. A better explanation would be highly appreciated. ![enter image description here](http://i.stack.imgur.com/DR6ED.png) Note: Three factor loadings are compulsory with a latent variable. In this case, I cannot reduce pp_4 or PP_5 or PP_3


title: Implementing Pettitt test in R

text: I'm trying to implement Pettitt test in R following papers like this pdf (pp. 5 & 6), or this pdf. But, I'm misunderstanding something, because having tested it with some data, I think that output is not correct. Here is the code:               pettitt <- function(x, alpha=0.99) {     # Pettitt AN. 1979 A non-parametric approach to the change point detection.     # x is a vector     # alpha, integer, level of significance     x <- na.omit(x)     o <- rank(x)     s <- c()     L <- length(x)     for (i in 1:(L-1)) {           s <- c(s, 2*(colSums(as.matrix(o[1:i]))) - (i*(L+1)) )     }     vc <- sqrt((-1) * log(alpha) * (L^3 + L^2)/6)     output <- list(abs(s), vc)     return(output)     }      Testing with `larain` and `tempdub` dataset from `TSA package`:               library(TSA)     data(larain)     data(tempdub)     pettitt(larain)     [[1]]       [1]  78 118 180  76  30  30 144  90 124 148 224 334 314 298 362 444 356 334      [19] 300 302 194 121  83  55  45  57  25  95 175 195 193 287 181 231 175 213      [37] 301 331 421 345 392 322 282 354 372 274 194 130 188 248 175  97  85 153      [55] 105 171 181 189 245 297 401 375 449 557 467 551 594 576 602 490 406 354      [73] 262 266 362 248 244 214 208 200 247 147  89  13   9  15  97   5   9  83      [91]   3  95 123  63  31  12  44   6  48  34  72 108 208 164 170 282 214 148     [109] 202 140 104   6 102  86          [[2]]     [1] 50.69224          > max(pettitt(larain)[[1]])     [1] 602          pettitt(tempdub)     [[1]]       [1]  83 161 226 235 164  60  80 169 220 219 188  74  57 177 266 281 228 147      [19]  19  82 125 140 102  41 100 197 235 254 233 141   1  97 144 153 112  26      [37]  73 206 255 258 235 137  28  49  98 101  46  29 149 252 281 274 247 160      [55]  43  70 115 126  79  22 157 248 317 328 287 224  96  27  86  79  27  82      [73] 225 348 407 406 351 256 125  10  58  77  32  61 200 314 381 386 353 216      [91] 124  40  35  70  35  36 173 302 365 386 321 242 131  10  51  38  19 146     [109] 241 319 342 359 330 223  89  45 113 144 111   2 123 228 280 275 250 177     [127]  34  50  89 102  59  22 131 248 334 359 302 198  73  46  83 100  73          [[2]]     [1] 70.96777          > max(pettitt(tempdub)[[1]])     [1] 407      I don't know if I lost something in pettitt test or there are error in my code.


title: Regression discontinuity and regression to the mean

text: Regression discontinuity designs aren't susceptible to regression to the mean. I'm not quite sure why, though. The best I can come up with is that regression to the mean is a phenomenon that occurs because you're dealing with a sample measured at different points, which have different underlying distributions at each point. With RD designs, however, you're measuring the difference in intercept along a column separating two smooth regression functions within a single bivariate distribution. Can anyone shed further light on this?


title: Variable with high number of categories

text: I would like some help with my analysis. I have an outcome variable that represents average number of errors somebody experiences per hour. It has about 26% zeros. 85% of values are smaller than 1. There are some outliers going up to >700. My explanatory variables are only categorical variables. One has 3 categories, one has 2 categories, one has six categories and one has 91 categories (area within a country). My question is whether anyone has some tips analyzing this dataset. I want to know the effects of the variables and potential interaction effects on the outcome. I would like to know as well if there are some locations that experience significantly high numbers of errors. Can anybody help me going in the right direction? I am just starting in statistics.


title: is it possible to use kernel function to calculate each instance of covariance matrix? If Yes why?

text: I saw a paper that uses Gaussian kernel for calculation at covariance matrix of given variables. Is it mathematically Correct or not? if it is okay,what is the intuition behind ? What about using any distance function instead of Gaussian kernel?


title: Significance difference test between two simple slopes within three level categorical predictor

text: I've come across related forums, but can't seem to fit them to my problem -- I have significant interaction in a multi-level framework between a three-level categorical predictor (condition assignment - one of which is a control condition and the other two experimental conditions) and a two-level categorical predictor. I would like to compare the simple slopes of the two experimental conditions, while ignoring the control. I've already run three separate simple slopes tests on each of the three conditions to see if they significantly differed from zero using the following formula (repeated three times), where each condition was dummy-coded while acting as the reference group for its particular simple slopes test:               omnibus.cc <- lmer(trust~(1+group|participant) + DC_cc1*group + DC_cc2*group,      data=study2_dummycodes_simp.effects)      Where I'm running into trouble though is coming up with the most appropriate coding system and also which significance test to use to compare the two simple slopes of my experimental conditions, as mentioned above. Any help would be much appreciated.


title: Calibrating a model using different timesteps, mean/skew

text: Suppose I have 1 year of hourly streamflow values - approx. 8760 values. I am developing a model that outputs an estimated streamflow value for each of the 8760 values. I examine the sum of squared errors (for different predicted and observed outputs) and try to minimize it during calibration. What differences may result from these three different focuses during calibration:   1. I calibrate the model using the differences between the observed and predicted values at the hourly timestep, examining every timestep individually.   2. I calibrate the model using the differences between the observed and predicted values at the daily timestep, summing the hourly predicted and observed values for each day, examining every day individually.   3. I calibrate the model using the differences between the observed and predicted values at the weekly timestep, summing the hourly predicted and observed values for each weekly, examining every week individually.   4. I calibrate to the mean of the 8760 values only.   5. I calibrate to the mean and skew of the 8760 values only.


How can we use Chi-Square test for categorical variables? For example, when we have {True,False} values for the variable under test. Is Chi-Square test suitable in this case? is there any other alternatives? Here is an example of my data. I have two adjacency matrices {observed, expected} for the same nodes, there is "1" if there is an edge between two nodes (i,j) and "0" otherwise. I want to apply Chi-Square test for these two matrices.


title: Interaction in multiple regression

text: I am a psychology undergrad student and I'm not the best at stats. My thesis is due soon but I'm having trouble with my analysis despite extensively searching for an answer. It would be helpful to me if any answers are given in simple terms! My dependent variable is victim blame and is continuous scale data. My independent variables are gender (male vs female) and scenario (scenarioA vs scenarioB vs scenarioC) as well as hostile sexism (continuous scale), benevolent sexism (continuous scale) and frequency of pornography use. It has been suggested that I carry out a hierarchical regression but I am having a great deal of trouble with this. My hypothesis is that victim blame scores will differ depending on what scenario was read and that this difference will be moderated by benevolent sexism scores. Does this involve creating an interaction term of scenario x benevolent sexism? Is it necessary to first dummy code the scenario variable? Or is regression not the right approach at all? I'd reeeeeaaaally appreciate some help with this! Many thanks.


title: Reshaping a distribution

text: Not sure what the exact term is for what I'm trying to do. I have a data set with random variable x with values X1, X2, ..., XN that has a standard deviation sigma and a mean m. I want to perturb the values of X1, X2, ..., XN so that the new distribution has a standard deviation of sigma' and mean m'. I would like an algorithm (or an implementation in Java) that does that. Thanks.


title: Correct terminology for describing relative confidence interval

text: My team manages a survey that has been reporting for many years as "relative sample error" a range of figures which are in fact the Xs in "+/- X percent" of an estimate where X is chosen to make up a 95 percent confidence interval. So in fact they are about 1.96 times the relative standard error. The users - who are not statisticians - appreciate these numbers because they can easily turn them into a confidence interval (which would take them one more step than they wanted if we quoted actual standard errors). My question is, what is the correct way to describe these figures?   * "Relative standard error" is certainly wrong;    * "1.96 * RSE" is clumsy and not quite correct (because they might be created with a bootstrap method which does not calculate the confidence interval that way)   * "Relative half confidence interval" is probably correct but not very plain English   * "95 percent confidence interval, expressed as plus or minus percentage of the estimate" is correct but verbose Putting aside any questions about how these are calculated, any views on good terminology?


title: How to generate more than two correlated series in R?

text: I'd like to generate more than two (i.e 3) correlated series, take an example when the series follows an IMA(1,1) process, at first I want to generate three correlated random errors `e[i]` and then I use the following equation to build the series: `d[i,] <\- mu + d[i-1,] - theta*(e[i-1,])+e[i,]` this is what we need to do for two series, how can we develop it for more than two series?               library(mvtnorm)     rho <- 0.5     mu <- c(10,10)     phi <- c(0.2,0.8)     theta <- c(0.3,-0.7)     d <- ts(matrix(0,ncol=2,nrow=50))     e <- ts(rmvnorm(50,sigma=cbind(c(1,rho),c(rho,1))))     for(i in 2:50)       d[i,] <- mu + phi*d[i-1,] - theta*(e[i-1,]+e[i,])     plot(d)


title: Is PCA appropriate where n<p?

text: This question is an extension of one I asked a few weeks back, here is the link: PCA and small number of observations I will restate that I am interested in the use of PCA in situations where _n = < p_; and generally am only interested in using the first few PC axes for descriptive purposes or as "synthetic" variables that reduce several dimensions into one. My question today revolves around a text "Numerical Ecology", 3rd edition by Legendre & Legendre. On page 450, they state "A full-rank dispersion matrix **S** [variance-covariance] cannot be estimated using a number of observations n smaller than or equal to the number of descriptors p. When _n = < p_, since there are n-1 DF in total, the rank of the resulting **S** matrix of order p is _(n-1)_. In such a case, the eigen-decomp. of **S** produces (n-1) real and _[p-(n-1)]_ null eigenvalues. Positioning n objects while respecting their distances requires _(n-1)_ dimensions only. A PCA where _n = < p_ produces (n-1) eigenvalues larger than 0 and the _(n-1)_ corresponding eigenvectors and principle components." In other words, I believe they are implying that it is OK to use PCA on a dataset where n =< p as long as you are only interested in using (n-1) or fewer of the PCs (as I am). I am interested in your opinion regarding this (their claim and my interpretation) if you have one; and would appreciate any additional literature that might corroborate this claim. Thanks ahead of time for your time and expertise.


title: Hypothesis test for comparing two F-measures

text: I have a classifier that classifies each sample point as positive or negative. Suppose I am evaluating it on a limited data set where I have ground truth and then computing the F-measure (F1-measure). In fact, I have two classifiers, and one has a higher score under the F-measure: the first classifier has a higher F-measure than the second. Is there a way to test whether this higher score is statistically significant? Is there a hypothesis test to test for statistical significance? I suppose the null hypothesis is roughly that the difference in F-measures can be explained by random chance, or something, but I'm not sure exactly how to formulate it. Is there a way to formulate this as a hypothesis testing problem and then apply some hypothesis test to determine whether the apparent difference appears to represent a real difference in the underlying abilities of the classifiers? For purposes of this question, I suppose we might as well assume that the F-measure is a reasonable way to evaluate the accuracy of the classifiers (maybe I don't have good data on the underlying base rate of positives and negatives, or on the relative cost of a false positive vs a false negative). If we take that for granted, is there any good way to compare the performance of two classifiers using the F-measure, or is that fundamentally impossible? I probably should also mention that both classifiers were evaluated on the same set of samples, so we effectively have paired data.


title: Partial, adjusted or Type-III SS of factorial experiment in R

text: The partial SS for two factor (N and P factors) factorial experiment with interaction can be calculated as: \begin{eqnarray*} \textrm{SS}\left(N_{i}|\mu,P_{j},\left(NP\right)_{ij}\right) & = & \textrm{SS}\left(\mu,N_{i},P_{j},\left(NP\right)_{ij}\right)-\textrm{SS}\left(\mu,P_{j},\left(NP\right)_{ij}\right)\\ \textrm{SS}\left(P_{j}|\mu,N_{i},\left(NP\right)_{ij}\right) & = & \textrm{SS}\left(\mu,N_{i},P_{j},\left(NP\right)_{ij}\right)-\textrm{SS}\left(\mu,N_{i},\left(NP\right)_{ij}\right)\\ \textrm{SS}\left(\left(NP\right)_{ij}|\mu,N_{i},P_{j}\right) & = & \textrm{SS}\left(\mu,N_{i},P_{j},\left(NP\right)_{ij}\right)-\textrm{SS}\left(\mu,N_{i},P_{j}\right) \end{eqnarray*} I coded these partial SS for two factor factorial experiment with interaction in R. I'm not getting the right answer and would appreciate if you point out my mistake. Thanks               # Partial SS     y <- c(55, 56, 57, 53, 54, 55, 51, 52, 53, 61, 62, 63)     Trt <- gl(  n = 4, k = 3, length = 4 * 3               , labels = c("N0P0", "N0P1", "N1P0", "N1P1")               , ordered = FALSE)     N <- gl(n = 2, k = 6, length = 2 * 6             , labels = c("Low", "High")             , ordered = FALSE)     P <- gl(n = 2, k = 3, length = 2 * 6             , labels = c("Low", "High")             , ordered = FALSE)     Data <- data.frame(y, Trt, N, P)     Fit1 <- aov(formula = y ~ N * P, data = Data)     ANOVASummary1 <- summary(Fit1)     print(ANOVASummary1)          library(MASS)          X <-        cbind(           model.matrix(object = y ~ 1,        data = Data)         , model.matrix(object = y ~ -1 + N,   data = Data)         , model.matrix(object = y ~ -1 + P,   data = Data)         , model.matrix(object = y ~ -1 + N:P, data = Data)         )          H  <- X  %*% ginv(t(X) %*% X) %*% t(X)          X1 <-        cbind(           model.matrix(object = y ~ -1 + N ,  data = Data)         , model.matrix(object = y ~ -1 + P,   data = Data)         , model.matrix(object = y ~ -1 + N:P, data = Data)         )     H1  <- X1  %*% ginv(t(X1) %*% X1)  %*% t(X1)     SSM1  <- t(y) %*% (H - H1) %*% y     SSM1                    X2 <-        cbind(           model.matrix(object = y ~  1,       data = Data)         , model.matrix(object = y ~ -1 + P,   data = Data)         , model.matrix(object = y ~ -1 + N:P, data = Data)         )     H2  <- X2  %*% ginv(t(X2) %*% X2)  %*% t(X2)     SSM2  <- t(y) %*% (H - H2) %*% y     SSM2          X3 <-        cbind(           model.matrix(object = y ~  1,       data = Data)         , model.matrix(object = y ~ -1 + N,   data = Data)         , model.matrix(object = y ~ -1 + N:P, data = Data)         )          H3  <- X3  %*% ginv(t(X3) %*% X3)  %*% t(X3)     SSM3  <- t(y) %*% (H - H3)  %*% y     SSM3          X4 <-        cbind(           model.matrix(object = y ~  1,     data = Data)         , model.matrix(object = y ~ -1 + N, data = Data)         , model.matrix(object = y ~ -1 + P, data = Data)         )          H4  <- X4  %*% ginv(t(X4) %*% X4)  %*% t(X4)     SSM4  <- t(y) %*% (H - H4)  %*% y     SSM4          print(ANOVASummary1)      **Edit** The `print(ANOVASummary1)` gives the ANOVA. I want to get the SS in ANOVA using Linear Model approach. Using the Linear Model approach I got the correct Sequential SS but could not figure out how to get the Partial, Adjusted or Type-III SS ( See here and here).


title: Does isomap give possible new bases with respect to importance order as PCA?

text: I am running isomap for feature reduction and I want to see what number of features are best for classification with experiments at different number of features. Do I need to run isomap for each try to get specified number of features or just run the algo. once and chose again and again the number of features I need? Is it like PCA? PCA gives possible new bases with their relative importance (eigenvalues) and you might choose the bases incrementally and test classification performance so you do not run PCA again and again for each experiment.


title: 95% CI for an estimated X given Y in a simple linear regression model

text: Considering a simple linear regression model Y= beta0+beta1 x X, with beta0 and beta1 computed, I have to estimate the expected X given a new Y and 95% confidence intervals. I used the formula X=(Y-beta0)/slope. How can I compute in R the 95% interval for the calculated value of ind, given a new value of the height? > head(df)                  ind height     1 4.27    174     2 3.60    159     3 3.61    175          summary(lm(df$ind~df$height))          Call:     lm(formula = df$ind ~ df$height)          Residuals:      Min       1Q   Median       3Q      Max           -0.56263 -0.27596  0.03866  0.26632  0.55440           Coefficients:              Estimate Std. Error t value Pr(>|t|)        (Intercept) -0.968895   1.512371  -0.641  0.52903             df$height    0.027871   0.008985   3.102  0.00562 **          Residual standard error: 0.3287 on 20 degrees of freedom     Multiple R-squared: 0.3248,     Adjusted R-squared: 0.2911      F-statistic: 9.622 on 1 and 20 DF,  p-value: 0.005621       # I tried:               pred.frame <- data.frame(ind=seq(3,5,0.25))      predict(bclm,int="c",level=0.95,data=pred.frame)         fit      lwr      upr     1  174.5780 169.3146 179.8414     2  166.7696 163.6419 169.8973     3  166.8862 163.7806 169.9917     ...............


title: Need to analyze: Unequal variances, Factorial Completer Randomized Design

text: I'm trying to analyze data that violates the assumptions of an ANOVA (homogeneity of variances), thus I think I have to conduct a non-parametric test. My study follows a factorial completely randomized design: 3 sources 4 treatments 4 replicates per treatment*source combination (replicates are equal in size) Dependent variable is seed germination. The Kruskal-Wallase test was suggested to me, but is there a way to do it and test for both source, treatment, and source*treatment interaction? Furthermore, can I draw pair-wise comparisons between source*treatment combinations much like Tukeys LSD would do? Could I use PROC GLIMMIX in SAS instead? If so, what are the assumptions that need to be met for that model and how do I make sure that I'm meeting them? Thanks!


title: Is it necessary to correct alpha in repeated measures ANOVA before any post-hoc comparisons?

text: I have a cross-over design where a sample of subjects undergoes 7 different interventions, and within each intervention there are 3 time points (pre- intervention, post-intervention1, post-intervention2). I firstly did an overall ANOVA with factors time x intervention, to check whether there were overall differences between interventions at any time point tested. For the within-intervention analysis, I first did a repeated measures ANOVA (alpha = .05) with factor time (3), with post-hoc multiple comparisons (pairwise t-tests) Bonferroni corrected. However, during review in a peer-reviewed journal, a reviewer claimed that my "time" effect for a certain target intervention was not significant (p = 0.028; < 0.05) because it wasn't corrected for multiple comparisons (which I only used in the post-hoc comparison). Is the reviewer's point valid? Do I need to correct alpha in the RM ANOVA? I've been trying to find a relevant explanation in "the interwebs", but I can't find anything conclusive. So far, the best I've come upon is this extract of a 2001 paper by Bender & Lange: > "Methods to adjust for multiple testing in studies collecting repeated > measurements are rare. Despite much recent work on mixed models [38,39] with > random subject effects to allow for correlation of data, there are only few > multiple comparison procedures for special situations. It is difficult to > develop a general adjustment method for multiple comparisons in the case of > repeated measurements since these comparisons occur for between-subject > factors (e.g., groups), within-subject factors (e.g., time), or both. The > specific correlation structure has to be taken into account, involving many > difficulties. If only comparisons for between-subject factors are of > interest, one possibility is to consider the repeated measurements as > multiple endpoints and use one of the methods mentioned in the previous > section. However, if the repeated measurements are ordered, this information > is lost by using such an approach."


title: Item information in IRT with item covariates (linear logistic test model)

text: Short question: How does the calculation and interpretation of IRT item information and test information change in the presence of item properties? Long question: There's a variation on IRT called the linear logistic test model (LLTM): $logit(Y_{p,j}) = \theta_p + \sum_{k=1}^K q_{j,k} \alpha_k$ For persons, $\theta_p$ is a random effect across persons $p$, just as in 1PL IRT. But unlike 1PL IRT, items have properties and each property is treated as a covariate. There are $K$ possible properties, and each item $j$ is coded with property values in the vector $q_{j,k}$. The effect of each property is the weight $\alpha_k$. For example, if your test has math problems and reading problems, one item property may be an indicator of whether the item is a math problem or a reading problem. Suppose the properties include indicator variables for each item $j$, i.e., that the LLTM item properties are a superset of the 1PL IRT item properties. That means we have a per-item effect, a.k.a. "difficulty" in IRT parlance. Knowing the difficulty of each item allows us to compute item information, and summing up information across item tells us the information of the full test. In the presence of item properties, can we still talk about some questions being more informative than others? How does the calculation and interpretation of IRT item information and test information change in the presence of item properties?


title: Predict binary outcome with R

text: I have a table includes the following data:               Status     | Age | Tenure | Function | Gender | Race | Grade Level     -----------------------------------------------------------------     Active     |  34 |    1.2 | PSO      | Female | White| 26     Terminated |  24 |    0.2 | Finance  | Male   | Asian| 32     Active     |  50 |    4.0 | HR       | Female | Black| 28     Terminated |  23 |    2.9 | Sales    | Male   | Hispa| 20     :     :     :      I want find out how are the different variables have impact on the binary outcome - Status, and possibly create a predictive modeling to predict the outcome of Status based on the other variables. Where do I start? I am thinking of doing a factor analysis to find the P-value for each factor. Is that a right path?


title: what can I say about the following kernel estimation?

text: I have a discrete random variable $X$ and two random vectors $\vec{Y}$ and $\vec{Z}$. For a given $x$ (i.e. an instance of $X$), I am interested in estimating: $$E[\vec{Y} \vec{Z}^{\top} | X = x]$$ based on samples $(x_1,y_1,z_1),\ldots,(x_n,y_n,z_n)$. I thought of using a kernel estimation method, and estimating the above as: $$E[\vec{Y} \vec{Z}^{\top} | X = x] \approx \sum_{i=1}^n K(x_i, x) y_i z_i^{\top}.$$ My question is: is there any mention of that in the literature? I am looking for conditions under which this estimation is "consistent" (i.e. it converges to the true $E[\vec{Y} \vec{Z}^{\top} | X = x],$ and maybe some anchor that I can look for something similar being done. The closest thing I could find is kernel density estimation, but it is not quite the same.


title: Most important journals in data mining/ML, NLP and IR?

text: Can you please provide with me with the names of the most important journals in data mining, machine learning, natural language processing and information retrieval?


title: What are the must-read papers on data mining and machine learning?

text: What do you think are the papers that each graduate student in machine learning and data mining should be aware of?


title: Star Coordinates vs. principal component analysis

text: I currently preparing a presentation for a university course in "Visual Data analysis". And one of my topics is the "Star Coordinate" visualization. Star Coordinates As Star Coordinates perform a transformation of high-dimensional data, and the well known PCA technique does so too, I wonder whether PCA can be mimicked by Star Coordinates? I think of rearranging the coordinate axes in a way that they represent a linear combination of the original variables? But this is only an idea. Can someone confirm or disprove this?


title: Difference between logistic regression and neural networks

text: How do we explain the difference between logistic regression and neural network to an audience that have no background in statistics ? Thank you.


title: How to find the rows that meet some conditions in a sequence data set

text: Here is the summary of my sequence data generated in the SPELL format. All sequences are supposed to have the same length of 1440, but the summary tells me that they are not the same (see "min/max sequence length:91/1440"). I want to find the rows:   1. that do not have the 1440 sequence length    2. that have missing states (*)    3. that have void elements (%) **Do you have any suggestions?**               > summary(diary.seq)        [>] sequence object created with TraMineR version 1.8-3        [>] 4971 sequences in the data set, 4160 unique        [>] min/max sequence length: 91/1440        [>] alphabet (state labels):            1=Cr (Cr)            2=Ec (Ec)            3=Eo (Eo)            4=Hb (Hb)            5=Hc (Hc)            6=He (He)            7=Hr (Hr)            8=Hs (Hs)            9=Ir (Ir)            10=Lt (Lt)            11=Or (Or)            12=Pb (Pb)             ...        [>] dimensionality of the sequence space: 27360        [>] colors: 1= 2= 3= 4= 5= 6= 7= 8= 9= 10= 11= 12= 13= 14= 15= 16= 17= 18= 19= 20=        [>] symbol for missing state: *        [>] symbol for void element: %


title: how to use global weights to update local ones

text: I have a population of data points. Each data point expresses a different name. Each such name belongs to a different ethnicity, say {A, B, C, D, E}, with a probability. The sum of probabilities for a given name equals to 1. We also have some global weights; the probability a random picked name belong to each of the five ethnicities. Let's assume that these probabilities are: 0.10478, 0.12615, 0.23906, 0.25961, 0.27040. What I want to do is use the global weights so as to update \ re-weight \ re-balance the local ones. For instance, if Mike is more common in the US (with prob 40%) but our population is mainly consisted of US citizens (with prob 50%), then all local US probabilities (for all names) should go down a bit and all the others should go a bit up. The idea sounds pretty easy, I believe, but I am really struggling in finding out how to do that (keep in mind that the updated probabilities should also sum up to 1). Any ideas?


title: Random walk with momentum

text: Consider an integer random walk starting at 0 with the following conditions:   * The first step is plus or minus 1, with equal probability.    * Every future step is: 60% likely to be in the same direction as the previous step, 40% likely to be in the opposite direction  What sort of distribution does this yield? I know that a non-momentum random walk yields a normal distribution. Does the momentum just change the variance, or change the nature of the distribution entirely? I'm looking for a generic answer, so by 60% and 40% above, I really mean **p** and **1-p**


title: How to interpret a predictor with a positive structure coefficient and a negative standardised coefficient in discriminant function analysis?

text: I am doing a discriminant function analysis and I have four continous independent variables and one categorical dependent variable (that has 3 groups). I have chosen to do this analysis to see how these groups can be predicted by the four independent variables.   * Firstly, when we are interpreting the results, do we look at standardised coefficients or the structure matrix coefficients? From my reading of the literature, I have found contradictory opinions. Some say that structure coefficients are better since they are more stable while others say that both standarized and structure coefficients should be looked at.   * Secondly, when I do look at both structure and standardized coefficients, two of my dependent variables have the opposite signs for their standardized and structure coefficients (for both of them the structure coefficient is positive and more than .30 whereas the standardized coefficient is negative and has a very small value). I have read up some stuff on it and I think these are suppressor variables (they are related to the other IV's but do not directly predict the DV). How are we supposed to handle these suppressor variables?


title: How to figure out the behind concept of statistical-look problems?

text: The following is my experience doing some researches linked to statistical concepts. In the most situations I had no idea on `how to organize the questions`, `how to ask for answers`, and `how to provide feedback` showing my acceptance or rejection. I'm certainly not talking about the language issue but pointing on that it was not clear to me `what the concept behind my question is`. Sometime it comes out that there is a misunderstanding regarding the question. Some people do not tolerate the lack of formality (~standards) in the question and so their replies are not useful but discouraging. I believe that there is some sort of kindness and helpfulness trying to answer questioner with the hope being satisfied.   Provided this background, I would appreciate it if you would provide me a list of: `tools / tips-and-tricks / methods / skills / techniques / strategies / ...` to find out how to seek for `solution / help / ...` in facing with statistical problems.   For example this site (Stats.StackExchange.com) is a tool providing a forum based technology to ask statistical questions. Some people probably recommend reading books etc which I consider them not a shortcut solution although for young scholars it could be a good idea. The reason is the questioner in my question has lots of thing to do and also has other profession.   My experiences showed that asking a clear/complicated question from an academic (unknown/famous) who does not know you has no way to be succeed. For a keen researcher it is however disappointing.   For this question I tried to be clear as I could if you see it not clear enough, you got the problem. Any ideas? **Please read this before commenting:**   I accept that the main question (the idea behind this post) might be a bit hard to be caught up by immediate-commentators. Briefly, in the first time facing with a problem we would have some guesses about the related topics/subject/methods in statistics that could help us to understand the problem and to propose solutions. This post discusses **the problem solving strategies** `to connect a proper statistical analysis to the problem being investigated` even roughly. Although productivity could be developed by the help of professional statistician however this is not the case to be addressed in this post. The main purpose of this post is to help self-learning practitioners who have enough background about the main concepts being investigated however do not have a clear idea if there is any modern/classic developed concept regarding the problem or not. For statistics as for other parts of sciences there are always `tips-and-tricks` that could help the researcher to speed up the progress of productivity, at least I have a bunch of them. I would recommend commentators to give a fair opportunity for discussions even if it is not their favorite. Please think positive and deep you'll see good points in this post. Nevertheless, to me what you comment or mark or vote will be useful even if those are not directly addressing the question and regardless the tone/language all are appreciated.   **Comments-on-comments**   @cardinal: Thanks. the provided link is an interesting post.


title: Random Forests / adaboost in panel regression setting

text: Tools such as random forests or adaboost are powerful at solving cross- sectional binary logistic problems or prediction problems where there are many weak learners. But can these tools be adapted to solve panel regression problems? One could naively introduce a time index as an independent variable but all this does is to provide an additional degree of freedom to the fitting algorithm. What we would like is a solution that allows information from period T-1 to have bearing on period T. If there is not a straightforward way to do this using these algorithms, is there an alternative algorithm that can perform a panel regression making use of the information in both the cross-section and time-series?


title: Species Richness, Dominance and Diversity Differences

text: I have collected 70 organisms from 4 different sites; two sites of treatment 1 and two sites of treatment 2. I also have a continuous explanatory variable (average temperature) which is different for each site. How can I test if measures of richness and diversity differ between sites or by temperature? **Richness** is the number of species in each sample **Diversity** is a weighted average of the proportions of each species present. In this case we are using the Shannon Index which is: $^1\!D= \exp\left(-\sum_{i=1}^R p_i \ln p_i\right)$ where $p_i$ is the proportion of each species $i$ at that site, and $R$ is the richness of the site. What sort of model can I use when diversity and richness are my response variables? * * * _This question had been abandonned by the OP without giving enough info to answer it properly, the above is an attempt to provide an answerable question in the spirit of what was asked._ Original Question: > I have collected ca. 70 species of organism from 4 sites. 2 sites of > Treatment 1 and 2 sites of Treatment 2. How do I test using R, whether the > richness, dominance, abundance and diversity is different between the two if > I have 6 explanatory variables?


title: Is it possible to obtain Risk Ratios after clogit?

text: Zhang and Yu (1998) proposed a formula to convert Odds Ratios (OR) to Risk Ratios (RR) for cohort studies: $$RR = \frac{OR}{(1 - P_u) + (P_u\cdot OR) }$$ ($P_u$ = incidence rate of the unexposed group) Is there a similar strategy for a clogit model (matched case control study)? Or is it even possible to use Zhang and Yu's formula also in this case? * * * Reference: J. Zhang and K. Yu, 1998. What's the Relative Risk, JAMA, Vol 280, No 19, pp 1690-1691.


title: Use of two Likert scale for correlations

text: I am struggling with two questionnaires. One comprises of 61 items explaining several components of emotional intelligence with 6 point scale while on the other hand 10 items making empathetic listening with 7 point scale. Can I correlate the factors obtained from these two questionnaires? If any scale modification is required,how to proceed the same using SPSS?


title: Box-constrained optimization using optim() in R

text: For my research I am currently working in R and I have created a function which gives me the loglikelihood. Now, I would like to optimize it for multiple parameters. So I decided to use the optim() function in R. Since I have 1 vector parameter and another parameter which are constrained to have values between 0 and 1, since they are proportions, I used the option method="L-BFGS-B". Applying this I get results but one of my restricted parameters is estimated as exactly 1.0000. If I do not constrain it, it gives a sensible answer. But, the standard error is large (still smaller than under the constrained). Is there a way to remain with the restriction but that I can solve the fact that the program ends with 1.000? I give some piece of my code.               Ps<-1-c(0.02083333, 0.03067485, 0.11290323, 0.13684211, 0.23629490)          loglik.negbin.perageclass<-function(param.qphi){     #q.param<-rep(param.qphi[1],5)     #phi.param<-rep(param.qphi[2],5)     #size.param<-param.qphi[3]     q.param<-param.qphi[1:5]     phi.param<-rep(param.qphi[6],5)     size.param<-param.qphi[7]     shape.param<-param.qphi[8]     scale.param<-param.qphi[9]     #I_days[1,]<-param.qphi[8:12]     delta.touse<-dgamma(0:4,shape.param,scale.param)/sum(dgamma(0:4,shape.param,scale.param))     delta<-function(t_current){         diff<-rep(0,t_current)         res<-matrix(0,nrow=t_current,ncol=5)         for(k in 1:t_current){             diff[k]=t_current-k             if(diff[k]==0){                 res[k,]<-delta.touse[1]             }             else if(diff[k]==1){                 res[k,]<-delta.touse[2]             }             else if(diff[k]==2){                 res[k,]<-delta.touse[3]             }             else if(diff[k]==3){                 res[k,]<-delta.touse[4]             }             else if(diff[k]==4){                 res[k,]<-delta.touse[5]             }         }         return(res)     }     week<-1     for(k in  2:35){         I_a_days[k,]<-phi.param*(q.param*C_a[,])%*%(apply(delta(k-1)*Ps*N*I_days[1:k-1,],2,sum))         I_s_days[k,]<-(1-phi.param)*(q.param*C_s[,])%*%(apply(delta(k-1)*Ps*N*I_days[1:k-1,],2,sum))         I_days[k,]<-I_a_days[k,]+I_s_days[k,]         if(k%in%newweek){             beginningweek<-k-6             for(l in beginningweek:k){                 I_week[week,]<-I_week[week,]+I_days[l,]                 I_s_week[week,]<-I_s_week[week,]+I_s_days[l,]                 I_a_week[week,]<-I_a_week[week,]+I_a_days[l,]             }             week<-week+1         }            }     print(I_week)         # Plotting         plot(GP.data$Weekno[24:29],log(GP.data$X0.4[24:29]),pch=19,lwd=2,main="GP data",ylim=range(0,log(700)),col=1,xlab="week",ylab="number of cases")         points(GP.data$Weekno[24:29],log(GP.data$X5.14[24:29]),pch=19,col=2)         points(GP.data$Weekno[24:29],log(GP.data$X15.44[24:29]),pch=19,col=3)         points(GP.data$Weekno[24:29],log(GP.data$X45.64[24:29]),pch=19,col=4)         points(GP.data$Weekno[24:29],log(GP.data$X65.[24:29]),pch=19,col=5)         for (i in 1:5){lines(c(200925:200929),log(I_week[,i]),col=i,lwd=2)}         for (i in 1:5){lines(c(200925:200929),log(I_week[,i]),col=i,lwd=2)}         # End plotting     we<-5/t(t(GP.matr[25:29,13:17])/N)/apply(1/t(t(GP.matr[25:29,13:17])/N),1,sum)     ll<-rep(0,5)     mu.param<-rep(0,5)     for(a in 1:5){         gp.a<-a+4         mu.param<-I_week[,a]         #print(c("debug",N,N[a]*GP.matr[25:29,gp.a]))         ll[a]<-sum(we[,a],log(1E-7+dnbinom(GP.matr[25:29,gp.a],mu=mu.param,size=size.param)))     }     ll.report<-sum(ll)     return(-ll.report)     }          q_start<-rep(0.05,5)     phi_start<-0.10     size_start<-1.2     shape_start<-2.9     scale_start<-0.5          starting.values<-c(q_start,phi_start,size_start,shape_start,scale_start)     loglik.negbin.perageclass(starting.values)          res<-optim(starting.values,loglik.negbin.perageclass,method="L-BFGS-    B",lower=c(0,0,0,0,0),upper=c(1,1,Inf,Inf,Inf),     control=list(trace=TRUE,maxit=500),hessian=T)      Does anyone has any idea if there is a solution for this? Thank you very much! Kim


title: Understanding a double peaked distribution

text: I am running a Monte Carlo Simulation and am sampling randomly from about 65 Normal Distributions, each with a different $\mu$ and $\sigma$. I end up with the Mixture Distribution graph shown below amd trying to evaluate the cause of the 2nd peak. The plot appears to be on the way to becoming bi-modal, with increased iterations _possibly_. Will it be worth trying to evaluate the cause of the modes as is mentioned here? Image: ![enter image description here](http://i.stack.imgur.com/1rOAJ.jpg) Thanks a lot


title: A statistically valid method of comparing values lying outside of a distribution with values within that distribution

text: I recently interviewed farmers about their memories of annual crop production and compared this with diaries they had kept recording actual quantities harvested. I would like to know where interviewees’ recollections sit within the range of production values and how much confidence I can place in their recollections, hence I asked them to recall good, average and poor crop yields over the period they had kept records. In analysing my data, I’ve made the following choices so far:   1. I cannot use absolute values as I am dealing with different types of crops (the only way to get my sample size up); hence I need to use some sort of relative measure to compare between farmers.    2. I’ve discounted comparing percentage difference between recalled values and recorded median, maximum and minimum values, as this fails to take the overall distribution of each data set into account.    3. I’ve settled on percentile ranking of recollections as my chosen way forward as this will illustrate each recollection within the context of the full range of actual production.  However, as some of you may have realised already, point three has thrown up a problem; some farmers recalled crop yields that were either higher than the maximum yield recorded, or lower than the minimum yield recorded. This places them outside of the percentile ranking distribution and thus presents me with a dilemma as to how to include them. They certainly need to be included within the analysis, but to derive confidence intervals I need to present them using a unit that is comparable to percentile rank. I am considering comparing any recalled data falling outside of the recorded values as a percentage difference from the maximum or minimum value. Because the data is (by and large) right skewed this has the effect of increasing the weight of outlying data for resulting confidence intervals. I’d rather underestimate levels of confidence than overestimate, but the problem is, I’m not sure whether this approach (comparing percentile rank and percentage) is statistically robust? Is there any way that I can incorporate these outliers using comparable units to percentile rank that will enable me to construct valid confidence intervals for the position of a poor/typical/good recalled yield within the range of actual yields? As stated before, I don’t just want to compare the relative difference between an arbitrary value i.e. max production/median production and a farmer’s recall, as this does not take the full distribution of the data in to account. I’ve reviewed quite a few papers that have conducted similar research, but they were either able to use absolute values or their interviewees didn’t seem to have such bad memories! Any help you can provide – either comments on the abovementioned approach, or additional thoughts on how I could intertwine the whole set of data without committing a statistical no-no, would be very gratefully received. Thank you very much, RHT


title: Eigenvalue cleansing methods in R?

text: Are there R packages that offer eigenvalue de-noising methods grounded in Random Matrix Theory? Various cleansing methods include the Power Law and Krazanowski filter.


title: How to deal with an error such as "Coefficients: 14 not defined because of singularities" in R?

text: When doing a GLM and you get the "not defined because of singularities" error in the anova output, how does one counteract this error from happening? Some have suggested that it is due to collinearity between covariates or that one of the levels is not present in the dataset (see: interpreting "not defined because of singularities" in lm) If I wanted to see which "particular treatment" is driving the model and I have 4 levels of treatment: `Treat 1`, `Treat 2`, `Treat 3` & `Treat 4`, which are recorded in my spreadsheet as: when `Treat 1` is 1 the rest are zero, when `Treat 2` is 1 the rest are zero, etc., what would I have to do?


title: If you convert factors into indicator variables, do you treat them as continuous predictors?

text: Let's say I have a data matrix X where one feature is a factor with 8 levels. If I change this to be 7 indicator variables of 1's and 0's, do I need to make these columns factors as well? Or if I am training a model on this converted dataset should I just treat these new 1 and 0 values as continuous and scale accordingly?


title: Different AIC values for the same model using step()?

text: I'm working with a GLM to try and optimize the model, and there are 152 predictive variables. A LOT of these are not significant, so I'm trying to figure out which ones to remove through use of the step() function. However, I'm having trouble interpreting the AIC output because it seems as if the models chosen have higher AICs than the original models in almost every case. Let's take a look at a garbage dataset I invented just for example:               > x = c(1,2,3,4,5)     > x2 = c(234,2134,432,1,41)     > x3 = c(123,456,789,123,456)     > y = c(1234,4321,1234,3456,7654)     > bull = glm(y ~ x + x2 + x3,family=Gamma(link=log))     > step(bull,direction="backward")     Start:  AIC=89.13     y ~ x + x2 + x3                 Df Deviance    AIC     - x3    1  0.65042 **88.561**     <none>     0.27488 89.127     - x2    1  1.34701 91.222     - x     1  2.43317 95.370          Step:  AIC=91.5     y ~ x + x2               Call:  glm(formula = y ~ x + x2, family = Gamma(link = log))          Coefficients:     (Intercept)            x           x2         6.3315474    0.4664784    0.0004815            Degrees of Freedom: 4 Total (i.e. Null);  2 Residual     Null Deviance:      2.434      Residual Deviance: 0.6504       **AIC: 91.5**      Notice the starred AIC values are different, and the bottom-most one is higher than the original? Yet, the model is chosen based on an AIC (also starred) that is lower. I'm wondering how to reconcile those 2 values and what exactly this model is doing.


title: Identifiability for time invariant state space models

text: Kevin Murphy's Kalman Filter toolbox (for Matlab) contains an example where it's the fact that the state space system in not identifiable causes problems. I include the example in it's entirety but you won't actually be able to run the code without installing the toolbox first. See the comments between the code blocks where he notes the systems lack of identifiability. Q1.) Can someone please explain why the system below is not identifiable. Q2.) Can someone please given general rules for determining whether a time invariant state space system is identifiable or not.               % X(t+1) = F X(t) + noise(Q)     % Y(t) = H X(t) + noise(R)          ss = 4; % state size     os = 2; % observation size     F = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1];     H = [1 0 0 0; 0 1 0 0];     Q = 0.1*eye(ss);     R = 1*eye(os);     initx = [10 10 1 0]';     initV = 10*eye(ss);          seed = 1;     rand('state', seed);     randn('state', seed);     T = 10000;     [x,y] = sample_lds(F, H, Q, R, initx, T);          % Initializing the params to sensible values is crucial.     % Here, we use the true values for everything except F and H,     % which we initialize randomly (bad idea!)     % Lack of identifiability means the learned params. are often far from the true ones.     % All that EM guarantees is that the likelihood will increase.     F1 = randn(ss,ss);     H1 = randn(os,ss);     Q1 = Q;     R1 = R;     initx1 = initx;     initV1 = initV;     max_iter = 100;     [F2, H2, Q2, R2, initx2, initV2, LL] =  learn_kalman(y, F1, H1, Q1, R1, initx1, initV1, max_iter,1,1);      Thanks Baz


title: Can ICC values be negative?

text: I am calculating the ICC as per the method of Shrout-Fleiss (1979), using the (3, 1) model in which the judges are fixed as are the targets. Some of my ICC values are negative. Is this to be expected? How to interpret this? I understand that the maximum ICC is 1.0, but is there a minimum value? How are negative values interpreted?


title: How to interpret autocorrelation

text: I have calculated autocorrelation on time series data on the patterns of movement of a fish based on its positions: _X_ (`x.ts`) and _Y_ (`y.ts`). By using R, I ran the following functions and produced the following plots:               acf(x.ts,100)      > ![enter image description here](http://i.stack.imgur.com/Jmh60.jpg)               acf(y.ts,100)      > ![enter image description here](http://i.stack.imgur.com/s3CET.jpg) My question is, how do I interpret these plots? What information is needed to report any sort of pattern? I have been surfing the internet and have yet to find a concise way that effectively explains it. Also, how do you decide the correct amount of lag to use? I used 100, but I am not sure if that is too much.


title: Inference about the outcomes of two logistic regressions

text: I have run two separate logistic regressions and would like to assess which model fits the data better. Each model has 1 predictor. Here's the output for both: Model 1:               Deviance Residuals:      Min      1Q  Median      3Q     Max       -1.280  -1.046  -1.046   1.078   1.315            Coefficients:                  Estimate Std. Error z value Pr(>|z|)       (Intercept)  -0.3185     0.1643  -1.938   0.0526 .     L1            0.5564     0.2317   2.402   0.0163 *          Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)          Null deviance: 421.32  on 303  degrees of freedom     Residual deviance: 415.49  on 302  degrees of freedom     AIC: 419.49      Model 2:               glm(formula = Response ~ L2, family = binomial, data = log)          Deviance Residuals:      Min       1Q   Median       3Q      Max       -1.6651  -0.7235  -0.7235   0.7585   1.7138            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -1.2068     0.1927  -6.264 3.75e-10 ***      L2            2.3054     0.2687   8.580  < 2e-16 ***          Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)          Null deviance: 421.32  on 303  degrees of freedom     Residual deviance: 334.99  on 302  degrees of freedom     AIC: 338.99          Number of Fisher Scoring iterations: 4      Is there something similar to an R square for logistic regression that I could use to assess how well it fits, and also that I can compare to the other model? Similarly if I run the lrm command from the rms package, which discrimination index is best to look at and how can I compare the two models?


title: Adaptive simulated annealing

text: I have been doing some reading about adaptive simulated annealing and as far as I know it is an algorithm that is really useful when it comes to finding the global maxima/minima of some functions, which is useful for calibration purposes. However, the explanation on Wikipedia for both ASA and SA are from the "physics" aspect, which involves heavily on the energy state, etc. I am just wondering if anyone could somehow explain to me, in the simplest math or non-technical ways, how this algorithm actually works and why it is better than the other algorithms when it comes to finding global extremes.


title: When will a Kalman filter give better results than a simple moving average?

text: I recently implemented a Kalman filter on the simple example of measuring a particles position with a random velocity and acceleration. I found that Kalman filter worked well, but I then asked myself what's the difference between this and just doing a moving average? I found that if I used a window of about 10 samples that the moving average outperformed the Kalman filter and I'm trying to find an example of when using a Kalman filter has an advantage to just using the moving average. I feel like a moving average is far more intuitive than the Kalman filter and you can apply it blindly to the signal without worrying about the state-space mechanism. I feel like I am missing something fundamental here, and would appreciate any help someone could offer.


title: Image reconstruction using compressed sensing

text: This question deals with an example of image reconstruction related to this other question on signal reconstruction. I have different issues in both examples but there could be an underlying factor. I tried to reconstruct an image using compressed sensing just as in Coursera's course "Computational Methods for Data Analysis" and described in some detail in this pdf (page 414.) By the way, this is the low resolution image: ![enter image description here](http://i.stack.imgur.com/oDZc4.png) However, I have a question regarding the reconstruction of the image using $f \approx \Psi x$ where $f$ is the original image, $\Psi$ is a basis (in this case, a discrete cosine transform) to express the image and $x$ its coefficients. This is what I have done, but I'm not entirely sure it's correct.               from PIL import Image     from scipy.fftpack import dct, idct     from sklearn.linear_model import Lasso          # Loading image in grayscale and obtaining its dimensions     im = Image.open('coffee-cup.jpg').convert('L')     nx, ny = im.shape     # Number of sample points used to reconstruct image     k = 1000          # Create a permutation from 1 to nx*ny and choose the first k numbers     Atest = zeros((nx, ny)).reshape(1, nx*ny)     r1 = permutation(arange(1, nx*ny))     r1k = r1[1:k+1]          # Suppose the first number in the permutation is 42. Then choose the 42th     # elements in a matrix full of zeros and change it to 1. Apply a discrete     # cosine transformation to this matrix, and add the result to other matrix     # Adelta as a new row.     Adelta = zeros((k, nx*ny))     for i, j in enumerate(r1k):         Atest[0,j] = 1         Adelta[i, :] = dct(Atest)         Atest[0, j] = 0          # Use the same permutation to sample the image to be reconstructed     image1 = im.reshape(nx*ny,1)     b = image1[r1k]          # Solve the optimization problem Adelta * x = b     lasso = Lasso(alpha=0.001)     lasso.fit(Adelta,b)      After this, lasso.coef_ contains the coefficients $x$. This is the part I'm not sure about. I transformed the coefficients using the discrete cosine transform. However, the construction of Adelta seems to suggest other more elaborated $\Psi$. Yet, when I plot the reconstruction using the inverse discrete cosine transform (IDCT) and the discrete cosine transform (DCT), this is what I get:               recons2 = dct(lasso.coef_).reshape((nx,ny))     recons = idct(lasso.coef_).reshape((nx,ny))          fig = figure()     ax = fig.add_subplot(121)     ax2 = fig.add_subplot(122)     ax.imshow(recons)     ax2.imshow(recons2)      ![enter image description here](http://i.stack.imgur.com/mfHAn.png) It worked reasonably well, but notice that the first image (in which IDCT was used) performed better. I would have assumed that, if something, DCT was more appropriate choice since at least we used the result of a DCT in the construction of Adelta. What is the correct procedure here?


title: How to "explain away" an effect size using an independent variable?

text: I have many measures of Y for two groups A and B. Following the more modern approach of research in education, I calculated the effect size of the difference of Y for both groups (Cohen's d), and the significant interval for the effect size. The fact that the effect size was positive and significant (both in statistical sense - the interval did not include the 0, and in the relevance sense - the value was large) was surprising. So I want to explain away the effect size by using a new independent variable X. The argument goes like this - the difference between A and B is due to the fact that A's have on average larger X's than B's, and it is known that larger X "causes" or is correlated to larger Y's. If I wanted to explain away Y by using X, I would do the regression:               Y ~ C1 + C2*X + C3*AB + C4*AB*X      where AB is the categorical variable that indicates whether the group is A or B, and then show that the coefficients C3 and C4 are not significant (in the statistical sense). Thus, the coefficients that include the variable AB are not really different than 0 and thus the Y can be fully explained by the X alone. In fact I did that, and the coefficients are significant, but I cannot interpret if they are large of not because they are in units of Y, and I want in effect size "units" (effect size do not have units and that is what makes them so useful!). The variable X is an integer with a small range of values, so I could calculate the effect size of AB for each of the values of X, and than perform a regression on those values and X, but I find this solution odd. In summary, given an effect size E, based on a categorical variable AB, how do I measure how much it is explained away by an independent variable X? EDIT - after @Peter Flom's answer I was not clear to specify that I do not know how to calculate the Cohen's d if there is a covariate variable X. Cohen d = ($\mu(Y_a) - \mu(Y_b))/\sigma_{ab}$ where $\sigma_{ab}$ is a pooled standard deviation. With the regression model above I now have Y ~= C1 + C2*X + C3*AB + C3*AB*X and I no longer have groups to calculate the mean of, or the standard variation! I would very much like remain in the Cohen d domain. The first analysis - by calculating the d for the groups A and B was standard in the literature, and I could compare that d with others to claim that the effect was large. I would like now to show that when X is taken into account the d reduces to a very small (or hopefully not significant) number.


title: Binomial coefficient in JAGS

text: I'm trying to write such likelihood using JAGS in `R` (in order to estimate parameters $\xi$ and $\pi$, being $m$ fixed equal 7) using the "ones trick" and got stuck in defining the binomial coefficient (`term3`). How would you define the binomial coefficient? $$ P(Y=y)=\left\{ \pi\binom{m-1}{y-1} (1-\xi)^{y-1}\xi^{m-r}+(1-\pi){1\over m}\right\} $$ This is my model:               model {       for (i in 1:N) {             term1[i] <- pai[i]*pow(1-csi[i],(y[i]-1))*pow(csi[i],(m-y[i]))         term2[i] <- (1-pai[i])*(1/m)         term3[i] <- logfact(m-1)/logfact(m-y[i]-1)*logfact(y[i])              L[i] <- term3[i]*term1[i]+term2[i]         phi[i] <- L[i]/C         ones[i] ~ dbern(phi[i])               csi[i] ~ dbeta (alfa, beta)         pai[i] ~ dunif (0,1)       }     }


title: Estimating the reliability of a test section

text: Consider me a "naive researcher." We have a test of 25 items grouped into 5 item types. Each item is correct or incorrect. All items are math questions, but we expected 4 item types to have different language confounds. We have a population of 28 group A and 32 group B. We expected group A to under perform group B in the 4 item types with language confounds. Overall the test scores are negatively skewed but the SDs are similar between groups. In t-tests for each section we found a significant difference in those 4 item types in the expected direction and no significant difference for the fifth type. Cohen's D shows medium effect sizes in the sections with significant differences. Does it make sense to estimate and report reliability for each of the item types? If so what would be the best estimate?


title: For a continuous random variable, why does $P(a < Z < b) = P(a \leq Z < b) = P(a < Z \leq b) = P(a \leq Z \leq b)$

text: My textbook puts this in a sidebox with the heading "Note" and doesn't explain why. Could you tell me why this statement holds? $P(a < Z < b) = P(a \leq Z < b) = P(a < Z \leq b) = P(a \leq Z \leq b)$


title: Which machine learning methods can I use to study/analyze this ordinal data? (generated by a hypothetical turing machine)

text: I have a set of ordinal data, generated by a hypothetical turing machine. The data consists of a sequence of symbols such as:               LB1, LB2, UU, UB1 ...      The labels have a natural ordering, and the machine "prints" out the symbols at an irregular pace. I would like to analyze, for the purposes of determining (with a degree of certainty), given the history of printed symbols:   1. How long before the next symbol is printed by the machine.   2. The most likely next symbol to be generated by the machine.


title: The reciprocal of $t$-distributed random variable

text: How is $1/T$ distributed if $T$ follows a Student's $t$-distribution?


title: Is this longitudinal data too complicated for GLMM or GEE?

text: After writing this post, I've realized that I am running around in circles, chasing my tail. Any help approaching this problem would be greatly appreciated, as I think I just need to bounce ideas around and don't have colleagues that can help with statistical methods at my new position. I am working on an incomplete longitudinal/repeated measure dataset, which is very new to me. After spending a week working with the data, I have come to the (possibly incorrect) conclusion that I should be approaching the data with a GLMM or GEE. I say incorrect, because I wonder if I have too many dimensions to model through these approaches (potentially leading to a PCA to reduce my dimensions?). I've been attempting to approach the data with lme4() in Program R, but am having a major issue in wrapping my head around how to account for the various variables of interest. I believe my data has a layer or two of complexity beyond the 'sleepstudy' data in Bates Chapter 3 (on the use of Mixed Effects Models in R), so would appreciate any advice on what I am missing and what direction I may pursue, as my brain has reached it's carrying capacity for selflearning at the moment My simulated dataset can be found here. A treatment (drug vs control) was administered to a set of subjects at TimePeriod 1. A Measurement of interest was taken at TimePeriod 0 (before treatment), and at TimePeriods 1, 2, 3, 4, 5 (spaced each 12hrs). These treatments are conducted on newborn babies born at different gestational ages (GA), but are conducted at similar (but not identical) HOL (hour of life). The dataset includes:   * 'GA' (gestational age at birth)   * the biological 'Measurement' of interest    * 'TimePeriod' (which timeperiod the measurement was taken at. 0=predrug/control administration, 1=post, 2=12hrs after, 3=24hrs after, etc)   * HOL (hour of life at eachmeasurement)   * 'Group' (term vs preterm birth, which is a categorization of GA that I don't think is needed)   * 'Treatment'(drug/control).  There are 6 measurements per subject each taken at a similar HOL between subjects. The main question is if treatment (drug/control) affects the measurement of interest. Another question is how GA and HOL relate to the measurement of interest. We can open this in R via               newData<-"[paste linked-data here]"     Data <- read.table(text=newData,header=TRUE)      We then can use ggplot2 to visualize the data, using my 4-dimensions of interest (GA, HOL, Treatment, and Measurement)               library(ggplot2)     library(grid)     s <- ggplot(Data, aes(x=HOL,y=GA, shape=Treatment,fill=Subject,color=Subject))     s <- s+geom_point(aes(size=Measurement))     s + opts(legend.key.height=unit(.3, "cm"))      .![enter image description here](http://i.stack.imgur.com/ergRC.png) What we see with this plot is GA (gestational age) on the y-axis, and HOL (hour of life on the xaxis). The size of the points give information on the Measurement of interest, while shape gives Treatment group. Clearly, the Measurement increases with increasing GA and HOL, and each subject has a unique GA with a ~15week spread of values, while the 6 measurements per individual are at similar HOLs. It is important to note that this is based on Biological Relevance. In our real data, we expect a strong effect of increasing HOL and GA on Measurement 1) The first question to analyze. Does Treatment (control/drug) affect the measurement of interest vs 'TimePeriod' (which tracks time of dosing, where solid black line is the dose/treatment. X-axis has been 'jittered' for clarity). ![enter image description here](http://i.stack.imgur.com/konzS.png) I envisioned a GLMM where TimePeriod is X, Measurement is Y (graph above), and then HOL and GA are the Random Effects. But then, how do I test for a Treatment effect between Drug/Control? I have created a plot of mean values and confidence intervals (normalized for within-subject variability via Morey 2008), which shows no effect of Treatment, but I wonder how to approach this with a GLMM (if that's appropriate), instead of just visualizing group means and CIs. I believe I understand how to make a GLMM for the Drug or Control group separately, but do not see how to compare the two through a statistical test (beyond plotting them on the same graph). I am clearly overlooking something in my brain-fog. 2) The first question is really my main interest. The second question is then, assuming there is no Treatment effect, can I model the relative importance of HOL and GA on the measurement of interest (again, we believe we will find no Treatment effect, but will find a HOL and GA effect)? This relationship is seen in the first graph I provide (measurement increases as you move up on both axises), but for some reason, I am having a complete mental block on approaching this question statistically. EDIT (6/26/14): I've been working through this, and have come to this approach.                lmer.full<-lmer(log(Measurement)~HOL*Treatment+GA+(1|Subject))      This allows a treatment of longitudinal data through the (1|Subject) random effect. It also simplifies the data since 'TimePeriod' and 'HOL' are so heavily correlated. Instead of asking how the Measurement changes per time- period, I can explore the interaction between Hour-of-Life and Treatment. Here is my new issue. Gestational Age is heavily correlated with the Subject effect (which makes sense since each subject has it's own Gestational Age). My questionis, if due to this high correlation, I should leave out GA, and just focus on the (1|Subject) random effect to control for this variability. EDIT 2: I have improved the model on advice of Bodo Winter, who has a great personal website with some nice beginner tutorials I highly recommend www.bodo-winter.net First, I zero-centered my Time parameter                Data$HOL<-Data$HOL-mean(Data$HOL,na.rm=TRUE)      I also removed all initial measurements before treatment administration                Data2<-subset(Data,Data[,6]>0)       This removed convergence issues I was having when adding a correlated random slope/intercept as such                lmer.full.slope<-lmer(log(Measurement)~HOL*Treatment+GA+(1+HOL|Subject))


title: Maximum likelihood estimation for state space models using BFGS

text: I have written some code that can do Kalman filtering (using a number of different Kalman-type filters [Information Filter et al.]) for Linear Gaussian State Space Analysis for an n-dimensional state vector. Let us consider the linear Gaussian State Space model $$y_t = \mathbf{Z}_{t}\alpha_{t} + \epsilon_{t},$$ $$\alpha_{t + 1} = \mathbf{T}_{t}\alpha_{t} + \mathbf{R}_{t}\eta_{t},$$ where $y_{t}$ is our observation vector, $\alpha_{t}$ our state vector at time step $t$. The quantities in bold are the transformation matrices of the state space model which are set according to the characteristics of the system under consideration. We also have $$\epsilon_{t} \sim NID(0, \mathbf{H}_{t}),$$ $$\eta_{t} \sim NID(0, \mathbf{Q}_{t}),$$ $$\alpha_{1} \sim NID(a_{1}, \mathbf{P}_{1}).$$ where $t = 1,\ldots, n$. Now, I have derived and implemented the recursion for the Kalman Filter for this generic state space model by guessing the initial parameters and variance matrices $\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$. where the points are the Nile River water levels for Jan over 100 years, the line is the Kalamn Estimated state, and the dashed lines are the 90% confidence levels. Now, in order to estimate the initial parameters of my model ($\mathbf{H}_{1}$ and $\mathbf{Q}_{1}$) I have adopted the Expectation Maximisation method to maximise the loglikelihood function $$\log L(Y_{n}|\psi) = -\frac{np}{2}\log(2\pi) - \frac{1}{2}\sum^{n}_{t = 1}(log|\mathbf{F}_{t}| + v^{\mathsf{T}}_{t}\mathbf{F}_{t}^{-1}v_{t}), \;\;\;\; (1)$$ with respect to parameter vector $\psi$. Where $v_{t}$ is the state error and $\mathbf{F}_{t}$ is the state error variance. The vector $\psi$ contains the unknown parameters (i.e., either the measurement and state disturbance variances themselves, or a re-parametrisation of these same variances, see below) A crucial role in the maximisation of (1) is played by the gradient or score vector $$\partial_{1}(\psi) = \partial \log L(Y_{n}|\psi)/\partial\psi,\;\;\;\; (2)$$ At least two algorithms can be used to maximise (1). The first is the so- called EM (Expectation-Maximisation) algorithm, and the second is the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm. I have successfully implemented the EM algorithm and this is providing correct results but is very slow; in light of this I now wish to implement the BFGS algorithm as this has better performance. So, from (2) using Taylor's theorem we can write $$\partial_{1}(\psi) \simeq \tilde{\partial}_{1}(\psi) + \tilde{\partial}_{2}(\psi)(\psi - \tilde{\psi}),$$ where $\partial_{2}(\psi)$ is the Hessian matrix or second derivative of the loglikelihood function wrt $\psi$ and $\tilde{\partial}_{N}(\psi) = \partial_{N}(\psi)|_{\psi = \tilde{\psi}}$. So now we can write our updated parameter vector as $$\bar{\psi} = \tilde{\psi} - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi) = \tilde{\psi} - s\tilde{\pi(\psi)}.\;\;\;\; (3)$$ where here $\tilde{\pi(\psi)} = - \tilde{\partial}_{2}(\psi)^{-1}\tilde{\partial}_{1}(\psi)$ and $s$ is our step size. So the basic outline of the BFGS algorithm to get our parameter vector and maximize (1) using BFGS is:   1. Initialise parameter vector $\psi = \psi^{*}$.   2. Then apply the Kalman and disturbance smoothing filters and thus for score vector (2) at $\psi = \psi^{*}$. When running the Kalman Filter evaluate the log-likelihood function.   3. Use (3) to obtain new values for $\psi$ given by $\psi^{+}$ and go back to step 2 until the value of loglikelihood function (1) no longer improves. This seems fine, but my problem is how to calculate the step size $s$ in (3) above. I know a line search is required, but remember $\psi$ is a vector and it is not clear what function I can use in a traditional line search technique in order to evaluate $s$. **Question: what line search method should I use in order to evaluate $s$?** * * * Ps. Assume that I would like to use the zoom line search algorithm (algorithm 3.5 of Numerical Optimization by Nocedal & Wright).


title: Machine-Learning algorithms for Forecasting

text: For work, I'm working on an app where you essentially forecast the failure rate of the overall machine through different factors such as the historical failure rates for the components used to build it or the failure rates of the factories that manufacture it, or even the historical rate for the machine itself. The idea is that for any machine you can make a solid prediction, so I need some algorithm to self-build a good model for each of the 1000s of machines. I've been able to implement this using ARIMAX models, but I just don't feel good about using auto.arima and then just cross-validating to see how many external regressors to add in. I've also tried SVM, but what seemed to happen was that the model was not good at dropping irrelevant factors, and therefore the prediction was a flat line. I feel like boosting would be a promising area, but I was wondering if anyone had other options and could more importantly, point me to examples of how the specific algorithm was implemented in R? I'm actually an undergrad intern majoring in statistics, so I'm not too strong in the actual programming side of things, so am not very good at implementing the theory I read about into R code. Also, would a normal GLM be good enough? I used ARIMAX because I wanted to correct for autocorrelation.


title: Cochran-Armitage test for trend in R

text: I have a 4x2 table that I want to do a Cochran-Armitage trend test for. Basically, I want to test for lack of independence as the dose (the rows) increase. I'm trying to figure out the R code but am having some trouble. As there's no clear documentation in the coin package. What's more my rows are not linearly spaced but are spaced at 0,10,30, 100. How would I need to modify:               independence_test(tumor ~ dose, data = lungtumor, teststat = "quad", scores = list(dose = c(0, 1, 2)))      if it's in table form (not long form). I'm assuming I would need to change the "scores" too? My table is called table1. Thank you!