title: Determining characteristics of sampling sets for EFA/CFA/SEM

text: Dividing sample data into several sets seems to be a common approach in statistics. This is especially evident in predictive modeling, where samples are traditionally divided into two sets, usually called _training set_ and _test set_ ( **two-way partition** ). Training sets are supposed to be used for _developing models_ and test sets for _evaluating models_. Recently, I became aware that a **three-way partitioning** scheme is also used, where data set for developing models is called _training set_ , for evaluating models - _validation set_ , and for evaluating the selected model - _test set_. I have also read about specific implementations of two-way partitioning, such as _multi-fold cross-validation_ , as well as _bootstrapping_ , which is popular in PLS-SEM. My **question** is threefold: 1) What **sampling approach** for data analysis is **correct/best/optimal** , considering that I plan to perform EFA and CFA, with following SEM? 2) What are the **guidelines (heuristics)** for training/test and training/validation/test **ratios** in terms of the _size_ of each subgroup of a sample? Any specific **strategy**? 3) Do I need to care about **specific partitioning implementations** (cross- validation, bootstrapping, etc.), or just good old **simple random sampling** will be sufficient? **Please include references, if possible. Thank you much in advance!** UPDATE: After giving these ideas some additional thought, I came to the following realization. It seems to me that a combination EFA/CFA basically represents a two-way partitioning approach (use training set for EFA and test set for CFA), while SEM can benefit from a three-way partitioning, when using alternative models variant of SEM. Please comment on this realization as well.


title: What is the Net Promoter Score Data Type

text: As many people know Net Promoter Score is a customer feedback measure calculated from the response to the question: "How likely would you be to recommend [us] to a friend or colleague". The response is a discrete value from 0 [very unlikely] to 10 [very likely]. The responses are then used to calculate NPS: NPS = % of 0-6 responses less % of 9-10 responses. In an effort to understand what drives NPS there is a desire to identify high correlations and models that have NPS as the dependent variable and other business variables as independent variables. Now, although NPS looks like a continuous variable it is calculated from Ordinal values (ranges of discrete values). So can it be treated as continuous for the purposes of these calculations? If not, why not and what are the alternatives.


title: Kendall correlation between $X$ and $X^2$

text: Let $X \sim N(0, 1)$. Can you confirm that the Kendall's correlation coefficient between $X$ and $X^2$ is equal to $0$. The way I interpret this, is that Kendall's tau only measures monotonic dependence.


title: How can I specify cross-correlations between differential equations?

text: for my Master thesis I want to model co-movements between two commodity forward curve. First, I specify the model for a single forward curve $F(t,T)= \bar{F}_0(t)e^{s(T_i)-y(t,T-t)(T-t)}$ where $F(t,T)$ is the forward curve at time $t$ for maturity $T$, $\bar{F}_0(t)$ is the forward curve level (is my first factor), $s(T_i)$ is the seasonal premium per month for $i^{th}$ year of maturity (is deterministic), $y(t,T-t)$ is the convenience yield (second factor). Now I want to simulate and describe the dynamics of the two factors, $d(\log \bar{F}_{0,j}(t))=dX_i(t)=k_i (\mu -X_i (t))dt+ \sigma dW_i^1 (t) \mbox{ with } d W_i^1dW_j^1=\rho^1_{i,j}dt$ where $dW_i^1$ is a wiener process for commodity i and $dW_j^1$ for commodity j, The first factor of commodity i should be correlated to the first factor of commodity j. The dynamics of the second factor are, $dy_i(t, T-t)= -a^{T-t} y(t,T-t)dt + \eta _i dW_i^2(t) \mbox{ with } dW_i^2dW_j^2=\rho^2_{i,j}dt$ again, from my historical dataset the convenience yield of both commodities should be correlated with each other.. Therefore, $dW_i^2dW_1^2=\rho^2_{i,j}dt$. So far, I do not think I made any mistakes, but now I want to specify the correlation relationship for the same commodity but between the two factors? How can I specify the correlation relationship for one commodity, i, but for different factors? I thought may be something like: $dW_i^1dW_i^2=\rho^{1,2}_{i}dt$ Thanks in advance for your reply


title: Find the UMVUE for $\Phi(\mu)$

text: > Question: Suppose $X_1, \cdots, X_n$ are $iid$ normal random variables with > unknown mean $\mu$ and known variance $\sigma^2$. Find the UMVUE for > $\Phi(\mu)$, where $\Phi(\cdot)$ is the cdf of a standard normal random > variable. I used to guess the desired UMVUE is or similar to $\Phi(\bar X)$ since it is a function of complete and sufficient statistic. However, it is not really unbiased (see here). **Could anyone shed light on how to find the UMVUE, please? Thank you!**


title: SPSS CHAID/CRT Query

text: Hello I am using SPSS trees - CART functionality for classification. I have a number of classifying variables such as Business_Size, Location, Previous_Record etc. I want to know how do I set configuration so that a particular classification variable is used at only one level of the tree. E.g say at the first level the records were classified using Location then at no subsequent level of the tree, Location should be used. Currently, at the first level sub-trees are getting formed based on location, then on the next level sub-trees are getting formed using Business_Size but on the next level again Location variable is getting used for further classification. Any help on this would be appreciated. My colleague told me that in SAS such a setting is possible. I do not have SAS. So want to know how it can be done in SPSS.


title: Finite mixture models with bounded data

text: I am trying to fit a finite mixture model to a dependent variable which is bounded (practically) between -0.594 and 1 (theoretically, the latent variable is bounded between -Inf - 1). The data are also bimodal, with a large number of values at '1'. The objective of the analysis is prediction of the dependent variable. My current approach has been to fit a mixture of normal distributions using the `flexmix` package in R, but I'd really like to account for the bounded nature of the data, as a recent study found this to be important (I also choose k=3 components based on this study). Using `flexmix` for truncated data appears non-trivial, as suggested here. Is there an R package that will permit mixture models with bounded data? I've noticed that actually predicted values do not seem to fall outside the bounded range; i.e. predicted values are not in practice greater than 1. Is this just a fluke of my data, or is it a feature of the methods I've used? Is the bounding even a problem in this context? As an alternative, I've tried transforming the data by simply taking 1-the dependent variable, thereby giving me a (zero-inflated) variable bounded by 0 and Inf which I have tried to model as a mixture of zero-inflated poisson models but I get the error:               Error in FLXfit(model = model, concomitant = concomitant, control = control,       : 1 Log-    likelihood: NaN      Is it possible to model non-integers with the poisson family in this context? Any suggestions or thoughts would be greatly appreciated, I'm very new to mixture modelling and indeed GLMs etc. Here's some simulated data: https://dl.dropbox.com/u/65336009/mydata.csv Here's my code:               require(flexmix)     require(ggplot2)     mydata <- data.frame(read.csv("mydata.csv", head=T))     attach(mydata)          #Plot of y var     summary(y)     ggplot(mydata, aes(y)) + geom_histogram(binwidth = .1)          #Simplified example of my current 'best' approach####     m1 <- flexmix(y ~ x1 + x2 + x3,                   data = mydata,                   k = 3)          #Predict cluster membership     clusters <- data.frame(clusters(m1, newdata = mydata))          #Predict y     a <- data.frame(predict(m1, newdata = mydata))          #Select prediction based on predicted cluster membership     mydata$flexmix.norm <- ifelse(clusters[,1]==1, a[,1],                                        ifelse(clusters[,1] == 2,                                               a[,2], a[,3]))         print(max(mydata$flexmix.norm))          #Plot predicted values     ggplot(mydata, aes(flexmix.norm)) + geom_histogram(binwidth = .1)          #Maybe it's more natural to model as 1 - y, which is bounded (0,Inf) ####     y.d <- 1 - y     ggplot(mydata, aes(y.d)) + geom_histogram(binwidth = .1)          #Error here ***     m2 <- flexmix(y.d ~ x1 + x2 + x3,                   data = mydata,                   k = 3,                   model=FLXMRziglm(family="poisson"))     rm2 = refit(m2)          #Predict cluster membership     clusters <- NULL     clusters <- data.frame(clusters(m2, newdata = mydata))          #Predict y (note back on original scale of y)     b <- 1 - data.frame(predict(m2, newdata = mydata))          #Select prediction based on predicted cluster membership     preds$flexmix.pois <- ifelse(clusters[,1]==1, b[,1],                                   ifelse(clusters[,1] == 2,                                          b[,2], b[,3]))          ggplot(mydata, aes(flexmix.pois)) + geom_histogram(binwidth = .1)      Thanks


title: Dummies instead of the Chow test

text: I have found somewhere a mention to the possibility of using dummies variables instead of the Chow test to test whether the coefficients in two linear regressions on different data sets are equal. Could someone please direct me to a reference where I could study the details?


title: logistic regression with independent categorical variables with more than two possible values using stata

text: I have to do a logistic regression with independent categorical variables with more than two possible values. Which is the best way to deal with such variables using Stata or spss? I need to have Odds ratios in the output. Thanks.


title: Prove F test is equal to T test squared

text: I need to show that F test is equal to T test squared, when the T test is for 2 independent groups and assuming variances are equal. I know that $F=\frac{MSB}{MSW}=\frac{SSB/k-1}{SSW/N-K}$ and I know that $T=\frac{X-Y}{S_p \sqrt{\frac{1}{n}+\frac{1}{m}}}$, so $T^2=\frac{(X-Y)^2}{S_p^2 ({\frac{1}{n}+\frac{1}{m}})}$ I've seen this proof in Regression but here we're not using MSE and MSR, so i'm not sure how to connect between the two. Any help will be much appreciated! :)


title: How helpful is a quant job in Goldman Sachs for later PhD in Machine Learning?

text: I am a masters in Computer Science and am interested in pursuing a career in Machine Learning, possibly academic, in the long run. I have been offered a position related to Financial Modelling at Goldman Sachs that mostly would provide experience in optimization and stochastic methods, which are also required for developing Machine Learning algorithms. How helpful or detrimental would the given job experience be in   1. chances of getting into PhD on Machine Learning in a top/moderate US/Europe Grad School   2. research during PhD after working few (around 2-3) years in the given job ? Any examples of similar shifts, possibly with ensuing results, would be extremely helpful


What is the significant difference between the Fligner Wolfe test and 2-sample Wilcoxon test?   How can I compare the treatment groups with the control groups of different studies after pooling?


title: How do we predict rare events?

text: I am working on developing an insurance risk predictive model. These models are of "rare events" like airline no-show prediction, hardware fault detection, etc. As I prepared my data set, I tried to apply classification, but I couldn't obtain useful classifiers because of the high proportion of negative cases. I don't have a lot of experience in statistics and modeling data beyond a high school statistics course so I'm kinda confused. As first thought, I have been thinking of using an inhomogeneous Poisson process model. I classified it based on event data (date, lat, lon) to get a good estimate of the chance of a risk at a particular time on a particular day in particular place. I'd like to know, what are the methodologies/algorithms to predict rare events?   What do you recommend as an approach to tackle this problem?


title: Reporting interaction effect for a 2x3 ANOVA

text: I have used a 2x3 mixed between-within subjects ANOVA to assess the impact of three feedback conditions (Grade, Comment and No Feedback) on participants’ scores on the Specific-Self Efficacy Scale on two occasions, pre- and post- task. I did not get a significant main effect but I did get a significant interaction effect. My tutor has told me that I can interpret the results from the interaction graph alone and that there is no need to run tests to see if there are significant differences between pre and post task scores in each of the feedback conditions. However, when discussing these results in my manuscript, I do not feel confident in relying on a graph alone. Is it appropriate to only use the interaction graph to discuss results or we you recommend running some paired samples t-tests? Thank you for taking the time to read my query!


title: SVAR - Time series econometrics

text: How do I determine the Cholesky ordering for a set of endogenous variables to obtain the variance decompositions in eViews? Is there anyway we can cross check if our order is right apart from intuition.


title: Testing whether two regression coefficients are significantly different (in R ideally)

text: If this is a duplicate question, please point to the right way, but the similar questions I've found here haven't been sufficiently similar. Suppose I estimate the model $$Y=\alpha + \beta X + u$$ and find that $\beta>0$. However, it turns out that $X=X_1+X_2$, and I suspect $\partial Y/\partial X_1 \ne \partial Y / \partial X_2$, and in particular, that $\partial Y/\partial X_1 > \partial Y / \partial X_2$. So I estimate the model $$Y=\alpha + \beta_1 X_1 + \beta_2 X_2 +u$$and find significant evidence for $\beta_1,\beta_2>0$. How can I then test whether $\beta_1> \beta_2$? I considered running another regression $$Y=\alpha +\gamma(X_1 - X_2) + u$$ And testing whether $\gamma>0$. Is this the best way? Also, I need to generalize the answer to many variables, i.e. suppose we have $$Y=\alpha + \beta^1 X^1 + \beta ^2 X^2 + \dots + \beta^nX^n + u$$ where for each $j=1,\dots,n$, $X^j=X_1^j+X_2^j$, and I would like to test for each $j$ whether $\partial Y/\partial X_1^j \ne \partial Y / \partial X_2^j$. By the way, I am primarily working in R.


title: Timeline of machine learning and data mining breakthroughs

text: Is there any timeline or historical overview of the most important breakthroughs in machine learning and data mining?


title: What to do With Curvilinear Relationships?

text: If I have determined a curvilinear relationship between my dichotomous y and continuous x, what should I do before running a logit regression? Should I log transform my x variable, etc?


title: Multiple hypothesis testing and F tests?

text: Does the multiple hypothesis testing problem apply to the calculation of an F statistic for joint significance? It seems to me that the more variables you are including in your test for joint significance, the more you are accepting that any one of your tests will produce a false positive, right? If so:   1. For all intents and purposes, is it really that big of a deal?   2. Is there a way to get around this? * * * EDIT: The responses make me think that I am understanding either the multiple hypothesis testing problem or F tests wrong. Here's an explanation of the conflict that is in my head, which may be incorrect. My understanding of the multiple hypothesis testing problem is this: Our alpha level (ex. $\alpha=0.05)$ is our accepted Type I error level, which is sort of a theoretical concept. If we are testing multiple hypotheses simultaneously, like $H_1: x_1 = 0$ and $H_2: x_2 = 0$, then we are implicitly testing $H_0:(x_1 = 0 |x_2=0)$, right? In which case we would add $0.05+0.05$ to get the probability of making a Type I error for the joint test, right? And it's my assumption that theoretically, this is what an F test does. So if you are running an F test on 20 variables, for example, then you are guaranteed, theoretically, to get a Type I error, right? Or, now that I think about it, perhaps I am understanding Type I error incorrectly. Any help would be appreciated.


title: Correlations with a linear combination means correlation with individual variables?

text: I was having this conversation with someone, and I feel like I have an explanation, but would greatly appreciate if someone could verify whether what I'm thinking is just completely wrong. He made the argument that if you know that a variable [a] is correlated with a linear combination of variables [C], then you know something about the correlations between [a] and the constituent variables of [C]. And vice versa. I feel like this is wrong, but I don't have a great explanation of why. There are two reasonings I've come up with, but both feel like I just pulled them out of nowhere, and I could use some feedback:   1. Creating a linear combination of variables, variables which potentially are uncorrelated with each other, could be seen as equivalent to adding random noise to any one of those variables, which would mean that you can't really say anything about the correlations with some outside variables after the combination.   2. I had this explanation that is very vague in my head and involves vectors and conic sections, but I'm not sure how to make it make sense. Basically, if you know the correlations between [a] and the constituent variables of [C], you can construct a particular vector space within which the correlation of [a] and [C] must lie, but I don't think that works in the opposite direction if you start with a correlation between [a] and [C]. Any feedback would be greatly appreciated.


title: A distributional test based on entropy and self-information

text: Say that I have a real-valued discrete distribution $p(x)$ and $N$ samples, $x_1, \ldots, x_N$, and I want to test whether the samples came from the distribution without making any further assumptions whatsoever. Note that there are very few samples; in the application that motivated me to make this post, we have $N = 5$. Thus Kolmogorov-Smirnov and Chi-squared tests are not expected to have much power. I had a simple idea for doing this under the assumption that one can efficiently sample from $p(x)$. Being a bit statistically naive, I'm having difficulty figuring out whether it exists in the literature or not, and hoping someone can point me to the right resource. The idea in a nutshell is to compare the self-information of the sample, $\hat{I}_N$, to the distribution of the self-information $I_N$ of $N$ random samples from $p(x)$. Formally, recall that the self-information of a random variable $X$ having distribution $p$ is given by $$ I = - \log p(X), $$ and the self-information of $N$ iid random variables $X_1,\ldots X_N$ each having distribution $p(x)$ is given by $$ I_N = -\sum_{i=1}^N \log p(X_k). $$ The self-information of the sample data we have, $x_1,\ldots,x_N$ is denoted $$ \hat{I}_N = -\sum_{k=1}^N \log p(x_k). $$ As a test statistic one might consider $\hat{C} = C(\hat{I}_N)$, where $C(s) = \mathbb{P}(I_N < s)$ is the cumulative distribution function. If $\hat{C}$ takes a value very close to 0 or 1, then it has an extreme value compared to the distribution of $I_N$, and is unlikely to come from $p(x)$. For concreteness, we could use the standard thresholds for extreme value tests, such as 0.95 and 0.05 for the high and low ends respectively. As an example application of this test, say that $p(x)$ is some strange multimodal distribution with multiple humps, and the samples $x_1,\ldots,x_N$ lie in the valleys between the humps. It is not clear that there is any test suitable from the literature for such a problem, but intuitively we can see that the samples are unlikely to have come from $p(x)$ because the values of $p(x_i)$ are so small, or equivalently, the self-information $-\log p(x_i)$ is much larger than typical. In terms of the above discussion, $\hat{C}$ will be very close to 1, and the hypothesis will be rejected. The main problem I see with this test is that the distribution of $I_N$ may be difficult to compute. That may be, but in many cases I would imagine that a few million (or billion) Monte Carlo samples would suffice to get a good approximation of the distribution. Analytical/asymptotic approximations could be used to speed things up, get theoretical results, etc. For example, the first moment of $I_N$ is $N$ times the Shannon entropy, and higher moments could be computed without great difficulty.


title: Selecting the best model using cross-validation on coefficient of determination and/or mean squared error

text: Sorry if this question is oft repeated. Let's say I'm doing regression and I want to know whether I should use Linear Regression or Random Forests. I do 10-fold cross validation on each model to obtain scores of R^2 and Mean Squared Error. How do I then compare these lists of numbers to see which is better? Do I choose the one that has the highest mean R^2? Or what about the lowest mean MSE? Or maybe the highest minimum R^2, or the lowest maximum MSE, or some weighted average of these quantities? What if the R^2 is better on one but the MSE is better on another? Is it arbitrary? Does it depend entirely on the particular problem I'm trying to solve? Is there a generally accepted silver bullet?


title: In R, does "glmnet" fit an intercept?

text: I am fitting a linear model in R using `glmnet`. The original (non- regularized) model was fitted using `lm` and did not have a constant term (i.e. it was in the form `lm(y~0+x1+x2,data)`). `glmnet` takes a matrix of predictors and a vector of responses. I've been reading `glmnet` documentation, and can find no mention of the constant term. So, is there a way to ask `glmnet` to force the linear fit through the origin?


title: Are there any recommended approaches for analysing data from genetic algorithms?

text: After running a study based on interactive genetic algorithms, I have a univariate data file containing multiple participants each doing multiple generations (blocks) of multiple trials. Is there an established approach used to analyse this? I've looked at using a grand mean on each of the component variables. I've also looked at the modal values of these variables after filtering out all but the final generation. Neither approach seems to properly capture the richness of the data, so I'd like to know how other people would handle this.


title: How does one do a Type-III SS ANOVA in R with contrast codes?

text: Please provide R code which allows one to conduct a between-subjects ANOVA with -3, -1, 1, 3 contrasts. I understand there is a debate regarding the appropriate Sum of Squares (SS) type for such an analysis. However, as the default type of SS used in SAS and SPSS (Type III) is considered the standard in my area. Thus I would like the results of this analysis to match perfectly what is generated by those statistics programs. To be accepted an answer must directly call aov(), but other answers may be voted up (espeically if they are easy to understand/use).               sample.data <- data.frame(IV=rep(1:4,each=20),DV=rep(c(-3,-3,1,3),each=20)+rnorm(80))      **Edit:** Please note, the contrast I am requesting is not a simple linear or polynomial contrast but is a contrast derived by a theoretical prediction, i.e. the type of contrasts discussed by Rosenthal and Rosnow.


title: How to refer to AIC model-averaged parameters and confidence intervals

text: I am writing up results from regression analysis where I used AICc model averaging to arrive at my final parameter estimates. I am wondering how best to refer to these parameters and their 95% confidence intervals. It seems like "significantly different" is taboo in the AIC world, but writing out "the parameter was x.x and its CI does not cross zero" seems much more laborious to me and the reader than saying "x.x was significantly different from zero." This seems like it might be an issue that would not come up if I had just selected the lowest AICc as my best model, which is what many folks do (against Burnham and Anderson repeatedly stating otherwise). Selecting the best model let's you say "the parameter is important b/c it is in the final model." Also, I'm wondering if there is an AIC model averaged equivalent to "marginally significant." I have parameters that have the predicted sign, indicate a fairly sizeable effect, but whose CI creeps over 0.0. Philosophically I like model averaging, and I also have many good models that often only differ by an extra covariate or an interaction. **EDIT:** This inquiry can probably be summarized by asking "In an AICc model averaging framework how does one interpret parameters whose confidence intervals span zero by only a small amount?"


title: How to transform an outcome which censored from both sides?

text: I'm trying to fit a `tobit` model for a response outcome where both sides were censored and the distribution is heavily skewed, even for the none-censored observations. Here is an example:               set.seed(123)     x<-10:100     a<-rnorm(length(x), sd=5)     b<-(x+a)^(-1/1.5)     y<-c(rep(min(b), 20),b,rep(max(b),100))     hist(y)      To fit a `tobit` model, `y` is assumed to be normal distributed. My question is how to transform this `y` so as to be normal-like distributed. I believe the below model is incorrect without a proper transformation of `y`.               age<-abs(as.integer(50+rnorm(length(y), sd=2)))     require(censReg)     fit<-censReg(y~age, left=min(y), right=max(y))     summary(fit)


title: Question about PCA data recovery equation

text: In PCA , consider a 4 x 3 data matrix ( 4 examples each with 3 features ). After getting the 3 eigenvectors (a/b/c) and projecting data on the first 2 vectors, the equation looks like this : [ first 2 eigen vectors transposed ] x [ data set transposed ] = [ data projected ] so the first example in the data after being projected would consist of: $$[a_1x_1 + a_2x_2 + a_3x_3]\\ [b_1x_1 + b_2x_2 + b_3x_3] $$ On recovering data the equation looks like this: [first 2 eigen vectors] x [projected data] = [original data set] so that the first example consists of: $$ a_1[a_1x_1 + a_2x_2 + a_3x_3] + b_1[b_1x_1 + b_2x_2 + b_3x_3]\\ a_2[a_1x_1 + a_2x_2 + a_3x_3] + b_2[b_1x_1 + b_2x_2 + b_3x_3]\\ a_3[a_1x_1 + a_2x_2 + a_3x_3] + b_3[b_1x_1 + b_2x_2 + b_3x_3] $$ My question is how is this equivalent to what the original example looked like $[x_1 x_2 x_3]'$? I know that unit vectors will cancel each other but the equation will be like $$ a_1[a_1x_1] + b_1[b_1x_1]\\ a_2[a_2x_2] + b_2[b_2x_2]\\ a_3[a_3x_3] + b_3[b_3x_3] $$ giving : $$ x_1 + x_1\\ x_2 + x_2\\ x_3 + x_3 $$


title: Interpret the effect of factor in GAMM

text: this will probably look like a trivial question but I am having trouble explaining the difference between the output plots of two models I built. My data are disease counts on algae. I would like to test the relationship between disease counts and percentage cover (continuous covariable) including the effect of Site (fixed categorical) and Zone (random, nested within Site).               mod1<-gamm(Nblesions~s(Cover),random =list(Site=~1,Zone=~1),family=poisson,data=df)      with the resulting plot: ![enter image description here](http://i.stack.imgur.com/U4RtO.png)               mod2<-gamm(NbLesions~s(Cover)+Site,random=list(Site=~1,Zone=~1),     family=poisson,data=df)      with the resulting plot: ![enter image description here](http://i.stack.imgur.com/6Tuby.png) The difference between the two models is the absence of factor Site in mod2. In both models there is NO significant effect of cover. Could anyone explain to me what is happening between the two and how to interpret the change in the relationship when Site is added? Any advice would be highly appreciated, thanks!


title: Evolution strategies in libsvm

text: I'm working on protein multi-classification problem. I'm using libsvm and the edit distance kernel. This kernel depends from a parameter (gamma). I'm able to get the best parameters (gamma and C) through grid search. If i use a kernel that depends from 3 or more parameters, the grid search is computationally heavy, so i'm thinking to approach the problem with an evolution strategy (for instance CMA-ES). Is there a way to do it in libsvm?


title: Appropriate test for difference in trials with varying calibration

text: Let's say I collect data from an experiment where for each measurement device, 5 trials are conducted at different temperatures. Each device is calibrated separately, and the 5 trial temperatures are constant across devices. In other words, I have values at Temperature 1, 2, 3, 4, and 5 from each device. **What is an appropriate test to see whether the values from a given measurement device, adjusted for calibration, are significantly different from the other measurement devices?** My guess is that the test must be conducted on the differences between the trial values to compensate for the different calibrations.


title: How would you fit ARIMA model with lots of autocorrelations?

text: This is really puzzling... I have this data which has a lot autocorrelations... The data is about 60000 data points of 15min data. I tried fitting it to ARIMA(6, 0, 6) and even GARCH(1, 1) with mean model ARMA(6, 6), still there are lots of autocorrelations in the residuals. ![enter image description here](http://i.stack.imgur.com/HqAHJ.jpg) I almost wanted to try ARIMA(100,0,100), but I think even that is not enough... I am doing these in R. How to I get out this swamp? Please shed some lights on me. Thanks so much Update: I have fitted a gjrGARCH model with distribution "sstd" in R using "rugarch" package... QQ-Plot of Standardized Residuals looks great - it's a straight line... However the ACF plots of the residuals and the residuals squared are a bit strange, especially the big spike in the residual squared ACF plot. Could anybody please shed some lights on us, esp about that big spike? Thank you! ![enter image description here](http://i.stack.imgur.com/MZxZ2.jpg) ![enter image description here](http://i.stack.imgur.com/KAMGE.jpg)


title: Cluster data into categories; train one classifier per category

text: Let me first give you a hypothetical example to make things more clear: Let's say your task is to classify art as either professional or amateur work, based on image data. You extracted 100 features from each artwork using image processing. The artworks fall into different genres (like sketch, abstract paintings, portraits, still-life, oil painting, watercolor painting etc.). These genres are responsible for most of the variance in the features. The genres (or even their number) are not known to you _a priori_. Indeed, there may be different ways to cluster this data into genres/categories (e.g. based on subject, painting technique, colors used). The class distribution (professional, amateur) varies per genre. For instance most amateur works in your dataset could be still-life. So if you would use a linear classifier it would basically just classify everything as still-life or not still-life. * * * Is there a model that is known to work well for a dataset like this? The obvious approach would be to throw a non-linear classifier (like kernel svm, decision forest) at it and hope that it can deal with the problem. However, that may work less well then a model that can make use of this property of your dataset explicitly. Another approach would be to cluster the data first, then train one classifier per cluster. However, that may lead to clusters that are not optimal for classification. So is there a model that is know to work well for problems like this on (bigish) real world datasets? Maybe modeling the problem probabilistically as $$P(\text{professional}) = \frac{P(\text{genre}) \times P(\text{professional}|\text{genre})} { P(\text{class} | \text{professional})}$$ and applying EM, starting from with a random class alignment. Or modeling it as a neural network, that has this structure hardwired and applying backpropagation.


title: Does the probability of rolling a '6' in dice change over time?

text: Suppose I have a die and I throw it a large number of times (a billion times) and have proven that it lands on each side an equal number of times so it is proven to be fair. Next suppose I started a second round of throwing of the same die and it lands on 1 through 5 say a 10,000 times each and the count for 6 is 100. Suppose someone ask you to bet on the next billion throws. Would you bet on 6?


title: Testing for statistically significant difference in time series?

text: I have the time series of the prices of two securities, A and B, over the same period of time and sampled at the same frequency. I would like to test whether there is any statistically significant difference over time between the two prices (my null hypothesis would be that the difference is null). Specifically, I am using price differences as a proxy for market efficiency. Imagine A and B are a security and its synthetic equivalent (i.e. both are claims to exactly the same cash flows). If the market is efficient, both should have the exact same price (barring different transaction costs, etc.), or a zero price difference. This is what I would like to test for. What is the best way to do so? I might have intuitively run a two-sided t-test on the "difference" time series, i.e. on the A-B time series, and tested for $\mu_0$=0. However, I have the suspicion that there might be more robust tests, that take into account things like potential homoskedastic errors or the presence of outliers. In general, are there things to watch out for when working with the prices of securities?


title: Correlation between dichotomous and continuous variable

text: I am trying to find the correlation between a dichotomous and a continuous variable. From my ground work on this I found that I have to use independent t-test and the precondition for it is that the distribution of the variable has to be normal. I performed Kolmogorov-Smirnov test for testing the normality and found that the continuous variable is non-normal and is skewed (for about 4,000 data points). I did the Kolmogorov-Smirnov test for the entire range of variables. Should I split them into groups and do the test? I.e., say if I have `risk level` (`0` = Not risky, `1` = Risky) and cholesterol levels, then should I:   * Divide them into two groups, like              Risk level =0 (Cholestrol level) -> Apply KS     Risk level =1 (Cholestrol level) -> Apply KS        * Take them together and apply the test? (I performed it on the whole dataset only.) After that, what test should I do if it is still non-normal? **EDIT:** The above scenario was just a description I tried to provide for my problem. I have a dataset which contains more than 1000 variables and about 4000 samples. They are either continuous or categorical in nature. My task is to predict a dichotomous variable based on these variables (maybe come up with a logistic regression model). So I thought the initial investigation would involve finding the correlation between dichotomous and a continuous variable. I was trying to see how the distribution of the variables are and hence tried to go to t-test. Here I found the normality as an issue. The Kolmogorov- Smirnov test gave a significance value of 0.00 in most of these variables. Should I assume normality here? The skewness and kurtosis of these variables also show that the data is skewed (>0) in almost all cases. As per the note given below I will investigate the point-biserial correlation further. But about the distribution of variables I am still unsure.


title: Spatial cluster analysis

text: Let's say I have a structure like this : ![enter image description here](http://i.stack.imgur.com/UID7o.png) This is a spatial region with measurement of plant population in each site. Black and red represent two regions with different intensities.The question is that I want to ask about the degree of spatial dependency. I've tried Moran's I statistic, but I think Moran's I can only tell us if the spread is spatially independent or not. Is there a way to test about spatial dependency but also to tell if we have different clusters on it ? Like in this example, I know that I have 2 different clusters (black and red). Thanks...


title: Are gridlines and grey backgrounds chartjunk and should they be used only on an exception basis?

text: It seems that most authorities agree that dark or otherwise prominent gridlines in plots are "chartjunk" by any reasonable definition and distract the viewer from the message in the main body of the chart. So I won't bother to give references on that point. Equally, we can all agree that there will be times that _pale_ gridlines to create a reference for viewers _will_ be necessary. Tufte argued the need for(and used) gridlines occasionally, as pointed out for example in this article. And I agree with the Hadley Wickham's approach in ggplot2 to make such gridlines white on a pale grey background, when you need to use them. What I am uncertain about however is whether such gridlines and grey background should be the _default_ , as they are in ggplot2. For example, there seems to be no reason for the grey background other than to case the white gridlines into relief - which further begs the question of whether either is needed. I have recently started using ggplot2 for most of my graphics needs and think it is awesome, but it has challenged my "no box, no background, no gridlines" approach to graphics that I used before. I used to think that gridelines=OFF should be my default unless there is a particular reason for adding them in - basically the approach recommended in this article, for example. Of course, it is straightforward to define a theme in ggplot2 to avoid the gridlines and background shading (and in fact we have done this at my work), but ggplot2's approach is so awesome and generally its default aesthetic choices are good that I wonder if I am missing something. So - I would be grateful for any references on this point. I am sure it has been well thought through (for example by Hadley Wickham in setting the ggplot2 defaults) and I am very open to being pointed in the right direction. The best I've been able to find is a couple of links on the ggplot2 google group but the most helpful looking reference by Cleveland isn't available at the link given.


title: Difference between sample variance and statistic variance

text: Conceptually we learn that the sample variance formula, we have to divide by n-1 such that it gives the correct estimation for the true variance. However if we have this example, let a sample $X_i = 1$ with probability $p$, and $0$ with probability $1-p$. so if we define a statistic $\hat{p} = \frac{1}{N} \sum X_i$, we find that $E(\hat{p}) = \frac{1}{N} N p = p$; and working from first principle we can find that variance of $\hat{p}$ $= \textrm{var}(\frac{1}{N} * \sum X_i)$ $= \frac{1}{N^2} Np(1-p)$ $= \frac{p(1-p)}{N}$ So my question is that is this variance obtained the same as the sample variance formula (since there is no divide by n-1 term)? If not which should be the correct answer? Thank you. Edit: Is it because the N-1 formula is talking about the sample variance, and my example we are actually looking at the variance of the sample mean (i.e. sample mean variance)?


title: How to interpret logistic regression with days of the week and months of the year as predictors together with continuous and binary predictors?

text: I am using logistic regression to predict electricity spike prices (price that exceeds a certain threshold). I directly use the following variables as my independent variables together: (a) demand/load - a continuous IV   (b) is a binary IV   (c) days of the week   (d) months of the year. I used the reference category for (c) and (d) and end up `sept` and `tue` for the reference. How do I interpret the results using marginal effects for all my predictors? I am using Stata 9.


title: Tried to overfit a Bayes net, but mean prediction error is worse than learned network?

text: I have variables A, B, C, D, and E. I am interested in building a classifier for A. I learned a Bayes net structure from the data using greedy search and BIC as a score. Call this network 1. Using cross validation, I got the mean prediction error for node A in network 1. I then, thought to create a network structure where all the arcs are incoming into A (B->A, C->A, D->A, E->A). Call this network 2. I also learned a third network that was constrained by having all the arcs in network 2, but could allow additional arcs between B, C, D, and E. Call this network 3. It seemed to me that the network 2 and network 3 would be better at predicting A because I was forcing all the info in the data as a direct input into predicting A -- i.e. overfitting. So I expected the mean prediction errors (MPE) for these networks to be less than the first network. But instead they were higher. Also, when I repeated the cross validation several times, the mean prediction error for network 1 was always the same. But the value varied for networks 2 and 3. Why would networks 2 and 3 have a higher MPE, and why would MPE for those networks be variable while not variable for network 1?


title: IRT how to graph the Item Information Function

text: Following the Stata Manual - Item–response theory (IRT) models' example 28g — One-parameter logistic IRT (Rasch) model: use http://www.stata-press.com/data/r13/gsem_cfa , clear gsem (MathAb -> (q1-q8)@b), logit var(MathAb@1) it is explained how to obtain the Item–characteristic curves graph the conditional probability of a particular response given the latent trait, as follows: predict pr*, pr predict ability, latent(MathAb) and thus we can obtain the item–characteristic curves for all eight questions by typing: twoway line pr1 pr2 pr3 pr4 pr5 pr6 pr7 pr8 ability, sort xlabel(-1.5(.5)1.5) or, a less busy graph might show merely the most difficult and least difficult questions: twoway line pr2 pr3 ability, sort xlabel(-1.5(.5)1.5) But, my question, how can we obtain the data to draw the so-called Item Information Function for each item as well as for the Latent Trait(s)? If anyone knows how to do this, much appreciated!


title: Fit of negative binomial regression model

text: I have ran a negative binomial regression. I'm guessing the use of a negative binomial regression is not ideal given my design, but I'm hoping I can 'get away with it', as it seems to be working fairly well. Here's my design:   * Count the the total number of geese in a flock   * Count the number of these that are vigilant.    * Then divide number vigilant/number in flock, to get a percentage that are vigilant    * Predictor variables are flock size (count), water height in metres, date, subspecific identity of the flock and the part of my study site that the flock were in (ssp.zone) Below is a plot of my response variable (`alert.integer`) ![enter image description here](http://i.stack.imgur.com/RbEJU.jpg) Below is a summary of the model               Call:     glm.nb(formula = alert.integer ~ ssp * water.height + ssp * ssp.zone +          ssp * count + ssp * I(as.Date(datetime)), data = dataScanSampling_sub,          init.theta = 1.168496598, link = log)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.0035  -0.9556  -0.3382   0.3886   2.0942            Coefficients:                                                         Estimate Std. Error z value Pr(>|z|)         (Intercept)                                        3.971e+02  6.239e+01   6.365 1.95e-10 ***     sspLight-bellied Brent Goose                      -4.369e+02  7.294e+01  -5.990 2.10e-09 ***     water.height                                       2.978e-01  8.719e-02   3.415 0.000637 ***     ssp.zonelight zone                                -1.214e+00  2.103e-01  -5.775 7.70e-09 ***     count                                             -5.226e-03  8.506e-04  -6.143 8.09e-10 ***     I(as.Date(datetime))                              -2.461e-02  3.904e-03  -6.303 2.93e-10 ***     sspLight-bellied Brent Goose:water.height          2.743e-03  1.015e-01   0.027 0.978438         sspLight-bellied Brent Goose:ssp.zonelight zone    1.344e+00  2.635e-01   5.102 3.36e-07 ***     sspLight-bellied Brent Goose:count                 4.174e-03  8.766e-04   4.761 1.93e-06 ***     sspLight-bellied Brent Goose:I(as.Date(datetime))  2.726e-02  4.564e-03   5.973 2.33e-09 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Negative Binomial(1.1685) family taken to be 1)              Null deviance: 558.21  on 367  degrees of freedom     Residual deviance: 433.87  on 358  degrees of freedom     AIC: 3236.4          Number of Fisher Scoring iterations: 1                             Theta:  1.1685                Std. Err.:  0.0894            2 x log-likelihood:  -3214.4300       And here are residuals from model: ![enter image description here](http://i.stack.imgur.com/hnBEW.jpg) So taking into account the description of my design, the plot of the response variable (`alert.integer`), the residuals plot and the output from `summary`(), can I get away with a negative binomial regression model here?


title: Inverse of arcsin

text: So this might be not the smartest question. But I defined this function:               trans.arcsine <- function(x){                                                              asin(sign(x) * sqrt(abs(x)))                                           }       which I got from here. My question is how to invert it? **Edit:** My domain of interest is [0,1], which was not stated in the beginning. @Glen, however, provided a very good answer for [-1,1] as well.


title: How to algorithmically determine the best order of fit?

text: I am doing a least squares polynomial interpolation for 10,000 data sets that look mostly like one period of a sine curve, but whose values are not evenly spaced in the time domain, and can sometimes be quite noisy. A lot of them do just fine with a 3rd order fit (a + bx + cx^2 + dx^3) but some oscillate wildly (Runge's phenomenon) and need a 2nd order instead, and some are much better fitted with a 5th or even 8th order fit. What methods exist for determining which fit is best with respect to avoiding both underfitting and overfitting? One thing I think that I can do is to choose the largest order $m$ such that $m<2\sqrt{n}$ where $n$ is the number of points in my data set, and also check to make sure that there are some number of points in $2\sqrt{n}$ equally sized bins spanning the space. On the other hand, however, I would also like to penalize for having more parameters than necessary. Maybe I can also minimize $SSE(m)/(n-m-1)$. What better methods are out there for doing this?


title: Estimating how many individuals are infected

text: Suppose you are simulating a population of 1000 individuals, each of which have a probability p of being infected with some condition. How do you determine how many individuals are infected in a individual round. To get a number for each round you could do: total = 0; for i = 1 to 1000 { if random number in [0,1] < p then total = total + 1 } It feels like there's a simple way to do this without the for loop, but I can't figure it out...


title: Cox Regression Prediction of Event Count For Given Period

text: I have a standard event history data set, which I used to fit a Cox proportional hazards model in Stata and R. The time variable is measured in days. There are no time-varying covariates on the RHS of the equation. I have recovered the usual output (coefficients, hazard ratios, SEs, etc.) as well as the predicted baseline survival for each observation. What I need to do now is estimate the predicted number of events in the next 30 days. Any ideas on how to do this? Thanks very much!


title: Moving average filter for outlier removal

text: I am using a moving average filter to smooth data for outlier removal. By changing the number of average points, I am getting different result. My data are multi-dimensional feature vectors. I applied the moving average to the entire matrix and then on individual variables. They give different results. So, how to choose/guess the number of points to average over and should it be applied on the entire matrix or on a one by one basis?


title: Should I report credible intervals instead of confidence intervals?

text: After stumbling upon the concept in a statistics textbook, I tried to wrap my head about it, and finally came to a conclusion which seems to fit all the explanations which I have seen so far: A credible interval is what non- statisticians think a confidence interval is. * * * _Digression for those like me-from-an-hour-ago who don't know the difference_ If we observed data and predicted some parameter from it, let's say the mean $\mu$, the credible interval is the interval $[\mu_{\text{min}},\ \mu_{\text{max}}]$ for which we are 95% sure that mu falls inside (or some number other than 95%, if we used another level). The confidence interval taught in introductory statistics classes _can_ overlap with the credible interval, but will not _always_ overlap well. If you want to brave the explanation, try reading this and this question on Cross Validated; what helped me finally understand, after much head-scratching, was this answer. * * * Does it mean that it would be scientifically preferable to use a credible interval over a confidence interval in my results? If yes, why haven't I seen any publications which use it?   * Is it because the concept _should_ be used, but the measuring scientists have not yet caught up with the correct statistic methods?    * Or is the meaning of the original confidence interval better suited to explaining results from empirical studies?   * Or is it that in practice, they are so often overlapping that it doesn't matter at all?    * Does the choice depend on the statistical distribution we are assuming for our data? Maybe with a Gaussian distribution, they always overlap numerically, so nobody outside of pure statistic cares about the difference (many studies I have read don't even bother to calculate _any_ kind of interval, and maybe around 1% ever give space to the thought that their data might not be normally distributed).    * Does it depend on our scientific-theory standing? For example, it feels like the confidence interval should be used in positivist work and the credible interval in interpretivist work, but I am not sure that this feeling is correct.


title: Price elasticity of demand and time lags

text: This is probably a very basic question so apologies in advance, but we are stuck. We currently change prices by varying amounts (both positively and negatively)every 3 months,and would like to calculate the price elasticity of demand. Questions: 1) Should we eliminate a period of time immediately after the price change to allow time for people to adjust behavior? 2)Should we include the entire time between quarters (i.e. the three months before the price change (so everything right after the previous price change all the way up to the current price change) and the 3 months after the price change (so everything after the current price change all the way up to the next price change)? Any help would be greatly appreciated!!!!


title: Minimizing the Training data

text: I have a grey-box model of the form Y= a + b X1 + c X2. Where a, b and c are the coefficients based on regression. The regression variables X1 and X2 are determined based on experimental/simulation data. I can experimentally run for a large number of cases and train the model (but it is time consuming) or try to figure out the optimum training data required and only run for those particular experiments. So, my objective here is to reduce the training data as much as possible without compromising much on the predictive capability of the model. I am very new to this statistics field and I have no idea as how this could be done. So can anyone please give me suggestions about this.?


title: How adapt MCMC when (constant) weights for oberservations are introduced?

text: I have the following problem: I have already set up a model and MCMC sampler for a mixed model without weights, i.e., every observation contributes the same amount of information. Now I would like to implement weights, thus I change the likelihood to $f(y|\theta) = |2\pi R\sigma_e^2|^{-1/2} \times \text{exp}\left(-\frac{1}{2\sigma_e^2} (y - X\beta - Z\alpha)^T R^{-1} (y - X\beta - Z\alpha) \right)$, where the matrix $R$ contains constant, known weights. For example, for an observation that would be effectively the mean of two previous observations, I enter 1/2 at the corresponding diagonal position in $R$, and so on. Is it possible to use the same sampler without weights by just doint a transformation of my data like $w = \sqrt{R^{-1}}$, $y^* = w y$, $X^* = w X$, $Z^* = w Z$ So we can rewrite the likelihood using the transforms $f(y|\theta) = |2\pi R\sigma_e^2|^{-1/2} \times \text{exp}\left(-\frac{1}{2\sigma_e^2} (y^* - X^*\beta - Z^*\alpha)^T (y^* - X^*\beta - Z^*\alpha) \right) $ And since $R$ is constant $f(y|\theta) \propto (\sigma_e^2)^{-k/2} \times \text{exp}\left(-\frac{1}{2\sigma_e^2} (y^* - X^*\beta - Z^*\alpha)^T (y^* - X^*\beta - Z^*\alpha) \right) $ It seems straightforwards, but also a little bit too simple to be correct. Is it possible to do it like that?


title: Why is an ARMA model a parsimonous approximation of an AR model?

text: I am reading a book on time series and I came across the following: "In addition to being a parsimonous approximation to a high-order AR(p) model, ARMA models...". Why is an ARMA model a (parsimonous) approximation to an AR model? Also, somewhat later on I read that the coefficients of an MA model can be estimated by writing the model in an 'AR form'. Could anyone also please explain that statement to me? Thanks in advance.


title: A statistic to describe a correlation with progression

text: I want to see if the `proportion of 'I don't know'` responses correlates with the progression through the survey. Generally we expect this proportion to increase through the survey due to a fatigue effect. However, the following chart suggests this effect is not present in my data: ![enter image description here](http://i.stack.imgur.com/Y2HQs.png) I would however like to back this up with a statistic. I figure I need a `correlation` measure between the question order and proportion of DK. Does anyone know about such measure, and maybe how to produce it in `SPSS` (R could also be an option)? PS: I eventually also want to control for question and respondent characteristics.


title: Algorithm for winning casino roulette

text: I would like to try the following algorithm in order to win in the roulette:   1. Be an observer until there are 3 same parity numbers in a row ($0$ has no defined parity in this context)   2. Once there were achieved 3 numbers with the same parity in a row: init $factor$ to be $1$   3. Bet the next number will be the opposite parity if you were wrong double $factor$ and return `3` else `goto` `1`. Here is a python code for it - which bring profit all the time:               import random          nums = range(37)     bank = 10**10     games = 3000     factor=1     bet_amount=100     next_bet = None     parity_in_row ={"odd":0, "even":0}          for i in xrange(games):         num = nums[random.randint(0,36)]              if next_bet == "odd":             if num % 2 == 1:                 bank += factor*bet_amount                 factor = 1                 next_bet = None             else:                 bank -= factor*bet_amount                 factor *= 2         elif next_bet == "even":             if num % 2 == 0:                 bank += factor*bet_amount                 factor = 1                 next_bet = None             else:                 bank -= factor*bet_amount                 factor *= 2                   if num > 0 and num % 2 == 0:             parity_in_row["even"] += 1             parity_in_row["odd"] = 0         elif num % 2 == 1:             parity_in_row["odd"] += 1             parity_in_row["even"] = 0         else:             parity_in_row ={"odd":0, "even":0}              if parity_in_row["odd"] > 2:             next_bet = "even"         elif parity_in_row["even"] > 2:             next_bet = "odd"         else:             next_bet = None          print bank        * If I did the calculation correctly the probability to have 4 numbers in a row with the same parity $ < (1/2-\epsilon)^4 $ [where $\epsilon < 1/36$ is a compensation for the probability to achieve $0$].   * Keep doubling $factor$ ensures you will have positive expectation, right?   * Isn't the probability of getting odd number is $0.5$ ? since $Prob(odd | even, even, even) = \frac{1}{2}$ ? **Please try to supply a rigorous proof why it is wrong, and try to explain why my python code always return with a positive outcomes**


title: Is there any statistical reason for diagonal lines in scatterplot on a log scale?

text: I'm perplexed by some vertical lines that show up in these scatter plots on a log scale. Population is on the y-axis and the proportion of the neighborhood with the attribute mentioned in the panel label on the x-axis. Is this just an artifact of the transformation? (I thought perhaps this was due to digit preference or something like that, but I can't think how this might be possible in this dataset from census data.) ![enter image description here](http://i.stack.imgur.com/I4sao.png)


title: Simultaneous heteroscedasticity and heavy tails in a regression model

text: I'm trying to create a prediction model using regression. This is the diagnostic plot for the model that I get from using lm() in R: ![diagnostic plots from R](http://i.stack.imgur.com/JRNM0.png) What I read from the Q-Q plot is that the residuals have a heavy-tailed distribution, and the Residuals vs Fitted plot seems to suggest that the variance of the residuals is not constant. I can tame the heavy tails of the residuals by using a robust model:               fitRobust = rlm(formula, method = "MM", data = myData)      But that's where things come to a stop. The robust model weighs several points 0. After I remove those points, this is how the residuals and the fitted values of the robust model look like: ![Residuals vs Fitted for the robust model](http://i.stack.imgur.com/4RGSu.png) The heteroscedasticity seems to be still there. Using               logtrans(model, alpha)       from the MASS package, I tried to find an $\alpha$ such that               rlm(formula, method = "MM")       with formula being $\log(Y + \alpha) \sim X_1+\cdots+X_n$ has residuals with constant variance. Once I find the $\alpha$, the resulting robust model obtained for the above formula has the following Residuals vs Fitted plot: ![Residuals vs Fitted for log-transformed response](http://i.stack.imgur.com/Lg2WB.png) It looks to me as if the residuals still do not have constant variance. I've tried other transformations of response (including Box-Cox), but they don't seem like an improvement either. I am not even sure that the second stage of what I'm doing (i.e. finding a transformation of the response in a robust model) is supported by any theory. I'd very much appreciate any comments, thoughts, or suggestions.


title: Alternative to SAS Enterprise Guide

text: Is there any freeware or commercial alternative to SAS EnterpriseGuide? I'm looking mostly for it's ability to save complex structured queries so the statistical libraries are of secondary importance.


title: Decision boundary plot for a perceptron

text: I am trying to plot the decision boundary of a perceptron algorithm and am really confused about a few things. My input instances are in the form [(x1,x2),target_Value], basically a 2-d input instance and a 2 class target_value [1 or 0]. My weight vector hence is in the form: [w1,w2] Now I have to incorporate an additional bias parameter w0 and hence my weight vector becomes a 3x1 vector? is it 1x3 vector? I think it should be 1x3 since a vector has only 1 row and n columns. Now let's say I instantiate [w0,w1,w2] to random values, how would I plot the decision boundary for this? Meaning what does w0 signify here? Is w0/norm(w) the distance of the decision region from the origin? If so how do I capture this and plot it in python using matplotlib.pyplot or its matlab equivalent? I would really appreciate even a little help regarding this matter.


title: How to determine if survey response samples are comparable?

text: I'm looking for a method to determine if the number of survey response received from a regional office are substantial enough to compare to each other year over year. For example, one year we might get 29 survey responses in a month from one site but only 5 responses for the same month the next year. Can I say those are truly comparable? I guess what I am thinking of is similar to a Gallup poll or something like that, where they test to see if their response make up approximately the same mix of people between two time frames. The reason we are trying to do this is to avoid a survey bias, as we know some regional sites consistently perform better than others, so if our sample contains more responses from the better performing sites, our results will be skewed favorably when they are really just due to having a better response rate (more observations) from certain cites.


title: Are posterior probabilities from a Naive Bayes classifier reliable?

text: I have read that the posterior probabilities of Naive Bayes classifiers are unreliable. Is this true? and if so, in what sense, and why? Specifically, I am interested to know if the probabilities can be used as some measure of confidence of the predictions. For example, can I say that higher posteriors are more likely to be correct, or use a threshold to say that if the posterior is below some value then it is low-confidence.


title: Estimating a first order plus dead time model

text: The data generating process is given by the following differential equation:   $y(t) = a + b u(t - \theta) + c \frac {dy} {dt}$   Now imagine having as data a long time series for both $y$ and $u$. If $\theta$ is known, this is an easy linear regression. But is there a "correct" procedure when the delay $\theta$ is unknown? Right now I treat it as a model-selection problem, running multiple regressions and choosing the one that predicts with the smallest error. That is clunky though, so if a better method exists, great! To clarify:   * $u(\cdot)$ is a function. It is an arbitrary stream of inputs   * $a$,$b$,$c$ and $\theta$ are parameters. Fixed and unknown.   * You can assume a gaussian $\epsilon$ measurement error whenever $y(t)$ is observed.   * The only data I have is a long series of daily observations of the pair $(u,y)$, where $u$ is observed perfectly.    * I feed $y_{t}-y_{t-1}$ for the derivative, which is appropriate since the original model is really a difference equation.


title: How can you tell if samples are from the same population

text: Given the 7 different samples with mean and variance for each group. What can I use to tell if these samples from the same population with 0.05 and 0.01 significance level?               Sample   Mean   Variance        A        50      5     B        47      4     C        52      5     D        51      6     E        49      5     F        52      4     G        49      6      Based on the give data, I calculated F statistics of value 53.72. Given the the F data, how would you conclude that these samples are drawn from the same population? Is the conclusion the same for the both for 0.05 and 0.01 levels of significant? I am trying to follow this example: http://en.wikipedia.org/wiki/F-test#One- way_ANOVA_example but the problem is that I don't have the individual observation to calculate the within-group sum squares. I only have the mean, variance and sample size(10) data. Any ideas how I could approach this?


title: Effect sizes from regression analysis and factor loadings

text: To do a meta analysis, I try to calculate effect sizes. In an article by Shah and Ward, the authors do a regression analysis after combining various factors by a CPA to four factors. The factor loadings are provided in the article. The beta coefficients for these four factors regarding operational performance (a factor of dependent variables) are given too. Is it possible to calculate effect sizes for the independent variables with this information? Thanks!


title: LASSO to identify important variables in ordered logistic regression?

text: I have spent two days grappling with this question, and the range of ambiguous answers online has driven me to ask. I am working with R. I have a dataset where my dependent variable is an ordered categorical variable with 5 levels ("dislike very much" to "like very much"). I intend to use the `ordinal` package for the actual regression analysis, but have been trying to decide on the best method to specify the model. From my reading it appears that stepwise methods are not a good option, and that a LASSO regression technique is a better method of selecting important variables. Is it acceptable to use a LASSO method to choose the variables to include and then to use these variables in a separate proportional odds regression? My primary interest is in the significance of the terms, rather than the size of the coefficient. I want to know which variables have a significant effect, in which direction, and in what order of importance. For this reason I would rather do the final modelling using a proportional odds glm than with a LASSO model. Many thanks in advance.


title: 2x2 ANOVA - assess violations of homoscedasticity & normality

text: I have a 2x2 factorial unbalanced between-subject design, **_n_ = 355**. My DV is a subjective probability estimate (i.e., a number between 0 and 100). **My ANOVA model:**               aov1 <- aov(prob ~ utility * imagination, data = e3nndata)      To **check for homoscedasticity** I used Levene's Test, which was highly **significant** :               library(car)     leveneTest(prob ~ utility * imagination, data = e3nndata)          >        Df F value    Pr(>F)         > group   3   7.639 5.838e-05 ***     >       351                            To **check for normality of residuals** I used the Shapiro-Wilk-Test, which was equally **significant** :               shapiro.test(aov1$residuals)          > W = 0.897, p-value = 9.099e-15      In order to tackle both problems, I decided to **log-transform the DV** :               aov1.log <- aov(log(prob) ~ utility * imagination, data = e3nndata)      In fact, now Levene's test is **no longer significant** :               leveneTest(log(prob) ~ utility * imagination, data = e3nndata)          >        Df F value Pr(>F)     > group   3  0.6385 0.5907     >       351                     However, the Shapiro-Wilk-Test on the residuals **remains significant** :               shapiro.test(aov1.log$residuals)          > W = 0.9661, p-value = 2.415e-07      **Is it correct, though, to run the Shapiro-Wilk-Test on the entire sample or should it be calculated for each group individually?** From posts like this one I learnt that with larger sample sizes significant Shapiro-Wilk results are often rather meaningless and more focus should be put on graphical representations. Therefore I had a look at **two diagnostic plots**.               library(car)     qqPlot(aov1.log)      ![qq-plot of residuals from anova model with log-transformed DV](http://i.stack.imgur.com/ckCvo.png)               library(ggplot2)     ggplot() + geom_histogram(aes(x=aov1.log$residuals), binwidth=0.25)      ![histogram of residuals from anova model with log-transformed DV](http://i.stack.imgur.com/kUaGH.png) I don't have any experience with violations of model assumptions, so I cannot tell whether the two plots suggest such a strong violation of the normality assumption that ANOVA is no longer a suitable method. ## Questions:   1. What is the right conclusion from these plots?    2. Are there other plots/statistics to take into account?   3. As a consequence, how should the data be analyzed?


title: A/B Test with factors other than conversion rate

text: I've read "AB Testing other factors besides conversion rate", but I don't know if it's relevant with my question, even if it seems similar. Here is my problem. Let's assume I am creating a new tools designed to help user do their task. I A/B test this tools, where some user get the tools, and other doesn't. I want to know if the tools really helping user with the task. My null hypothesis: users without tools finish their task as fast as than users without tools. **1. Is my hypothesis valid?** Can we statisitically prove this? Should I use no less faster instead of slower? Assuming valid (or it's going the right way), **2. What is the data I need to prove this tools help users?** While on experiment, the user do the task more than one time. I don't know what the effect of that in the experiment, but user with tools will consistently have ability to use the tools, and user without will consistently not using the tools. Let's say I assume the time users needed to complete experiment is what I needed to collect. I have the data for each user, number of task they did, and how long they finish the task. The task is identical for all user in the experiment. And then I gathered following data:                            | sample |  (Average time / User) |     without tools|  184   |          17            |     with tools   |  176   |          12            |      **2. I this kind of data correct?** I mean, am I going the right way? Should I averages it? If not, how should I give the data? Should I only use the first time the user use it? But I want to prove the tools really helping people, so the first time uses might made people need time to learn. From the data, statistically intuitive, the tools helped users, because users can finish their work averagely faster, and there is smaller action to do. **3. How to prove it mathematically?** Is it even my intuitive correct? We can calculate the standar deviation for each set (with tools and without), and then calculate their standar error. From what I've search, we can calculate the standar error for each data as               se(withouttools-withtools) = sqrt(sqr(se(withtools)) + sqr(se(withouttools)) )      Then we can calculate the z-score as               z-score = (mean(withouttools) - mean(withtools)) / se(withouttools-withtools)      Based on the data-distribution, we can compare z-score with significance value with p = 0.05. Null hypothesis states average time of people withouttools will be as fast as with tools (no effect of the tools). We are only interested whether using tools will be faster than without, so we are using one-tailed test. That means if z-score of data > z-score of significant point, i.e. if z-score > 1.645, The null hypothesis can be rejected confidently, that is, using tools will make people finish their task faster than without. Here I assumed the data is normal distribution, because of Central Limit Theorem. But I don't know if it's valid assumption. **4. Am I calculating it correctly?** How should I actually calculate it? **5. If I'm still not confident, or have not enough data, how much data I need? How long I should run the test?** Thank you.


title: Cluster robust standard errors in plm package

text: I am a beginner in R. I am using the plm package and I want to know the easiest way to get cluster robust standard errors for the fixed effects estimator. My data (called Mydata) is on poverty, GDP, investment, and education in 130 countries over 33 years, and my index in plm is country and time. My model is something like this:               Fixedeffect<-plm(Y~ X_(i,t) + W_(i,t) + Z_(i,t) , data=Mydata, model= "within")      Where Y, X, W, and Z are all panel data variables. I read this post: Standard error clustering in R (either manually or in plm) , but it didn't help. Someone there suggested to use this:               coeftest(model.plm,vcov=vcovHC(model.plm,type="HC0",cluster="group"))      but it didn't work for me. I got the following error:               could not find function "coeftest"      Can anyone help, Please?


title: Complete machine learning library for Java/Scala

text: Python is plenty of ML libraries (like the great scikit-learn). Is there any good for java/scala, containing many algos (regression, classification, clustering, cross-validation, feature processing), stable & maintained and able to deal with massive dataset? I've just found Mahout, Breeze/Nak, and Weka, but they're not looking as great as Python ones. Additionally, if there's no equivalent, how can I efficiently connect java code with Python?


title: Telecommunication Billing DataSet

text: Has anyone seen any telco billing dataset available for analysis? I looked at Kaggle or http://www.d4d.orange.com/en/presentation/data but couldn't find anything.


title: What are speed differences beetwen ML implementations in different languages?

text: I am trying to write my own ML library. For speed reasons I started out writing things in C using BLAS, but then I learned that NumPy and Theano also use BLAS. I am wondering if there are huge speed differences between implementations of ML algorithms in C/Python/Matlab/Octave. Does anybody have some experience or can provide some data for the comparison? If there is no really good reason to write in pure C, I would rather not.


title: Logistic regression with an log transformed variable, how to determine economic significance

text: I am using a logistic regression model with continuous independent variables and two log transformed size variables (total assets and total deposits). My question is how to interpret the results and measure the economic impact? The normal steps in interpretation I take is as follows: I firstly run stata test on summary statistics (sum), then a logistic regression (logit) and than the marginal effects for the economic significance (mfx). These results can be explained by sentences as "a one standard deviation increase in $X$ (an independent variable) increases the probability of participation by $X\%$ I will multiply the $dy/dx$ of the marginal effects by the SD to get the % impact. So far so good, but how do I interpreted the results with the two ln(size) variables? To clarify my problem I will show you a part of the results:               sum stat:                       obs       mean       stdD        min         max     lnTCD|      2755   -1.469624    2.992985  -15.01948   4.827369          lnTA |      2767   -1.131045    2.796627  -16.51857   6.731258          logit res:                     Coeff       stERR       Z      P>Z             95%     lnTCD|   1.088141   .1424768     0.65   0.519     .8418451    1.406494     lnTA |   .7880959    .105272    -1.78   0.075     .6065653    1.023954          Margin eff:                      dy/dx       stERR      Z      P>Z           95%          X     lnTCD|    .0104326      .01619    0.64   0.519  -.021308  .042174  -1.09399     lnTA |   -.0294112       .0165   -1.78   0.075  -.061749  .002926   -.7886      In a normal situation I would multiply lnTA of the marginal effects (-0.0294112) by the Standard deviation of the summary statistic (2.796627) which results in -8.22% Although this holds for the other non log variables, intuitively this does not sound correct for these (ln) variables.


title: poisson vs logistic regression

text: I have a cohort of patients with different length of follow-up. So far I´m disregarding the time aspect and just need to model a binary outcome- disease/no disease. I usually do logistic regression in these studies, but another collegue of mine asked if Poisson regression would be just as appropriate. I´m not that into poisson and was left uncertain as to what the benefits and disadvantages of doing poisson in this setting would be compared logistic regression. I read Poisson regression to estimate relative risk for binary outcomes and I am still uncertain as to the merits of poisson regression in this situation.


title: Least Squares Estimators' Forms

text: Given a regression model of the form: $$y_{i} = \beta_{0} + \beta_{1}x_{i}$$ The least squares estimators for $\beta_{1}$ and $\beta_{0}$ are given by: $$\hat{\beta}_{1} = \frac{\sum x_{i}y_{i} - n\bar{x}\bar{y}}{\sum x_{i}^{2} - n \bar{x}^{2}}$$ and $$\hat{\beta}_{0} = \bar{y} -\hat{\beta}_{1}\bar{x}$$ In addition, $\hat{\beta}_{1}$ and $\hat{\beta}_{0}$ can be expressed in the form: $$ \hat{\beta}_{1} = \sum_{i=1}^{n}k_{i}y_{i}$$ and $$ \hat{\beta}_{0} = \sum_{i=1}^{n}k_{i}'y_{i}$$ where $k_{i} = \frac{(x_{i}-\bar{x})}{S_{XX}}$, $k_{i}' = \frac{1}{n} - \frac{(x_{i}-\bar{x})\bar{x}}{S_{XX}}$ and $S_{XX} = \sum x_{i}(x_{i}-\bar{x})$. My question is: how do you get from the first set of equations to the second set (functions of $k$ and $k'$) of equations? Although this may look very like a homework exercise, I can assure you that it's not. **EDIT:** I should probably provide a bit of context here. Properties such as the expected value and variance of the estimators appear to be more readily obtainable in the forms above containing the $k$ and $k'$ terms. In the event that I want to obtain such properties from alternative regression models, it would seem that it is advantageous to manipulate the estimators' equations until they are in a form comparable to the second set of equations given above.


title: Operations research versus statistical analysis?

text: What is the difference between operations research and statistical analysis?


title: Is it possible to use machine learning as a method for learning stats, rather than vice-versa?

text: During every machine learning tutorial you'll find, there is the common "You will need to know x amount of stats before starting this tutorial". As such, using your knowledge of stats, you will learn about machine learning. My question is whether this can be reversed. Can a computer science student learn statistics through studying machine learning algorithms? Has this been tested, at all? Are there examples where this is the case already?


title: Are Arma models on levels and first differences related?

text: Imagine I find a good fit of an `Arma 4,4` model on the levels, and on the first difference I find `Arma 3,4`. Are both models related somehow?


title: Forecasting estimate Demand

text: I have the following table               t    Demand    Level    Trend    Forecast     .     .     .     8    121,63    119,04    12.52    118,4     9                                 131.56      I want to calculate the Demand, Level and Trend for t=9 with the following formulas $$ L_{t} = \alpha D_{t} + (1-\alpha)(L_{t-1}-T_{t-1})$$ $$ T_{t} = \beta(L_{t}-L_{t-1})+(1-\beta)T_{t-1} $$ $$ F_{t+1} = L_{t} + T_{t} $$ $$ \alpha = 0.2, \beta = 0.4 $$ but I have no formula to calculate $D_{t}$


title: ARIMA model identification

text: Having a bit of difficulty identifying the appropriate ARIMA model by looking at ACF/PACFs. I know that AR(1) models, the ACF has a geometric progression from its highest value at lag 1 and the PACF has a spike at lag 1 and then then cuts off afterwards. Could someone possibly explain how I identify the rest? AR and MA please?!?


title: Normal distribution and monotonic transformations

text: I've heard that a lot of quantities that occur in nature are normally distributed. This is typically justified using the central limit theorem, which says that when you average a large number of iid random variables, you get a normal distribution. So, for instance, a trait that is determined by the additive effect of a large number of genes may be approximately normally distributed since the gene values may behave roughly like iid random variables. Now, what confuses me is that the property of being normally distributed is clearly not invariant under monotonic transformations. So, if there are two ways of measuring something that are related by a monotonic transformation, they are unlikely to both be normally distributed (unless that monotonic transformation is linear). For instance, we can measure the sizes of raindrops by diameter, by surface area, or by volume. Assuming similar shapes for all raindrops, the surface area is proportional to the square of the diameter, and the volume is proportional to the cube of the diameter. So all of these ways of measuring cannot be normally distributed. So my question is whether the particular way of scaling (i.e., the particular choice of monotonic transformation) under which the distribution does become normal, must carry a physical significance. For instance, should heights be normally distributed or the square of height, or the logarithm of height, or the square root of height? Is there a way of answering that question by understanding the processes that affect height?


title: How should I explain the lack of ERROR in a repeated measures ANOVA table using REML?

text: I have conducted a repeated measures ANOVA using a Linear mixed-effects model fit with REML in R. When I report the results I include the ANOVA table output by `anova(model_object)`               anova(O2.rmtest)                   numDF denDF   F-value p-value     (Intercept)       1   131 168.08616  <.0001     mmolO2.L          1   131   8.35259  0.0045     Lake              2   131  13.83037  <.0001     mmolO2.L:Lake     2   131   1.00465  0.3690      The editor of the journal wants to know why this table lacks an `Error` term. What is the proper explanation for the lack of an error term in this output and is there any other information that I should be including so that readers can properly evaluate these results?


title: What are the requirments for the sequence term in a power series distribution? Why is the binomial distribution one?

text: Given that a power series distribution is a distribution that can be expressed as $$f_{\theta}(x) = {{a\theta^x}\over{C(\theta)}}$$ where $a$ is a sequence of non-negative real numbers, I've also seen it expressed as $a(x)$. I understand that you can re-parametrize a function so that there are separate terms for $\theta^x$ and some function of $\theta$ but what are the limits on what $a$ could be? For example, according to this site, the binomial distribution is a type of power series distribution, such that $a(x) = {n \choose x}$ but $n$ is a parameter, and the sequence $a(x)$ is a function of a parameter and the random variable. What makes it okay to treat $n\choose x$ as part of $a$ but not, for example, $\theta^x$?


title: What R packages do you find most useful in your daily work?

text: **Duplicate thread:** I just installed the latest version of R. What packages should I obtain? What are the R packages you couldn't imagine your daily work with data? Please list both general and specific tools. UPDATE: As for 24.10.10 `ggplot2` seems to be the winer with 7 votes. Other packages mentioned more than one are:   * `plyr` - 4   * `RODBC`, `RMySQL` - 4   * `sqldf` - 3   * `lattice` - 2   * `zoo` - 2   * `Hmisc/rms` - 2   * `Rcurl` - 2   * `XML` - 2 Thanks all for your answers!


title: Which normalization method for different units (zero, negative, etc.)

text: Commonly, **normalization** of ratings seems to mean adjusting values measured on different scales to a notionally **common scale**. I have different metrics that measure different **performance aspects** of a server. I want to **aggregate these values**. The problem is:   * The values are on **different scale types** (interval and ratio scale)   * The values have different minimum and maximum values Most of them are between $0.0$ and $1.0$, where $0$ means **complete absence** of the measured attribute and $1$ means the attribute is maximized to the best characteristic it can gain.   * One of the ratio scale values can be $[0, \infty]$. Theoretically, speaking, practically something around $500$. More is here better $0$ is absence, too.   * The interval scale can be **negative, positive and zero**. Zero does not mean, absence. The less, the better. The more, the worse. How do I approach this? I need to aggregate the results and feel a bit lost. I see two options here.   1. Use a normalization, but which one? Ideally I would all map to $[0,1]$ so I can present them as percentage, although that would be unprecise for the one ratio and the one interval. I know **min-max scaling** (feature scaling?), but is that appropriate? What about **z-score?** With z-score however, it seems I will land again on a negative-positive scale.   2. Use a tendency to aggregate, arithmetic mean (needs normalization), geometric mean (false semantic on zero results), **adjusted geometric mean** , etc. How should I approach this, can someone advise me?


title: Model multilevel proportion data

text: Suppose that one has data on a proportion $P$ measured on a continuous scale. Further suppose that this data has a grouped structure -- some proportions are clustered according to a certain variable $g$ (location, say). Further say that you've got a covariate $X$ that your think is correlated with $P$. You want to know the effect of $X$ on $P$ in order to get conditional distributions of $P$ given $X$. One could simply run a linear model with random effects, but predictions may go beyond the 0,1 bounds, and intervals will be wrong. One could fit a binomial(logit) GLM, but this may be a poor model for the data, as it would lead to overly upward (downward) skewed intervals near the boundaries, and a perhaps-artificial sigmoidal shape to the predictions. Could a Beta regression work for this sort of problem? If so, how would one incorporate the grouped structure of the data? Any ideas for implementation in R?


title: Is there a formal name for this data normalization formula?

text: I am using a generalized formula for normalizing one data range to another but am having difficulty finding its formal name, if it even exists (sorry if my notation is strange): $$ x_b = min_b + \frac{(x_a - min_a)(max_b-min_b)}{max_a - min_a} $$ This converts an $x$ in range $a$ to its respective value in range $b$. For example, if you are mapping the value $2$ from range $[0,\ldots,10]$ into $[100,\ldots,200]$ then $$ 120 = 100 + \frac{(2 - 0)(200-100)}{10 - 0} $$ I chose these values because it is easy to validate in my head; it works for general ranges.


title: Team diversity measure as time invariant predictor of latent growth curve

text: I am specifying a latent growth curve model in MPLUS with a time-invariant predictor of the growh curve. The time-invariant predictor is a team diversity measure (standard deviation) which is based on a 7 item survey scale. Each team has 5 team members. My question is whether the correct approach would be to calculate the standard deviation on the basis of each individual's average score on the 7 items and include this in the model as an observed variable, or calculate the standard deviation per item across the 5 members, and include the SD per item as indicators of a latent variable (assuming that scale reliabilities and Rwg/ICC - within group agreement and intraclass correlation - are all satisfactory).


title: Inter quartile range and outliers in boxplot with logarithmic y-axis

text: Let's say I have a set of data which are lognormal-distributed and I want to create boxplots. The part of the values which is symbolized by the _box_ may lie somewhere around 1. In that case, the inter quartile range will also be single digit. With the lower values lying somewhere below 0.1, even the most extreme values will be covered by 1.5 times the IQR, and thus by the lower whisker. On the contrary, the upper whisker will be quite short, making every value of 10 or higher a possible outlier. As I have up to 9000 observations in each plot, ranging from 10^-7 to 1, I always get plots with very long lower whiskers and a huge field of outliers to the top. This makes me wonder: 1. Is there a rule (of thumb) for what should be considered an outlier among logarithmically distributed data? 2. Is there a more elegant way to configure the whiskers for boxplots with logarithmic y-axis? For R-users, this may serve as a minimal example: data <- rlnorm(50) # random- generated lognormal data IQR(data) # interquartile range boxplot(data, log='y')


title: How to visualize three different data sets on the same graph?

text: First of all, I would just like to point out I have next to nothing knowledge about data representation, so I would be more then thankful if someone could answer this question on as basic level as it could be. I have three different data sets which I need to plot on the very same graph. The first data set is calculated for each year (each graph represents one year) and its purpuse is to show borders between data (e.g. below that value discard all data, between this and this value use weight X, etc). The other data sets represent non-cloudy days wihtin the year and the third set of data are actual values measured every hour for the observed day (14 hours per day). As you can see, in theory, for one year I could have no data at all (all days in the year are cloudy) and for the next year I could have 365 days with each of them carrying 14 different values. Plain and simple, are there any known graphs for as simple as possible representing this kind of data? To explain a little bit I'm pasting one homemade graph for this problem. However, I fear that this kind of graph could turn to be pretty messy if there are lots of sunny days in the year. ![homemade graph](http://i.stack.imgur.com/RfDOl.png) * * * P.S. Attached dummy .csv. In the first column are values for given year (in the picture red line); in the second column are dates (black lines) for a gives year; in the rest of columns are values for a particular day in a year (blue curve). Link to .csv hosted on Google Disk


title: Comparing ROC-curves

text: I would like to find if there is a significant difference between two ROC- curves. I've found the roc.test in the pROC package. However, I cannot seem to find any information on how this test is actually performed. Can anyone explain it to me or refer to a webpage that has information on this? From the documentation for `roc.test()`, I've found that the used test is DeLong's test (based on the article 'Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach'). However, I don't quite understand the article. My ROC-curves have been found by comparing two logistic regressions, where one has a subset of attributes of the other seeing as I want to examine whether these attributes are useful. The probabilities are found using 5-fold cross- validation.


title: Formula for one-sided Hodrick-Prescott filter

text: I am not very familiar with filters. The Hodrick-Prescott filter as one can find it e.g. in wikipedia is two-sided. I also found an R implementation for this in the R package mFilter. There the filter is given as: find $(\tau_t)_{t=1}^T$ such that $$ \left(\sum_{t=1}^T (y_t - \tau_t)^2 + \lambda \sum_{t=2}^{T-1} (\tau_{t+1}-2 \tau_{t} + \tau_{t-1} )^2\right) \rightarrow Min. $$ If I try to formulate a one-sided version of it myself, then I would take backward looking second order differences. I.e. find $(\tau_t)_{t=1}^T$ such that $$ \left(\sum_{t=1}^T (y_t - \tau_t)^2 + \lambda \sum_{t=3}^{T} (\tau_{t-2}-2 \tau_{t-1} + \tau_t )^2\right) \rightarrow Min. $$ How is the usual formulation of a one-sided Hodrick-Prescott filter and does there exist a robust R implementation?


title: Bivariate CDF estimation

text: If I have a sample $(X_i,Y_i)_{i=1}^n$ I can estimate the joint CDF by: $$F_{X,Y}(x,y) = \frac{1}{n}\sum_{i=1}^nI[X_i\leq x,Y_i \leq y]$$ Assume now I observe only $(X_i)_{i=1}^n$, but I know that $Y\sim U[0,1]$. How does one estimate the bivariate CDF in this case? I would be thankfull for any reference or suggestion.


title: Hypothesis testing of significance of least squares predictor: Two methods?

text: I've read about two ways to test the significance of just one predictor when fitting a linear model with least squares. Suppose I fit the following linear model using least squares using standard assumptions. $$y = X \beta + \varepsilon \tag{1}$$ where $\beta = (\beta_1, \cdots, \beta_i, \cdots, \beta_p)^T$. The first method I've seen is the most basic way, where we compute the test statistic $$t = \frac{\hat{\beta}_i}{se(\hat{\beta}_i)} . \tag{2}$$ Then assuming that we've already shown that $\hat{\beta}_i \sim N(X\beta, \sigma^2 I)$ under the assumption that $X$ is non-stochastic, the test statistic $t$ follows a $T$ distribution with $\nu$ degrees of freedom, $$t \sim \frac{Z}{\sqrt{\chi^2/\nu}} = T_\nu .$$ The second approach I've seen is a likelihood ratio testing approach, where $\Omega$ is a linear model with $p$ parameters and $\omega$ is a linear model with $q$ parameters formed by imposing linear constraints on the model $\Omega$. The likelihood ratio test then gives us the following test statistic with an F distribution. $$f_{general} = \frac{(RSS_\omega - RSS_\Omega)/(p-q)}{RSS_\Omega / (n-p)} \sim F_{p-q,n-p}$$. If we let $\Omega$ be the model in Equation (1), and let $\omega$ be that model but with $\beta_i$ removed (or imposing the constraint that $\beta_i=0$) for some particular $i$, then we get $$f = \frac{(RSS_\omega - RSS_\Omega)}{RSS_\Omega / (n-p)} \sim F_{1,n-p} . \tag{3}$$ We know that the square of a T-distributed random variable is an F-distribution random variable, $$T^2_\nu = F_{1,\nu}.$$ So what's the difference between testing with the test statistic in Equation (2) versus with the test statistic in Equation (3)? I think that the Neyman-Pearson Lemma tells us that testing with Equation (3) would give us the most powerful test for $H_0: \beta_i=0$, but can the same be said of Equation (2)? Also, I don't see why it can't be that $t^2 \neq f$, so is it possible for (2) to fail to reject $H_0$ but then for (3) to actually reject $H_0$ for the same model and data?


title: Solving feature bias issues in Learning to Rank with implicit feedback

text: I have a learning to rank system where implicit feedback (from user clicks) is used to determine +ve and -ve examples for the training. The problem is that (obviously) the learner sees only the top few ranked results (+ve or more relevant examples will appear at the beginning). This introduces some kind of a bias in the learning problem. Is there a way to solve this?


title: Imputing missing values in a count time series with variable effort with the goal of trend estimation

text: I have a time series monitoring data set that looks like below: The response is a count.                  Site Year Response Effort     1     1    1        8      3     2     1    2       NA      0     3     1    3        9      3     4     1    4        7      2     5     2    1        4      2     6     2    2       NA      0     7     2    3        3      2     8     2    4        5      1     9     3    1        0      4     10    3    2        1      4     11    3    3        1      2     12    3    4        2      6      (note: full data frame has 20 sites and years) I want to estimate the trend (over time) of the response, and the standard error of the trend, but I first need to account for that fact there are missing values in the counts. A problem is that each data point also has an associated measure of effort (to be included as an offset term). I can make some imputation of the missing data if I assume year to year changes in counts are similar across all sites (not unreasonable assumption) and then use log- linear regression to estimate a trend - but how do I get at the standard error of the trend, accounting for some of the data being imputed. Thanks for any pointers. Thanks


title: Estimating the mean of an asymmetric distribution

text: The problem is how to estimate, by sampling, the mean of a very right asymmetrical distribution (even with right outliers)


title: Define Priors for Dirichlet Distribution parameters in JAGS

text: I'm defining a Multinomial-Dirichlet model in JAGS and want to assign some hyperpriors to the parameters of the Dirichlet distribution. In the WinBugs manual I read that "the parameters of Dirichlet and Wishart distributions and the order (N) of the multinomial distribution must be specified and cannot be given prior distributions. (There is, however, a trick to avoid this constraint for the Dirichlet distribution). What's this trick?


title: Dataset of repeated measures with unbalanced sample size and unequal variace - What to do?

text: I have a dataset of three groups of cells treated with 10 different compounds and am not sure how to check for significant differences between those treatments. Within each group the data is also paired. I tried doing an two-way ANOVA; however, that doesn't work because of the unbalanced sample size. I then used GLM in SPSS, which at least gives me the correct means and standard deviations for each treatment, but it doesn't calculate correct p-values since it assumes equal variances. Could somebody suggest the correct test in this case or knows how to modify the GLM accordingly? Edit: So here is how my dataset looks like in priciple, I only have 8 more treatments: ![enter image description here](http://i.stack.imgur.com/XYNnj.png) I now would like to compare the treatment conditions to each other, so for example 'is there a significant difference between the number of puncta in control and starvation conditions?'. So far I used univariate GLM with the Puncta as dependent variable and group and treatment as independent factors.