title: What type of Research Design?

text: In an experiment, 25 subjects (Group A) are assigned to play an adventure game (with audio and non-audio) while another 25 subjects (Group B) are assigned to play a strategy game (with audio and non-audio). Each game has two levels (audio and non-audio). The independent var are: (i) groups (A and B) (ii) game types (with two levels: adventure and strategy) and (iii) game versions (with 2 levels: audio and non-audio). Group A and Group B are independent of each other.   * One within-subjects factor is the game version (audio and non-audio)   * One between-subjects factor is the game types (adventure and strategy)    * One between-subjects factor is the group (A and B) Can this research design be considered a 3-Way mixed-ANOVA, with one within- subjects factors and two in-between subjects factors?   Kindly advise


title: Approximating distributions in expectation propagation

text: Can the approximating distributions for various factors in Expectation propagation be different distributions but still from the exponential family. For example, I have the following posterior form: $$ p(w, \lambda, \phi) = p(\phi)p(\lambda)p(w|\lambda) \prod_{i}p(y_{i}|w, \phi, \lambda) $$ So, I need to approximate each of these factors. My question is can, for example, $p(\phi)$ be approximated as a gamma distribution and $p(w|\lambda)$ be a multivariate Gaussian and $p(\lambda)$ be some thing else from the exponential family.


title: discrete survival analysis reorganise data

text: I cannot figure out how I will reorganize and run my data for discrete survival analysis. I have 1169 farms repeated from 2001 to 2010 and I measure years in low income but poverty is not death. So spss measures 10*1169= 11690 cases and no 1169.


title: Bootstrapped p-value

text: I have a p-value that I generate via resampling. **Resamples = 5000** **Positive findings = 1000 positive findings** **P-value = 1000/5000 = 0.2** How can I compute the 95% confidence interval for this p-value? I would assume it's a function of the number of positive findings and resamples 1000 and 5000 in this case, respectively. Is the answer **p-value +- 1.96*sqrt(1000/5000 * (1 - 1000/5000) / 1000)**? Why this matters: Resampling is very expensive in my code and I'd like to stop resampling as soon as the 95% confidence interval for the boostrapped p-value doesn't include 0.05. Right now I'm doing a million resamplings to estimate every p-value and it is very slow.


title: Accounting for the nonlinear behaviour of a binary predictor

text: I have a continuous response variable $Y$ that I have broken down into four categories, {Low, Medium, High, Very High}. I'd like to model $Y$ against a binary variable $X$. Calculating the proportion of $X=1$ for each of the four categories, I notice that the proportion of $X=1$ is low in the first and last category (Low and Very High), but higher in the middle categories (Medium and High). ![proportion of x=1 for each response category](http://i.stack.imgur.com/w5Cys.jpg) So clearly $X=1$ is more predictive of a Medium or High response. I'd like to somehow account for this in any models I run, as I feel it is important to be able to adjust for the effect of this variable when adding further predictors into the model. I understand the polr function in R assumes that the proportion of $X=1$ moves linearly as I move linearly between the response categories, so I'm uncertain as to what type of model I can use instead. My question is: what can be done to account for the nonlinear behaviour of x? And what kind of packages in R handle these kind of models?


title: Proximal SVM and classification algorithms

text: I was trying to understand the parallel machine learning using Proximal SVM. The implementation I'm following is from atbrox.com. It was a no-brainer to execute the python script and get the output. But the problem is, I am unable to do the actual classification, since the implementation was not there. So, after training is done, the output ((omega, gamma) = (I/mu + E.T*E).I*(E.T*D*e)) we get is,               "outputkey" "[[-0.05385383942136976], [-0.05354571802246949], [-0.03507846021345456], [-0.010885702939765561], [-0.0024890810132077723], [-0.09428323024469554], [-0.03059508326981175], [-0.03772250284905853], [-0.0011720079823984084], [-1.3891341572668838]]"      So to do the next step - (x.T*omega - gamma), we need to figure out the omega and gamma from the above output. Can any one show me how to do x.T*omega - gamma with the above output? (Let x be one from the input sample itself).


title: Test of Normal Distribution

text: Before I start be easy on me I am not a stats person - I am just on the quest for information :-) I was tasked to generate many graphs for someone who wanted to see if the data presented with a bell shape curve vs a graph with > 1 peak - this was no biggie but silly me I thought maybe it would be easier to run the data thru some stats package like R that would flag the populations that need to be looked at with a graph. The quest has led me to reading Shapiro-Wilk, Anderson–Darling test etc but various people with more knowledge than I still pointed to the graphing of data as the only way to verify Normal Distribution. So my question is is there any Statistics method I can run the following data thru to flag the data for further investigation instead of graphs ??   Instead of Counts of Patients by days should I just look at full list of DaysInHospital Thanks > Data is plotted with a graph per disease with the data plotted   >  x-axis - Days in Hosiptal   >  y-axis - Cnt of Patients DiseasePopulation #1               ( DaysinHospital,Cnt of Patients) = {(1,1),(2,1),(3,3),(4,6),(5,2),(6,1),(7,1),(8,1)}      DiseasePopulation #2               ( DaysinHospital,Cnt of Patients) = {(1,1),(2,1),(3,1),(4,2),(5,6),(6,3),(7,1),(8,1)}


title: Forecasting ar(p) for several counties

text: I have a data set of prices, these prices vary across time and across area. I have 18 areas with 32 time periods. What i want to do is forecast these prices, i have found that a AR(3) process fits the data for each area using first difference. \begin{equation} \Delta \hat{y}_t = \hat{\beta}_0 + \hat{\beta}_1 \Delta y_{t-1} + \hat{\beta}_1 \Delta y_{t-2}+ \hat{\beta}_1 \Delta y_{t-3} + \epsilon_t \end{equation} I have been told to investigate panel data methods, however the information i find about dynamic panel data methods requires "large" amount of time series observations (32 not really "large"). My concern is for the latent heterogeneity, assuming prices in one area affects other/neighboring areas, over time. How should my estimation be done when the ultimate goal is forecasting the individual area prices. I have considered several estimation methods, but lack specific knowledge about how to select the best. Estimation methods considered:   * DPD (dynamic panel data)   * SUR (Seemingly unrelated regression)   * SAR (sptial autoregressive models)  Best regards


title: Vector autoregression with gradient descent

text: I am no expert in statistics, but I have been asked to implement a VAR model using gradient descent in R. I have written a code that, from what I have been told, it makes sense. However, the estimated parameters are completely different from the ones calculated with gretl. This is more or less what I have been told to write:               for (k in 1:20) {         J[k] <- 0         costs <- c(0,0,0,0,0,0,0,0,0)              for (i in 1:m) {             if (i > 1) {                 h <- theta[2]*x[["a"]][i] + theta[3]*x[["a"]][i-1] + theta[4]*x[["b"]][i] + theta[5]*x[["b"]][i-1] + theta[6]*x[["c"]][i] + theta[7]*x[["c"]][i-1] + theta[8]*x[["d"]][i] + theta[9]*x[["d"]][i-1]                 y <- x[["y"]][i] - x[["y"]][i-1]         } else {             h <- theta[2]*x[["a"]][i] + theta[4]*x[["b"]][i] + theta[6]*x[["c"]][i] + theta[8]*x[["d"]][i]             y <- x[["y"]][i]         }              difference <- (h - y)         costs[2] <- costs[2] + difference * x[["a"]][i]         costs[4] <- costs[4] + difference * x[["b"]][i]         costs[6] <- costs[6] + difference * x[["c"]][i]         costs[8] <- costs[8] + difference * x[["d"]][i]              if (i > 1) {             costs[3] <- costs[3] + difference * x[["a"]][i-1]             costs[9] <- costs[9] + difference * x[["b"]][i-1]             costs[5] <- costs[5] + difference * x[["c"]][i-1]             costs[7] <- costs[7] + difference * x[["d"]][i-1]         }              J[k] <- J[k] + difference^2     }              for (i in 2:9) {             theta[i] <- theta[i] - costs[i]*learning_rate/m             print(theta[i])         }              J[k] <- J[k]/(2*m)     }      Could anyone tell me what is wrong, or if gradient descent cannot be used to estimate a VAR model? The dataset is very big, that is why we need gradient descent.


title: Comparing samples across time periods

text: I have about 80,000 observations that I collect each day, which include a numeric variable called amount (it's a bid). I often stuff with determining whether a decrease or increase in the amount, be it by hour, day, or week, is meaningful enough that it warrants further analysis. So let's say I have a week's worth of data and I notice that Thursday had a lower average bid amount than the other days of that week. I want to use a 'test' to determine if the values from Thursday are 'significantly' different from the rest of the week. My initial thought was to use a one-sample t-test as it allows me to assess difference between a sample mean (Thurs's mean bid amount) and the population mean (weekly bid amount) to verify if the null hypothesis is correct. However, I'm a little confused how to proceed given that I'm dealing with time series data. How then can I go about comparing bid amounts from different time slots or days to determine whether there is a significant increase or decrease in the bid amount? PS: Question tagged as 'R' because I'll be using R, although the question itself isn't R specific.


title: Appropriate model for proportional response variable

text: I'm researching the behaviour of a bird called Brent Goose. I plan to collect data on 100 individuals of geese. My response variable is divided into four categories: feeding, walking, aggression and other. I plan to watch a bird for 5 minutes, and record the total proportion of that 5 mins that is spent feeding, walking, in aggression and doing other. For example, a bird watched for 5 mins might spend a total of 2 mins feeding, 1 min walking, 1 min aggression and 1 min doing other. Predictor variables are all categorical or integers and include: habitat, age and date. Could someone suggest an appropriate statistical model to analyse these data? In particular, a model that can handle the proportional response variable.


title: Mixed effects model or mixed design ANOVA in R

text: I'm wondering if you could quickly advise on the best statistical analysis for my data. I have a designed experiment with 12 plots within one larger area. This area was blocked into 6 to overcome issues of natural gradients present there, so there are 2 plots in each block. Treatment 1 is a two level treatment split between these 12 plots i.e. 1 plot in every block receives one level of the treatment ("High") and the other receives another level ("Low"). Nested within this, each plot is split into two (24 subplots), and given a second two level Treatment 2 (levels "A" and "B"). Every single plot receiving Treatment 1 has a Treatment Two Subplot receiving level "A", and a subplot receiving level "B". As such, I have a two-way factorial design, with 4 types of treatment combination: High/A, High/B, Low/A and Low/B. This experiment was run for three years, and species richness measured as a response (a simple count of the number of species present in each Treatment 2 subplot) each year. I am interested to see what the effects of the two levels of Treatment Two are on Richness change with time, in the context of the level of Treatment 1 received. I am not sure if I need:   * A two-way mixed ANOVA, with within subjects effects (Year?) and between subjects effects (Treatment 1 and Treatment 2?)   * A mixed effect model, with Treatment 1 and Treatment 2 as fixed effects, and Year and Plot as random effects. (Is Year a random effect? I am unsure...). Should subplot also be a random effect here?   * A mixed effect model with a time series.    * I am also unsure as to what to do with block. Block doesn't have levels of treatment, it was just an attempt to reduce the incidence of similarity between plots as a result of the natural gradients in the study area...   * In addition, because the richness data is a count of species, i'm not sure if this affects my choice of analysis... If you have any thoughts on any of these points, I would be glad to hear them. Thanks!


title: Monte Carlo simulation of low probability events with impact

text: I wonder if I'm lacking the terminology to phrase this question: I would like to model, in **Excel** , the risk of a fairly **rare event** , E, happening over time, where time is divided into chunks called periods, each with n trials. So the event has a low probability (p < 0.1%) of occurring but a **high variable impact** , hence the client's interest. If the event impact were constant I would use the **inverse binomial** function to model the occurance of events in each period as `BINOM.INV(n, p, RAND())`. A possible approach is then to layer the impact distribution (assume Normal) on top of this: NORM.INV( RAND(), $\bar{x}$, $s$ ). By repeating this many times I can build up a **monte-carlo** simulation model (into which I need to feed other parameters). It feels contrived and unrobust, probably inevitably: I don't have much data on how often E occurs, or its impact impact, so clearly I can't be confident about when these events will occur in the future. _Which doesn't stop people asking._ Here's some data: In a sample of 7000 "trials", the event happened 5 times with impacts 6.3, 6.3, 6.7, 6.8 and 7.6. There are 100 trials per period. My specific questions, which may not be the right ones, are:   1. Given the need to model this event, how appropriate is this combined simulation model?   2. What better ways are there to model these kinds of events - that have both a probability of occurring and an variable impact? Thanks in advance for your patience and ideas.


title: Cronbach's $\alpha$ - two instances

text: From a survey I have the following data: 44 questions are supposed to measure 11 variables (underlying theory, 4 questions per variable). These 11 variables are supposed to fit into two categories (underlying theory, 5 variables category A, 6 variables category B). I am unsure how to test this. I was advised to use Cronbach's $\alpha$ instead of conducting a factor analysis, I have some doubts though. Further I may have a situation like this: I use Cronbach's $\alpha$ and find let's say the reliability of the scales is the highest if always 3 questions are used for each of the 11 variables. Out of these 33 questions I create the 11 variables. And now the bit that confuses me (even more): **To determine the reliability of the categorization, shall I test $\alpha$ of the 5/6 variables, or the 15/18 questions? Or would it be the same?** ![Model](http://i.stack.imgur.com/T5rbK.jpg) I am quite open to other approaches. I was thinking about looking into McDonald's $\omega_{h}$ or other FA approaches.


title: How to exclude covariates of a Cox model?

text: I am using the Cox regression model to forecast failures of machines. I have several possible covariates. I ran the analysis in SPSS. The method I used to add covariates was backward stepwise (Wald). My question is that after I got the results, which is the table 'Variables in the Equation' (see below), I can see that some covariates are not significant: for example, Cate and PCPT are not significant. Can I still use this solution? I mean excluding the ones which are not significant, but keep the rest, and use the coefficients $\exp(B)$ for the other covariates in the Cox regression equation $h(t)=h_0(t)\exp(zx)$? Or I should exclude these two non-significant covariates when I add them to SPSS and get a new solution in which all the covariates are significant? Variables in the Equation                         B     SE   Wald df  Sig. Exp(B)     ACT     -.004 .001 31.814  1 .000   .996     TPCPT    .002 .001  6.287  1 .012  1.002     ALT      .008 .001 31.046  1 .000  1.008     MALTCPT 1.315 .204 41.465  1 .000  3.723     LTCPT   2.101 .592 12.601  1 .000  8.176     Cate     .076 .046  2.743  1 .098  1.079     PCPT     .004 .002  3.611  1 .057  1.004


title: How do I compare 2 models?

text: I would like to show that a new biomarker (let's say percent tumor volume reduction at 3 weeks of treatment estimated with automatic methods) performs better than an old one (same measure, but based on manual volume measurements) to distinguish groups or predict outcome. I have these different settings:   1. outcome: progression free survival (PFS) (I have all data, no censoring) method: a) log-rank test for trends in stratified Kaplan Meier b) cox regression (is that a good idea if I have all data?) c) logistic regression (binary oucome PFS<6month)   2. outcome: treatment response (3 groups) method: pairwise test on means (Mann-Whitney-U) I'd like to know if there are better approaches for these settings, and how can I compare the resulting p-values (for the new and for old biomarker). Can I compare them directly, should I run bootstrapping etc. thanks.


title: Odds ratio for continuous variable

text: I have a question > How would you interpret the odds ratios for a continuous variable?


title: Bounded in probability ($o_{P}$)

text: In Larry Wasserman's lecture notes on $o_{P}$ and $O_{P}$, I am not able to follow the derivation of the following example in page 9. Consider $m$ coins with probabilities $p_{1}, \ldots ,p_{m}$. Then \begin{align*} \mathbb{P}(\max_{j} | \hat p_j - p_{j}| > \epsilon) & \le \sum_{j=1}^{m} \mathbb{P}(\hat p_{j} - p_{j}) \quad \text{(Union bound)} \\ & \le \sum_{j=1}^{m} 2 e^{-2 n \epsilon^{2}} \quad \text{(Hoeffding's inequality)} \\ & = 2 m e^{-2 n \epsilon^{2}} \end{align*} I thought of concluding as $n \rightarrow \infty$ we get $\mathbb{P}(\max_{j} | \hat p_{j} - p_{j}| > \epsilon) \rightarrow 0$ and thus $$ \max_{j} | \hat p_{j} - p_{j}| = o_{P} (1) $$ But the author bounds $m$ in terms of $n$ as follows: Suppose $m \le e^{n^{\gamma}}$ where $0 \le \gamma \le 1$. Then \begin{align*} \mathbb{P}(\max_{j} | \hat p_{j} - p_{j}| > \epsilon) & \le 2 m e^{-2 n \epsilon^{2}} \\ & = 2 \exp(-(2 n \epsilon^{2} - \log m)) \\ & \le 2 \exp(-(2 n \epsilon^{2} - n^{\gamma})) \rightarrow 0 \end{align*} Then he concludes $$ \max_{j} | \hat p_{j} - p_{j}| = o_{P} (1) $$   * The necessity of bounding $m \le e^{n^{\gamma}}$ is to avoid the cases where $m$ is large. Is my understanding correct?    * What happens if the $m \le e^{n^{\gamma}}$ is not satisfied? Can we prove it is not convergent?


title: How to calculate SD from pre SD and post SD?

text: I have two values of SD calculated from paired samples of size $n=10$. Results are:   1. Pre test - mean=29.1 and SD=19.05   2. Post test - mean=34.9 and SD=16.74 How do I calculate the SD?


title: Dynamic time warping between more than two sequences

text: Dynamic time warping is a technique to find an optimal alignment between two given sequences. Is there a technique to find an optimal alignment between more than two sequences? Does it as an implementation in python or R? Thanks.


title: Confidence interval for a small number of iid Poisson

text: I want to calculate the confidence interval for $\lambda$ from a small ($n=10$) set of repeated observations from a Poisson distribution. That is, I have $X_1, \dots, X_{10}$ which I believe are i.i.d. and I want a 95% confidence interval given these observations for $\lambda$. I know I can use $\bar x \pm 1.96\frac{s}{\sqrt{n}}$ when $n$ is large but that doesn't seem to be a good assumption here. I also know for a single instance $y$ you can use               qchisq(0.025,2*y)/2, qchisq(0.975,2*(y+1))/2      so my problem appears to be somewhere between the two.


title: Questions about G*Power

text: Since I found an insignificant result, I performed a power calculation in G*Power. But I am not entirely sure if I input the right number in some fields. My design is a 2(between Ss) x 2(between Ss) design, so I have 2 variables, 4 different conditions (20 participants in each one). I ran an univariate test in SPSS. ![enter image description here](http://i.stack.imgur.com/ciO3l.png)   1. Is the **SD sigma within each group** the averaged SD of the 4 conditions?   2. is the **number of measurement** 4?   3. I am not sure what to out in **Corr among rep measures** , so I ran a bivariate correlation test in SPSS (Pearson, 1-tailed). I found the correlation between the 2 variables is 0. So I put 0 there. Am I correct? by the way, have anyone read an experiment in which researchers found an insignificant result because the power is not enough (due to small sample size). I want to read an article like this since I am having the exactly same situation now. Thank you so much! :)


title: Geweke diagnostic of a Markov chain: why does the first window have to cover the burn-in?

text: I read the following statement in this document$^{[1]}$, at the bottom of page 11: > Too wide A will some times “hide” the burn in part within the converged part > of the chain and the difference in means may not be large enough for the > diagnostic to give an indication of non convergence. where A is the first window in Geweke diagnostic. Geweke statistic is defined as: $$Z_n=\frac{(\bar \theta_A - \bar \theta_B)}{\sqrt{\hat S_\theta^A(0)/n_A+\hat S_\theta^B(0)/n_B}}$$ where $A$ and $B$ refers to the first window and second window respectively, $\bar \theta$ is the mean value, $\hat S_\theta(0)$ is spectrum density, $n$ is the number of data points. What I don't understand here is: since this quantity should converge to $N(0,1)$ if the two intervals have the same distribution, why does interval A have to cover the burn-in to give an indication of non-convergence? [1] Sahlin, Kristoffer (2011),   ' _Estimating convergence of Markov chain Monte Carlo simulations_ ',   Masters thesis, Stockholms universitet, Stockholm


title: What are the proper ways how to add noise to a dataset? (algorithm data sensitivity analysis)

text: Consider: **dataset** defined as `n` datapoints `x_i` in `m`-dimensional space. And there is a **label** `y_i` defining one of the classes belonging to `x_i`. There are let's say **5 classes** 1,2,3,4,5 (and there is _total order_ among the classes, i.e. 1<2<3<4<5). What I want to do is to **analyse the sensitivity of the algorithm to noise** in the dataset. It means that I will sequentially add more noise to the dataset and check how good the classifier will be when learned on the noisy data. The question: **What is the proper way of adding (generating) the noise?** My personal guess is that I will need to _normalize_ the values and somehow add noise based on gaussian distribution. But I am not sure about the particular proper way. I don't want to make any statistical or mathematical mistake.


title: Minimum sample size for a dichotomous outcome

text: I have two questions.   1. I am running an experiment where I am interested in determining the **sample size** required for a certain CI and error, where values range between $<1$ and $>-1$. However, I am only interested in the minimum sample size required as per the sign + or - not as per the specific values. Could this problem be framed as a **binary** problem (+ or -) as to calculate the minimum sample size (or should cathegorical values be applied instead and how) and what would be the correct method?    2. I have a model where I split the data into **train** and **test** samples. In this case and considering the above (interest is in + and - sign only not on the specific values) I would like to find out what the optimal breakdown sample size would be for train / test. In other words, at what **training set size** the test error becomes significant when compared with the null hypothesis of a random prediction. In this respect I have reviewed power analysis, however I have doubts if this can be applied to a **binary test** (or for the same token if values are cathegorical) given that the size effect calculation involves a difference in means that could be different from the one calculated applying the actual values or result in 0 and also there is the issue of on what values to calculate the standard deviation. In this situation how would I determine the sample size for training.


title: What is a consistency check?

text: I was asked such a question as "Did you do any consistency check in your daily work?" during a phone interview for a Biostatistician position. I don't know what to answer. Any information is appreciated.


title: How to exclude interaction effect in linear regression modelling?

text: I am working with biological dataset. I am looking for significantly differentially expressed genes that stand out under effect of 3 different factors (e.x. gut microbiota level has/not, tissue, Lcells has/not). The data is not good so I cannot perform normal ANOVA test to measure effect of all factors at once. So I have done linear regression modeling for each factor separately. Then I compare all one-way linear regression and found those genes that were significantly expressed only in one factor comparison (e.x. only based on having or not having gut microbiota) were selected. Now I am curious to know if I can clam these are significantly expressed by effect of this factor (e.g. gut microbiota)? Because I know in ANOVA test there is also one test for interaction which I might don’t have it here!


title: What kind of test to use for measuring a factor during time?

text: I am working with microbiological data. I want to test how gut microbiota is developing in my samples. I have one base line, which doesn’t have any gut microbiota, and then by day 1 they receive microbiota and I am monitoring the development of microbiota till day 7. The biological question is: How gut microbiota is developing and this is measuring by finding significantly differentially expressed genes at each step. I have data from day 1, day 3 and day7. Since the number of replicates in each group is low (only 3 per each) I cannot use normal ttest. So I went for linear regression modeling. Up to now I tried to do linear regression modeling between each day and the base line. But I was wondering if this is a good approach for this kind of analysis or I should consider other tests? Thanks in advance! Rozita


title: Calculating Orwin's (1983) modified Fail-safe N in a meta-analysis with Odds Ratio as summary statistic?

text: I am attempting to calculate Orwin's (1983) modification of Rosenthal's (1979) Fail-safe N for my meta-analysis of Odds Ratios. However, all the equations I am finding are using Cohen's d, which I cannot calculate (I don't have two groups). I have STATA, SPSS, and MedCalc at my disposal.


title: From logistic regression coefficients to marginal distributions of the dependent variable

text: I have the coefficients for a logistic regression model with $25$ independent variables of mixed continuous and categorical types. Considering one particular categorical independent variable with $2$ levels so that in the linear model we have a term $x_{i1}$ taking a value of either $0$ or $1$ and an associated coefficient $\beta_1$, is it possible to infer from the value of $\beta_1$ what proportion of $y_i$ in the data used to train the model had a value of $1$?


title: How to calculate sample size for comparing the area under the curve of two models?

text: Because I would like to calculate the sample size for comparing the area under the curve (AUC) of 2 models (cross-sectional study, predictor = continuous variable). Can you point me which function in R solves it.


title: How to get LogLiklihood value from logistic regression in R

text: ### Question:   * Is it possible to get log liklihood values for my stepwise glms? ### Context: I am able to get a logliklihood value using lmer with the following model. My study involves unbalanced repeated females, two sites (females don't exchange between sites), 8 predictors, and a response.                   (glmfit1 <- lmer(Response ~ 1 + SITE + (1|SITE:FEMALE) + Variable A + Variable B +                             Variable C, data = data))      Removing the repeated females and taking out the Site factor gives me this formula.               (lrfit1 <- glm(Response ~ 1+Variable A + Variable B + Variable C, data = data))     summary (lrfit1) #just gives me p values for my variables and an AIC but NO LOGLIK.


title: Free and concise statistics tutorials to give out as a refresher

text: Is anyone aware of a statistics resource (preferably 20 to 40 pages maximum) that reviews basic stats for people who took statistics classes already? This resource could be handed out as a refresher to those who need it. The reason why I'm not looking for a book is that I find that people will more likely read a 20/40 pager than a 500 page book that goes into too much detail for the intended scope. The ideal resource will explain the statistical methods, hypotheses, various methods, etc. It has to cover methods like Chi-squared. It has to be written in an easy to read and digest manner.


title: A critical proof or counterexample regarding independence

text: Does independence of $X^2$ and $Y^2$ imply independence of $X$ and $Y$?


title: How to calculate confidence intervals for pooled odd ratios in meta-analysis?

text: I have two datasets from genome-wide association studies. The only information available are the odd ratios and their confidence intervals (95%) for each genotyped SNP. My want to generate a forest plot comparing these two odds ratios, but I can't find the way to calculate the combined confidence intervals to visualize the summary effects. I used the program PLINK to perform the meta-analysis using fixed effects, but the program did not show these confidence intervals.   * How can I calculate such confidence intervals? The data available is:   * Odd ratios for each study,    * 95% confidence intervals and    * Standard errors.


title: How do you choose MAX or MIN values from a range of incomplete data?

text: I am currently working with large datasets that have missing data points. My data ranges from the years 1979-2009. I want to choose, by year, the maximum value from that year. Missing values are replaced with a "." I am currently using the formula: {=MAX((E13:E2863=1979)*(H13:H2863))} This array returns the max value where the appropriate data is available. However when i use the array over the entire dataset, a #VALUE error occurs within the box. My question is how do you overcome this issue? Is there a way to develop a formula that will not take into account for the missing data points i.e. cells with "." Below i will post a screenshot representing my spreadsheet/issue.               B       C    D     E     F =MAX((B6:B18=1979)*(E6:E18))      1979    0    53    67    38         1979    0    57    72    41         1979    0    63    .     47         1979    0    64    79    49         1979    0    65    78    52         1979    0    57    71    43         1980    0.37 43    52    33         1980    0    47    60    34         1980    0    55    70    39         1980    0    64    79    48         1980    0    66    82    50         1980    0    69    87    51         1981    0    71    85    57      Any suggestions?


title: Manually calculated chi-squared statistic differs from that given by chisq.test in R

text: I'm using R to do a chi-squared test on two sets of data to check if they follow a distribution. In the first case the chi-squared statistic given by `chisq.test` matches the one calculated by hand, but in the other case it does not. First case               > e1 = c(18, 38, 78, 107, 107, 72, 35, 15)     > o1 = c(10 ,30 ,80 ,180 ,60 ,40 ,40 ,30)     > sum((o1-e1)**2/e1)     [1] 105.6762     > chisq.test(o1,p=e1,rescale.p=T)$statistic     X-squared      105.6762       Both methods give `105.67` as the statistic. But in the other case:               > o2=c(401, 1235, 2989, 5682, 8489, 9966, 9196, 6668, 3800, 1701, 599, 165)     > e2=c(385, 1204, 2945, 5643, 8472, 9966, 9186, 6635, 3754, 1665, 578, 157)     > sum((o2-e2)**2/e2)     [1] 5.111822     > chisq.test(o2,p=e2,rescale.p=T)$statistic     X-squared      3.301292       Now the statistics are different and I don't understand why. Am I doing something wrong or missing something?


title: How to code and analyze skiped questions with SPSS

text: sorry if my question may seem rather silly, but I'm not well up on statistics and need a piece of advice. I have some cases in my dataset where respondents have abandoned the survey half-way, so I coded the missing values as 9999 "Refused to answer". In addition, for likert-scale type of items I applied the code 8888 for "Don't know" in the data set. Finally, I have some skip patterns: some variables of these skip patterns are categorical and some are scale and ordinal. For the legitimately skipped answers, I have coded them 7777 "Legitimate skip" in the data set. I've set up all these codes in the "Values code" column in SPSS. My questions are:   * For the codes 8888 and 9999, should I define them as "Discrete missing values" in the "Missing" column in SPSS as well?   * For the code 7777, how should I consider it, as a missing value as well? (I don't think "Legitimate skips" are actually missing values in the proper sense of the word, but if I don't code it as Missing, then for scale variables it will be a problem, as it will be added up in means, etc. right?) Thank you to anyone who can help!


title: MATLAB: leave-one-out-cross-validation for PCA

text: I'm trying to write my own function for PCA (of course there's a lot already written but I'm just interested in implementing stuff by myself). The main problem I encountered is the cross-validation step and calculating PRESS value (actually doesn't matter which cross-validation I use. It's a question mainly about the theory behind). From the theory I found out that in order to perform LOOCV you need to:   1. delete an object    2. scale the rest    3. perform PCA with some number of components    4. scale the deleted object according to parameters obtained in 2)    5. predict the object according to the PCA model    6. calculate PRESS for this object    7. re-perform the same algorithm to other objects    8. sum up all the PRESS values    9. profit  Because I'm very new in the field, in order to be sure that I'm right, I compare the results with the output from some software I have (also in order to write some code I follow the instructions in the software). I get completely the same results calculating the RESS and R^2 but calculating PRESS is a problem. Could you please tell to me if what I implement in the cross-validation step is right or not:               case 'loocv'          % n - number of objects     % p - number of variables     % vComponents - the number of components used it CV     dataSets = divideData(n,n);               % it's just a variable responsible for creating datasets for CV               %  (for LOOCV datasets will be equal to [1, 2, 3, ... , n]);'     tempPRESS = zeros(n,vComponents);          for j = 1:n       Xmodel1 = X; % X - n x p original matrix       Xmodel1(dataSets{j},:) = []; % delete the object to be predicted       [Xmodel1,Xmodel1shift,Xmodel1div] = skScale(Xmodel1, 'Center', vCenter,                                                    'Scaling', vScaling);                % scale the data and extract the shift and scaling factor       Xmodel2 = X(dataSets{j},:); % the object to be predicted       Xmodel2 = bsxfun(@minus,Xmodel2,Xmodel1shift); % shift and scale the object       Xmodel2 = bsxfun(@rdivide,Xmodel2,Xmodel1div);       [Xscores2,Xloadings2] = myNipals(Xmodel1,0.00000001,vComponents);                % the way to calculate the scores and loadings                     % Xscores2 - n x vComponents matrix                     % Xloadings2 - vComponents x p matrix       for i = 1:vComponents         tempPRESS(j,i) = sum(sum((Xmodel2*(eye(p) -                                             Xloadings2(1:i,:)'*Xloadings2(1:i,:))).^2));       end     end     PRESS = sum(tempPRESS,1);


title: Least Squares Question to show proportion

text: I have a question. Let's say data has been collected for people with jobs and people without jobs. People with Jobs are represented by 1, people without jobs are represented by 0. How do you use a least square method to show the proportion of people with jobs or the proportion of people without jobs? Thanks


title: Adjusted vs. unadjusted effects in regression

text: Is "unadjusted" basically just simple linear regression whereas "adjusted" is multiple regression? For example, looking at the effect of x on y adjusting for other variables like a, b and c versus not adjusting for them.


title: When to use ridge estimator / naive Bayes

text: I used the Logistic function in weka, to predict a binary class. I have used SimpleLogistic before, but Logistic also seem to give me good results. I did want to clarify if I understand some things well on this classifier, pecially about the ridge estimator Logistic does multinomial logistic regression with a ridge estimator.   * multinomial means for more then binary classes as well. So it would appear I didn't really need this function   * ridge estimator: I have been reading about this, but it is hard to summarise in one sentence. When attributes are highly correlated this might seem a good technique because XXt becomes close the zero otherwise. Is this correct? Do you have a better explanation in a few sentences?   * If my interpretation of ridge is correct. It seems to be valid in exactly the opposite case as when you should use Naive Bayes, because this assumes independent variables. Thanks for validating and providing a short explanation of ridge (which I don't fully get)


title: Why is it possible to have a non-significant MANOVA but multiple significant univariate ANOVAs?

text: The only source I have for this being possible is http://en.wikiversity.org/wiki/Multivariate_analysis_of_variance, but this page does not explain why it is possible.


title: Non-parametric tests to compare numeric profiles

text: I am trying to compare two different numerical profiles and determine whether they are the same or not, computing a p-value. These profiles are composed of 200 values each (each of which corresponds to a particular position in a DNA strand). They do not have any particular distribution so the first thing i tried was a paired Wilcoxon test in which the first position of the first profile is compared to the first position of the second profile and so on. This however proved to be way too sensitive and even randomly generated profiles (which I would use to simulate a background distribution) are all significantly different. Here's an example of two randomly generated profiles: 1-               [1] 0.05410946 0.05405054 0.05399161 0.05393269 0.05387377 0.05381484 0.05375470 0.05369350 0.05363223 0.05357187 0.05351339     [12] 0.05345778 0.05340602 0.05335908 0.05331795 0.05328360 0.05323382 0.05317267 0.05310136 0.05305059 0.05297459 0.05289146     [23] 0.05278836 0.05269685 0.05259318 0.05247068 0.05235686 0.05227757 0.05225131 0.05230355 0.05233573 0.05237550 0.05245256     [34] 0.05258071 0.05273305 0.05280477 0.05280422 0.05281129 0.05281803 0.05269052 0.05247504 0.05226766 0.05219812 0.05219966     [45] 0.05228203 0.05238135 0.05265886 0.05308829 0.05343443 0.05376170 0.05415894 0.05451981 0.05468740 0.05465393 0.05445623     [56] 0.05408778 0.05361698 0.05306193 0.05241602 0.05184365 0.05128075 0.05081200 0.05039226 0.05011328 0.04988996 0.04971936     [67] 0.04959749 0.04950916 0.04934904 0.04918447 0.04901725 0.04885474 0.04864097 0.04846153 0.04834330 0.04825480 0.04817072     [78] 0.04809311 0.04804392 0.04805756 0.04809029 0.04809125 0.04808202 0.04805034 0.04799174 0.04789605 0.04778036 0.04766744     [89] 0.04756477 0.04745339 0.04735238 0.04728329 0.04724986 0.04722612 0.04718636 0.04713367 0.04707984 0.04701450 0.04693417     [100] 0.04682481 0.04671401 0.04658613 0.04643983 0.04628397 0.04609772 0.04590572 0.04571484 0.04554022 0.04539498 0.04526802     [111] 0.04516986 0.04507842 0.04499786 0.04490481 0.04480493 0.04469124 0.04460219 0.04450903 0.04444782 0.04440287 0.04436476     [122] 0.04431326 0.04428163 0.04426604 0.04425380 0.04421425 0.04420969 0.04421636 0.04423634 0.04423125 0.04425017 0.04424474     [133] 0.04425620 0.04423391 0.04416950 0.04408473 0.04402678 0.04395341 0.04383386 0.04371769 0.04355977 0.04340431 0.04327231     [144] 0.04310986 0.04292362 0.04272375 0.04256255 0.04236743 0.04221407 0.04210954 0.04205756 0.04202356 0.04198690 0.04191072     [155] 0.04185391 0.04178650 0.04175152 0.04179372 0.04189675 0.04198836 0.04208730 0.04214621 0.04222863 0.04226568 0.04228232     [166] 0.04230990 0.04241991 0.04255647 0.04260372 0.04263997 0.04266844 0.04263244 0.04257633 0.04251360 0.04244118 0.04244753     [177] 0.04244423 0.04239680 0.04233157 0.04226786 0.04219283 0.04210845 0.04206977 0.04204039 0.04201743 0.04200761 0.04197215     [188] 0.04193485 0.04189571 0.04185473 0.04181238 0.04177013 0.04172815 0.04168658 0.04164562 0.04160541 0.04156525 0.04152508     [199] 0.04148491 0.04144474 0.04140457      2-               [1] 0.05351192 0.05322555 0.05293917 0.05265280 0.05236643 0.05208005 0.05179246 0.05150552 0.05122072 0.05093953 0.05066343     [12] 0.05039391 0.05014104 0.04988936 0.04965494 0.04943778 0.04923900 0.04905929 0.04892465 0.04879003 0.04869032 0.04860967     [23] 0.04854307 0.04849056 0.04846746 0.04845587 0.04845305 0.04846037 0.04846948 0.04846119 0.04847131 0.04847203 0.04847388     [34] 0.04849018 0.04851020 0.04851629 0.04852339 0.04853276 0.04853965 0.04855213 0.04856085 0.04857185 0.04858035 0.04858784     [45] 0.04858446 0.04857930 0.04856101 0.04852577 0.04850898 0.04849094 0.04847929 0.04846053 0.04843542 0.04840443 0.04837794     [56] 0.04834002 0.04829689 0.04827229 0.04824362 0.04817578 0.04810853 0.04801922 0.04794492 0.04786671 0.04777944 0.04765683     [67] 0.04752808 0.04740640 0.04728304 0.04720053 0.04715555 0.04714327 0.04714215 0.04719154 0.04721849 0.04727323 0.04733321     [78] 0.04741521 0.04750358 0.04758516 0.04759874 0.04754620 0.04747679 0.04740342 0.04735974 0.04742875 0.04755039 0.04778289     [89] 0.04813493 0.04859605 0.04926145 0.04997775 0.05067520 0.05147505 0.05233079 0.05298947 0.05340390 0.05362887 0.05380001     [100] 0.05369953 0.05334439 0.05268226 0.05213193 0.05170225 0.05093835 0.04997205 0.04910467 0.04843989 0.04785735 0.04722082     [111] 0.04668176 0.04621078 0.04588176 0.04535853 0.04488548 0.04472865 0.04490304 0.04508753 0.04531480 0.04567468 0.04595480     [122] 0.04619730 0.04663426 0.04696277 0.04734696 0.04752347 0.04754857 0.04731406 0.04696973 0.04642326 0.04572888 0.04515736     [133] 0.04468316 0.04414105 0.04369353 0.04325270 0.04294296 0.04274710 0.04266621 0.04269275 0.04284180 0.04304411 0.04319943     [144] 0.04330029 0.04331458 0.04325282 0.04322687 0.04322829 0.04317383 0.04307845 0.04297811 0.04284321 0.04272740 0.04262550     [155] 0.04258994 0.04260521 0.04265048 0.04270060 0.04271161 0.04269050 0.04266664 0.04264355 0.04264804 0.04263363 0.04261875     [166] 0.04258121 0.04252872 0.04249338 0.04243412 0.04236218 0.04232400 0.04228812 0.04224938 0.04221131 0.04217026 0.04213377     [177] 0.04209305 0.04204009 0.04198693 0.04192882 0.04186326 0.04178475 0.04170866 0.04163548 0.04156215 0.04148637 0.04141216     [188] 0.04133743 0.04126257 0.04118770 0.04111297 0.04103852 0.04096447 0.04089096 0.04081814 0.04074613 0.04067564 0.04060515     [199] 0.04053466 0.04046417 0.04039368      My question is whether there is a better/more correct way to compare profiles and determine their similarity. Keep in mind that this comparison is part of a larger plan to determine a background distribution by simulating random profiles via Monte Carlo sampling so any advice on how to execute that operation is welcome.


title: How to do correlation in MATLAB

text: I'm trying to do correlation between two .wav files or two .pcm files using MATLAB. I just can't figure out how to do that! If anyone can help me I would appreciate it!


title: Diagonal elements of the projection matrix

text: I am having some problem trying to prove that the diagonal elements of the hat matrix $h_{ii}$ are between $1/n$ and $1$. Suppose that $Range(X_{n,k})=K $ the number of columns of our matrix of data with a constant.⇒$H_{k,k}$ $H=X(X' X)^{-1}X' ⇒ H'=H ;H^{2}=H $ If $y = \beta x + \epsilon ⇒ HY =\hat Y; (I-H)=\epsilon $ $\mathbf{H}=\begin{bmatrix}h_{11} &... &h_{1n}\\⋮ & ⋱ &⋮\\ h_{n1} & ... & h_{nn}\end{bmatrix}$ If $\boldsymbol 1=(1, \ldots, 1) ∈ X ⇒ 1H^2=1H=1 ⇒ ∑ h_{i1}^2=∑ h_{i1}=1 , h_{ii}≤1$ So $∑h_{i1}^2∑h_{i2}^2...∑h_{in}^2=1 ⇒h_{11}^2h_{22}^2...h_{nn}^2<=1 $ How to prove that $h_{ii} >=(1/n)$? (exercise 3.4 from Meyer "Classical and modern regression with applications") Let $h_{ii}$ be the ith diagonal of the Hat matrix H. (a) prove that for a multiple regression model with a constant term, $h_{ii} >=(1/n)$ ; (b) show that $h_{ii}≤1$ (Hint: make the use of the fact that H is idempotent)


title: Combining evidence using Dempster-Shafer theory

text: Can someone post a simple explanation of Dempster-Shafer theory? There are lot of links available but the reading material in those sites is academic in nature and time consuming to read and comprehend. What I'm trying to do is calculate similarity scores between two words. I am using 3 different similarity measures. Next, I need to combine these scores using Dempster-Shafer theory.


title: Choosing non-informative priors

text: I am working on a model relying on an ugly parametrized function acting as a calibration function on a part of the model. Using a Bayesian setting, I need to get non-informative priors for the parameters describing my function. I know that ideally, I should derivate reference or at least Jeffreys priors but the function is very ugly, have many parameters and I am pessimistic on the possibility to obtain actually a result. So I decided to drop this possibility and to empirically choose my priors prying for them to be quite non- informative. Here are my two questions.   1. Can I make more than prying and give insights of their non-informativness from inference results ? **Edit :** I guess that plotting posterior Vs prior would be a first point. Maybe comparing MAP and ML estimations could be a second argument ?   2. Moreover, is that make sense to justify some aspect of the choice from a "dimensional analysis"? As an example, if I consider a likelihood structure of the form (in a simple regression setting): $$ Y | a,b,x = a.x+b.e^{-x} + \epsilon $$ Do you think that I can guess any "structure" for the prior on $a$ and $b$ based on the fact that one weighs $x$ while the other weighs $e^x$ ?


title: Derive the conditional pdf of data on prior parameters

text: In Bayesian statistics I see this derivation often. Given the likelihood function $f(X|\theta)$ and the prior $f( \theta |a, b)$, the author will derive $f(X|a,b)$. The steps in between are considered trivial but I cannot derive $f(X|a,b)$ myself. Can anyone please provide some hints, perhaps how to "introduce" $\theta$ into $f(X|a,b)$?


title: Which regression analysis should I use for ranked dependent variables and proportional independent variables?

text: I am analysing the effect of deprivation on breastfeeding and am wondering which type of regression analysis I should use. It is area level data. Deprivation data is available as a score from 0 - 100 (but no scores actually reach either extreme value) or alternatively as ranked data i.e. 1 = most deprived, 151 = least deprived (there are 151 areas). Breastfeeding data is available in percentages e.g. percentage of maternities where bf was initiated. Again, values do not tend to hit the extremes of 0 and 100. Most of the confounders are also expressed as percentages. I am confused as to how to incorporate the proportional nature of the dependent variables with the ranked/ordinal nature of the independent variables. Analyses will need to be performed in Stata. Thank you for any help!


title: How does a frequentist calculate the chance that group A beats group B regarding binary response

text: ... (optional) within the context of Google Web Optimizer. Suppose you have two groups and a binary response variable. Now you get the following outcome:   * _Original_ : 401 trials, 125 successful trials    * _Combination16_ : 441 trials, 141 successful trials The difference is not statistically significant, however one can calculate a probability that Combination16 will beat Original. To calculate "Chance to beat Original" I have used an bayesian approach, i.e. performing a two-dimensional monte carlo integration over the bayesian-style confidence intervals (beta-distribution, (0,0) prior). Here is the code:               trials <- 10000     resDat<-data.frame("orig"=rbeta(trials,125+1,401-125+1),                         "opt"=rbeta(trials,144+1,441-144+1))     length(which(resDat$opt>resDat$orig))/trials      This results in 0.6764. Which technique would a frequentist use to calculate "Chance to beat ..." ? Maybe the power function of Fisher's exact test ? **Optional: Context of Google Web Optimizer** Google Web Optimizer is a tool for controlling multivariate Testing or A/B-Testing. This only as an introduction since this should not matter for the question itself. The example presented above was taken from the explanation page of Google Web Optimizer (GWO), which you can find here (please scroll down to the section " _Estimated Conversion Rate Ranges_ "), specifically from figure 2. Here GWO delivers 67.8% for "Chance to beat Original", which slightly differs from my result. I guess Google uses a more frequentist-like approach and I wondered: What could it be ? **EDIT:** Since this question was close to disappear (I guess because of its too specific nature), I have rephrased it to be of general interest.


title: How to model logistic regression with continuous predictor in Bugs?

text: Is it possible to use a continuous predictor in Bugs? The simplest way of doing this would be turning the size variable in alligators example from discrete to continuous. Both Winbugs and JAGS examples use combination of values of covariates as indices as in               X[i,j,] ~ dmulti( p[i,j,] , n[i,j]  );      where `i` is the lake index (4 possible values), and `j` is the size index (2 possible values). With this approach, a continuous size variable would mean an infinite amount of indices. There must be something I'm missing here.


title: How to create n time series characterised by a defined average and correlation?

text: I have a time series (X) representing a natural phenomenon (wind speed, measured every 15 minutes) and I have to create similar time series (up to 20, Xdi, i=1,...,20) with the same structure (same average, same standard deviation, same percentiles distribution...) but with a predetermined correlation (about 0.7) between each other. Is there any defined method for this operation? Can you provide a link to a book, a paper, a page, or just the name of the method? Thank you very much, Andrew * * * A couple of clarifications: for "pre-defined correlation" I mean that if I have a seed time series X and I want to create three derived time series Xd1, Xd2, Xd3, then the correlation among two of the time series must be equal (almost) to a value chosen (e.g.: 0.7). For example Correlation(Xdi, Xdj) = 0.7 The comment of charles.y.zheng answers the question (thanks!), but the resulting time series (AZ in the example) does not have necessarily the same autocorrelation of the original seed. (in the original time series 0 tends to be followed by 0 and 1 tends to be followed by 1, the values are clustered). Sorry to haven't specified also this requirement, but I noticed after trying the solution proposed. I guess that I can fix the problem by manipulating the values in Z. This adds a second section of the question: Is there a defined method to create a time series with a pre determined autocorrelation? p.s.: sorry if the language is not correct, I'm not a statistician but a programmer, I try to do my best but if there is something not clear, just ask for details and clarifications.


title: Clustering spatial data in R

text: I have a set of sea surface temperature (SST) monthly data and I want to apply some cluster methodology to detect regions with similar SST patterns. I have a set of monthly data files running from 1985 to 2009 and want to apply clustering to each month as a first step. Each file contains gridded data for 358416 points where approximately 50% are land and are marked with a 99.99 value that will be NA. Data format is:                  lon     lat   sst     -10.042  44.979  12.38      -9.998  44.979  12.69      -9.954  44.979  12.90      -9.910  44.979  12.90      -9.866  44.979  12.54      -9.822  44.979  12.37      -9.778  44.979  12.37      -9.734  44.979  12.51      -9.690  44.979  12.39      -9.646  44.979  12.36      I have tried CLARA clustering method and got some apparently nice results but it also seems to me that is just smoothing (grouping) isolines. Then I am not sure this is the best clustering method to analyse spatial data. Is there any other clustering method devoted to this type of datasets? Some reference would be good to start reading. Thanks in advance.


title: Problem building multinomial logit model formula on huge data in R

text: I am attempting to build a Multinomial Logit model with dummy variables of the following form:   * The dependent variable represents 0-8 discrete choices.   * Dummy Variable 1: 965 dummy vars   * Dummy Variable 2: 805 dummy vars The data set I am using has the dummy columns pre-created, so it's a table of 72,381 rows and 1770 columns. The first 965 columns represent the dummy columns for Variable 1; the next 805 columns represent the dummy columns for Variable 2. I'm on a Sun Grid Machine at my university, so memory won't be an issue... I have been able to generate the factors and generate `mlogit` data using code:               mldata<-mlogit.data(mydata, varying=NULL, choice="pitch_type_1", shape="wide")      my `mlogit` data looks like:               "dependent_var","A variable","B Var","chid","alt"     FALSE,"110","19",1,"0"     FALSE,"110","19",1,"1"     FALSE,"110","19",1,"2"     FALSE,"110","19",1,"3"     FALSE,"110","19",1,"4"     TRUE,"110","19",1,"5"     FALSE,"110","19",1,"6"     FALSE,"110","19",1,"7"     FALSE,"110","19",1,"8"     FALSE,"110","19",2,"0"     FALSE,"110","19",2,"1"     FALSE,"110","19",2,"2"     FALSE,"110","19",2,"3"     FALSE,"110","19",2,"4"     FALSE,"110","19",2,"5"     TRUE,"110","19",2,"6"     FALSE,"110","19",2,"7"     FALSE,"110","19",2,"8"     TRUE,"110","561",3,"0"          ...      The `mldata` contains 651,431 rows. If I try to run this full data set I get the following error:               > mlogit.model<- mlogit(dependent_var~0|A+B, data = mldata, reflevel="0")     Error in model.matrix.default(formula, data) :     allocMatrix: too many elements specified     Calls: mlogit ... model.matrix.mFormula -> model.matrix -> model.matrix.default     Execution halted      Smaller datasets (`mldata` with only 595 rows) and `mlogit` works fine and generates the expected regression output. Is there a problem with `mlogit` and huge datasets? I suppose this is perhaps not the best way to assess this kind of data, but I am trying to replicate a previous analysis that was completed on a similar amount of similar data.


title: Using time series analysis to analyze/predict violent behavior

text: This is a bit of a flippant question, but I have a serious interest in the answer. I work in a psychiatric hospital and I have three years' of data, collected every day across each ward regarding the level of violence on that ward. Clearly the model which fits these data is a time series model. I had to difference the scores in order to make them more normal. I fit an ARMA model with the differenced data, and the best fit I think was a model with one degree of differencing and first order auto-correlation at lag 2. My question is, what on earth can I use this model for? Time series always seems so useful in the textbooks when it's about hare populations and oil prices, but now I've done my own the result seems so abstract as to be completely opaque. The differenced scores correlate with each other at lag two, but I can't really advise everyone to be on high alert two days after a serious incident in all seriousness. Or can I?


title: Bootstrapping in ANCOVA AND MANOVA in R

text: I have a difficult one for strong solvers Can anyone provide us please an example with coding in R of the bootstrapping method in ANCOVA and one in MANOVA Thank you very much T


title: Covariance of a random vector after a linear transformation

text: If $\mathbf {Z}$ is random vector and $A$ is a fixed matrix, could someone explain why $$\mathrm{cov}[A \mathbf {Z}]= A \mathrm{cov}[\mathbf {Z}]A^\top.$$


title: How do I extract a keyword plus its corresponding number out of a text

text: I am finding my way around in the different text mining tools, but I can't find the technique that does the following: Extract out of the follwing advertisement texts               "apartment with bedrooms: 2, bathrooms: 2"      "appartment with 2 bathrooms"      "appartment with 2 beaufitul batrooms"      that these three appartments all have 2 bathrooms. The text of the advertisement is already tokenized. What are my options now? Do I have to start writing algorotimhs like _"find the word bathroom and a number just before it"_? Or are there better techniques?


title: Comparing and visualizing densities of two samples in R (2)

text: Comparing the density function of two specific samples in R, both areas are clearly different and cannot be equal to one (see graph linked below). I would like both areas equal one on the same graph as to determine easily the set of abscissa values for which the pdf of a sample is larger than the pdf of the other. With my script, one of the pdfs largely overlaps the other, which is almost invisible. How can I fix it? Thanks Here are the pdfs I get when using plot() and lines(): ![Figure](http://i.stack.imgur.com/X82lw.jpg) Here are the x-values of the sample #1 (9 values):               0.043, 0.022, -0.000015,  12.98, 0.316,  -1.785, 10.410, 0.633, -0.0058      Here are the x-values of the sample #2 (448 values): Google shared text file I use a similar script to this one (provided by ashkan):               x <- rnorm(1000, mean=3, sd=2)      y <- rnorm(500, mean=3.5, sd=3)     dx <- density(x)     dy <- density(y)     plot.new()      plot.window(xlim=range(c(dx$x, dy$x)),      ylim=range(c(dx$y, dy$y)))     with(dx, lines(x, y))     with(dy, lines(x, y, lty=2))


title: Clustering with probabilities

text: I will ask my question beginning by an example (I am novice in stats) I have a set of probabilities from a given observation A={0.3,0.2,0.001, 0.02, ...} I want to partition A on subsets or clusters so as each one will contain the closest probabilities. And after, i will be interested in calculating the distances between the clusters. Is there any simple algorithm I can use for this? Thanks


title: What's the deal with autocorrelation?

text: To preface this, I have a pretty deep mathematical background, but I've never really dealt with time series, or statistical modeling. So you don't have to be very gentle with me :) I'm reading this paper about modeling energy use in commercial buildings, and the author makes this claim: > [The presence of autocorrelation arises] because the model has been > developed from time series data of energy use, which is inherently > autocorrelated. Any purely deterministic model for time series data will > have autocorrelation. Autocorrelation is found to reduce if [more Fourier > coefficients] are included in the model. However, in most of the cases the > Fourier model has low C.V. The model may, therefore, be aceptable for > practical purposes that does (sic) not demand high precision. 0.) What does "any purely deterministic model for time series data will have autocorrelation" mean? I can vaguely understand what this means--for example, how would you expect to predict the next point in your time series if you had 0 autocorrelation? This isn't a mathematical argument, to be sure, which is why this is 0 :) 1.) I was under the impression that autocorrelation basically killed your model, but thinking about it, I can't understand why this should be the case. So why is autocorrelation a bad (or good) thing? 2.) The solution I've heard for dealing with autocorrelation is to diff the time series. Without trying to read the author's mind, why would one _not_ do a diff if non-negligible autocorrelation exists? 3.) What limitations do non-negligible autocorrelations place on a model? Is this an assumption somewhere (i.e., normally distributed residuals when modeling with simple linear regression)? Anyway, sorry if these are basic questions, and thanks in advance for helping.


title: Using a time series model to forecast future values in R

text: I would like to use a set of weather-related historical data to fit a time series (let's say 1970-2000, Fourier terms plus ARIMA terms), but then use the fit on recent data (i.e., the last week/month of data) to forecast the upcoming day/week/month. All of the functions that I've found forecast from the endpoint of the dataset used to fit. Can someone point me in the right direction? Or let me know it doesn't exist and I have to write it out the long way (i.e. `Tomorrow = fit$coef[1]*Yesterday + ...`)?


title: Statistical Approach to restricting predictions between 0 and 1

text: I am working with tobacco prevalence data where I am smoking rates at the country level. I and am in the process of adjusting estimates across different frequencies of use(ie, daily versus occasional smokers). There are some surveys which report both of these values (ie, daily smoking = 30%, occasional smoking = 40%), so I evaluated the relationship between occasional and daily smoker via linear regression:               daily_smokers ~ B + B1*occasional_smokers + B*contols_variables      However, I am prediciting prevalence (restricted to 0-1), so I put both the outcome variable of interest and the predictor variable into logit space (putting them both into logit because I believe the relationship is linear). Yielding:               logit(daily_smokers) ~ B + B1*logit(occasional_smokers) + B*contols_variables      However, when I adjust the data using the betas, the adjustment is inconsistent (ie, not all points shifted the same direction in real space): ![enter image description here](http://i.stack.imgur.com/0hSuG.png) What is an appropriate mechanism for restricting my outcome variable range while applying a consistent adjustment? Thanks, Simulation Code in STATA               clear     set obs 10000     gen prev = rnormal(.5,.1)     gen logit_prev = logit(prev)     gen adjusted_logit = logit_prev*.78     gen adjusted_prev = invlogit(adjusted_logit)     gen diff_norm = prev - adjusted_prev     tw scatter adjusted_prev prev || function x=x, range(prev) xtitle("Unadjusted") ytitle("Adjusted") legend(off)


title: Standard error of estimates of covariance parameterized in tems of cholesky

text: In unconstrained optimization, covariance matrix $C$ is parameterized in terms of its Cholesky ($C=LL^\prime)$. In other words, the parameter vector $\theta$ involves elements of the lower triangular matrix $L$. Derivative of $-2LL$ with respect to $\theta$ is computed exactly. Hessian at the final-iteration is computed numerically (central difference). Analytic derivatives of $C$ with respect to $L$ are also available (Jacobian?). The goal is to compute variance of elements of computed covariance $var(C_{ii})$ and $var(C_{ij})$.


title: how to impute missing monthly climate data using AMELIA package?

text: I have been trying to impute future monthly climate data based on sets of downscaled data using the AMELIA package. I have monthly precipitation data from 1960-2099 for different geographic regions and I want to impute precipitation data from 2099-2150 using the past trend of precipitation. I want to use bound feature (i.e., max and min values) of different months of a year specific to the geographic region to impute monthly data from 2099-2150 but I can't find ways to do it. I have figured out how to set one bound (Max and min values) for a particular geographic region but can't figure out to set multiple bounds so that data are imputed using that bound for a given month. Is there any way to perform such task? I know there are several limitations behind using AMELIA package to impute climate data, but I was wondering if people have attempted to impute climate data using AMELIA package? I find none. I will appreciate any suggestions/comments on this post. Thanks!


title: Prove that $\limsup_{n\to\infty} |X_n|/n \le1 $ almost surely

text: Suppose {$X_n$} a sequence of random variables. If $\sum_{n=1}^{\infty}P(|X_n|>n)<{\infty}$ Prove that $$\limsup_{n\to\infty}\frac{ |X_n|}{n} \le1 $$ almost surely **What i have done so far:** I thought using the Borel-Cantelli lemma could lead me somewhere, but i didn't have any luck. From Borel-Cantelli lemma we know that if $\sum_{n=1}^{\infty}P(|X_n|>n)<{\infty}$ then $P(|X_n|>n)=0$ How could I proceed? I would appreciate any help, advice. Thank you all very much in advance for your time and concern.


title: How to compare observed distribution to sampled distributions, using cumulative distibution function?

text: I have a distribution of observed measurements and I want to compare it to sampled distributions, using R. I have a program that samples distributions, according to a certain low of probability / constraint. Suppose I have:   * $m$ observed measurements   * two sets of $n$ samples of $m$ measurements I'd like to compute:   * the average cumulative distribution function (a-cdf) of the first set of samples (its value at each evaluation points is the mean value of the cdf of each sample at this evaluation point, and the evaluation points are the pooled $n\cdot m$ measurements)   * the 95% envelope of the first set of samples   * the spatial distribution index (sdi) of my observed measurements using the second set of $n$-samples. To compute the sdi, we compute all the maximum differences between the cdf of each new sample and the a-cdf, and the maximum difference between the observed cdf and a-cdf. That makes $n+1$ points, and the sdi is the rank of the observed difference among the $n+1$ differences.  I'd like to know if thoses function exists in R, and if not, what would be the smartest way to implement them. Hope i'm clear enough.


title: Good classifiers for small training sets

text: What are good classification methods in the case of continuous independent variables (features) and a small training set (particularly where the number of training examples is approximately equal to the number of independent variables)? Here, small means about 50. I am particularly interested in being able to know which variables are “significant”. Ideally, I'm looking for a method for which the training step is computationally efficient; I care much less about the computational cost for the actual classification task.


title: Determining Population Normalcy and Proper Analysis Techniques

text: Good Day - I am working with a dataset that is not normal (general ledger information). I am attempting to determine the following:   1. Given it is not normal, what is(are) the best method(s) to determine variance (e.g., from median or mean) in a meaningful way?   2. How can I measure how valid various statistical tests are for such a population?   3. What tests should I use, if those applying to a Normal population do not provide valid, usable descriptive statistics? For example, if I calculate mean, median, STDDEV and variance on the dollar amount of the transaction, and I find a small set of results that is 27 standard deviations from the mean (or median, for that matter), can I assert that this is truly as much of an outlier as it would be in a Normal population? I strongly suspect not, but don't know how to quantify the 'how much less valid is it' or how to analyze the population to get more better results - in this context, better being defined as 'the same level of significance in our population as one would get from 27 standard deviations in a Normal population'. Ultimately, we are doing this analysis to identify transactions to test as part of our internal audit work - we want to 'stack the deck' of our samples with transactions that are unusual by some defined Risk Indicator. For example, the frequency of expense reports submitted with a total amount of \$24.99 - $.01 less than our 'requires a receipt' threshold. There are many such defined Risk Indicators, the above is only an example - several others do in fact lend themselves to a more 'statistical' approach, as opposed to this simple stratification example. I am trying to come up with a way of presenting statistical analysis around this population to people who understand enough statistics to be dangerous - and my only claim to 'not being dangerous' is that I am questioning the validity of the findings, as in the above example. Any assistance is greatly appreciated.


title: What's the safest way to resample to ensure equal class frequencies in training data?

text: I am working on a number of EEG data sets for binary classification. A good example of one is publicly available here. If you look at what Matthias Kaper did to classify that set, one thing that jumps out is that he made sure to train with an equal number of samples from each class. Say you have such a data set, with, say, 300 trials from class 1, and 700 trials from class 2, with class frequencies determined by experiment (not inherent to the process generating the signals). If you sample, say, 70% of the trials for training uniformly, you end up with roughly the same frequencies in the training set. Let's say you still want 70% of the **total trial count** for training, but you also want equal class frequencies. If you want 50% of 70% of 1000 trials for each class, you need 350 trials from each class, but that's more than you have for class 1! Is there a safe way to resample up to this fraction of the available data while keeping class frequencies equal, or should you give up and use something like 70% of the smaller class trials (210 from each)?


title: Accounting for discrete or binary parameters in Bayesian information criterion

text: BIC penalizes based on the number of parameters. What if some of the parameters are some sort of binary indicator variables? Do these count as full parameters? But I can combine $m$ binary parameters into one discrete variable that takes values in $\{0,1,...,2^m-1\}$. Are these to be counted as $m$ parameters or one parameter?


title: Normal prior for Binomial likelihood

text: Pardon my ignorance, i am new to Bayesian Analysis. I am trying to use Normal prior for a binomial likelihood, which of these are most likely candidates ( $\bar{x} $, $ \mu $, $ \sigma $ ) for prior and most likely candidates for constants ? Please advise.


title: What is a main difference between RBF neural networks and SVM with a RBF kernel?

text: What is a main difference between RBF neural networks and SVM with a RBF kernel, from a practical point of view?


title: Test for Statistical Significance in the Accuracy of a Machine Learning System

text: I have what I imagine is an elementary question about evaluating statistical significance, but while I know a lot about probability I can't t-test my way out of a paper bag. From here I'm hoping to get a pointer to where I should look for the answer. I have a machine learning system and a test set of _N_ questions. I run the test and the system gets _r_ questions right and _w_ questions wrong where _r + w = N_. I measure the performance of my system with its _accuracy_ , which I define to be _r/N_. I make a change to my machine learning system and rerun the same test. Now I get a different accuracy. Is is possible to ask whether or not the change in accuracy is statistically significant? I suppose the null hypothesis is that the change in accuracy is just due to chance. The system is deterministic, so in theory there is no chance involved. However, certain questions will be hard ones that lie on its decision boundary, making them so sensitive to the system's configuration their correctness is essentially random. Is it meaningful to talk about statistical significance in this instance? If so, what test should I use? * * * Following up on the cross-validation suggestions in the comments below, would the following work? I have a test data set of size _N_. It is divided into _M_ disjoint partitions. To test my system, I calculate its accuracy on each of the _M_ partitions. Then I can take the mean and standard deviation of this accuracy. To compare the performance of two systems, I run them both on the same set of _M_ partitions, and then see if the difference between the mean accuracies is statistically significant. Would I use Welch's Test for this? Here the "randomness" in the accuracies arises from overfitting.


title: How to find local peaks/valleys in a series of data?

text: Here is my experiment: I am using the `findPeaks` function in the quantmod package: I want to detect "local" peaks within a tolerance 5, i.e. the first locations after the time series drops from the local peaks by 5:               aa=100:1     bb=sin(aa/3)     cc=aa*bb     plot(cc, type="l")     p=findPeaks(cc, 5)     points(p, cc[p])     p      The output is               [1] 3 22 41      It seems wrong, as I am expecting more "local peaks" than 3... Any thoughts?


title: How to find the distribution of the weighted sum of independent Bernoulli random variables for positive non-integer weights

text: How do I find the distribution of the weighted sum of independent Bernoulli random variables if the weights are non-negative real numbers? I have N number of independent Bernoulli distributed random variables lets say X1, X2, X3...XN and suppose I have a set of weights W1, W2, W3,...,WN which are non-negative real numbers and the sum of all the weights is equal to N then I need to find the distribution of the random variable Z which is the weighted sum of the N independent Bernoulli random variables i.e, Z =sum(Wi*Xi), where i=1,2,...,N please help! Thanks in advance!


title: Distribution of the t-statistic

text: Suppose we have a simple regression model $y=\beta_0+\beta_1x+u$ under Gauss- Markov assumptions, and $u$ has a uniform distribution on $[0,1]$. Find sampling distribution of the $t_{\beta_1}$-statistic. $t_{\hat{\beta_1}}=\frac{\hat{\beta_1}}{se(\hat{\beta_1})}$, where $se(\hat{\beta_1})=\frac{[(\sum\limits_{i=1}^n{\hat{u_i}}^2)/(n-2)]^{1/2}}{ (\sum\limits_{j=1}^n(x_i-\bar{x})^2)^{1/2}}$ We have that ${\hat\beta_1}=\beta_1+\sum\limits_{i=1}^nd_iu_i$, where $d_i=\frac{x_i-\bar{x}}{\sum\limits_{j=1}^n(x_i-\bar{x})^2}$. Conditionally on the sample $d_i$ can be treated as nonrandom variables, so $\hat{\beta_1}$ is constant plus the linear combination of independent uniform variables and it seems possible to find its distribution. Even if it is found, there still exists denominator of the $t$-statistic. Probably another approach should be taken here.


title: How to transform time-varying covariate measure of response time in a multi-level model of longitudinal data?

text: I am trying to fit a multi-level model to some longitudinal data that I have. As an example, let's pretend participants had to make 10 basketball free throws, and I measured how long it took them to make each one (in seconds). In this case, the time variable is severely positively skewed because although there is a minimum amount of time it takes to make 10 throws, some participants could take much longer than others. When trying to fit the model with the original data, I ran into all sorts of convergence errors using lmer in R. After transforming the time variable (in seconds) by applying a fourth-root transformation, the model fits fine, and the coefficients are interpretable. I'm just wondering if (a) this is a valid thing to do, and (b) what I should be aware of/careful about when doing this. edit: If someone could direct me to further reading on why the non-transformed model might fail to converge, that would be great!


title: Trying to reproduce predictive likelihood

text: I'm trying to reproduce the average predictive likelihood from page 13 (Eq 4.1) from the following paper and I'm struggling to find the function f_M,k(...). It details more about it on page 14 but I don't know how to actually produce numbers for it. Any advice or help would be awesome Paper http://www.atl-res.com/finance/MAHEU.pdf


title: Multiple measurements within groups test

text: I am measuring the X, Y, Z variables of i subjects in k states, and for each subject and state i am taking j measurements. What is the proper way to test for differences between the k states? It seems likes a one way repeated measure ANOVA problem, however how should I deal with the fact that in each subject and within each state i take j measurements of X, Y, Z? Should I take the mean of the Xj, Yj, Zj within each subject and state and with this mean do a one way repeated measures ANOVA? Many thanks?


title: First difference or log first difference?

text: I am evaluating the effect of covariances between series on returns. That is I run the following regression: $$ r_t = \beta_0 + \beta_1\text{Cov}(Y_t,r_t) + ...$$ I have conducted my analysis with both first difference and log(first difference) on the series. That is I can take either $r_t = P_{t+1}- P_{t}$ or $\ln(P_{t+1}/P_{t})$. (and similarly for $Y_t$) However, the level of significance of my coefficients is considerably reduced (sometimes they are just not significant anymore) when I use the log rather than the simple first difference to **compute the covariance.** Is there something incorrect in using **only** the first difference to compute the covariances rather than a log(first difference)? What am I missing here?


title: More specific test to see where it is significant rather then whether it is significant at all? Degree courses against levels of motivation to study

text: I've been doing chi-square tests to see whether there is a correlation between e.g. motivation levels when studying and degree course. With chi-square, it only reveals whether there is a significant difference between degree course and motivation, but not _where_ i.e. it looks like architects are significantly more motivated than politics students, but is there a significant difference between language students and psychology students? Or between psychology and politics students etc. ![Percentages of motivation](http://i.stack.imgur.com/DL7lQ.png) ![Stats of motivation](http://i.stack.imgur.com/rfsql.png) Which test would this require? Would it be run for every combination of degree course pair? Thank you. **EDIT:** Thanks I've been having trouble with ANOVAs though - the degree is nominal and I've put the motivation as scale, but it won't run an ANOVA as only one variable shows up. Could you expand on what GLM is? In the post hoc section I couldn't find it. ![enter image description here](http://i.stack.imgur.com/OlXM3.png)


title: How to know whether a vector is a row or a column?

text: When the text simply says that $X_i$ is a J-dimensional random variable, how do I know that it's supposed to have $1 \times j$ or $j \times 1$ dimension? An example where this matters is Cramer-Wald device, taken from my professor's notes: > Let $\{X_i\}_{i=1}^n$ be a sample of J-dimensional random variables and let > $X$ be a J-dimensional random variable. Then, $X_n \xrightarrow d X \iff > \lambda'X_n \xrightarrow{d} \lambda'X, \forall\lambda \in R^J $ I am taught that the $\Leftarrow$ direction of this device is used to prove multivariate convergence based on univariate convergence. Thus, I am able to infer that $\lambda'X_n \in R$, i.e. $(1 \times j)(j \times 1)$. But is there a way to tell the dimension of $X_n$ based solely on notation, without context?


title: Overfitting a logistic regression model

text: Is it possible to overfit a logistic regression model? I saw a video saying that if my area under the ROC curve is higher than 95%, then its very likely to be over fitted, but is it possible to overfit a logistic regression model?


title: Training set as donor for test set in binary classification problem

text: I am wondering if there exists such a method in machine learning that: Given a binary classification problem, for each person in the test set the person that most closely resembles this person in the training set acts as a donor for the value; either 0 or 1. Could someone show point me in a direction? **Note:** For instance, you have the variables age, sex and income. Sex has less spread, so we might standardize this variable to make it comparable with standardized age and income. It might then be calculated which person most closely relates to each test person.


title: Hidden Markov Model to predict the next state

text: I am learning to use HMM and I am trying to solve the following problem. There is a robot moving around the nodes in graph. The robot can move to adjacent nodes with certain probabilities. Each time the robot steps into a new "node", a (noisy) information about the node is generated. That is, I do not know the exact node. I have the following data:   * Each node node is a hidden state (finite number)    * A transition matrix that defines the probabilies for the transitions between nodes. A   * Emission probabilities for each hidden node Using HMM standard functions I should be able to predict the "hidden state" at time t if I have "t observations" (observations from 1..t) P(X_t|O_1::t). Is there a way to predict the next move (hidden state at t+1)? If I have "t-observations" (all observations at time t), is it possible to predict the most probable "hidden state" at time "t+1"? Which HMM principle should I use?


title: consensus scoring model that maximizes ranking agreement with real data

text: Consider the following dataset: OBJECT: 1 2 3 4 ... 15 SFA: 3.5 2.2 4.0 3.6 ... 2.8 rankingA: 3 9 1 2 ... 7 SFB: 3.1 3.1 3.7 3.4 ... 2.7 rankingB: 4.5 4.5 2 3 ... 9 SFC: 3.3 3.0 3.7 3.5 ... 2.1 rankingC: 4 7 2 3 ... 11 real_ranking: 2 8 1 3 ... 9 I have used 3 scoring functions (SF[A-C]) to score 15 objects. From these scores I have estimated their rankings. I also have their real ranking which I can use to compare with the predicted rankings (ranking[A-C]) using a statistic like Kendall's W. What I seek is that linear combination of the 3 scoring functions that yields the highest Kendall's W. In other words I seek the a,b,c,d values in the following equation: consensus_score = a*SFA + b*SFB + c*SFC + d The consensus_score should yield a ranking that agrees the most with the real ranking. What statistical modeling method should I use to achieve that?


title: Getting the prediction standard error from a natural spline fit

text: I'm fitting a natural spline fit to some data points. I'd like to estimate the prediction error for the predicted value. In linear regression (I agree that natural spline is also a linear regression with a specific type of design matrix), we know: $\hat{\beta} = (X^T X)^{-1}X^TY \rightarrow \text{ assuming var(Y) = } \sigma^2 I \text{ then : }var(\hat{\beta}) = (X^TX)^{-1} \sigma^2 $ Now consider $\hat{Y} = {X_i}^T \hat{\beta} + \epsilon_i$. We can then write: $Var(\hat{Y}) = {X_i}^T ((X^TX)^{-1} \sigma^2) (X_i) + \sigma^2$ This is easy to calculate for linear regression. How should I do it with natural spline? I can get the design matrix for natural spline. I can get $(X^TX)^{-1} \sigma^2$ but how can I get the rest of it: Here is an example in R:               set.seed(12345)     x <- c(1:100)     y <- sin(pi*x/50)     epsilon <- rnorm(100, 0, 3)     knots <- c(10, 20, 30, 40, 50, 60, 70, 80, 90)     myFit <- lm(y ~ ns(x, knots = knots))      Now consider x = 32.5 . How can I get the variance for the $\hat{Y}$ corresponding to x = 32.5 ? I know we can use the predict function. however, what I do really want is to get calculate it similar to linear regression by getting the design matrix and multiplying them together. I really appreciate your help.


title: Feature Selection - Mutual Information with response variable that takes three values

text: I am trying to calculate Mutual Information scores for Feature Selection. I have successfully implemented the Mutual Information to test each feature against the binary response variable. Each feature in my case is an n-gram with values 1(appears), 0(does not appear), and the binary response class takes values 0 or 1. However, I now want to apply the same test, but this time the response variable takes three values(0,1,2). I have searched for a similar example but I haven't found anything online. Is this possible to perform? If yes, is this a simple expansion to the existing formula, or does this introduce a more complicated problem?


title: Bayesian inferencing: how iterative parameter updates work?

text: I have been struggling with this for a while. A typical optimisation problem can be viewed as optimising some cost function which is a combination of a data term and a penalty term which encourages certain solutions. And normally there is a weighting term between the two. In the bayesian setting, this can be interpreted with the usual prior and the likelihood function. In the current problem that I am trying to understand, I model the prior as a multivariate normal with zero mean and the precision matrix equal to $\lambda \Lambda$ where $\lambda$ can be thought of this regularisation weighting and $\Lambda$ is some appropriate precision matrix structure that encodes the plausible solutions somehow. In my particular example, the precision matrix encodes some smoothness constraints on the estimated parameters i.e. the prior encourages smooth solutions. In this case, $\lambda$ denotes the strength of this smoothness penalty term. A $\lambda$ of zero would mean the ML estimate where we only optimise the cost function i.e. the likelihood function. This is because as $\lambda$ decreases, the precision decreases and hence the variance of each of the parameter in the prior increases. So, low values of $\lambda$ will move towards the unregularized solution. Now, a typical thing I have seen is that there is some sort of an iterative scheme, where we first start with an approximation to $\lambda$ and compute the distribution over the other parameter of interest using some approximate scheme like variational Bayes or Expectation Propagation and then use this approximation to update our estimate of $\lambda$ (assuming priors over $\lambda$ are of the conjugate form, usually done with a Gamma distribution which also keeps it positive). Now, my question is that if I start with a very low value for $\lambda$ as my approximation, then the prior term would hardly have any effect. Would this not push the estimated distribution towards solutions that are less plausible i.e. basically give high probabilities to unregularized solutions? I am having a lot of trouble understanding how this update scheme can actually find good values for $\lambda$ i.e. finding the value of $\lambda$ that is optimal with respect to the observed data. So, basically what I have trouble understanding it is what is stopping the inference to drive this value of $\lambda$ down to zero or close to zero to prefer the unregularized maximum likelihood estimate? I really do not see how this value of $\lambda$ is being driven by the data or the evidence term.


title: Decision Trees on training data

text: Wouldn't any decision tree trained on a training data set have no errors in classification? In other words, wouldn't every data point be classified correctly in the training data set? How would this tie in with the misclassification rate?


title: How to interpret coefficients of $x$ and $x^2$ in same regression

text: If I have the below functional form for an OLS regression, how do I interpret the $x$ and $x^2$? I cannot interpret them separately, correct? Do I interpret them as a summation of the two coefficients, so the effect of a one unit change of $x$ on $y$ is $0.5 + 0.3 = 0.8$? $$ y = a + 0.5x + 0.3x^2 + e $$ * * * _Update:_ Now with 2 squared terms, and one that is an interaction term with the other and a binary variable, is there no rule for how to interpret the coefficients for all values of $x$? Or does it depend on the value of $x$?   \begin{align} y = 0.04\; &+\; 0.05\times {\rm drop}\; +\; 0.008\times {\rm season}\; -\; 0.0004\times {\rm season}^2 \\ &+\; 0.002\times {\rm dropseason}\; -\; 0.0003\times {\rm dropseason}^2\; -\; 0.029\times p \\ &-\; 0.01\times r_1\; +\; 0.003\times r_2\; +\; 0.01\times r_3\; +\; 0.003\times v_5\; +\; e \end{align}


title: Neutral option on Likert Scale

text: I would like to know how do we analyze or regard the Neutral response (such as Don't Know) in a questionnaire using the following Likert Scale: 1:Strongly Disagree 2:Disagree 3:Don't Know 4:Agree 5: Strongly Agreee


title: spec.pgram in R: how to check for seasonality with longer phase?

text: I am using spec.pgram(tsobj,spans=6, plot=TRUE) to obtain a periodogram for my univariate time series of monthly observations which were sampled over 86 years (so I have 1032 observations in total). It's defined as an object of class "ts" with frequency=12. The periodogram nicely depicts the seasonality for 1 year and has further peaks on following years. However, it stops at 6 years. I have reason to believe that there might be a seasonality with longer phase, say 8 years. How can I change the command to obtain values for longer seasonalities? I am sorry, it sounds like a stupid question and possibly it is, but starring further on the help files, playing around with the available options and googling around didn't brought me any inside. Thanks alot for your help! * * * Here is the plot that I obtain: ![graph](http://i62.tinypic.com/7286kz.png)


title: How to enter/analyze likert data

text: My research involved a questionnaire with open-ended questions, multi-answer multiple choice (select all that apply; entered as "1" for circled answers, "0" for non-circled), and likert scale questions. I want to examine whether there are differences in responses to the survey (especially the likert questions) based on a variety of independent variables. Independent variables are things such as: age, ward number, intervention group (3 colour-coded groups - this was a randomized clinical trial), treating physician/dietitian (entered into SPSS as string variables), presence/absence of dysphagia (entered as binary data 1,2)...and many others. I would like to treat the likert data as ordinal data rather than interval. This is my first time conducting a research study and would appreciate pointers as to how to even begin the analysis. My limited understanding is that non-parametric tests are more suitable for my scenario. Thanks!


title: Dickey-Fuller for 2D signals?

text: I would like to run a stationarity assessment test on a 2D signal. Is there a suitable test, e.g. a 2D DF or ADF, etc? Thanks. Edit: A 2D signal is a signal defined on a lattice, i.e. instead of a time- series, `x[n]`, we have `x[m,n]`, where `m,n \in {0,1,...,N}`.


title: Standard Deviation After Subtracting One Mean From Another

text: I've got two sets of data from some fluorescent cells. The first set is when the cells don't have their fluorescence switched on, but they are still faintly glowing. The second set is when they do have their fluorescence switched on. I ran each of the two experiments three times, so n=3, and I therefore have a mean for each dataset and a SD. I need to subtract the non-switched on set from the switched on set in order to determine the amount of fluorescence that arises as a result of being switched on. How do I then calculate the SD of the final value? I know that for discrete random variables, $$ E[X+Y] = E[X] + E[Y] $$ So I assume that holds true for subtraction as well, but I can't find the rules for continuous random variables, as these are. Example data for one data point:                        ON            OFF     MEAN:  33956.6666    3835.66667     SD:    457.47301     38.0905


title: Which Statistical Technique (ANCOVA or Taguchi methods) for this problem?

text: I have one qualitative factor (Treatments) and two uncontrollable factors, temperature and moisture, with one response variable. I'm confused whether I should use ANCOVA with temperature and moisture as covariates or Taguchi methods can be used for this problem. Any advise and point to reference will be highly appreciated. Thanks


title: Applying the distribution of a large sample to estimate probability for a smaller sample

text: I have the sample size, mean and standard deviation for the number of births per year for my county, but I'd like to transfer that information to a smaller population within the county (ie. a school within the county). So in this situation, in the age group I'm concerned with, the population of women in the county is 1,000 and there are 83 births per year with a standard deviation of 5.6. How can I apply that distribution to a subsegment within the county with a sample size of 20?


title: The rule of the 3 sigmas ---- how many times multiply the sigma to get 85%?

text: $\sigma c=\sqrt{\overline{R}(\sigma S)^{2}+\overline{S^{2}}(\sigma R^{2})}$ where, $\sigma c$ --- the buffer stock --- needed for maintaining the risk of stockouts, is used by logistitians $\overline{R}$ --- the mean of the reserves $\overline{S^{2}}$ --- the mean of the production $(\sigma S)^{2}$ --- the square of the std. deviation of the production series $(\sigma R^{2})$ --- the square of the std. deviation of the reserves With this value of the stock buffer i am able to assure 84% of the stockouts. If i wanted 97.6% of the cases i would calculate for the value of 2 sigmas instead of one like in the equation. Now these percentages are given values......but what if i wanted to have the stock buffer to cover a custom percentage not only the cases where sigma is 1, 2 or 3....but 1.2 for example. The problem is i need to get how many times i have to multiply the sigma in order to get 85% of the stockouts covered by the possible stockouts.


title: Methods / approach to improve the predictive accuracy of a logistic regression model

text: **The situation:** I have a logistic model that should predict a defect (`1=defect, 0=no defect`). My model uses 4 out of 14 parameters, which are significant for my dependent variable (tested through `summary()` and the `anova()` chi-squared test). Furthermore, I used 80% (~5000, with ~1200 defects) of my data to train the model and 20% (~1200, ~300 defects) to test it. **My results are:** _If the cutoff I use is 50%_ *   * Sensitivity = 23% (percent of correctly predicted defects/total number of actual defects)   * Specificity = 98% (percent of correctly predicted non defects/total number of non defects)   * Accuracy = 80% (ratio of correctly predicted units/total number of units) _If the cutoff I use is the percent of defects in the training set [i.e., training defects / all training rows (~1200/~5000 = ~24%)]_   * Sensitivity = 55%    * Specificity = 81%   * Accuracy = 75% **My question is:** How can I train my model to get a better result for my sensitivity (defect rate)? (Or is there any fault in my approach, i.e. how can i check my model for scientific correctness??) I'm somehow lost at this point and appreciate any help, reference to a book or link that guides me in the right direction.


title: Joint pdf of functions of order statistics

text: Let $ Y_1 < Y_2 <\ldots <Y_{10}$ be the order statistics of a random sample from a continuous type distribution with cdf $F(x)$. How would I begin to show that the joint distribution of $V_1=F(Y_4)-F(Y_2)$ and $V_2=F(Y_{10})-F(Y_6)$ is given by: $$h(v_1,v_2)= \frac{10!}{1!3!4!} v_1 v_2^3 (1-v_1-v_2)^4 $$ Any hints? Thank you.


title: Compare variance and point statistics in a 3 way group

text: I need some help on this analysis. I have 3 groups of individuals, different the level and kind of physical activity. For each individual I have the heart rate, taken every minute for ~24h. (each person has a **slightly different number of measurement** due to differences in observation time) I need to compare the heart rate between the groups and the in-person variability among the groups. Any suggestions? The simplest approach I thought of was to generate one point statistic (mean or median) for each individual and than make a simple anova among the three groups. But I don't like this solution, mostly because I would loose all the information about the in-person variability, which is something I want to analyze. Some people suggested me a factorial repeated anova, but there is the problem that the number of persons per group and the number of in-person measurements are different. Another approach that has been suggested is mixed effects model (I should learn how to do it). Any suggestion for a good approach? (NB. I'm not an expert in statistics)


title: What are appropriate statistical test to compare historic Census variables for different locations?

text: I am wondering if anyone has suggestions on which statistical test would be appropriate to compare demographic data obtained from 1900 Census schedule sheets? I have a database of boarding school student information (age, gender, race, birthplace,...) for 25 different boarding schools across the country. I want to compare the Census data from each school to attendance/enrollment criteria to see whether or not each school enrolled students of the proper age, gender, and race. Any suggestions and comments would be highly appreciated. Thanks everyone.


title: How do we find concordance between two sets of non-linear, non-monotonic data

text: Is there a statistical procedure that can measure the concordance between two sets of data that have, say, a sinusoidal tendency?


title: Does the concordance correlation coefficient make linearity or monotone assumptions?

text: Does Lin's concordance correlation coefficient assume that the 2 data sets have linear or monotonic tendencies? Or, can I measure the concordance between 2 data sets that have a sinusoidal tendency? Thanks NG