title: Sample size requirements when studying total population

text: I work for a fundraising organization and we want to implement an A/B test to determine if a particular donor recovery technique is working. The population under study is our lapsed donors. These are donors that gave to us in the past but did not give during our last campaign. For the study we plan to solicit _all of our lapsed donors_. All lapsed donors will get the same solicitation package but the A group (test group) will also get some sort of handwritten personalized note. The B group will get no personalized note. Success of the technique will be determined by comparing the number (proportion) of 'recovered' donors in each group. My challenge is to determine how many personalized notes we have to write. For argument sake, I'll say we have 8000 lapsed donors so it is not realistic to simply divide the group in two and write 4000 personalized notes (we don't have the resources). What tests do I need to run to determine the size of the test group? I would have assumed power analysis but does it matter that I am not doing inferential stats here? We will be studying all lapsed donors. I've done some reading and I understand that there is a lot of nuance here. For example, one could argue that, in a way, I will be doing an inferential test because I want to know that the test will inform our approach to, not only the population under study, but all _future_ lapsed donors as well (the hypothetical super- population?). But even if we set the inferential bit aside, I would still want a robust enough number in the test group to be confident the post test results are reliable. Please advise or suggest resources. Also note that I am fairly proficient in R so feel free to add R code. I've also explored the 'pwr' package a bit. Many thanks in advance.


title: Double Seasonal Holt winter method using dshw

text: I have a two weeks data set which have intraday and intraweek cycles so I decided to use dshw in r. Although it gave me a pretty good MAE and RMSE, when I wanted to see SSE, it showed me a null value. Does it mean that Im using a wrong model for my dataset?I only used auto.arima command to see the best fitted model but the recommended model did not provide any seasonality values while the time series plot clearly shows seasonality in hours of a day and days of a week. Hence I switched to dshw. On the other hand, what is the commend to find predictions for the next cycles? I tried command "predict" but it did not work. Thanks


title: Difference between log-normal distribution and logging variables, fitting normal

text: Context: I have a set of data that is bimodal, so I used the mixtools package in R to fit a bimodal normal distribution to it. It looked as if the normal did not fit very well, and given other similar sets of data that I have (that are not bimodal) and are log-normally distributed, I figured log-normal would probably make more sense. However, mixtools does not have a way to fit a bimodal log-normal distribution, so I took the log of the data and refit. This image is below, and it fits to my satisfaction. Really, this is kind of a software issue (not knowing what package could explicitly fit a bimodal log-normal) but it's made me think about something I have wondered often: is fitting a normal distribution to logged data equivalent to fitting a log-normal distribution to the original data? I suspect not, but I am not sure why? (Ignoring cases of having zeros or negative data. Assuming all data are fairly large positives, but skewed). Also knowing that of course you'd need to back transform to get to the original value estimates. I tried to test this out with some toy data and realized I don't even know why the meanlog associated with a log-normal distribution is NOT what you get when you take the mean of the logged normal distribution. So maybe my understanding is already broken when it comes to the log-normal's parameters.               library(fitdistrplus)          set.seed(1)     test <- rnorm(1000, mean=100)     test[test<=0] <- NA #Unnecessary since no values <= 0, but just to prove     test<-na.omit(test)          log.test <- log10(test)     mean(log.test)     sd(log.test)     #1.999926 is mean for log.test     #0.004496153 is sd for log.test          fitdist(log.test, dist="lnorm", method="mle")     #However, "meanlog" is 0.693107737 and "sdlog" is 0.002247176     #The means are so different, not sure why?      ![enter image description here](http://i.stack.imgur.com/Y14OB.png)


title: What are best practices in identifying interaction effects?

text: Other than literally testing each possible combination of variable(s) in a model (`x1:x2` or `x1*x2 ... xn-1 * xn`). How do you identify if an interaction SHOULD or COULD exist between your independent (hopefully) variables? What are best practices in attempting to identify interactions? Is there a graphical technique that you could or do use?


title: Method to reliably determine abnormal statistical values

text: I'm searching for a statistical method to determine if a player is cheating in an online game. The game is a Quake3 like game (ego-shooter). Given a number of positive points and a number of negative points per player (score) and given n players (n<=64). The score comes together like this (positive/negative seen from anti-cheat- perspective):   positive = number of time the player died himself   negative = number of opposite team's players killed         Additional available values:   t = Time in which this score has been accomplished   c = number of current number of players   m = number of times a player killed someone from the same team (teamkiller, very negative if higher than 1 or 2 over a 15 minutes period)      Each player can cheat independent of each other. I used the standard deviation (on the value v = (negative+1)/(postive+1) per player, assigning the exceeding of 3 times standard deviation with absolute certainty (100%) of cheating - although my experience has been that exceeding stddev alone would be sufficient in most cases), which works great on a single player but fails miserably as soon as several players are cheating at the same time. Also the practise of defining 3 times stddev as 100% allows for > 100% probabilities. Unfortunately, the percentage of players cheating in an unfiltered game is roughly 70 to 80 percent or even higher, if accounting for more subtle cheats. One more thing to consider is, if there is a group of people cheating, and one that isn't, the distribution is no longer a normal distribution, but a camel- bump like distribution. Is there any useful algorithm or formula for this problem ?


title: Aggregation of association rules

text: ## Problem: A customer "Cu" is purchasing two items : (A,B). I'd like to know which item should I recommand to him between items C and D. The natural idea I have is to analyze the basket of previous customers and see if (A, B, C) is more frequent than (A, B, D)               Example :     (A,B,C) : 4 customers     (A,B,D) : 6 customers      It's seems better to recommand D to customers who are purchasing (A,B) The problem is I don't have enough users (A,B,C) or (A,B,D) to extrapolate. I think I should also analyze behaviours of customers who purchase (A,C), (A,D) and (B,C), (B,D) because there is a common item with the customer "Cu".               Example:     (A,C) : 60 customers     (A,D) : 40 customers     (B,C) : 90 customers     (B,D) : 10 customers      Based on these users, it's seems better to recommand C, but the relation between them and "Cu" is weaker (only one item in common) ## Question : Does it exist a standard way to aggregate association rules and give more weight to stronger association (with more common items) ?


title: Is Monte Carlo simulation more appropriate than parametric tests for constructing confidence intervals for weighted means?

text: A colleague suggested that Monte Carlo simulation should be preferred for constructing confidence intervals for weighted means calculated from a sample. How and why exactly might Monte Carlo perform better than the traditional t-test in this use-case?


title: Correlation between sensors

text: Background: A home wired with multiple sensors, measuring attributes like temperature, light, motion etc. In addition, a multitude of actuators that can carry out an action like opening a door, switch on a light etc. Problem statement: If a user were to use this home like a normal home, I'd like to be able to find recurrent patterns in his interaction with his home. Leading to an automation of these activities. In such a situation where the deployment and usage of the equipment varies but the premise is people tend to repeat things they do at home pretty reliably. What sort of a learning algorithm do you see using in this scenario? Essentially input data is going to be a collection of time series streams. Some are discrete "events" like light switch on, but others are numeric values like temperature.


title: Compare Means of Groups with Unequal Sample Size

text: I have to compare two groups of data and see whether there is a significant difference between the means. I wanted to perform a Paired-Samples T-Test with SPSS, however, it will only focus on a part of the data. ![Overall Means and N](http://i.stack.imgur.com/G3gkS.png) ![The new N after performing the test](http://i.stack.imgur.com/6Dgbn.png) As you can see, the N drops from 10148 to 1166 in the Paired Samples Statistics. Is it possible to compare means for the WHOLE group?! I have not the best knowledge of statistics, but I really need to compare the two groups... Should I choose another test? Am I using the wrong one? I really hope someone can help. Thanks in advance! Best, David I will try to give some more background. I have two columns in excel with stock reactions on a recommendation. Let us say that in the first column we have 1166 recommendation reactions from 2005 and in the second column we have 10148 recommendation reactions from 2006. The means of the two groups are totally different (0.79 and 0.36). However, I want to compare these means in a statistical way. I hope this helps..


title: Tobit with proportion as dependent variable

text: I want to run a censored regression with a dependent variable that is a percentage (a bunch of values are 0 because of the way the data was measured, and the rest between 0 and 1). Tobit regressions don't work with such variables because the log of a variable smaller than one is negative, right? What should I do?


title: Can pivot be used for testing

text: A pivotal quantity $Q(X, \theta)$ can be used to construct a confidence interval. I was wondering if it can be used to construct a test statistic and rejection region? In simpler cases involving a simple null hypothesis it is obvious how to use a pivot to construct a test statistic and rejection region. What about composite null hypotheses? Thanks and regards! > **Pivotal quantities are fundamental to the construction of test > statistics** , as they allow the statistic to not depend on parameters – for > example, Student's t-statistic is for a normal distribution with unknown > variance (and mean).


title: Compute sum of vectors drawn from multivariate normal, subject to a linear constraint

text: I want to compute $S = \sum_{i=1}^n x_i$ where $w^t x_i>-1, \; \forall i$ and $x_i \tilde{} \mathcal{N}(\mu, \Sigma)$ for known $w$, $\mu$ and $\Sigma$. I know $S$ can be approximated by sampling techniques but I was wondering if S can be computed analytically.


title: Calculate a 95% confidence interval and p-value for the change in C-statistic using bootstrap with R

text: I am using `boot()` and `boot.ci()` to furnish confidence intervals for the difference in the $c$-statistic (AUC) between models with and without a novel biomarker. My method is below. I am not sure if this is the correct way to go about it. Here is example data (n = 500, 3 columns: disease status (1/0), covariate (1/0) and biomarker (1-6)):               data <- read.table("https://notendur.hi.is/eyb8/dataframe.txt", header=TRUE)      Here is a function that returns the difference in $c$-statistic between the models, for use in `boot()`:               library(rms)     Cindexdiff <- function(data, indices){       data <- data[indices,] # select obs. in bootstrap sample         # C-statistic without the biomarker:       C1   <- lrm(disease ~ covariate, data=data, x=TRUE, y=TRUE)$stats["C"]          # C-statistic with the biomarker:       C2   <- lrm(disease ~ biomarker + covariate, data=data, x=TRUE, y=TRUE)$stats["C"]        as.numeric(C2-C1) # returns the difference     }      And now the 95% CI for the difference in $c$-statistics with bootstrapping (repeated 999 times):               library(boot)     set.seed(1)     b <- boot(data, Cindexdiff, 999)     boot.ci(b)      Now, I have **a few questions:**   1. Is this a sound method for calculating bootstrap confidence interval for the change in C-statistic?   2. How should I choose which type of confidence interval from `boot.ci()` to use? "Normal" or bias-corrected (BCa)?


title: Test for comparing "mean of means"

text: I have a treatment and a control group. Both groups consist of several teams, each with several individuals. When I want to compare the performance between the groups, I could run a t-test or Wilcoxon ranksum test on the individual values within both groups, and that's fine. 1) However, because those teams are very different in size, the "mean of means" looks more like what I get in parametric regressions later. With "mean of mean" I refer to first taking the mean of individual values within the team, then take the mean of team means within the group. Denote by $t$ team, $i$ individual within the group, then $$\text{mean of means}=1/T\sum_{t=1}^{T} (1/n_t \sum_{i=1}^{n_{t}} x_{it}).$$ The question is: are there tests that allow me to test the hypothesis that the two "mean of means" are equal? 2) If not, how would you compute the standard deviation of the mean of means? Clearly, just using the SD of the team means wastes information. In wikipedia I came across pooled standard deviations, which seem to be suitable in this context. Would it be valid to construct standard errors based on these pooled SDs and use them for the hypothesis test? Thank you very much!


title: Log-linear analysis vs Breslow-day test

text: I would like to test independence of several categorical variables for a few datasets. I believe Breslow-day is possible for some of the analyses I want to do, and log-linear may be possible for all, but besides applicability conditions, how are these two analyses different and which is more appropriate for the cases where both are possible? Experimental participants' responses to a comprehension question were scored as strong or weak in each of two respects: correctness and generality. Participants came from one of two training conditions and were also classified into one of two previous knowledge groups. So I have altogether 4 binary factors: respect, score, training, and group. Score and group were observed, while respect and training were manipulated. Respect is within-subjects, while the others are all between-subjects. In a later experiment I add an additional binary experimental factor, pretraining. I initially didn't take account for group, so I did a Breslow-day test with respect, score, and training. I found a significant effect, which I interpret as "the effect of training on score differs depending on respect". In particular, one training condition promoted correctness better, while the other promoted abstractness better. I later wanted to see whether adding group made any difference. As far as I know, I cannot add another factor in Breslow-day. Thus, I tried a log-linear analysis. I first did this with the same 3 factors used in the Breslow-day test, so I could compare the results. Here's R code:               T = table( D$score, D$training, D$respect, dnn=c('Score','Training','Respect') )     loglm(~ Score + Training + Respect, data=T)      This result is NOT significant, which I interpret as meaning the data can be explained completely by main effects of the factors, without any interaction term. That seems to contradict the result of the Breslow-day test, so I'm now not sure which one is more accurate. Next, I'd like to add the other factors (group, pretraining). Am I right that I MUST use the log-linear analysis in this case?


title: Best transformation for sinuous data sets?

text: I am analysing annual behavioural patterns and am currently scratching my head over how to best transform data to test for correlation where my independent variable is the moon phase (given as perecentual moon illumination, i.e. 0% new moon, 100% full moon). May be a rather trivial question, but I haven't come across this problem before so am a bit at a loss. Cheers.


title: Identify points significantly outside a regression model

text: I would like to update my previous question "Methods to determine reliability of measurements using median and median absolute deviation". I have data from biological experiments, that look like this:               v1 2 1.8 1.5 1.9 2.1 1.78 1.95 2.0 2.1       v2 2 100 -5.2       v3 1 -1.3 -2 2.3       v4 1 1.5 1.6 1.9 2.1 2.0 2.4 -1.1 2.3 1.5 1.6 1.9 1.8 1.6      These represent gene expressions. Now, I would expect that all values of each variable(genes) are more or less similar, since the values are ripetute measurement of the same gene. Having a variable with such huge difference, as v2, doesn't have sense. That means that it does not make sense to have a gene with such a huge variation in expression. Therefore, it has to come from a methodological error. I was looking for a method (possible a statistical test) which could identify the "average variability" among my samples and report me which variables(genes) have a variability significantly greater. This means that for these gens my data are not good enough to estimate the expression, and I have to discard them. ** _In the first place I would really appreciate any suggestion on test I could use for my purpose._ ** I thought at something like that: I consider for each gene the ratio between **median absolute deviation(mad)** and **median(med)**. I would expect this ratio has a normal-like distribution, assuming that the variability comes from error in the measurement occurred by chance. Therefore, I could identify then all variables(gene) that have a **mad/med** ration significantly higher that the average distribution, or if you want that are "positive outliers". For me it sounds quite intuitive but I'm not sure how to perform that in an appropriate way. I would really appreciate your opinions about my method and also your suggestions on how I could perform the analysis.


title: Knapsack problem with multiple groups where items belong to more than one group

text: As follow is a Knapsack Problem where there are multiple groups and each item belongs to one group. The goal is to maximize the profits subject to the constraints. In this case, only one item from each group can be chosen with the exception of the special group. Here is the problem :   Maximize: $$\sum _{j=1}^{N}\sum _{j=1}^{\left | G_{i} \right |}p_{ij}x_{ij} +\sum _{j=1}^{\left | G_{s} \right |}p_{sj}x_{sj}$$ Subject to: $$ \begin{align} (\sum _{i=1}^{n}\sum _{j=1}^{\left | G_{i} \right |}w_{ij}x_{ij} +\sum _{j=1}^{\left | G_{s} \right |}w_{sj}x_{sj})&\leq c \\\ \sum _{j=1}^{\left | G_{i} \right |}x_{ij} &= 1 \\\ \sum _{j=1}^{\left | G_{s} \right |}x_{sj} &= k \\\ x_{ij} \in \\{0,1\\}\\\ x_{sj} \in \\{0,1\\} \end{align} $$ Notation:   _N_ the number of items   _n_ the number of groups   _c_ the capacity of a single knapsack   _G_ = {$G_{1},...,G_{n}$} the set of groups   ${\left | G_{i} \right |}$ the number of items in group $G_{i}$   ${\left | G_{s} \right |}$ the number of items in group $G_{s}$ (special group)   $p_{ij}$ the profit of item _j_ of group $G_{i}$   $w_{ij}$ the weight of item _j_ of group $G_{i}$   $p_{sj}$ the profit of item _j_ of group $G_{s}$   $w_{sj}$ the weight of item _j_ of group $G_{s}$   _k_ restraint of group $G_{s}$   This is fairly easy to solve. However, I am having trouble describing a situation where the items may be members of more than one group but can only be selected once. For example, lets say you have a fantasy sports roster. You want to optimize your selections based upon your estimated point contributions. The weights are the player's salaries. The previously provided formula will work where player to position (item to group) is in a one to one relationship. But lets say a player can fill two or more positions. The question becomes: how can this be written as a constraint in the above? Is this an entirely different problem altogether?


title: Issues with estimating the sparse inverse covariance matrix with Glasso

text: I am trying to estimate the sparse inverse covariance matrix of my gaussian graphical model. I installed the `glasso` package in R and tried out some examples. After that I ran the `glasso` software on my own data. So I fed it my empirical covariance matrix. However, it seems to get stuck and doesn't give me the results. So I was wondering if I was running the software incorrectly. I followed this manual: http://cran.r-project.org/web/packages/glasso/glasso.pdf So if `S` is my empirical covariance matrix. I just ran the following:               a <- glasso(S, rho=.01, trace=TRUE)     # outer loop, m = 1     # outer loop, m = 2      After the above, it gets stuck. I am attaching my empirical covariance matrix as well which I feed to the software. http://rapidshare.com/files/141198734/emp_covariance_matrix.mat Any guidance will be much appreciated.


title: Combinatorics problem applied to GAM logistic regression

text: I'm doing biomedical research and I need to set a GAM Logistic Model which get the maximum AUC score as possible. I have 4 disease markers; $Y_1, Y_2, Y_3, Y_4$ with different data in each one, and the model must have the following specification: $\text{Model} = \log(\text{Odds Ratio} (\text{Illness}~|~Y_j,Y_i)) = \alpha + \sum\limits_{i,j}^n f_k(Y_k) $ The problem is that I need to code in R all the possible combinations between this 4 markers without repetitions, in order to do multiple comparisons, such as: $\text{Model}_1 = \alpha + f_1(Y_1) + f_3(Y_3)$ or $\text{Model}_2 = \alpha + f_1(Y_1) + f_2(Y_2) + f_3(Y_3)$ etc. When I code each model manually I use sentences like these. This is an example for a model with only one marker:               model=gam(group~s(y1),family='binomial')     pred=predict(model,type='response')     r=ROCEmpiric(group,1-pred)      Where `group` is always in, as indicator of disease, and `ROCEmpiric` is a function which calculates the ROC value. Does anyone know the way to automate the calculations? It's seems ridiculous to write each function after doing the summation of permutations.


title: Cronbach's Alpha on true and false

text: I have a 10 item knowledge scale (true/false). Can I do Cronbach alpha test?


title: R - What does "incomplete block on file" mean?

text: I am using RStudio, and am trying to use packrat with my current project. I click the checkbox "Use packrat with this project" and press OK, where I get the following output from the console:                   > packrat::init()     Initializing packrat project in directory:     - "/Users/Ash/Dropbox/Uni/2014/Thesis/Code/R"          Adding these packages to packrat:                 _                  packrat   0.4.0.12          Fetching sources for packrat (0.4.0.12) ... Error in snapshotSources(project, activeRepos(project), allRecordsFlat) :        Errors occurred when fetching source files:     Error in untar2(tarfile, files, list, exdir, restore_times) :        incomplete block on file      I can't find much on this error, but based on this code (searching for the error string, incomplete block on file), it seems this relates to an unexpected file length/size. I assume the download may be corrupt? But I have tried multiple times. Or perhaps I don't have the correct permissions? But the project file should have no issues here. Has anybody else had this issue?


title: Appropriate analysis for this experimental design - non-normal data

text: I've run an experiment with 8 subjects, which had to perform [a series of 4] navigational tasks under 4 different vision conditions each. The experiment was repeated weekly, 16 times with the same subjects. The order of the conditions was randomised at each repetition. As a result, I have now various measures of navigation performance (dependent variables) such as speed, no. of errors, distance travelled etc., each for 4 tasks x 4 vision conditions x 16 repetitions x 8 subjects. For analysis, I'm currently only looking at one dependent variable and one task at a time. My main interest is whether overall performance varies between the 4 vision conditions, and if possible, how this relates to repetitions. (Post-hoc?) pairwise comparisons are a very central goal. Could you please help me in choosing an appropriate statistical test for this? Originally I had been using a two-way ANOVA, but I've looked at the data, and it is not always normally distributed. Mainly, one or two of the vision conditions deviate from normality. Also depending on the dependent variables, there are large deviations from normality. Can I justify using two-way ANOVA, and if not, what would be a valid non- parametric alternative? I starting trying Friedman test - just not sure if all the data for the subjects' repetitions can be appended and fed into the test to find differences between the conditions. I'm using SPSS and Matlab. Thank you!


title: How to compare fit of discrete process with discrete underlying process?

text: I am basically looking for an equivalent to something like an $\mathbb{R}^2$ for a model on a dataset that is itself simply a collection of points. That is, if my data set is (trivial example):               x <- 1:10          y <- c(1,6,3,7,7,10,4, 7, 8, 10)      and my 6-point fit grid is:               x.fit <- c(1,2,4,6,8,10)     y.fit <- c(1,4,6,6,7,10)      how might I compare this fit, with 6 coordinate pairs, with one of fewer or more? Basically I am looking to minimize the number of points I need to represent the underlying process with the fit, while ensuring the fit still models the data well. This 'well' is what I am looking to statistically analyze, so that I can plot my model with $N$ points from 1 to `length(process$x)`, and find a balance. I am uncertain how to process because I am used to fitting polynomials or predictor point-invariate models. Any help would be greatly appreciated.


title: Comparison of close data sets

text: I'm studying around **100 sets** of temperature ( **$N_{sample}=500$** ), which depends $4$ explicative variables such as power or speed. The dependency is always the same in each set, but sometimes the mean and the variance is different.   * I would like to **group similar sets together** and study each group separately: find a multiple regression model for each group.   * Then I would like to **classify future sets** I will discover in one of those groups. I don't know if it concerns clusterisation or pattern recognition or maybe something else. But I have no idea of how to do classifiy/compare my sets, except comparing the mean and the variance one by one. Does anybody have advices or suggestions? Here is an example of 3 sets: ![Temperature vs power](http://i.stack.imgur.com/WaOzk.jpg) ![Temperature vs speed](http://i.stack.imgur.com/mF5ub.jpg) Thanks for helping.


title: False alarm in Neyman-Pearson Test. A doubt in text?

text: Neyman-Pearson Lemma says that the most powerful test, $\phi(x)$, of size $\alpha$ (probability of false alarm), for testing $H_0:\theta=\theta_0$ versus $H_1:\theta=\theta_1$ is the likelihood ratio test of the form: $$ \phi(x)=\begin{cases} 1 & l(x)>k \\\ \color{#C00}{p} & l(x)=k \\\ 0 & l(x)<k \\\ \end{cases} $$ where, $p$ is the probability of accepting $H_1$ at $k.$ $k$ is threshold. $l(x)$ is the likelihood ratio, $l(x)=\frac{f_{\theta_1}(x)}{f_{\theta_0}(x)}.$ Probability of false alarm, $P_f=\color{#C00}{p} \cdot P[l(x)=k \hspace{2 mm}| \hspace{2 mm} H_0 ]\hspace{2 mm}+\hspace{2 mm} P[l(x)>k \hspace{2 mm}| \hspace{2 mm} H_0 ].$ Probability of detection, $P_d=p \cdot P[l(x)=k \hspace{2 mm}| \hspace{2 mm} H_1 ]\hspace{2 mm}+\hspace{2 mm} P[l(x)>k \hspace{2 mm}| \hspace{2 mm} H_1 ].$ **Question** Why Probability of false alarm, $P_f$, is $\mathbf{not}$ defined as: $P_f=\color{#C00}{(1-p)} \cdot P[l(x)=k \hspace{2 mm}| \hspace{2 mm} H_0 ]\hspace{2 mm}+\hspace{2 mm} P[l(x)>k \hspace{2 mm}| \hspace{2 mm} H_0 ].$ My understanding for saying above equation is that: if the probability of occurring $H_1$ is $p$, then the probability of occurring $H_0$ is $(1-p).$ And therefore, the probability of selecting $H_0$ at $k$ is $(1-p).P[l(x)=k \hspace{2 mm}| \hspace{2 mm} H_0].$ Where my understanding is wrong? Any comments please.


title: deeper lags in dynamic panel regressions using xtabond2 in STATA

text: I know this kind of question may be asked a number of times in this community or in other web sites. I am kind of knowing what to do. But just to make sure I did not misunderstand anythings, I pose my question. Currently I am using XTABOND2 to run difference and system GMM to model annual change of TFPt = TFPt-1 + some explanatory or predetermined variables + time dummies. The code is like that for difference GMM: xtabond2 dtfp ltfp x time2-time12, gmm(ltfp) iv(x time2-time12) nol robust for system GMM: xtabond2 dtfp ltfp x time2-time12, gmm(ltfp) iv(x time2-time12) robust dtfp = annual change of TFP at time t ltfp = TFP level at time t - 1 x = just some explanatory variables time = time dummies The number of panel is 123,800 and time is from 2000 to 2011. It is a extreme large unbalanced panel with average observation per group is around 4. My questions are here   1. Arellano-Bond test for AR(2) for both difference and system GMM are statistically significant, p-values are 0. That means lag 2 serial correlation exist and the current instrumental variable is not sufficient. I tried to use a deeper lag, simply change gmm(ltfp) to gmm(ltfp, lag(2 .)). When I did that, AR(2) is statistically significant but AR(3) and AR(4) are not, p-values between 0.6 and 0.7. From these results, can I say gmm(ltfp, lag(2 .)) is sufficient instruments?   2. Some people do suggest add lag dependent variable to this equation, so xtabond2 dtfp l.dtfp ltfp x time2-time12, gmm(ltfp) gmm(ldtfp) iv(x time2-time12) nol robust. I just wonder is it the right thing to do? If so, how do I know how many lags of dependent variables are sufficiently enough? still using arellano-bond correlation tests in AR2 AR3?   3. In XTABOND2, there is an option to do a difference transformation or orthogonal transformation. From David's how to use xtabond2, he seems to suggest using orthogonal transformation. One obvious benefit for that is to save observations at the first time period (I know there are few more benefits). From my observations, they return quite different coefficients - the coefficients of ltfp in difference and orthogonal transformation in difference-GMM are -0.63 and - 0.19. Are there anyway to know which transformation is doing a better job?   4. XTABOND2 offers difference-in-hansen tests for both GMM and IV instrments. Just to check I get it correctly. If both Hansen test excluding group and difference in GMM instruments for levels have very small p-values, it suggests system-GMM is invalid. If both tests in iv have small p-values, it suggests some of variables in IV are strictly exogenous and they should be treated endogenous or predetermined? There is my questions. Thanks in advance for your helps. Cheers, Fish


title: How to evaluate a trained model using parameters other than AUC in RapidMiner?

text: I am using RapidMiner to build predictive models trained and cross-validated by a set of medical data(65 cases. 18 attributes), I am now running trials by trying different combinations of learners and parameters, by this way I hope to explore the features of the operators and how they work, also see how they can fit with my dataset to create better models. I know that merely recording AUC for each situation is not enough because the issue of overfitting. So in the results, how can I see overfitting when it happens, or which numbers should I write down in the trial notes for further analysis of overfitting? I know my dataset is kind of small, I am doing a short-term research project, so that's all I got in hand. I really don't know if anything meaningful will come out when I even can't believe my own validation results. Any suggestions about the way to detect overfitting or any other steps I can take to research on this dataset is welcome. I am pretty new in this area, so anything could be helpful. Thank you.


title: Too big (?) histogram values when using normed histo options in SciPy and matplotlib

text: I have been trying to create a normed histogram using either SciPy or matplotlib (or anything for Python). When I create my histogram with 'normed' option disabled, it looks like below (this example is for 10 bins, but the same happens for a larger number of bins): (The first number represents the start of a bin, the second the bin's height)               -2.83785600931e-17   1182     5.6688145554e-15   1137     1.13660076709e-14   1031      1.70632007864e-14   950     2.27603939019e-14   912     2.84575870174e-14   802     3.41547801329e-14   853     3.98519732484e-14   948     4.55491663639e-14   1315     5.12463594794e-14   870      Which is absolutely fine, and what I was expecting. However, I later need to fit this histogram to another histogram, and for that I prefer to have a normed version of this histogram so that fitting the height of those histograms is easier. Strangely, when I use the option `density=True` (for `scipy.histogram` version) or `normed=True` (for `matplotlib.pyplot.plt` version) my histogram bin heights get very large values, like below:               -1.44880082614e-17   2.00318764844e+13     5.71138595513e-15   1.98921598219e+13     1.14372599185e-14   1.8040914044e+13     1.71631338819e-14   1.52465807942e+13     2.28890078453e-14   1.56133370332e+13     2.86148818087e-14   1.4617855813e+13     3.43407557721e-14   1.50020766348e+13     4.00666297355e-14   1.74296536456e+13     4.57925036989e-14   2.3769297206e+13     5.15183776622e-14   1.50020766348e+13      I hardly know anything about statistics, but I expected "normed" to mean "sums up to one". Am I incorrect in my thinking, or is the output normalized histogram wrong after all?


title: Choosing levels of random factor in a mixed model

text: My apologies in advance if this question is as much about experimental design as it is about statistics. I'd like to know if there is a "best" way to choose levels of a random factor in a mixed model. Let's say I'm interested in how people differ in performance on a task either blindfolded vs. not. Suppose also that I can vary the difficulty of the task on a scale of 1-100 arbitrary units. I am not interested in the effects of task difficulty - it is a "nuisance" factor that I'd like to model as a random effect. If I can choose only 5 levels for my random factor, should I aim for levels that are closer together (say, task difficulties of 30, 32, 34,... etc), or for levels that are as far apart as possible (e.g., 20, 40, 60 etc.), and why? How might this choice affect my statistical model? Thank you!


title: Regression by multiple dependent variables with constraints & feature selection

text: I have a data set of 1000 records. Each record has three dependent variables $y_1, y_2, y_3$ and 100 independent variables $x_1,...,x_{100}$, where the dependent variable $y_i$ satisfies:   1. $0\le y_i \le1$   2. $y_1 + y_2 + y_3 =1$ I.e. $y_i$ represent the probability of a observation belonging to one of the three classes. Q1: How can I build a (multivariate linear) model using $x_1,...,x_{100}$ to predict ($y_1,y_2,y_3$)? Is any R package available? How can I implement it? Q2: Since there are so many independent variables $x_1,...,x_{100}$ (features), is it possible to do feature selection using LASSO, SCAD, or elastic net for this multivariate linear model using tools in `glmnet` package?


title: In OLS is it methodologically correct to use the variance of a variable as an explanatory variable?

text: Are some OLS assumptions not satisfied if I use the variance of a variable as a proxy of uncertainty in a regression? For instance, would it be methodologically correct if I use moving averages of rainfall and moving averages of the variance of rainfall in a regression were the dependent variable is production at the farm level? I have not found examples of regressions containing the variance of a variable in the X's so I wonder if that is because there are some modeling issues with that approach. Thank you for your help.


title: plot: straight regression line and prediction bands

text: I have realized a multiple linear regression. To know how well my model does in terms of prediction, I can compute prediction intervals bands and decide if they are narrow enough to be of use. If they are too wide, then they probably are not useful. Using R, I have written these R codes.               LinearModel.1 <- lm(GDP.per.head ~ Competitivness.score + Quality.score,                          data=Dataset)     summary(LinearModel.1)     predict(LinearModel.1, se.fit = FALSE, scale = NULL, df = Inf,              interval = c("none", "confidence", "prediction"), level = 0.95,              type = c("response", "terms"), terms = NULL)      Here is an example code:                # generate random numbers for the domain      x = rnorm(10)      y = rnorm(10)           # make up some sample data to be fit.      z = 3*x+2*y - 5 + rnorm(10,sd=3.5)           # Perform the regression      fit = lm(z~x+y)           # Find some points to use to plot the confidence interval      # and then use the predict command to get the bounds on      # the confidence interval      x1 = (-20:20)/10      y1 = (-20:20)/10      p = predict(fit,data.frame(x=x1,y=y1),interval='confidence')           # Now plot the original points, then plot the predicted values and the      # upper and lower bounds of the confidence interval.      plot(x,z)      points(x1,p[,1],col=2,type='l')      points(x1,p[,2],col=3,type='l')      points(x1,p[,3],col=4,type='l')      I try to do the same with my LinearModel.1. I mean I would like to plot the original points, then plot the predicted values and the upper and lower bounds of the confidence interval. I don't get the R code to do that with my LinearModel.1. Thanks for your help


title: What is the optimal selling price of an item based on its historical price/volume data?

text: I'm trying to figure out how to price items on a market that appears to move randomly but within certain bounds. For example, here is a brief look at the time series data:               DateTime               Price       Volume     2013-08-27 00:00:00 1.9299999476    4     2013-08-27 01:00:00 1.6849999428    2     2013-08-27 02:00:00 1.7024297714    2     2013-08-27 03:00:00 1.2264721394    2     2013-08-27 04:00:00 1.75            7     2013-08-27 05:00:00 1.5             2     2013-08-27 06:00:00 1.9299999476    3     2013-08-27 07:00:00 1.75            5     2013-08-27 08:00:00 1.7000000477    5     2013-08-27 09:00:00 1.8600000143    4     2013-08-27 10:00:00 1.8236290216    2     2013-08-27 11:00:00 1.8600000143    7     2013-08-27 12:00:00 1.8046891689    6     2013-08-27 13:00:00 1.9107458592    2     2013-08-27 14:00:00 1.9414916039    6     2013-08-27 15:00:00 1.9605424404    6     2013-08-27 16:00:00 1.8037977219    6     2013-08-27 17:00:00 1.871355176     4     2013-08-27 18:00:00 1.727257967     5     2013-08-27 19:00:00 1.6200000048    5     2013-08-27 20:00:00 1.8400000334    9     2013-08-27 21:00:00 1.8099999428    1     2013-08-27 22:00:00 1.8099999428    10     2013-08-27 23:00:00 1.7849999666    2     2013-08-28 00:00:00 1.8500000238    2     2013-08-28 01:00:00 1.8799999952    11     2013-08-28 03:00:00 1.9900000095    7     2013-08-28 04:00:00 1.8511325121    3     2013-08-28 05:00:00 1.8049999476    6     2013-08-28 06:00:00 1.8419520855    4     2013-08-28 08:00:00 1.7999999523    1     2013-08-28 09:00:00 1.9550000429    8     2013-08-28 10:00:00 1.9500000477    5     2013-08-28 11:00:00 1.9510477781    4     2013-08-28 12:00:00 1.9403985739    2     2013-08-28 13:00:00 1.75            3     2013-08-28 14:00:00 1.8400759697    6      So lets say I purchase one of these items at 1.50. The question is: at what selling price will I be able to maximize profit? Also, how can I account for any cyclical price movement? (as a side note, I am somewhat familiar with R and Rapidminer)


title: What is a good voting system that is based on a value and number/frequency of votes?

text: So let's suppose there are _n_ candidates running for office. A bunch of random people begin to selectively vote the candidates, each vote being a number from _1_ to _n_. What is a good way to order the candidates based on their vote number and also the number of votes/frequency of votes? So for example: If candidate A was voted by a bunch of voters, but was the votes were fairly spaced out, he would probably be considered "not as great" as a candidate who was voted less/same amount/more but the votes were shortly after one another. I was thinking of something like Reddit's Karma system, but that only takes into accounts +1's and -1's (and I'm not even sure if Reddit every gave information on how their Karma system works!) * * * So for example: Candidates = [ Bob, Sally, Joe, Sam, Rob] Person A ranks Sally as #1 (best choice) and [Bob, Rob] as #5 (worst choice) and leaves (he doesn't vote for others). Person B ranks Sam as #2, Joe as #4, and leaves. Person C ranks Sally as #1, Joe as #3, and Bob as #5. So in this case it may be compelling to say Sally is close to being #1, while Bob is probably closer to being #4 or #5. But Sam was only voted once as #2 and as such there's not a whole lot of grounds to say "yes he's definitely #2". Similarly, Rob was only voted once as being #5 but there isn't sufficient votes to definitively say he's worse than Bob. As such, the raking may become something like: Sally > Sam > Joe > Rob > Bob You can probably twist this example to have even more voters and create a few outliar, but I hope his gives you a general idea of what I'm looking for.


title: In a neural network with N covariates, are >N hidden units useless?

text: I'm fitting a neural network to this example data I found online: Machine Learning Repository I am cross validating 1 to 10 hidden units (in only 1 layer), and I have the minimum error with 10 hidden units. However, I'm somehow thinking of linearly dependent design matrices when introducing 10 hidden units for only 3 input variables (the level of Red, Green and Blue). Is this concern justified, or can I just use 10 hidden units here? Maybe the (sigmoid) transformation does something to avoid the linear dependency?


title: The statement of homoscedasticity of variance when describing the OLS model

text: In an applied econometrics paper, the author states the model to be estimated as: ![enter image description here](http://i.stack.imgur.com/1NSxJ.png) Why does the author claim homoscedasticity? This isn't making sense to me; can't the _population_ variance-covariance have elements of the diagonal that are conditional upon $t$ if there's something like omitted variables?


title: Can logistic regression's predicted probability be interpreted as the confidence in the classification

text: Can we interpret posterior probability obtained from a classifier that outputs a predicted class value and a probability (for example, logistic regression or Naive Bayes) as some kind of a confidence score that is assigned to that predicted class value?


title: Analysing a balanced mixed models design?

text: I am designing an experiment and I'm in need of some help. Participants are presented with two videos on a website. In one period, a feature (a table of contents for the video) is enabled, while in the other period the feature is disabled. Cognitive load (the DV) is measured during both videos. I know that it's sensible to counterbalance the order of treatments. However, since the participant must watch two videos, the video order needs to be counterbalanced also (even though it is not of interest). Therefore I need 4 groups to represent all of these combinations:               F=Feature     N=No Feature     V1=Video 1     V2=Video 2          Group   Period 1   Period 2     1         F&v1       N&v2     2         N&v1       F&v2     3         F&v2       N&v1     4         N&v2       F&v1      According to this website about counterbalancing (last section) I can use two between subject variables "order" and "combination" to represent the above table of group combinations. So my data looks like this:               Order   Combination cogload1    cogload2     1       1           5.00        3.00     2       1           4.00        6.00     1       2           5.00        7.00     2       2           8.00        5.00      This is where I get stuck. Given that feature/video are encoded in order/combination now, how can I determine the main effect of feature or the interaction between feature and video? Also, I'm still getting to grips with general linear modelling, interactions and interpreting the results from SPSS. Does anyone have any resources that they could recommend?


title: Presenting the indexation of a cross-sectional model ($i = 1,...,N$)

text: Many papers in applied econometrics present the indexation of a cross sectional OLS model as $i = 1,...N$. For example: $(1) Y_i = a + b X_i + e_i$ $E[e_i] = 0$ $i = 1,...,N$. Would it be okay to present it as follows? $(1) Y_i = a + b X_i + e_i$ $E[e_i] = 0$ $i \in \\{1,...,N\\}$. Perhaps researchers do the former because the latter doesn't necessitate that $(1)$ is estimated $\forall i$ (only that $i$ is a member of the set $\\{1,...,N\\}$, and $\\{i\\}$ could be some proper subset of $\\{1,...,N\\}$). Want to know best practice. Thanks.


title: Dimensionally weighted distance between two points in n-dimensional space

text: I am not sure if this question makes sense or is even ideal, but here is my problem: Many proposed distance formulas only give a distance in terms of 1-dimensions. A common distance formula is the Euclidean distance (the most obvious one). Euclidean distance $d_{ecdn}$ of two points located at $(x_1, x_2, x_3, ... , x_n)$ and $(y_1, y_2, y_3, ... , y_n)$, respectively, is given by $d_{ecdn} = \sqrt{\displaystyle\sum_{i=1}^{n} (x_i-y_i)^2}$. Geometrically, in any number of dimensions (but imagine 3 dimensions as an example), if you were to draw a line between the two points, whose length expresses the Euclidean distance, the line itself exists in 1 dimension. Why is this a problem? I am doing a project involve data mining and information retrieval. Basically, each document in a corpus is mapped to a location in an $n$-dimensional space where $n$ is the total number of terms within the whole set of documents. They say that if document A is geometrically closer to document B than document C, then document A is more similar to B than C. Euclidean distance gives the wrong impression of similarity since each dimension matters. I am wondering if anyone has any good algorithms that will calculate the "distance" between two points, but with a weight on each dimensions. An example would be the Pearson Product-Moment Correlation Coefficient. This formula is given by ![pearson](http://upload.wikimedia.org/wikipedia/en/math/6/8/6/68654488b517714870216d44f1ce8459.png)


title: How to understand the notation behind the standarded residuals?

text: I don't understand my lecture notes. It states that standarded residuals is by definition $$R_i=\dfrac{y_i-\hat{y}_i}{\sqrt{\dfrac{1}{n-2}\sum_{i=1}^n (y_i-\hat{y}_i)^2 }}$$ I understand that I can compute $\hat{y}$ by the ordinary least-squares method. But as we have $i$ used as a subscript of $R$ and as an index of summation, I don't understand whether it defines real numbers $R_1,\ldots,R_n$ or what. I'm familiar with notation $f_n=\sum_{i=1}^n$ but $f_i=\sum_{i=1}^n$ is new for me. So, can anyone explain how I should understand the definition? If I just plug $i=2$ to the equation I can see that $R_2$ looks like a reasonable notation but what is $\sum_{2=1}^n$?


title: Statistical analysis of mutant proteins sequences

text: Let's say that there is a table like the one below, where for each mutant (A, B, C, etc.) protein or peptide we have a binding and a functional information (true/false) column, and for each amino acid in that protein, we have some kind of amino acid property like hydrophobicity:               Protein Binding Functional     1  2  3  4  5  6  7  8  9 10       A          0         0        13 96 39 77 70 94 96 29 22 82         B          0         1        94 45  2  2 11 46 50 77  7 99         C          0         1        66 71 97 37 14 77 89 92 12 72         D          1         1        11  8 94 73 16 53  2 27 54 97         E          1         1        31 62 49 51  2 86 91 49 61  7         F          1         0         2 42 65 42 54 41 45  9 71 20         G          0         0        26 44 56 65 61 43 56 90 70 86         H          0         1        54 99 68 64 94 81 85  0 50 84         I          1         1        27 52 76 12 46 38 24 74 11 90         J          1         1         1 58 77 50 72 51 87 99 47 67        What statistical tests would you recommend to:   1. See if the difference between means for each mutant is statistically significant.    2. See if the the aminoacid property can somewhere be related to binding and functional properties. Thanks in advance.


title: Poisson deviance - and what about zero observed values?

text: According to Crawley *), the Poisson deviance is computed as $2 \sum{O \cdot \log ({O \over E})}$ But what if observed value $O$ is zero? *) _Statistical Computing - An Introduction to Data Analysis using S-Plus, Michael J. Crawley, printed by Wiley 2004, page 539._


title: Response variable with exponential distribution - how to analyse?

text: I would like to analyse the following data:   * Number of observations: 430 trees   * 1 response variable: diameter growth of trees   * 2 categorical factors: tree species (9 levels) and treatment (2 levels) The response variable seems to have an exponential distribution (density, density of log-transformed growth). This plot suggests that tree species has a strong influence on dbh growth, and this plot that treatment hasn't. For some species, there seems to be an interaction between species and treatment. I would like to analyze influence of species, treatment and interaction between species ~ treatment. If the response variable was normally distributed, I guess an ANOVA would be adequate. However, due to the exponential distribution, I am unsure about how to proceed further. I made an ANOVA (`aov(dbh.growth ~ species * treat, data=df)` that showed significant `Pr(>F)` values for both factors + interaction (see analysis plots), but are the results usable? What possibilities are there for factor analysis with exponentially distributed data? Are there any transformations I might use the transform the data to a normal distribution? Thanks for any help... I'm really stuck here. ;)


title: Efficient way to populate matrix in R?

text: Given a table comprising three columns,               1 2 0.05     1 3 0.04     2 3 0.001      and a matrix with defined dimensions,               mat <- matrix(NA, nrow = 3, ncol = 3)              [,1] [,2] [,3]     [1,]   NA   NA   NA     [2,]   NA   NA   NA     [3,]   NA   NA   NA      Is there an efficient way to populate the matrix with the entries in the third column of the table with R, without having to iterate over the table, isolate the indices and insert value in a for loop? Thanks.


title: Who said "Tell a story with the data"?

text: Does anyone know the exact reference of a quote with words to the effect: **"Tell a story with the data".** I want to use this quote in my presentation but need to acknowledge it. It is the right thing to do! (I am unsure which tag I should use!)


title: Difference between soft k-means and fuzzy c-means

text: I often see in various tutorials and webpages that explains fuzzy c-means clustering and soft k-means clustering. But i was not able to find any material which differentiates them. Is it that both fuzzy c-means and soft k-means clustering are same or different?


title: Asking tweet classification

text: I want to ask you the process to classify the tweet data. Now, I am working to Twitter data but i have confuse how to classify the tweet data using Mallet Tool. Example; I have 200,000 tweets. The tweets content of about apple company. I want to classify the tweets into "false positive" and "positive". False positive:( not related to apple company) nice fish and chip.. good pie apple :-) I ain’t seen uncle Mac in a while positive (company related ) @Apple please make a squid #emoji so I can put one next to @squidneyy22’s name Maybe #Apple will release #iTunes11 tomorrow. Maybe? Its possible to use Mallet tool? How can i do? Thank your very much.


title: Understanding Differential Expression Analysis in Microarray Experiments

text: Can someone please provide me a simple explanation of how differential gene expression analysis works? I know that this method is described in [1] and used in the R package limma. [1] Smyth, Gordon K. "Linear models and empirical bayes methods for assessing differential expression in microarray experiments." Statistical applications in genetics and molecular biology 3.1 (2004). Thank you


title: Combinef in R HTS package- constrain to keep forecasts positive?

text: When using the combinef function from Rob Hyndman's very useful hts package for forecasting hierarchical and grouped time series, there does not seem to be a way to constrain the optimally combined forecasts to be positive- the starting forecasts can be positive, but can go negative through the reconciliation process. The forecast.gts and forecast.hts functions have an argument to keep forecasts positive, but this does not seem to be an option when using combinef by itself with forecasts obtained by other methods. Am I correct in this understanding, and if so is there a decent workaround?


title: Is there a term for reducing the number of parameters of a distribution in the following way?

text: Let's say I have a general distribution that can be specified with 100 parameters $a_1,...,a_{100}$. Now, let's say that within this family of distributions, there is a subfamily that can be represented with just a few parameters, which impose constraints on the original parameters. A trivial example would be if $a_1=...=a_{100}$, meaning we could represent the distribution with a single parameter. I would like to know if there is a name for this process and if there is a name for the new reduced parameter set. It seems that something like _meta- parameterization_ and _meta-parameter_ might be appropriate, but I could not find anything on the subject.


title: What is the distribution of inverse Wishart times a Fisher-Von Mises?

text: Suppose $p$-dimensional vector $x$ takes a Fisher-Von Mises distribution with parameters $\kappa, \mu$, and $M$ is inverse-Wishart with $n$ degrees of freedom and matrix parameter $I_p$, the $p\times p$ identity matrix, with $x$ and $M$ independent. Is it the case that $Mx$, projected to the sphere, takes a Fisher-Von Mises distribution? _related_ : What transformations preserve the von Mises distribution?


title: How do I conduct moderation using composite variables

text: I am running a few moderation analyses using a hierarchical multiple regression. All my variables are composites of parent and child versions of each scale. I created the composites by adding the z-scores from each scale together. I think this is a problem though because when I create the interaction term by multiplying the scores of the IV and moderator, the negative values become a big issue when multiplied together. None of the textbooks are of much help with this, anyone have any ideas? Thanks in advance, Elliot


title: Show that Y/Z does not have finite expectation

text: The unit interval (0, 1) is divided into two sub-intervals by picking a point at random from inside the interval. Denoting by Y and Z the lengths of the longer and the shorter sub-intervals respectively, show that Y/Z does not have finite expectation. Let X~U(0,1) and U be the length of the shorter of the intervals (0,X) and (X,1); that is Z=$min$(X,1-X) and let Y=1-U.So,Y and Z be the lengths of the longer and the shorter sub-intervals respectively Then to find Expectation I need to find P.D.F of Y/Z but I can't do that. Please Help.


title: Time Series: correcting the standard errors in a huge panel time series data set

text: I have stock returns at every 5 minute interval of each trading day for over 2 years for 40 stocks. I want to run a Fama-Macbeth regression by time interval (5min intervals) and then correct the standard errors correlation using Newey- West in SAS. This is my code:               ods listing close;     ods output parameterestimates=pe;     proc reg data=dset;      by time_interval;      model depvar = indvars; run;     quit;     ods listing;          proc means data=pe mean std t probt;      var estimate; class variable;     run;          proc sort data=pe; by variable; run;          %let lags=1000;     ods output parameterestimates=nw;     ods listing close;     proc model data=pe;      by variable;      instruments / intonly;      estimate=a;      fit estimate / gmm kernel=(bart,%eval(&lags+1),0); run;     quit;     ods listing;      As you might have I guessed, this is very computationally intensive. I can't get any results. Should I run instead the Fama-Macbeth regression by stock rather than time interval? Would I get the same thing? Is it because I am using ods listings?


title: How to measure clustering in time duration?

text: I wish to measure clustering in the duration between stock trading. For example, a trade occurs at 1:59:19 and the next trade follows at 1:59:23 - the inter-trade duration is 4 seconds. I have roughly 50,000 trades per day for a particular stock. I was told that I can square the inter-trade duration in order to capture clustering in trades. Is this correct? What if I have a lot of trades that occur at the same time? Hence, inter-trade duration is 0 seconds. Can this have an impact on my clustering measure?


Let the distribution of $X$ be $U(0,1)$. Let U be the length of the shorter of the intervals $(0,X)$ and $(X,1)$; that is, $Z=min(X,1-X)$ and let $Y=1-Z$ be the length of the larger part. Show that, for $t>1$, $P[\frac{Y}{Z}\leq t]=\frac{t-1}{t+1}$;Find the PDF of $Y/Z$. I am unable to find the distribution function.I can find PDF of $Y/Z$. Please help me to find out $P[\frac{Y}{Z}\leq t]=\frac{t-1}{t+1}$


title: Transform log posteriors to the original posteriors

text: Using a two stage RJMCMC process I obtain a vector of log posteriors pertaining to each model. My goal is to convert these to probabilities after normalization (to use in the second stage of the RJMCMC process). The problem is the log posteriors are large (ex. 1100). To convert to probabilities (not on the log scale) I would like to do:               exp(log.post)/sum(exp(log.post))      However this is not possible due to the large values of the log posterior. Any suggestions on how I can convert? I would like the probabilities on the original scale, not on the log scale.


I have a job doing data analysis in finance. My current job is such that I don't have much exposure to things happening in the rest of my industry or other industries. I have a fair amount of knowledge about Bayesian statistics. I'd like to keep myself marketable, so I am curious what data and statistics skills are currently in high demand and where. The software world is awash in data, so I would expect them to need statisticians really badly, but my impression is that they are not in high demand. A friend of mine has suggested that the software industry primarily needs "big data" skills, not statistics skills per se. What data and statistics skills are currently in high demand and where are they in high demand?


title: How to calculate confidence and prediction bands for a linear regression using Mathematica?

text: I need to calculate confidence and prediction bands for a linear regression that is in the general form of $y=ax$, where $a$ is the multiplier and $x$ the variable. I found a code in http://demonstrations.wolfram.com/ConfidenceAndPredictionBands/ , though it gives the confidence and prediction bands for $y=ax+b$ which is out of my interest. * * * This is the code I wrote based on @shabbychef's answer. Is it correct?               bands[list_, x_, type_, gamma_] :=       Module[{a, nlm, z, xs, ys, alphahat, betahat, sigmahat2, n, xbar,         multiplier},       xs = list[[All, 1]];       ys = list[[All, 2]];       alphahat = Mean[ys];       nlm = NonlinearModelFit[list, a z, {a}, z,          ConfidenceLevel -> 1 - gamma];       betahat = a /. nlm["BestFitParameters"];       n = Length[list];       xbar = Mean[xs];       sigmahat2 = nlm["EstimatedVariance"];       multiplier =         If[type == 1,          Sqrt[n*sigmahat2/(n - 2)*(1/n + (x - xbar)^2/               Plus @@ ((xs - xbar)^2))],          Sqrt[n*sigmahat2/(n - 2)*(1 +              1/n + (x - xbar)^2/ Plus @@ ((xs - xbar)^2))]];       {betahat*(x),        betahat*(x) -          multiplier*Quantile[StudentTDistribution[n - 2], 1 - gamma/2],         betahat*(x) +          multiplier*Quantile[StudentTDistribution[n - 2], 1 - gamma/2]}       ]


title: Adjusting performance by opponent skill level in a simple contest

text: I've struggled with this question for several months and found this website, which is amazing. Imagine a hypothetical league composed of players competing in a simple contest where the object of each contest is to score more than an opponent. Player strategies are simple, existing on a range from focusing on scoring, to focusing on preventing the opponent from scoring, to anywhere in between. Assume that players stick with a certain strategy over time. For each player in each contest, we then know the number of points scored, and the number of points allowed. The challenge is to evaluate how a player performed in a certain contest. Obviously the number of points scored by the player will not be enough since it does not speak to the strategy of the opponent (for instance, an opponent who has a strategy focused on scoring will allow more points to every player). I have thought up an approach, and am looking for feedback from everyone here, especially if this approach has a name or (as usually happens) has been otherwise thought up. It has three inputs:   * The number of points the player scored in the contest   * The average number of points all players scored in all the league’s contests   * The average number of points the opponent allowed in all his previous contests, but not this one. The simple formula looks like:               AdjPointsScored = NumPointsScored – (MeanOppPointsAllowedToDate – MeanLeaguePtsScored)      That basically expresses the player’s performance as +/- what a typical player playing against this opponent might be expected to score. I’m also wondering if it might be more accurate to express NumPointsScored and MeanOppPointsAllowedToDate as Z-scores based on the distribution of MeanLeaguePtsScored, subtract the z-scores, and then turn the result back into units of MeanLeaguePtsScored, but I’m unsure if that will improve the analysis or just needlessly complicate things. Much Ablidged.


title: What are some classic examples of feature selection in classification?

text: Is there a classic example showing the importance of good feature selection in classification? The ideal example would be simple, and very easy to understand. I've been volunteered/instructed to put together a short presentation for work to explain some basics of classification techniques. Many members of my audience have never done classification before, so I'd like to keep the example as simple as possible while still using real data. Even if you can't point me to a data-set, if you could provide a example idea, I would greatly appreciate it. The best idea I've been able to come up with so far, is a really old vending machine that needs to sort out which coins are real and which ones are fake.


title: Simulated annealing and k-means

text: One of my problems http://stackoverflow.com/questions/7783933/clustering-data- outputs-irregular-plot-graph suffers from the curse of dimensionality, which also makes it infeasible for exhaustive search or analytical methods, metaheuristics for combinatorial problems like simulated annealing might help in my endeavour how would I combine simulated annealing with k-means in my example problem? Sorry for the very short question im still trying to get my head around simulated annealing and how it could possibly help. Other questions that could help is the comparison of simulated annealing to pca, vds or how it can be used in conjunction. (sorry abit late for me) Edits are welcome if you undertsand my lack of knowledge and what I might be trying to get at.


title: Algorithms to find most frequent data

text: I want to find most frequent words in some cases that some words stand out. I have two data for each case, the word and the number of times it appears. I want to know which algorithms, techniques can be used to find the most frequent data. At the moment, I'm using the discard for data normalization (Mean + Standard Deviation * 3), but I think there are another techniques. What I can use? **Edit:** I know what I want, but I don't know what is possible because I don't know much about the statistics area. I have a data frame in R which contain words and the number of times each words appears. This data is get from a large database (million rows). What I want to find are the words that are **much frequent** than others. Example, I have a list of commons words that came after commas. I want to differentiate which are the really common words than the words that appearing only because there is a lot of data. **Indication of frequency algorithms/formulas are what I want** (I'm a software developer, not an statistician, I don't know the right names).


title: Bolasso in R, or other model selection techniques for parametric models

text: I can't find any packages which allow me to implement bolasso in R, does anyone know of one? Otherwise, I am interested in model selection techniques which can be implemented in R for logistic models. My sole concern is predictive ability. The data has many covariates (or features) and I would like to investigate up to 3rd or even 4th level interactions and 2nd order terms, this is several thousand covariates (there are ~200k observations) so it is important that the techniques are fast.


title: Are suvey responses of ranking significantly different from mean ranking?

text: I'm a stats dummy, so please use small words. (-: This might not be a valid test, so I'm sure you'll let me know if that's the case. I have some survey responses in which the responses are rankings between 1 and 10. I have a frequency for each ranking (2 1's, 3 2's, 2 3's, 10 4's and so on). Intuitively, it looks like a lot of these responses fall in the middle and a few are on the tails, and I want to see if the reponses in the middle ranges are significantly different than the mean. Is this a valid thing to do? If so, whats the best (simple) test I can do? I tried chi-squred, but a few of my rating categories have fewer than 5 responses, so is that valid?


title: Maximum Likelihood Estimator (MLE)

text: Salmon sometimes carry a parasite anisakis simplex which they pick up when feeding on krill at sea. The number of parasites on each fish might be assumed a random variable $X$ having a probability distribution $p(x)$ below: $$p(x)= \frac{\theta^x e^{-\theta}}{x!}\,,\quad x= 0,1,2...$$ You may assume that $E(X) = \theta$ and $V(X) = \theta$. A random sample of $n$ salmon were examined and the number of parasites per fish, $x_1,$ $x_2,$ $\ldots,$ $x_n$ counted. Show that the Maximum Likelihood Estimator (MLE) of $\theta$ is the sample mean. Show that the MLE is unbiased for $\theta$. Find an expression for the variance of the MLE in terms of the parameter $\theta$ and hence show that the MLE is a consistent estimator of $\theta$. A random sample of $n = 20$ salmon taken from a particular Scottish river were examined and the number of parasites recorded as follows:               xi : 2  0  4  0  1  0  3  2  0  5  1  2  0  3  0  1  6  2  2  7      Using these data, estimate $\theta$ giving approximate 95% confidence limits for the true value. Show that the probability of one or more parasites being found on a fish is given by $$P(X \geq 1) = 1 - e^{-\theta}$$ A random sample of $n = 60$ fish were taken from a different river and examined. Suppose all that is known is that $43$ fish were found to have at least one parasite and $17$ salmon were found to have no parasites. Find an expression for the likelihood of this and hence obtain the MLE of $\theta$ given this information.


title: How to calculate error of percentage ratio?

text: I have the following problem: I am counting a subset of cells from a tissue in Drosophila (fruit fly) at different days after "birth". At Day0 I obtain a population of flies, and dissect some of these flies at each time-point (I.e. I have a population of 30 flies at Day0. I dissect 3 of them right away, 3 more after 1 day, another 3 after 3 days, and so on up to 15 days after birth). To obtain these results I need to dissect (thus killing) the fly, so I cannot assay the differences in a single individual fly changing overtime, but I can assay the differences in the whole population. What I obtain is a n number of around 50 for each timepoint. I could just plot the average for each timepoint on a chart and show the variation overtime this way. I would use standard error of the mean to build up the error bars. My problem is that I need to show these data not as raw data, but as a ratio over the first Day, in order to better show the decrease in the cell number overtime. To do this, I am showing my data as a percentage ratio over day0 [(Day1Mean/Day0Mean)*100]. To make it more clear: Day0 mean: 5 Day1 mean: 4 Day3 mean: 3.5 Day6 mean: 3.5 Day9 mean: 3 Day12 mean: 2.5 Day15 mean: 2 Instead of making a dispersion chart with those data (that are the average of the number of cells I counted) I would like to show them this way: Day0: (Day0/Day0)*100 = 100% Day1: (Day1/Day0)*100 = 80% Day3: (Day3/Day0)*100 = 70% Day6: (Day6/Day0)*100 = 70% Day9: (Day9/Day0)*100 = 60% Day12: (Day12/Day0)*100 = 50% Day15: (Day15/Day0)*100 = 40% Now my problem is this: If I am using the "raw" data, I can just plot them in and use the standard error calculated on the mean. I cannot use the standard error with the second kind of setup since my number is completely different. If it was possible I would have liked to "convert" the standard error to the new visualization. I thought it could be possible with some kind of mathematical transformation, in a similar way as reported in this document: http://www.census.gov/acs/www/Downloads/data_documentation/Accuracy/PercChg.pdf Many thanks for Your help


title: matlab gmdistribution.fit 'Regularize' - what regularization method?

text: I am wondering what is behind matlab 'Regularize' option for method gmdistribution.fit. If it is simply adding a 'little' value to diagonal elements of covariance matrix, so as to make covariance matrix inversible, is convergence of fit algorithm guaranteed? Is some more sophisticated method used here, like shrinkage covariance regularization? Doc for this function is here. In general, what do you think of adding small value (for instance 10E-4!) to diagonal elements? Is convergence guaranteed? Thanks


title: Pearson's r over time

text: I have a Pearson's r for the same variables for 4 different periods of time. (4 different years). How can I find the Pearson's r for all of them put together? Here are the years and the related coefficients: `1969: -.169, 1975: .427, 1982: .308, 2006: .632`. **Can I average these to the overall Pearson's r coefficient?**


title: What if a numerator term is zero in Naive Bayes?

text: I'm trying to predict the probability that a user will visit a particular website based on several factors (day of the week, time since last visit, etc). My question is what to do if one of the numerator terms goes to zero? For instance, suppose I visit www.google.com often, but I've never visited on a Monday. $p(Monday|google)$ is zero. Do I simply remove this term from the equation altogether?


title: How do you write up Tukey post-hoc findings?

text: What is the proper way to write a Tukey post-hoc result? There are several examples with different results? Say you have North, South, East and West.               North N=50 Mean=2.45  SD=3.9 std error=.577 LB=1.29 UB=3.62       South N=40 Mean=2.54  SD=3.8 std error=.576 LB=1.29 UB=3.63      East  N=55 Mean=3.45  SD=3.7 std error=.575 LB=1.29 UB=3.64      West  N=45 Mean=3.54  SD=3.6 std error=.574 LB=1.29 UB=3.65      North is statistically significant with East (sig=.009) and West(sig= .040) but not South (sig=.450). East is statistically significant with South(.049).


title: How to correct for inter-experimental variation in cell culture research?

text: I am testing cell cultures. I grow cells in a "control media" for 2 weeks. During these 2 weeks, the cells produce an antibody, and at the end of the 2 weeks, I measure the levels of the antibody. I also measure cell growth (the number of living cells). The problem is, because it is a biological living system, if I do the experiment 5 times (in exactly the same manner as described above), I will have substantially different numbers for cell growth and antibody production. For example, in 5 experiments I found cell growth: 46, 59, 50, 77 and 74, and antibody production: 2558, 3089, 3498, 2115 and 1958. These variations can even be higher, depending on my test sample. Can anyone suggest, how can I correct for this kind of inter-experimental variation. thanks.


title: Dichotomous DV and 3 points of measurement per condition within each subject

text: I have been reading a lot on this site and I think I am going to analyze my data using a mixed model ANOVA- random effect= subject and 2 fixed effects- comparison and pricing. All of my subjects were in all 4 conditions (a cross between the 2 factors) and answered 3 items in each condition. The catch is that my DV is dichotomous. I read that it is not a good idea to assign a mean (i.e. proportion in this case) of the 3 items for each subject and then do the ANOVA using that. Does this mean that I will have number of subjects * 3* 4 observations?


title: What are Pu and Pz in the Phillips & Ouliaris cointegration test?

text: I see the type parameter for Phillips & Ouliaris Cointegration Test in the `urca` package that should be Pz or Pu. I read the details, but I don't understand what type is useful for me. If I wanted to check the cointegration with that method for two series of stock prices, what type would be better? In other words, what does the following mean? > The test "Pz", compared to the test "Pu", has the advantage that it is > invariant to the normalization of the cointegration vector, i.e. it does not > matter which variable is on the left hand side of the equation. In case > convergence problems are encountered by matrix inversion, one can pass a > higher tolerance level via "tol=..." to the solve()-function.


title: What model to classify a discretized continuous variable

text: Consider a variable $y$ (typically, $y_i$ is something like _number of inhabitants of city $i$_ ) and some given features $X$. Let us assume that these features are continuous (eg. total city area, number of shops, roads). We want to use features $X$ to explain $y$. Although the variable $y$ should naturally be seen as a continuous one, it has been discretized in my training set to five classes: $0 \leq y < 300$, then $300 \leq y < 1.000$, and so on. Thus, even for my training set, I do not know the continuous value of $y$, only its discrete class $C(y)$. I need to use this training data to predict $C(y)$ for the other cities (ie. cities that are not part of my training set). Please note that I only wish to predict to which discrete class each other city belongs (I am not seeking for an accurate prediction of $y$ itself). I am wandering what could be a theoretically satisfying approach to this kind of problem. Should I handle this as a classification problem even though there is an underlying hidden continuous variable? Can I identify borderline data points and give them less weight? Should I convert $C(y)$ to a continuous variable using the middle of the class as an estimate? **I would appreciate any idea or reference on this topic. Thanks to all.**


title: Is this a valid test?

text: I have 5000 cases vs 5000 controls with positive/negative outcome. `chi.test` shows no significance for this data. I get significance when subsetting cases on _VarX_ , 1000 cases vs 5000 controls. I suspect that doing 1000 cases vs 5000 controls makes up the significance. To eliminate this I am randomly picking 1000 controls and running `chi.test` 1000 times vs 1000 cases. Here is the histogram of pvalues and the red line is `Pvalue` for 1000 cases vs 5000 controls: ![enter image description here](http://i.stack.imgur.com/Nc0qA.jpg) From this can I conclude that there is no significance? If there is a better way of doing this please advise.


title: Correct time lag in regression

text: Let's say you have two time series. you know they are highly correlated with the right time lag. Is there a function that's used to calculate the correct time lag, or do you just have to run a bunch of different regressions to find it?


title: Contradicting p-values for Anova and Kruskal-Wallis on same data: Which is right?

text: I have a data file with task duration values for three groups, and I want to determine group effect on task duration (tasks were executed by individuals; each group had 7 different individuals; each individual executed the same three tasks; and the data for one individual in group B was not recorded because of a setup problem during the experiment). I have created from the data file the following box plot (red dots are the means, and "n" is the number of time values in each group): ![enter image description here](http://i.stack.imgur.com/tQu8J.png) and also the following histogram (duration given in "min:sec"): ![enter image description here](http://i.stack.imgur.com/G0o8L.png) My data sample per group is small, and "Shapiro-Wilk normality test" tells me that group A does not come from a normal distribution, and that groups B and C come from a normal distribution. Because groups are small and one group is non-normal, I decided to run Kruskal-Wallis one-way analysis of variance (non- parametric) and its result was:               Kruskal-Wallis rank sum test     data:  Duration by Group      Kruskal-Wallis chi-squared = 4.2811, df = 2, p-value = 0.1176      so I should accept that the effect of the groups was not significant (p-value > 0.05). However, when I run one-way Anova (sanity check just in case Kruskal-Wallis was not the correct choice), Anova's result was:               ANOVA Duration ~ Group                  Df    Sum Sq   Mean Sq F value  Pr(>F)        Group        2 0.0003692 1.846e-04   6.473 0.00293 **     Residuals   57 0.0016257 2.852e-05                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Tukey multiple comparisons of means        95% family-wise confidence level                       diff          lwr           upr     p adj     B-A -0.0055763154 -0.009704328 -0.0014483027 0.0054132     C-A -0.0048032407 -0.008769307 -0.0008371744 0.0138750     C-B  0.0007730747 -0.003354938  0.0049010874 0.8943085      That is, Anova returns p-value < 0.05, that is, it is telling that the group effect is significant (in this case, group A was significantly different regarding B and C). However, because of a small number of samples per group and the fact that one group is not normally distributed, I tend to accept Kruskal-Wallis result, but I am not sure. So my questions are: Should I accept Kruskal-Wallis result because this test is more adequate than Anova for this particular case? How should I have used the box plot and the histogram to decide for the most adequate test? Thanks


title: Basic Multivariate Regression Analysis with R

text: I'm trying to look for simple patterns in weather data. Here's a simplified version of what I'm working with.               day <- c(4, 5, 6, 8)     temp <- c(97, 100, 98, 80)     humidity <- c(62, 46, 50, 55)      Suppose I wanted to predict the temperature on day 7 and day 9 (neither are listed) using the following factors. Which R command would I use to do this? I know how do predictions using 2 variables.   _(updated using Greg Snow's answer)_               fit <- lm(temp ~ day)     newdat <- data.frame( day=c(7,9) )     predict(fit, newdat)      But how would I go about using more than 2? Your help would be greatly appreciated.


title: How to find the inflection point from my density plot curve?

text: I sampled the data qtl_500(50,000 points), and get the density curve of it as follow: > plot(density(qtl_500k)) ![enter image description > here](http://i.stack.imgur.com/hajiN.png) Now I hope to get the inflection point of this curve. Could anyone please give me some idea? Thank you very much!


title: Given an interaction how to find the region of var1 at which var2 has significant effect on response Y?

text: I have run a binomial `glm` in R using proportional data in the form `y <- cbind`. I have two continuous explanatory variables, from which `var1` is significant and `var2` is not, but the interaction is also significant. So I want to find which are the values of `var1` for which `var2` results in an effect on `y`. I have discovered a function called `getRegion()`, which I found in two tutorials on the internet, but I haven't found the package that contains it. Are there any other functions that could work in a data set like mine (i.e., proportional response to two continuous explanatory variables)?


title: auto.arima from Forecast package

text: I am trying to fit a time series using the function auto.arima and I face some strange results. As a first try, I use the command               auto.arima(data,d=0,D=1,max.p=2,max.q=2,max.P=2,max.Q=2,max.order=8, xreg=xreg_past,trace=TRUE,ic="aic")      The model I get is an ARIMA(2,0,2)(0,1,1)[12] with an AIC equal to -300.14. But since I know that this command will make use of the stepwise selection algorithm, I want to make a try with the tests of all possible models using the option stepwise=FALSE. I thus try the command               auto.arima(data,d=0,D=1,max.p=2,max.q=2,max.P=2,max.Q=2,max.order=8, xreg=xreg_past,stepwise=FALSE,trace=TRUE,ic="aic")      And now, the model I get is an ARIMA(0,0,2)(2,1,0)[12] with an AIC equal to -293.14. Since my second attempt takes all the models into account, this result is strange as the previous model had a lower AIC. Furthermore, If I take a look in the trace of the last function call, I see that the ARIMA(2,0,2)(0,1,1)[12] model has now an AIC of -245.13 which explains why it has been rejected. Why did the AIC value change ? At least, if I use the simple command               arima(data, order=c(2,0,2), seasonal= list(order=c(2,1,2), period=12), xreg=xreg_past)      I get an AIC value of -319.15, which is better that the two models provided before. I think I am missing something important but I am not able to see what. Can somebody help me ? Thanks in advance, Regards, Ludo


title: Do methods based on trimmed means require homoscedasticity?

text: I understand (from Wilcox' book [1]) that methods using the trimmed means are robust to distributional assumptions. Do the methods (e.g. `bwtrim`, from the `R` package WRS) assume that the distribution is the same (although not normal) across all conditions being modelled? In other words, what if the data have heteroscedasticity (more variance in some conditions than others) indicating that distibutions are not the same across conditions? Example:               library(reshape2)     library(ggplot2)     library(WRS)          data <- data.frame(     A1=rnorm(100),     A2=rnorm(100, 1, 2),     B1=rnorm(100, 5, 5),     B2=rnorm(100,10,10)     )          mdata <- melt(data)     qplot(variable, value, data=mdata)      ![qPlot](http://i.stack.imgur.com/wQGnK.png)               # Would it be appropriate to use:     bwtrim(2, 2, data)          $Qa     [1] 149.6878          $Qa.siglevel          [,1]     [1,]    0          $Qb     [1] 43.91871          $Qb.siglevel                  [,1]     [1,] 8.856752e-09          $Qab     [1] 16.68218          $Qab.siglevel                  [,1]     [1,] 0.0001269782          # or does that assume homoscedasticity?      Could anyone recommend an alternative for analyzing this kind of data (i.e. unequal variance)? [1] Wilcox, R. Modern Statistics for the social and behavioral sciences. CRC Press, 2012.


title: How does extreme random forest differ from random forest?

text: Are they more efficient implementation -- is the difference important from practical point of view, there is R package which implements them. Is it new algorithm which overcomes "generic" implementation (RandomForest package from R) not only in terms of efficiency or also in some other areas? Extreme Random Forest http://link.springer.com/article/10.1007%2Fs10994-006-6226-1


title: Preliminary estimates of ARIMA in R?

text: We know that dealing with model involving MA factors is not easy to estimate, since there are past values of errors to be computed recursively. And this recursive estimation requires preliminary (initial) estimates of the parameters. For example, an ARMA(1,2) $$z_t=\phi z_{t-1}-\theta_1 \varepsilon_{t-1}-\theta_2 \varepsilon_{t-2}+\varepsilon_t$$ To estimate the parameters, we need to compute first the values of $\varepsilon_{t-1}$ and $\varepsilon_{t-2}$, since these are not available yet. And they are computed using $$\varepsilon_t=z_t-\phi z_{t-1}+\theta_1 \varepsilon_{t-1}+\theta_2 \varepsilon_{t-2}$$ Procedures for obtaining preliminary estimates of the parameters is available in Box and Jenkins, Time Series: Forecasting and Control. And this estimation is already available in much statistical software. My question is, "Is there a function for obtaining a preliminary estimate of the parameters in R?" I need this to obtain a preliminary estimate for my Space-Time ARIMA. Another question is, "How does `arima` function of R compute preliminary estimates of the parameters when there are MA factors involved?" Thanks in advance!


title: Test with Student's t

text: What kind of null hypothesis can be tested using a Student's t? It's very likely to be one of the questions of the statistics exam. In order to test a null hypothesis using the Student's t, does it need to have particular features? Thank you.


title: online resources for spatial analysis

text: I am wondering about some good online resources for learning how to implement spatial statistical techniques, like Moran's I and Geary's C, or even K Nearest Neighbors...any help is appreciated. Thanks!


Is it possible to extract fixed effects in the fixed effects model in Stata? I know that in R I can do something like this with `fixef` command:               withinmodel <- plm(formula=y~x,data=Data, model="within")             summary(fixef(withinmodel))      But is it possible in Stata?


title: How do I obtain prediction intervals for fixed effects of Generalized Linear Mixed Model?

text: Suppose I use R to fit a Generalized Linear Mixed Model from the binomial family and with a logit link. How do I obtain the prediction intervals (as opposed to the confidence intervals) for the fixed effects (as opposed to the random effects) in a way that incorporates the variability described by the random effects? Thank you.


title: Combining categorical and continuous variables to calculate a factor

text: I have a categorical predictor (income) and three continuous predictors (area, no. of bed rooms and no of cars). How can I form a single factor from these variables? In other words I want to combine all these variables to get one single variable.


title: What is the proper way to compare two estimated densities using sample data?

text: Say if have a dataset $X \subset \mathbb R^d$. I have two candidate probabilistic models M1 and M2 (e.g., M1 is a mixture of 2 gaussians and M2 is a mixture of 3 gaussians). I want to know which model best describes my data. I am wondering if it is ok to do a cross-validation similar to what we would do in supervised learning, where each fold $i$ would involve the following:   * separate $X$ into $X_\mbox{train}^i$ and $X_\mbox{test}^i$   * fit M1 and M2 using the train sets: $$\theta_1 = \mbox{argmax}_\theta {\cal L}_1 ( \theta, X_\mbox{train}^i) $$ $$\theta_2 = \mbox{argmax}_\theta {\cal L}_2 ( \theta, X_\mbox{train}^i) $$   * compute the likelihoods on the test sets: $$ {\cal L}_1 ( \theta_1, X_\mbox{test}^i) $$ $$ {\cal L}_2 ( \theta_2, X_\mbox{test}^i) $$ In the end, the best model would be the one with the highest average likelihood over the test sets $$\max( \sum_i {\cal L}_1 ( \theta_1, X_\mbox{test}^i), \sum_i {\cal L}_2 ( \theta_2, X_\mbox{test}^i) ) $$ Is this valid?


title: What are some Bayesian alternatives to the Kolmogorov-Smirnov test?

text: What are some Bayesian ways to see if two empirical distributions come from the same data-generating process? Are any of the Bayesian alternatives implemented in Matlab/R/Python?


title: Principal component analysis in 3D Scatterplot

text: I am new to principal component analysis (PCA). I performed PCA for a dataset with 54 samples. When I project them in 3D scatterplot, I can see samples with similar characteristics are grouped together separately. The X,Y and Z axes in 3d scatterplot represent PC#1, PC#2 and PC#3 respectively. Along the axes positive and negative values are represented. What does these values convey, especially negative values ? If a sample is found along an axis with negative value, what does that imply? Also the overall variance all 3 PCs is 40% (PC#1-20%,PC#2-13% and PC#3=7%). What does that imply? Why it is not 80-90%? Is my data of good quality?


title: PCA, LDA (Linear Discriminant Analysis), proportion of variance

text: I have some basic questions regarding PCA and LDA (Linear Discriminant Analysis ) but I am a bit lost and I will appreciate the help.   1. In PCA there is a way to calculate the proportion of variance explained, is it also possible to do it for LDA?. If so, how?    2. Is the “Proportion of trace” output from the lda function (in R MASS library) is equivalent to “proportion of variance explained”?


title: Calculate $\mathbb{E}[Z_i]$ where $Z_i = \min(X_i, Y_{i-1})$, $X \sim Beta(\alpha,1)$

text: Let $X$ be a IID random variable with support in $[0,1]$ and CDF given by a Beta distribution, i.e. $X \sim \mathrm{Beta}(\alpha,1)$. Let $Z_i = \min(X_i,Y_{i-1})$, $\forall i >1$. I would like to compute $\mathbb{E}[Z_i]$, knowing that: \begin{align} Y_1 &= Z_1 = X_1, \; \mbox{when} \; i=1,\\\ Y_i &= \frac{1}{i} \sum \limits_{t=1}^{i} Z_t, \forall i >1\. \end{align} How can I compute $\mathbb{E}[Z_i]$ ?


title: Repeated measures with multiple time points for the predictor and dependent variable: Does Xt-1 predict Yt better than Yt-1 predict Xt?

text: I have a question on what type of analysis I should be looking into to analyze some data I have: Suppose you have 2 runners X and Y, and they take turns sprinting 100 meters, with runner X going first. Because the two runners are able to see each other's performance, their performance on a certain trial (Xt or Yt) might be influenced by the other's previous trial (Xt-1) or (Yt-1). If you collected data for 10 time points, how might you test to see who is more susceptible to the other's influence? I've been doing some reading on distributed lag models, but I'm not sure if this is the right approach.   * Is the distributed lag model generalizable to answer my question, where I have multiple dependent time points?   * Is this even the right approach? If not, what should I be reading up on? Thank you!


title: Obtaining factor scores from a factor analysis of mixed types (both ordinal and scale) variables

text: I have to analyze a data set having both ordinal and scale variables. What I was trying to do was to factor analyze the data. The data was in SPSS. But I think SPSS factor analyzes the data by calculating Pearson's correlation by default which might be misleading as my data set contains mixed types of variables. I could use `hetcor` in the R package `polycor` and then use the correlation matrix as an input for factor analysis using the commands in package `psych`. But most of the commands do not provide me with factor scores which is very important to obtain for the purpose of my analysis. This is because my input is the heterogeneous correlation matrix, not the data itself. I heard that using `fa.poly` would give me the factor scores. But whenever I run this command it produces the following:               > fa.poly(dat,nfactor=4,rotate="varimax")      Polychoric |.      How can I extract the factor scores? Do I have any alternative? For your convenience to answer this here is the data. I need the factor scores to be calculated. What do you think should be an appropriate procedure? Thanks and best regards.


title: R mgcv gam: variable can be used as a predictor but not as an offset - Why?

text: I am trying to use the R mgcv package to fit a gam. I have one variable that I want to use as an offset. It works fine as a predictor, but when I use it as an offset gam fails with an error message about non-conformable arrays.               > x <- gam(gbf_3 == "G" ~ s(pred_lo) + s(avm_prop_21_mult), family=binomial(), weights=round(1000*wt), data=ma, subset=(!is.na(score_2)))     > x <- gam(gbf_3 == "G" ~ pred_lo + s(avm_prop_21_mult), family=binomial(), weights=round(1000*wt), data=ma, subset=(!is.na(score_2)))     > x <- gam(gbf_3 == "G" ~ offset(pred_lo) + s(avm_prop_21_mult), family=binomial(), weights=round(1000*wt), data=ma, subset=(!is.na(score_2)))     Error in x %*% coef + offset : non-conformable arrays      Why would a variable be acceptable as a predictor but not as an offset and what would I have to do to use the variable as an offset? Thanks. Ross


title: Bootstrap sampling from empirical data, is this the best way to estimate the underliyng distribution?

text: I have observed an experiment where 500 individuals (one after another) touch surfaces in a room and then leave. So I have a data set of recorded number of surface contacts $n$ for each person. I need to make a pdf of this but not sure how best to do it: This is a histogram converted to pdf:![enter image description here](http://i.stack.imgur.com/8cmOo.png) Would you recommend using a `bootstrap` sample to calculate the pdf of this distribution by `kdensity`?


title: Why use control variables in differences-in-differences?

text: I have a question on the differences-in-differences approach with the following standard equation: $$ y= a + b_1\text{treat}+ b_2\text{post} + b_3\text{treat}\cdot\text{post} + u $$ where treat is a dummy variable for the treated group and post. Now, my question is simple: Why do most papers still use additional control variables? I thought that if the parallel trend assumption is correct, then we should not have to worry about additional controls. I could only think of 2 possible reasons for why to use control variables:   1. without them, trends would not be parallel    2. because the DnD specification attributes any differences in trends between treatment and control group at the time of treatment to the intervention (i.e. the interaction term treat*post) - when we don't control for other variables, the coefficient of the interaction may be over-/understated Could anyone shed some light on this issue? Do my reasons 1) or 2) make sense at all? I don't fully understand the use of control variables in DnD.