title: Linear relationship between explanatory variables in multiple regression

text: I was reading the multiple regression chapter of Data Analysis and Graphics Using R: An Example-Based Approach and was a bit confused to find out that it recommends checking for linear relationships between explanatory variables (using a scatterplot) and, in case there aren't any, transforming them so they **do** become more linearly related. Here are some excerpts of this: > 6.3 A strategy for fitting multiple regression models > > (...) > > **Examine the scatterplot matrix involving all the explanatory variables. > (Including the dependent variable is, at this point, optional.** ) Look > first for evidence of non-linearity in the plots of explanatory variables > against each other. > > (...) > > This point identifies a model search strategy - **seek models in which > regression relationships between explanatory variables follow a "simple" > linear form**. Thus, if some pairwise plots show evidence of non-linearity, > **consider use of transformation(s) to give more nearly linear > relationships**. While it may not necessarily prove possible, following this > strategy, to adequately model the regression relationship, this is a good > strategy, for the reasons given below, to follow in starting the search. > > (...) > > If relationships between explanatory variables are approximately linear, > perhaps after transformation, it is then possible to interpret plots of > predictor variables against the response variable with confidence. > > (...) > > It may not be possible to find transformations of one or more of the > explanatory variables that ensure the the (pairwise) relationships shown in > the panels appear linear. **This can create problems** both for the > **interpretation of the diagnostic plots** for any fitted regression > equation and for the **interpretation of the coefficients** in the fitted > equation. See Cook and Weisberg(1999). Shouldn't I be worried about linear relationships between dependent variables (because of the risk of multicollinearity) instead of actively pursuing them? What are the advantages of having approximately linearly related variables? The authors do address the issue of multicollinearity later in the chapter, but this recommendations seem to be at odds with avoiding multicollinearity.


title: Is statistics.com class worth the money?

text: I'm debating taking some classes on statistics.com and was curious if anyone has taken classes there or knows any one who has or has any suggestions for good ways to learn stats/machine learning on the cheap-ish and on-line. NOTE: I'm already signed up for Andrew Ng's ML class for spring '12. Thanks.


title: How to generalize precision, recall and F-score to non-classification problem?

text: I am working on a problem to predict a value between 0 and 1. The data is skewed, so that only 18% of the values are greater than 0. I implemented a machine learning system, and what is happening is that the lowest error rate on the validation set (18%) is found when all predictions are 0. So, I wanted to implement precision, recall and F-Scores, to find ways of more accurately measuring my results (instead of error rate and MSE). But precision, recall and F-Scores seem to apply to classification problems. What is a better approach to measuring the results in this case? Is there a good way to use precision, recall and F-Scores? EDIT: To clarify, I think it is a regression problem and should output a continuous value between 0 and 1. It is a recommendation engine using Collaborative Filtering. I am creating a matrix between users and items, and minimizing the MSE.


title: What problem do shrinkage methods solve?

text: The holiday season has given me the opportunity to curl up next to the fire with _The Elements of Statistical Learning_. Coming from a (frequentist) econometrics perspective, I'm having trouble grasping the uses of shrinkage methods like ridge regression, lasso, and least angle regression (LAR). Typically, I'm interested in the parameter estimates themselves and in achieving unbiasedness or at least consistency. Shrinkage methods don't do that. It seems to me that these methods are used when the statistician is worried that the regression function becomes too responsive to the predictors, that it considers the predictors to be more important (measured by the magnitude of the coefficients) than they actually are. In other words, overfitting. But, OLS typically provides unbiased and consistent estimates.(footnote) I've always viewed the problem of overfitting not of giving estimates that are too big, but rather confidence intervals that are too small because the selection process isn't taken into account (ESL mentions this latter point). Unbiased/consistent coefficient estimates lead to unbiased/consistent predictions of the outcome. Shrinkage methods push predictions closer to the mean outcome than OLS would, seemingly leaving information on the table. To reiterate, I don't see what problem the shrinkage methods are trying to solve. Am I missing something? Footnote: We need the full column rank condition for identification of the coefficients. The exogeneity/zero conditional mean assumption for the errors and the linear conditional expectation assumption determine the interpretation that we can give to the coefficients, but we get an unbiased or consistent estimate of something even if these assumptions aren't true.


title: When is half normal distribution useful?

text: Julian Faraway, in his book on linear model with R, used half normal plot to find hat(a value in hat matrix) with too large value. That's the only example where I saw half normal distribution is used. But other than that, I can't come up with real application of it though I guess it might be useful for approximating some positive numbers whose mean is near zero. What are real applications of it?


What are the **variable/feature selection that you prefer** for binary classification when there are many more variables/feature than observations in the learning set? The aim here is to discuss what is the feature selection procedure that reduces the best the classification error. We can **fix notations** for homogeneity: for i=0,1, let x1^i,...,xni^i be the learning set of observations from group i. So n0+n1=n is the size of the learning set. We set p the number of features (i.e. the dimension of the feature space). If x is a vector of R^p, x[i] is the ith coordinate. Please give full references if you cannot give the details. EDIT (updated continuously): Procedures proposed in the answers below   * **Greedy forward selection** Variable selection procedure for binary classification   * **Backward elimination** Variable selection procedure for binary classification   * **Metropolis scanning / MCMC** Variable selection procedure for binary classification   * **penalized logistic regression** Variable selection procedure for binary classification As this is community wiki there can be more discussion and update I have one remark: in a certain sense, you all give a procedure that permit ordering of variables but not variable selection (you are quite evasive on how to select the number of features, I guess you all use cross validation?) Can you improve the answers in this direction? (as this is community wiki you don't need to be the answer writter to add an information about how to select the number of variables? I have openned a question in this direction here Cross validation in very high dimension (to select the number of used variables in very high dimensional classification))


title: Question about hierarchical regression

text: I have 4 models in the hierarchical regression output. The models are as follows: Model 1: Control variable Model 2: Control variable and 11 independent variables. Model 3: Control variable, 11 independent variables and moderator. Model 4: Control variable, 11 independent variables, moderator and interaction terms. **My question is this:** **Which of the 4 models do I consider for final output interpretation?** Whether there is a significant relation between different predictors and the outcome variable changes from model to model, in that the relationships that are significant in Model 1 become insignificant in subsequent models. Please advise and thank you very much for your time. Kind regards Muzi


title: Measure the uniformity of distribution of points in a 2D square

text: I have a 2D square, and I have a set of points inside it, say, 1000 points. I need a way to see if the distribution of points inside the square are spread out (or more or less uniformly distributed) or are they tending to gather together in some spot inside the square. I need a mathematical/statistical (not programming) way to determine this. I googled, found something like goodness of fit, Kolmogorov, etc., and just wonder if there are other approaches to achieve this. Need this for class paper. Inputs: a 2D square, and 1000 points. Output: yes/no (yes = evenly spread out, no = gathering together in some spots).


title: Continually updating naive Bayes classifier

text: I am attempting to use a Naive Bayes classifier to classify text. To accomplish this I have created an Excel sheet with a binary distribution for three variables. The workbook can be found here. Assuming that my math is correct, my questions are:   1. Can my training set can be expanded as I classify new inputs? In other words, every time I check a classification that the model has produced, I can add it to the training model but then I might have an uneven number of examples for each class. Is this a problem?   2. How can I incorporate a prior distribution to the equation? If I for example know from prior data that Class A is twice as likely than Class B?   3. How can I incorporate tf–idf to the equation? I can analyze all the data sets a priori and the frequencies of each word in both the corpus and each document, but am unsure how to incorporate this into the Classifier. Thanks in advance for everyone help. AMAS


title: How to implement a regression with non-independent observations in the DV?

text: I have a design where a number of object pairs were rated for association strength. These pairwise association ratings along 5 criteria of interest are expected to be used as IVs in a regression. These pairs include a target object being paired with one of 4 accompanying objects to create the different types of association of interest. Therefore, I have 4 relation types. The DV represents the frequency each of these relation types is chosen over the others, in many 2AFC trials where 2 of the 4 possible accompanying objects are presented with the target object. Thus, the problem with independence of the observations. When I quantify the frequency of choosing A over B, C and D, this is not independent of the frequency of choosing B over A, C or D. An easy fix would be to quantify the frequency of choosing A over B, A over C and A over D independently but I think I am losing important variance this way and I really want to avoid it. One assumption is that the 5 IVs of interest might be significant predictors for some but not all of the 4 relation types. Therefore, I want to run a regression which predicts the DV (frequency of choice A, B, C and D) based on the 5 IVs. **Question 1** : My prefered way of quantifying the DV violates the assumption of independent observations, making the use of a multiple linear regression difficult. Would there be a different approach I can use that is not so sensitive to the independence of observations? E.g. SEM, PCA or some other type of regression (e.g. multi-model regression)? **Question 2** : What would be the effect of this violation of independence of measures? If I do use a multiple linear regression despite that, would it work against me and lead to a conservative fit of my model, potentially hiding significant effects, or would it overfit the data showing an inflated R square? **Addendum** Let's assume that A is chosen 70% of the time. The other 30% will be spread over B, C and D. Let's assume B is chosen 60% of time, the rest will be spread between A, C and D, and so on. Is there any case where I could theoretically say that I am not so concerned about this "minor" independence violation and still do a multiple linear regression on these, rather than split the variance in the way described in my 3rd paragraph above, where the accompanying variance loss would be detrimental to my aim?


title: How to define contrast coefficient matrix?

text: I have this data               data <- as.data.frame(cbind(y,x1,x2,pre))     data$x1 <- as.factor(data$x1)     data$x2 <- as.factor(data$x2)          data$x1 <- C(data$x1, contr.treatment, base=3)     data$x2 <- C(data$x2, contr.treatment, base=2)          y x1 x2 pre     1  16  1  1  14     2  15  1  1  13     3  14  1  2  14     4  13  1  2  13     5  12  2  1  12     6  11  2  1  12     7  11  2  2  13     8  13  2  2  13     9  10  3  1  10     10 11  3  1  11     11 11  3  2  11     12  9  3  2  10          'data.frame':   12 obs. of  4 variables:      $ y  : num  16 15 14 13 12 11 11 13 10 11 ...          $ x1 : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...       ..- attr(*, "contrasts")= num [1:3, 1:2] 1 0 0 0 1 0       .. ..- attr(*, "dimnames")=List of 2       .. .. ..$ : chr  "1" "2" "3"           .. .. ..$ : chr  "1" "2"      $ x2 : Factor w/ 2 levels "1","2": 1 1 2 2 1 1 2 2 1 1 ...           ..- attr(*, "contrasts")= num [1:2, 1] 1 0           .. ..- attr(*, "dimnames")=List of 2           .. .. ..$ : chr  "1" "2"       .. .. ..$ : chr "1"          $ pre: num  14 13 14 13 12 12 13 13 10 11 ...      And I fitted the following model               lm(y ~ x1 + x2 + x1*x2)      My design matrix 'x' is               x <- as.matrix(cbind(rep(1,12), data$pre, c(1,1,1,1,0,0,0,0,0,0,0,0),      c(0,0,0,0,1,1,1,1,0,0,0,0), c(1,1,0,0,1,1,0,0,1,1,0,0),      c(1,1,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,1,1,0,0,0,0,0,0)))          x       [,1] [,2] [,3] [,4] [,5] [,6] [,7]      [1,]    1   14    1    0    1    1    0      [2,]    1   13    1    0    1    1    0      [3,]    1   14    1    0    0    0    0      [4,]    1   13    1    0    0    0    0      [5,]    1   12    0    1    1    0    1      [6,]    1   12    0    1    1    0    1      [7,]    1   13    0    1    0    0    0      [8,]    1   13    0    1    0    0    0      [9,]    1   10    0    0    1    0    0     [10,]    1   11    0    0    1    0    0     [11,]    1   11    0    0    0    0    0     [12,]    1   10    0    0    0    0    0      I'm trying to use this design to reproduce the following table:               Source DF Squares Mean Square F Value Pr > F     Model 6 44.79166667 7.46527778 12.98 0.0064     Error 5 2.87500000 0.57500000     Corrected Total 11 47.66666667          Source DF Type III SS Mean Square F Value Pr > F     pre 1 3.12500000 3.12500000 5.43 0.0671     x1 2 4.58064516 2.29032258 3.98 0.0923     x2 1 3.01785714 3.01785714 5.25 0.0706     x1*x2 2 1.25000000 0.62500000 1.09 0.4055      The first part is fine               XtX <- t(x) %*% x     XtXinv <- solve(XtX)     betahat <- XtXinv %*% t(x) %*% y          H <- x %*% XtXinv %*% t(x)      IH <- (diag(1,12) - H)     yhat <- H %*% y      e <- IH %*% y     ybar <- mean(y)          MSS <- t(betahat) %*% t(x) %*% y - length(y)*(ybar^2)      ESS <- t(e) %*% e      TSS <- MSS + ESS           dfM <- sum(diag(H)) - 1      dfE <- sum(diag(IH))      dfT <- dfM + dfE           MSM <- MSS/dfM      MSE <- ESS/dfE           Ftest <- MSM / MSE     pr <- 1 - pf(Ftest, dfM, dfE)      The contrast coefficient matrix for 'pre' seems correct.               L <- matrix(c(0,1,0,0,0,0,0), 1, 7, byrow=T)     Lb <- L %*% betahat      LXtXinvLt <- round(L %*% XtXinv %*% t(L), digits=4)      SSpre <- t(Lb) %*% solve(LXtXinvLt) %*% (Lb)      MSpre <- SSpre / 1      Fpre <- MSpre / MSE      PRpre <- 1 - pf(Fpre, 1, 12-7)      But I can't understand how to define the contrast coefficient matrix for x1, x2, and x1*x2. What's the problem with the rest of my code? Below an example for how I think I should calculate for x1               L <- matrix(c(0,0,1,1,0,0,0), 1, 7, byrow=T)     Lb <- L %*% betahat      LXtXinvLt <- round(L %*% XtXinv %*% t(L), digits=4)      SSX1 <- t(Lb) %*% solve(LXtXinvLt) %*% (Lb)      MSX1 <- SSX1 / 1      FX1 <- MSX1 / MSE      PRX1 <- 1 - pf(FX1, 1, 12-7)       The result I get for the second part of the ANOVA table is               pre     1    3.1250 6.0000  3.6822  5.4348 0.06711     x1      2    1.0439 3.9189 -3.4291  0.9078 0.46098       x2      1    0.2500 3.1250 -4.1457  0.4348 0.53881       x1:x2   2    1.2500 4.1250 -2.8141  1.0870 0.40554       Thanks!


title: Model from $\hat{Y}$s or model from residuals?

text: This is for modeling revenue by looking at historical data. I am trying to estimate a curve where $x$ = Fiscal Year Quarter and $Y$ = % of Revenue for lifetime of a bid(sale/opportunity/whatever term makes the most sense to you...) The goal of this is so when a sales person puts in a bid that says "This is worth 50,000 over 3 years", we can then use a rough estimate for about how much % of that 50k we will see each quarter over the next 12 quarters ($3\times 4=12$). Now there are 5 different time spans that fit 90% of our sales that is 1/2/3/4/5 years. There are also 4 different market types we sell to. Currently I have been breaking down each market into 5 models, one model for each year and using Polynomial Linear Regression. Granted the $R^2$ isn't the best but it is the best we can hope for. In the end I need to take the models from all 4 markets and make a single model that can be applied to all sales overall. That way we have 5 models not 20 ($4\times 5=20$) models. My question is when I have all of the models from each market should I take all of the 1 year models and extract the $\widehat{Y}$s from the models, plot those then make a model off those? Or should I take all of the 1 year models extract the residuals and make a model off those? Then do the same for the 2/3/4/5 year models to end up with 5 models, or is there a far better way to model this that someone might point me in the direction of? Thank you for any help or collaborative thinking in advance. _Edit 1_ Example : ![MarketModels](http://i.stack.imgur.com/mNydx.jpg) Sorry for not labeling the Axis on the graph but the $X$ is Quarters. So 1 is the first Quarter Revenue from that Market. $Y$ is Percent of Revenue of Lifetime Revenue, so when (1,1) For the BlueLine = then that Market is expected to generate about 1% of the lifetime revenue that item is estimated to make. Here are 3 of the Market's for a 5 year run (Each model is made out of different data as things sold as "Market 1" can't be sold as "Market 2". The Fiscal Quarter they come from vary as "Market 1" doesn't sell at the same time as "Market 2". I need to figure out a single model that can as best as I can estimate how a 5 year (20 quarter) would look no matter the market. That way I could take a bid for $X$ dollars over 20 quarters no matter the market and estimate my revenue projections off a single model. Hopefully this is a bit easier to follow. I would have originally posted these graphs but I was still in the process of finding the base models. So would I do as I did below the 3 lines and model just off the $\widehat{Y}$s or would I do as offered by @Glen_b and model off both the Residuals and the $\widehat{Y}$s?


title: Multiple Chi-Squared Tests

text: I have cross classified data in a 2 x 2 x 6 table. Let's call the dimensions `response`, `A` and `B`. I fit a logistic regression to the data with the model `response ~ A * B`. An analysis of deviance of that model says that both terms and their interaction are significant. However, looking at the proportions of the data, it looks like only 2 or so levels of `B` are responsible for these significant effects. I would like to test to see which levels are the culprits. Right now, my approach is to perform 6 chi-squared tests on 2 x 2 tables of `response ~ A`, and then to adjust the p-values from those tests for multiple comparisons (using the Holm adjustment). My question is whether there is a better approach to this problem. Is there a more principled modeling approach, or multiple chi-squared test comparison approach?


title: Regression-multiple observations per subject

text: I have data for about 1 year, 100 observations, multiple observations per subject, transactions occur on weekly basis but have 6-12 subjects per week, there is no order to this. There is a policy change in latter half of year, I want to model change in dependent variable due to policy change as a dummy variable: time1=0, time2=1.   * Is this a case for fixed effects estimation?  The number of weeks per subject varies a lot and the number of weeks in time1 is greater than time2. Computed means for time1, time2 and percent change=large change in dependent variable, estimated linear model: pay=X1 X2 Time(dummy). Dummy variable is not statistically significant.   * Any suggestions as to how to model this? Can I treat it as panel data?


title: Why is the average of the highest value from 100 draws from a normal distribution different from the 98th percentile of the normal distribution?

text: Why is the average of the highest value from 100 draws from a normal distribution different from the 98% percentile of the normal distribution? It seems that by definition that they should be the same. But... Code in R:               NSIM <- 10000     x <- rep(NA,NSIM)     for (i in 1:NSIM)     {         x[i] <- max(rnorm(100))     }     qnorm(.98)     qnorm(.99)     mean(x)     median(x)     hist(x)      I imagine that I'm misunderstanding something about what the maximum of a 100 draws from the normal distribution should be. As is demonstrated by an unexpectedly asymetrical distribution of maximum values.


title: Principal component analysis - numerous subjects

text: I have bought a book on factor analysis (on PCA more precisely) because it seems to be the algorithm that fits the most to mine my data. Here is the context: I have a table of "subject x variables", each line is a subject (there are numerous subjects) and each column is a variable (here it is network variables but it shouldn't change anything). All variables are quantitative. My goal is to determine which variables are the most correlate each other and which one are not. I would like to find a typology of my variables (and why not after find clusters among subjects) and have a summary view of all variables giving a good idea if relations are hidden in the amount of data. Do you think PCA is good approach for doing it? In my book, they often propose two summary graphs to problems: One for subjects and one for variables, and they say, the two graphs are two point of views of the same problem and must be interpreted considering each other. Nevertheless they study a small amount of subjects on examples (for example 10 countries...). In my case I can have thousands of subjects. The number of variables can be up to 50. In this case, I don't see how I could interpret the "subject graph". Is it a problem? Should I only get the "variables graph"? Moreover, in the examples, the subjects have a meaning (countries, socioprofessionnal categories...), in my case, they don't really have meaning (every subject is a set-top-box for example).


title: about construction of simultaneous test envelope

text: According to Davison and Hinkley's book "Bootstrap Methods and their Application", page 154, based on an `RxN` array where `R=999` is the number of simulations and `N` is the number of ordered samples (smallest to largest) in each row, we could construct:   1. **pointwise** test envelope: by extracting the 2.5th and 97.5th quantile of each column, and then connect the `N` (ordered) dots   2. **simultaneous** test envelope: "we first compute **columnwise ranks**. Then we calculate the proportion of rows in which either the minimum rank is less than or equal to `k`, or the maximum rank is greater than or equal to `R + 1 — k`, or both" Here `k = p(R+1) = 0.05*(999+1) = 50` and `R + 1 - k = 950`. I am confused about the construction of #2 "simultaneous test envelope" above (the explanation is unclear to me...) -- can it be directly made with the `R x N` matrix? Could anyone help me in explaining the procedure in detail? Thanks! Here is the plot that shows both the pointwise (dashed) and simultaneous (solid) envelopes from the book: ![enter image description here](http://i.stack.imgur.com/Z6xhY.png)


title: Variance explained of a mixed effects model in a new data set

text: I'm using the following function to calculate Edwards R^2 (formula 19 in Edwards et al. 2008) of a mixed effects model (I hope the implementation is correct):               r2lmer <- function(model) {        require(aod) # need the aod package for wald.test function        if (class(model) != "mer") stop("mer object expected")        n <- model@dims[['n']] # the number of observations        p <- model@dims[['p']] # number of parameters        df.mod <- n - p             wald.model <- wald.test(b = fixef(model), Sigma = vcov(model), df = df.mod, Terms = 2:p)        wald.F <- as.numeric(wald.model$result$Ftest[1])        ( (p - 1) * 1 / df.mod * wald.F ) / (1 + (p - 1 ) * 1 / df.mod * wald.F )     }      In order to internally valdiate a mixed effects model, I would like to estimate the variance explained when the model is applied to a new data set, that is, when only the data changes, but not the parameter estimates. I'm not sure if it's even theoretically possible. Any help is highly appreciated. ### Edit If it's not a good idea to use an R^2 statistic for LMM, what other performance measures could one use to internally validate LMM? I have a LMM with just one random effects factor (varying intercept due to repeated measures) and several fixed effects factors. I'am primarily interested in the fixed effects and it is my understanding that Edward's R^2 is a good measure for the variance explained by the fixed effects. Edward's R^2 is also recommended for cross-validation in this paper: Cheng et al. Real longitudinal data analysis for real people: Building a good enough mixed model. Statistics in Medicine. 2010, 29 504-520


title: Importance sampling in mixed, discrete/continuous variables

text: Consider the following model $$ X \sim |\mathcal{N}(X;0,1)| $$ $$ Y|X \sim Q(Y;X) $$ where I define $Q(Y=-x|X=x)$ with probability mass $\int_{-\infty}^{-x}\mathcal{N}(x;0,1)dx$, $Q(Y=+x|X=x)$ with probability mass $1-\int_{x}^{\infty}\mathcal{N}(x;0,1)dx$, and the density of $Q(Y|X=x)$ equivalent to a truncated Normal distribution in $(-x,x)$. Assume now that I an unknown $x_{unk}$ is sampled from $|\mathcal{N}(X;0,1)|$ and that $y$ is sampled from $Q(Y|X=x_{unk})$. I am given $y$ and would like to approximate the posterior distribution, $P(X|Y=y)$. Using importance sampling, I would like to take $N$ sample values for $X$ then weight them according to the "probability" of $Y|X$ such that for sufficiently many samples, $\mathbb{E}[X|Y=y] \approx \frac{1}{\sum_{i}w_{i}} \sum_{i} x_{i} w_{i}$. Were $Q(Y|X)$ entirely continuous, I would use $w_{i} = Q(Y=y|X=x_{i})$. Were it entirely discrete, I would use the probability mass function instead. As $Y|X$ is distributed according to a sort of hybrid, what should my weights be?


title: Combining 2 sets of coefficients, weighting one of the sets

text: I have two sets of coefficients from similar data taken at different times. What I want to do is combine the two sets of coefficients giving greater weight to the more most recent set. The goal is building a predictive model. So say I have dataset A from 2009, and dataset B from 2010. My coefficients for A are:               param1: 0.33     param2: 1.224     param3: -0.119      My coefficients for B are:               param1: 0.42     param2: 1.309     param3: -0.011      If I wanted B to be considered twice as important, would it be sound to use a formula like this:               (2*B + A) / 3 = New Coeeficient      And do that for each parameter? Or am I suggesting something that is fundamentally flawed? In general could one combine coefficients effectively using the basic forumula:               (Weight * DatasetACoeffcient + DatasetBCoeffient) / (Weight + 1)      **Edit** This is a multivariate linear regression problem where the datasets may not be available when someone decides something like this needs to be done.


title: Help in drawing confidence ellipse

text: car packages's ellipse function asks for a `radius` parameter. In help says that is the "radius of circle generating the ellipse". Could you please tell me which circle is this? Thank you very much


title: How to generated datasets for estimating risk differences in multiple groups using R?

text: I would like to generate datasets which has 3 treatment groups: A, B, and C based on 3 covariates: $x_{1}$, $x_{2}$ and $x_{3}$. The $x_{1}$ and $x_{2}$ are standard normally distributed and $x_{3}$ is Bernoulli (0.5) distributed. Besides these, assume risk difference between A & B = -0.36, A & C = -0.26 and B & C = 0.10. Can someone give me a hint how to generate these datasets in R? many thanks in advance


title: What kind of test should I use for?

text: I am designing a study for my project. I wanted to test if music affects reading comprehension. This study will be a between group design. The independent variable is type of music and the dependent variable is the score on a reading comprehension task. Half of participants will be randomly assigned to a room with classical (soft) music, and the remaining half of participants will be assigned to the room with rock (hard) music. Would this be a t-test? Then what kind of t-test should I need to use? Or is this a one way ANOVA?   * Someone answered my question as to use independent t test and since i have two independent sample, i would need to conduct unpaired test.... I am still a basic statistic student and i am not aware of this unpaired test... Then how could i make simple design? Could anyone explain why do i need to use unpaired t test? and what will happen if i don't?


title: One-out-of-K coding and standardization

text: I heard from an university course that a normal method for turning categoric variables into binary variables was to simply create a binary variable for each category (one-out-of-K coding). However, these variables should be standardized so the weight of the original variable remains unchanged. How should these new variables be standardized? Independently with their own standard deviation or should they be scaled using a shared standard deviation?


title: Conditional Distribution over the unit disc

text: How can I show that, $U$ and $V$, two independent uniform $(-1,1)$ random variables have a conditional distribution, given that $U^2 + V^2 <1$, that takes the form: $$f_{U,V|U^2+V^2<1} (u,v|w<1) =1/{\pi}, \quad u^2+v^2<1$$ I tried using the CDF technique for $W=U^2+V^2$ $$ \frac{P\left[ U^2+V^2 \leq w , U^2+V^2 <1 \right]}{P \left[ U^2+V^2<1 \right]}=P \frac{\left[ U^2+V^2 \leq w \right]}{\frac{\pi}{4}}$$ because the event $[U^2+V^2 \leq w] \subset [U^2+V^2<1] $ But after that, I am stuck as even if I was to evaluate the double integral at the numerator, it does not seem that I would get the required distribution. Any hints? Thank you.


title: How does my questionnaire analysis plan sound?

text: I will soon have the results of a questionnaire that I have carried out, and I have made a plan on how I will analyse the data. I was hoping that someone could comment on my plan to advise me. About the questionnaire: there are 3 groups of people (let's call them `A`, `B`, and `C`), and they have all answered a questionnaire with 18 questions, all on a 5-point Likert scale. The questions are grouped into 4 subsets. Group `A` is the reference group, and my aim is to see if the the responses of group `B` or `C` are different to `A`. The primary outcome measure is the percentage of the total possible score for each subset of questions. The secondary outcome measure is the percentage of the total possible score for each question. The current outline in which I would present the data is as follows:   * Descriptive analysis: I would detail the demographics of each group of patients and use ANOVA or Kruskal–Wallis (most likely as the data will probably not be normally distributed) to determine if there are differences in the demographics of the groups   * Primary outcome measure: Use ANOVA or Kruskal–Wallis   * Secondary outcome measure: Use ANOVA or Kruskal–Wallis As I understand with ANOVA and Kruskal–Wallis, you can find out where a difference lies if the statistic comes out below the threshold. This is something I would plan on doing. My problem is that I am unsure how I could go about testing to see if any of the demographic variables are associated with outcome. Regression? Furthermore if a demographic variable is associated with outcome, how can I adjust for that in my statistical analysis? I will be using SPSS.


title: interpreting my own data (negative confidence interval in epidemiology)

text: Im trying to measure the absolute change between two means of systolic bloodpressure, and used a simple method for calculating the confidence interval between to means. My goal is to estimate the absolute change for the TREND between 1980-2000. So to accomplish this im comparing the first and last year (i.e 1980 and 2000, which i already have an point estimate and 95%Cl for, achieved with proc STDRATE in SAS 9.4, using a backgroundpopulation of holland). Now to the problem -> For the year 1980 i have a point estimate with 95%cl for systolic bloodpressure at 94.9(73.1 to 120.8) For the year 2000 i have a point estimate with 95%cl for systolic bloodpressure at 80.3(70.1 to 96.8) So the absolute change is -14,6, this means that the asbolute decline is 14.6 with the calculated confidence interval of +/- 22,7 (normal distribution). ### Here's the thing-> i cant make any sense of how to interpret an absolute decline of 14.6 (-8.1 to 37.3). Since im trying to describe trends and already know that systolic hypertension have a decline with 14,6 estimates then how can i describe a lower limit of minus (-)8.1 if its already a decrease. Hope that this makes any sense, its probably a more epidemiological question then a mathematical or programming one. In desperate need of help :)!


title: What is an F test

text: Would someone please tell me what an F test is and what is shows? Also what is an alpha? And how do I evaluate the p values? I posted two questions earlier and one was 'put on hold' because it was not in the right place, and the second one was marked as a duplicate because of the one put on hold. Neither of which answered my question.


title: Calculating distance metrics between a sample set and a point

text: i have a list of text files and i know that these texts belong to a group, by using this group of text files (i.e this is my sample set) i'd like to calculate Jaccard index and Edit distance for each text file that are not in the sample set. By doing so i'd like to generate two features of each text file which are the similarity/diversity of each file from the sample set of text files. After constructing my dataset, classification will be used. how can i achieve this? Problem arises when i try to calculate Jaccard index and Edit distance for each text file. i'm not sure i should calculate the Jaccard index and Edit distance for each text in sample set and take the average of them and use these as a feature, this can be computationally inefficient. P.S: my reputation is too low that i couldn't create the tags of jaccard- index, edit-distance,similarity and divergence


title: Regression model with heteroskedasticity in both variables

text: I've been learning (lurking) from this site for a while and I finally have a question I haven't seen answered yet. I'm doing a flight test and trying to fit the resulting data to linear line. From a combination of ~14 different sensors, I can calculate the 2 values of interest. From uncertainty propagation, I can find the uncertainty in these 2 values (call them x and y), based on the uncertainty in any of the 14 sensors. So, for any point, there will be error in both x and y, and these errors will not be the same in x and y (I think about it as an ellipse). Also, as x and y change, the error in x and y will change, and the ratio will not remain constant. My understanding is these two conditions are called error-in- variables and heteroscedasticity, respectively. So for any given point, I'll have an expected value and error in both x and y. I'm interested in fitting a curve to this data, and I'm having trouble finding regression models that can handle both error-in-variables and heteroscedasticity. I could use some advice on models or good books for engineers that might help. I'm also interested in figuring out how good each of the linear coefficients are. I don't think I will be able to get very good accuracy for any given point, but I will have 10,000 to 100,000 points. Is there a way to leverage this fact to get an accurate fit, even though any single point isn't accurate?


title: Model Trees: M5, Cubist - What about the assumptions of the linear models being used?

text: From what I understood, these models differ from CARTs for regression, mostly because they fit a linear model at the leaves of the tree instead of simply taking an average. They also "smooth" the tree by generating linear models in the intermediate steps of the tree growth process. I have been using the R implementation of them a bit for regression getting very good results. But I wonder about the assumptions of a usual linear model? Multicolinearity, Autocorrelation, Non-Normality being violated doesnt worry people in this case?


title: Building separate logistic regression models for each categorical variable

text: I am building a binary logistic regression model. I am not sure if using the variables as interactions is a better choice than building separate models for level of a categorical variable. Is there a way to determine this? For ex, there are 2 categorical variable with 2 levels each - 1.Regions - East and West 2.AGI - High and Low The response variable behaves differently for each level i.e. East has higher number of positives response and so does high AGI. here is my question - Is it good to build separate models for each level i.e. have 4 models East and high, East and Low, etc.? Or use the 2 as categorical variables in one model? How do you determine when to split the data into separate models? I have built separate models when i knew the behaviour of the data (independent variables) are different in each level of a categorical variable. Thanks!


title: Confusion related to Gaussian Markov Random field

text: I was reading this paper related to Gaussian Markov Random field. I didn't get how they derived this equation from the standard multivariate gaussian distribution equation ![enter image description here](http://i.stack.imgur.com/wT1MZ.png) The multivariate gaussian distribution is given by $P(X) = \frac{1}{{(2\pi)}^{N/2}|\Sigma|^{1/2}} exp^{-\frac{1}{2}x^T\Sigma^{-1}x}$. I didn't get how they derived the above equation. Clarifications guys?


title: Testing significance of groups formed by classification tree leafs

text: I used mvpart to produce a classification tree (response variables concerning travel behaviour, explain variables are socio-demographic). In the output I can see the averages of the travel behaviour variables for the different segments formed by the tree. Is there a way to test if these differences are statistical significant different from each other? Note: In the output I only see the summed SME, nog SE for every travel behaviour variable.


title: Relation between two ordinal variables

text: For a university course, a friend of mine had to test whether there is a relation between two ordinal variables. These variables were opinion about the European Union (positive, neutral, negative) and whether people felt being (a) citizen of their country, (b) primarily citizen of their country, but also European, (c) primarily European, but also citizen of their country and (d) European. The data is available in form of a contingency table. Now, what my friend did was to calculate a chi-squared test. I believe this is wrong, because it does not consider the fact that the groups are ordered, i.e. the variables are ordinal and not categorial. However, chi-squared test seems to be common for this kind of analysis (e.g. here). In my opinion, the idea behind this question is whether or not there is a significant correlation between those variables, and since Pearson’s correlation cannot be calculated from ordinal data, Spearman’s Rho should be used. An alternative to that would be Goodman’s/Kruskal’s Gamma test, as suggested here. Perhaps the Kruskal Wallis test would also be an option? But can it be used in this setting? So, my first question is: Is it valid to use the chi-squared test here? What do you think of the alternatives I described? Question two: How can the significance of Spearman’s correlation be calculated in R? It seems to me most commands reqire observations, not contingency tables, as input.


title: Bayesian mixed model regression with a between subjects factor

text: I'm trying to specify a model in JAGS/rjags with one between subjects factor ( **a** , with two levels - _Y_ , _N_ ) interacting with one repeated measures continuous variable **x** plus subject varying slopes and intercepts that correlate. I can specify this model simply enough with the lmer function:               lmer(y ~ a + x + a:x + (1 + a | id))      My JAGS/rjags is very rusty (or very fresh). The below seems to me to be fitting a model with subject varying intercepts and subject varying slopes while estimating the slope for both levels of **a** , but I'm not sure I'm doing what I think I'm doing. There's also no correlation specified between the two.               modelstring = "     model {       for ( i in 1:Ntotal ) {         y[i] ~ dnorm( mu[i] , tau )         mu[i] <- a1[aLvl[i]] + s1[sLvl[i]] + a2[aLvl[i]] * x[i] + s2[sLvl[i]] * x[i]       }       # Prior:       tau <- pow( sigma , -2 )       sigma ~ dunif(0,1000)       for ( j in 1:2 ) {            a1[j] ~ dnorm( 0.0 , aTau )            a2[j] ~ dnorm( 0.0 , aTau )        }       aTau <- 1 / pow( aSD , 2 )       aSD <- abs( aSDunabs ) + .1       aSDunabs ~ dt( 0 , 1.0E-7 , 2 )       #       for ( j in 1:NsLvl ) {          s1[j] ~ dnorm( 0.0 , sTau )          s2[j] ~ dnorm( 0.0 , sTau )       }       sTau <- 1 / pow( sSD , 2 )       sSD <- abs( sSDunabs ) + .1       sSDunabs ~ dt( 0 , 1.0E-7 , 2 )     }     "      The framework for this comes from Kruschke and this has been of some help too. I would appreciate some pointers or links to examples of similar analyses.


title: Random forest is overfitting?

text: I'm experimenting with random forests with scikit-learn and I'm getting great results of my training set, but relatively poor results on my test set... Here is the problem (inspired from poker) which I'm trying to solve: Given player A's hole cards, player B's hole cards and a flop (3 cards), which player has the best hand? Mathematically, this is 14 inputs (7 cards -- one rank and one suit for each) and one output (0 or 1). Here are some of my results so far:               Training set size: 600k, test set size: 120k, number of trees: 25     Success rate in training set: 99.975%     Success rate in testing set: 90.05%          Training set size: 400k, test set size: 80k, number of trees: 100     Success rate in training set: 100%     Success rate in testing set: 89.7%          Training set size: 600k, test set size: 120k, number of trees: 5     Success rate in training set: 98.685%     Success rate in testing set: 85.69%      Here is the relevant code used:               from sklearn.ensemble import RandomForestClassifier     Forest = RandomForestClassifier(n_estimators = 25) #n_estimator varies     Forest = Forest.fit(inputs[:trainingSetSize],outputs[:trainingSetSize])     trainingOutputs = Forest.predict(inputs[:trainingSetSize])     testOutputs = Forest.predict(inputs[trainingSetSize:])      It appears that regardless of the number of trees used, performance on training set is **much** better than on test set, despite a relatively large training set and a reasonably small number of features...


title: Is it possible to generate means and standard deviations for two groups based on a Cohen's d?

text: **Is it possible to generate means and standard deviation for two groups from a Cohen's d effect size?** Even if the means and standard deviations are trivial/meaningless values, I would like to keep the same effect size. The formula for Cohen's d is as follows:               d = (mean(x1) - mean(x2))  / (pooled sd of x1 and x2)


title: How would you rate the following books in terms of usefulness for experimental data analysis?

text: I am looking at following books:   * Statistical Models (Cambridge Series in Statistical and Probabilistic Mathematics) - A. C. Davison   * Essentials of Statistical Inference (Cambridge Series in Statistical and Probabilistic Mathematics) - G. A. Young   * Asymptotic Statistics (Cambridge Series in Statistical and Probabilistic Mathematics) - A. W. van der Vaart. They all cover different topics, I skimmed through the pages of both of them and it seems to me they are all useful for data analysis, but I wonder which one would be the most useful to buy / read first?


title: Formula for single linear regression for dataset that has uncertainties on both x and y

text: I'm going to teach classes on Physics Laboratory on First year of Bachelor studies. In most of the excercises during data analysis students will have to fit a line to measurements they have taken. I need formulas that fit data to the model: $$ y_i = Ax_i + B + \epsilon_i $$ Where $\epsilon_i$ is error for point $i$. Since in most cases they will have significant errors on both $x_i$ and $y_i$ I need to include: $\sigma_{x_i}$ and $\sigma_{y_i}$). But I have found only formulas that only include errors on $y_i$ by including weights: $$ w_i=\frac{1}{\sigma_{y_i}^2} $$ I would also like to have formulas for uncertainities on $A$ and $B$ that if needed take into account $\sigma_{x_i}$ and $\sigma_{y_i}$. Any references to sources that discusses this problems would be also appreciated. PS. While I would like to be able to use OLS for this --- it is not necessary, model just must be simple enough to be used (and understood) by first year students, and simple enough to be implemented in OOCalc on Excel sheed withhout using advanced stats functions. PS II. I really checked most obvious sources for that formulas, so please be kind :)


title: Rejection-Sampling of Exponential Distribution

text: Consider the following question. > Consider the generation of random numbers following an exponential > distribution with some known mean. Give three reasons why the rejection- > acceptance method would not be the preferred approach for generating random > numbers for this distribution. What is the preferred approach? > > Assume a program is available to generate uniformly distributed random > numbers between 0 and 1. I'm not looking for direct answers to the question, but rather advice on how to go about solving it. I could not find a reason as to why the rejection method is not well suited for the exponential distribution. Thanks


title: What is a "posterior median"?

text: I know what the term "posterior distribution" means (namely the conditional probability distribution of some parameter with respect to some observed data), but I'm reading an author that uses the qualifier "posterior" _all over the place_ , mostly in contexts that I can't make any sense of. For example, "posterior median". What could this possibly be? Edit: based on the comments/answers I've received it looks like my question should have been closer to "what is the median of a _distribution_ (as opposed to a sample)?" More generally, is "posterior" just shorthand for "conditional on some observed values"?


title: assess the interaction of two independent variables, on the dependent one

text: I have a set of laws/policies that overlap with each other through time (for example one started in 1990 and finished in 2001, and the other started in 1995 and finished in 1999). I want to examine how the interaction of those two influence my dependent variable. I am using spss, and I was thinking to set the independent var. (laws) as dummy variable (0 for the years not in force, 1 for the years in force) and run a regression analysis. Does this make sense? Any help/suggestions would be really appreciated.. Thank you


title: Testing the goodness of fit for a Poisson model

text: I need help in calculating the goodness of fit for a Poisson model. The online literature on it is quite scarce, so if you have any useful articles and links I would appreciate it very much. Here is how the table looks like from SPSS. The outcome variable is a count of patents a company has. ![model](http://img827.imageshack.us/img827/1655/screenshot20130606at124.png)


title: Co-occurrence statistics for sets

text: I am looking for help in the following situation: I have a set of numbers $A=\\{1, 2, \dots, 100\\}$, and I am drawing subsets of 10 of these numbers $\\{a_1, a_2, \dots, a_{10}\\}$ according to an underlying distribution that I do not know. After doing this procedure $N$ times (which can be very large), I would like to know if certain subsets of $A$ (of 9 or fewer numbers) are drawn as part of a single larger subset of 10 more often than would be expected by chance. Essentially, I am trying to find subsets of $A$ the co-occur together more often than would be expected by chance when A is sampled in sets of 10. Are there any well established statistics for this sort of thing?


title: What is the pooled standard deviation of paired samples?

text: I am trying to do a priori sample size calculation based on published results. However, I am unable to obtain a reasonable estimate of the published effect size (which is not reported) as I am unable to obtain an estimate of the pooled standard deviation or the standard deviation of the difference. The problem resides in the fact that it is a factorial experiment with two within-subjects factors ($2 \times 3$ levels). I only have the cell means and standard deviations (i.e., for the $2 \times 3$ table) but not the marginal means for the first factor (with 2 levels) which I need. I know that the formula for the pooled standard deviation for independent samples is the square root of ![enter image description here](http://i.stack.imgur.com/sTVlk.png) (taken from wikipedia) .   **But what is the formula for pooled standard deviation for dependent samples?** * * * Means:                    1A   1B     2a 3.24 3.01     2b 2.91 2.56     2c 3.01 3.05      Standard deviations:                    1A   1B     2a 0.65 0.70     2b 0.68 0.60     2c 0.46 0.53      I want to obtain the effect size between 1A and 1B (so pool over levels of factor 2). Sample size is 27.


title: Adjusting the classification threshold of Naive Bayes

text: I've been involved in a machine learning project recently and am now in the process of writing the project up for a paper submission. We used the naive bayes classifier on the project and developed a method for adjusting the classification of datapoints depending on how exactly we want the classifier to perform, ie. increase or decrease recall on one class in a binary classification problem. This is similar to adjusting the classification threshold but is sensitive the feature counts unlike just adjusting the classification threshold. My question is if anyone knows of any previous reserach on methods relating to adjusting the decision boundary in any other way than just changing the classification threshold. **update** (as per steffen's question) thanks steffen for your reply 1) by feature counts I mean feature counts, as in the counts of the features (or tokens/words if you will) in the unseen document. I don't mean the number of articles seen, but the number of features (words) in those articles. 2) I am aware that adjusting the threshold is a common way to tweak the classifier to some particular case, what I would like references to research that mention some other way of achieving this tweaking. By adjusting the classification threshold I mean that instead of treating an unseen document whose total sum of log probabilities (and the prior) is below 0 as a case from class A, we say that the score needs to be (say) below 0.5. Presto we've changed the classification of a bunch of documents. There are however other methods for tweaking the classifications, like for instance weighting the features according to their importance.


title: How to find differentially expressed genes in R?

text: Can you tell me some procedures in R, which are able to find differentially expressed genes for a given microarray data? Moreover,If you can give some ideas related to clustering of significant genes for microarray data, it will be quite helpful for me.


title: How to estimate shrinkage parameter in Lasso or ridge regression?

text: I want to use Lasso or ridge regression for a model with more than 50,000 variables. I want do so using software package in R. How can I estimate the shrinkage parameter ($\lambda$)? Edits: Here is the point I got up to:               set.seed (123)     Y <- runif (1000)     Xv <- sample(c(1,0), size= 1000*1000,  replace = T)     X <- matrix(Xv, nrow = 1000, ncol = 1000)          mydf <- data.frame(Y, X)          require(MASS)     lm.ridge(Y ~ ., mydf)          plot(lm.ridge(Y ~ ., mydf,                   lambda = seq(0,0.1,0.001)))      ![enter image description here](http://i.stack.imgur.com/9xPJ6.jpg) My question is: How do I know which $\lambda$ is best for my model?


title: Help with change scores and/or ANCOVA design

text: I am carrying out an evaluation for an intervention using secondary data. Unfortunately the design of the study is weak as participants cannot be randomised. I am looking to see whether the intervention has had any impact on patients quality of life? A measure of QoL has been taken before and after the intervention, i am also interested in a number of variables and their potential effect on reported QoL including: age, gender, pharmacological treatment, disease duration, mariatal status, and employment status(re literature). These IV's are both continuous and categorical, but can be all be converted to categorical if needs be. I would appreciate your advice concerning which design is best suited to this piece of research. I am considering using a ANCOVA but my experience with violations of assumptions may further weaken the merit of an already compromised design. Warm Regards Dara


title: Possibility of analyzing the hypothesis using Friedman's test

text: Lets assume I have collected some data about sold ice creams on the local market on the daily basis. Days( just a number of observations, might be classified otherwise, it's not variable ) present the day collection of data took place, var1 represents sold ice creams between 09:00-10:00, var2 between 12:00-13:00, and var3 between 16:00-17:00. My Q. is, **Is it possible to conduct Friedman's test on this data with hypothesis** H0: Time of the day isn't the cause factor for the size of the traffic. H1: At least one time of the day causes on the size of the traffic. Friedman's test is as we know used for 3+ dependent variables, however none of those variables depend on each other. So is it possible to conduct such a test on this type of data, or should I use more appropriate K-W, or perhaps Cochran Q test. > days 1 2 3 4 5 6 7 8 9 10 > > var1 42 16 98 67 53 23 45 72 6 5 > > var2 23 34 12 53 12 42 51 62 73 74 > > var3 32 12 53 25 64 43 23 73 23 34 Thanks in advance.


title: Concern about collinearity when adding gender and gender-specific comorbidity for prediction of disease risk

text: I am build a model to predict the risk of having disease X, let say I have a series of variables and I select the variables to be included in multivariate logistic regression model by: i) clinical input and/or ii) p-value<0.25 in univariate regression. I wonder if there will be an issue of collinearity if I include both the `gender` and `gender-specific comorbidity` (e.g. benign prostatic hyperplasia) into multivariate regression? As a further question, I have read from books that if one found the beta of the variable is too big while the p-value is big also, this may suggest the presence of collinearity. Is it correct?


title: Non-significant ANOVA but significant pairwise comparisons: how is it possible?

text: When I do a 2 (tasktype: slogan vs. product) x 2 (novelty: old vs. new) ANOVA, the interaction and the main effects do not seem to be significant. But when I look at the pairwise comparisons, I see that old and new are significantly different from each other both in the slogan and product condition (p`s <= .01). Given this, should not the main effect of novelty come out significant when I do ANOVA? It seems odd to me. Is there an explanation? Thank you!


title: How to obtain decision boundaries from linear SVM in R?

text: I am in need of a package that can give me the equation for a linear SVM model. Currently I'm using e1071 like so:               library(e1071)     m = svm(data, labels, type='C', kernel='linear', cost=cost, probability=FALSE, scale=scale)     w = t(m$coefs) %*% data[m$index,]  #Weight vector     b = -model$rho #Offset      However, I'm not sure how `e1071::svm()` selects positive and negative classes, so I think this might screw up with different data sets. Can anyone confirm how this function decides which class is positive and which one is negative? Also, is there a better package for this?


title: Goodness of fit between columns of a 10000x2 table (with low counts)

text: I have a $10000\times2$ table (see below) where column **A** is made out of _observed data_ and column **B** is the _output of a model_ that I have which should fit the observed data **A**. $n_{A}$ is the total number of points in column **A** and $n_{B}$ the total number of points in **B** and they are **not** necessarily equal.               A   B     -----     1   0     1   0     1   0     0   1     1   0     2   0     1   0     1   0     1   1     0   0     0   0     0   2     0   1     5   7     1   0     1   0     6   4     (...)      How should I characterize the _goodness of fit_ of this model output ( **B** ) with my observed data ( **A** )? The _**Chi-Squared**_ test: $\chi^{2}=\sum\limits_{i=1}^{N} \frac{(O_i-E_i)^2}{E_i}$ (where $O_i$ is the observed frequency and $E_i$ is the expected frequency) has the issue of what to do I do if $E_i=0$? As you can see in the table above, there's **lots** of bins where $E_i$ is zero. The same happens if I use a _**log likelihood Poisson**_ test. This is how I've seen it expressed (taken from here): $-2ln\lambda=2 \sum\limits_{i=1}^{N} (E_{i} - O_{i} + O_{i}ln\frac{O_{i}}{E_{i}})$ it will clearly have the same problem as the _**Chi-Squared**_ test above, whenever $E_i=0$. _**Fisher's exact test**_ was mentioned in the original question (here) but later retracted. So what can I do whit this table? * * * This question comes from this original post. I've tried to phrase it as simple as possible. A very similar question was asked here but it was never fully answered.


title: What other predictive models should I consider

text: I have a numeric dataset involving 34 input variables and one output variable. I'd like to fit a model to the data so to be able to make predictions (regression problem). So far I've tested `MARS (Multivariate Adaptive Regression Spline, earth package), simple random forest (Random Forest, randomForest package), neural networks (Model Averaged Neural Network, package nnet), support vector machines (kernlab package, polynomial, RBF, laplace kernels)`. I'm using R. The best model (both in train and test data sets) is the Laplace kernel (R2 = 89%) but the SVR polynomial (degree = 2), the RBF and the neural networks are close seconds (with 87-88%) and the random forest with 85% and MARS with 82%. Do you suggest any other model that I could try to fit to my data set, based on the results I've fgot so far and on the type (numeric) and size (700 values for each variable) of the data set? Thanks link to data I'd like to predict R (last column of data set)


title: Good books covering data preprocessing and outlier detection techniques

text: As the title goes, does anyone know of a good, up to date book that covers data preprocessing in general and especially outlier detection techniques? The book doesn't need to be focusing exclusively on that, but it should deal with the aforementioned topics exhaustively - I wouldn't be happy with something that's a starting point and quotes a list of papers, explanations of the various techniques must appear in the book itself. Techniques for dealing with missing data preferable, but not necessary...


title: How can I measure model performance with weighted logistic regression?

text: I am working with some survey data that uses probability weights. A number of sources explain that likelihood-based tests and fit statistics like likelihood-ratio, AIC, and BIC are not valid in the context of the weighted MLE. Are there other tests, statistics, or graphics that one can use in this context to get an idea of how one model performs relative to another model?


title: Standardized coefficients and IV method

text: In a multivariate regression, suppose we want to calculate the metric coefficients from the standardized ones. Is the method (standardized coeffcient times standard deviation of the dependent variable) valid in both OLS and IV or not?


title: Difference between Research Design and Experimental Design

text: What is the difference between Research Design and Experimental Design? I can't see any difference. Both of them need to establish _Causality_. Both of them are the arrangement for collection and analysis of data in a manner that aims to achieve the objective in an economical way. Then where is the difference between them ?


title: Solutions to exercises in Nelsen's “An Introduction to Copulas”

text: I am paving my way through Roger Nelsen's "An Introduction to Copulas". The book has exercises (quite good actually), but no solutions. Does anybody have a solution manual for (some of those) exercises? Also references to Copula related (intro level) problem sets with solutions would be appreciated.


title: Designing an experiment where participants see all levels of one factor, but not of the other

text: Let's say I have an experiment with two factors _A_ and _B_ , each with a number of different levels. I want participants to   * see all levels of factor _A_ exactly $n$ times,    * never see a level of _B_ twice How would I generate a test plan that assigns me $x$ users to the various combinations of _A_ and _B_ given the above constraints? More specifically, suppose I have 5 pairs of speakers (factor _A_ ) and 10 audio samples (factor _B_ ). I would like each speaker to be tested by 12 people. I therefore need to recruit 60 ($5 \times 12$) people to listen to all pairs of speakers. Since the audio samples are so long they couldn't possibly listen to all of them. How do I assign them to the audio samples? I tried starting with a smaller version.. 3 speakers, 6 audio sources, and 6 listeners (A through F), meaning I get two listeners per speaker. But how can I generalize this? ![](http://i.stack.imgur.com/1Dehs.png) Perhaps I'm missing the terminology to look it up somewhere, but I only got as far as creating the full test matrix, which isn't really helpful here.


title: sub-questions with likert scales choices

text: I want to determine compliance to a certain standard in ISMS and to determine that, I have the standards started and it has sub-questions (say 2,3,4 questions), with each having a 5 point likert type choices. I have a sample of about 32 questionnaires filled. My difficulty now is how to treat or analyse the sub-question and draw one chart to represent the 5 point scales (for each standard, not each question) and also determine the degree of compliance. For example: > A. Information Security Policy   >  A1. Is a published policy document, approved by management:   >  ("Not compliant", "Partly Compliant", "Currently Compliant","Largely > Compliant", "Compliant") > > A2. Is the published policy reviewed regularly, (at planned intervals) or if > significant changes occur to ensure its continuing suitability, adequacy and > effectiveness?   >  ("Not compliant", "Partly Compliant", "Currently Compliant", "Largely > Compliant", "Compliant")


title: Error Propagation Calculation

text: I have a few machines that are used to calibrate each other. Machine 1 has is accurate to 0.025% Machine 1 is used to calibrate Machine 2, which has an accuracy of 0.005% Machine 2 is used to calibrate Machine 3, which has an accuracy of 0.025% Machine 3 is used to calibrate Machine 4, which has an accuracy of 0.04% Using the root of the sum of the squares gives an error for Machine 4 of just over 0.052%, but I need it to be below 0.05%. Is there any games I can play (like bootstrapping, maybe) to get this error down? In other words, I can get all sorts of empirical calibration trials...Can I use that data somehow to bring down that error propagation?


title: On mutual correlation of multiple time series

text: The question is probbaly easy, but I don't want to fool myself. If there are three (univariate real) time series $A$, $B$ and $C$ such that $A$ and $B$ both highly correlated with $C$ at some arbitrary lags, does it implies that $A$ and $B$ also highly ( _near_ as high) correlated at some lag? I presume it does, because I can't imagin a counter example.


title: Should I log transform my volatility variable?

text: I'm wondering if my volatility factor is specified correctly. My data consists of log returns on the S&P 500 index, a measure of news sentiment, and a newscount variable (# of articles published on that day). Now I want to include a volatility variable, and I tried several but will explain my current issue. I used the square of log returns, but as you can imagine, the distribution of it is skew: 10.03 and kurt:132.48. If I take the log of this, it becomes skew:-0.96 and kurt: 5.36. So it is closer to normality now. Does it make any sense to take the log again? What can it mean for inference if I were to conduct a VAR analysis? (elasticities?) I did for both and included the table for granger causality, using both log volatility (left) and (unaltered) volatility (right) measures. I hope someone can help me. Cheers ![Granger Causality table](http://i.stack.imgur.com/RBAvk.png)


title: Approximate the distribution of the sum of ind. Beta r.v

text: If $X_i$ has a Beta distribution $\beta(1,K)$. What is the best approximation for the distribution of $ S=\sum_{i=1}^N X_i$, when the $X_{i}$ are independent and $N$ is finite. Thanks


title: Cluster analysis on ordinal data (Likert scale)

text: I want to do clustering of my data in R, using kmeans or hclust (I am a new R user). My data is ordinal, Likert scale, to measure the causes of cost escalation. I have 41 causes "variables" that scaled from 1 to 5 (1: no effect, 5: major effect). I have about 160 observations "who rank the causes". I would like to cluster the variables (the columns) in terms of similarity of occurrence in observations, but I don't know how to start. Do I have to convert the scale to percentage or z-score before clustering? My data is available and shared as a Google Drive spreadsheet.


title: Apparent fallacy in Fisher-Irwin test

text: For a hypothesis $H_0\equiv p_1=p_2$ (for two Bernoulli distributions), we take a measurement and find that $X_1=i,X_2=j$, and we denote $k=i+j$ Now, to test $H_0$, the Fisher-Irwin test tells us to compute $$P_{H_0}(X_1=i|X_1+X_2=k)=\frac{\tbinom{n_1}{i}\tbinom{n_2}{k-i}}{\tbinom{n_1+n_2}{k}}$$ As I understand this, this probability is _supposed_ to be "the probability that we would get the measured data $i,k$ in a measurement assuming that $H_0$ is true". We then compare this probability with the significance level, as in any other hypothesis testing exercise -- we always calculate the probability that we would obtain the measured dataset, and then compare with the significance level. However, the given expression is _not_ the probability that we would get the data $i,k$ in a measurement assuming that $H_0$ is true". That probability is $$P_{H_0}(X_1=i \:\cap X_1+X_2=k )=P_{H_0}(X_1=i|X_1+X_2=k)\times P_{H_0}(X_1+X_2=k)$$ The probability distribution suggested by the Fisher-Irwin test seems to miss out the fact that there are two independent variables and not one ($X_1$ and $X=X_1+X_2$). The expression gives me the impression that it is presupposing $X=k$ and then going on to deal with the probability for $X_1$. Why do we use $P_{H_0}(X_1=i|X_1+X_2=k)$ in the Fisher-Irwin test and not $P_{H_0}(X_1=i \:\cup X_1+X_2=k)$? The explanation1 2 given in Sheldon Ross' book isn't quite convincing.


title: Plotting log likelihood function of Pareto distribution inR

text: How can I plot log likelihood function with type "I" Pareto distribution function $f(x|\alpha,k)=\alpha k^\alpha/x^{\alpha+1}$ in R? The problem for me is there are two parameters, how can I plot the log likelihood with two unknown parameters? Do I need to give a value of $\alpha$ and then plot the log likelihood function with respect to single parameter $k$?


title: Is it acceptable to rotate factors with PCA for binary data?

text: What issues, if any, might there be in rotating factors in order to obtain factor/component loadings of binary data? Is it acceptable to rotate the factors when doing a traditional PCA? (Assuming I’m using a polychoric correlation matrix.) I am not a statistician, but an applied educational researcher making sure my approach is appropriate to the problem at hand. My doubt comes from a note in the SAS support site that advises to use `PROC PRINCOMP` or `PROC CORRESP` instead of `PROC FACTOR` (in all cases to be careful to factor a tetrachoric correlation matrix). One important difference between these two procedures and `PROC FACTOR` is that the former do not implement factor rotation, even though a PCA can be specified within `PROC FACTOR`, and would otherwise seem to be preferable for exploratory analysis. It occurred to me that perhaps rotation of factors is inappropriate for binary data, but I cannot find a discussion of this in literature or stat forums. * * * The problem, in case you care to know the context: We have data from several hundred colleges around the U.S. regarding what programmatic/curricular elements they implement in certain kinds of cohort educational programs. 24 variables, all dichotomous. We are looking to find patterns of how those elements tend to group together (in addition to a separate latent class analysis or cluster analysis to categorize the programs). Factor analysis is not appropriate since there is not a theoretical underlying latent structure that these elements are “measuring”.


title: Assumptions of GAM

text: I am looking to understand the assumptions of using a generalized additive model. 1) Are the assumptions the same as the assumptions for each equivalent link function in a generalized linear model - e.g. as linear regression, logistic regression etc. along with the assumption that the additive smooth structure is correct and the errors are independent? 2) If these are the same assumptions, do they only play a role when you are making hypothesis tests (e.g. making inference about linear and smooth components)? **ADD** : From Simon Wood (author of mgcv package in R), sounds like one indeed treats the need to check assumptions of GAM as the assumptions underlying the generalized linear model with the same link.   * Since a GAM is just a penalized GLM, residual plots should be checked, exactly as for a GLM. I The distribution of scaled residuals should be examined, marginally, and plotted against covariates and fitted values. residuals(model) extracts residuals.


title: Can I compute a confidence level for a z-test based on the size of the samples?

text: I'm doing a two sample z-test for the difference between means using pooled variance. Suppose I've run two experiments. In the first experiment, the sample size for the two populations is only 100, and I get a z-value of 5.0. In the second experiment, the sample size for the two populations is one million, and I also get a z-value of 5.0. Can I quantify somehow a different confidence level or confidence interval for z-values from the two experiments? I'm hoping there's some generally applicable quantitative formula, not just general guidelines about "sample sizes being sufficiently large".


title: Does independence imply conditional independence?

text: If two or more variables A, B, C, etc. are jointly mutually independent of one another, does this imply that that they are also conditionally independent given some set of conditioning variables X, Y, Z, etc.? If you test variables A, B, C,, etc., and the hypothesis of independence is not rejected, are you generally safe in assuming that these variables will also survive a test of conditional independence given X, Y, Z, etc.? Is it safe enough that you do not need to bother running the latter test? Why or why not?


title: Incoherent ANCOVA?

text: I am doing an ANCOVA model in order to explain the gap in a given distance. So I have a control group, and I have several quantitative variables, right now I am trying to evaluate the impact of the quantitative variables individually. But I think the results I get are incoherent, lest see:                                          Df Sum Sq Mean Sq F value   Pr(>F)         as.factor(groupe)           1  738.8   738.8  21.931 1.03e-05 ***     BASE0008                    1   36.6    36.6   1.087  0.29992         as.factor(groupe):BASE0008  1  270.0   270.0   8.015  0.00576 **      Residuals                  87 2930.9    33.7       These are the results from the ANOVA table, we could say that the group has a significant effect so does the interaction, but when I look at the results of the regression model, I find this:                                             Estimate Std. Error t value Pr(>|t|)        (Intercept)                    1.16666    3.38404   0.345  0.73111        as.factor(groupe)NAV           4.48191    3.79515   1.181  0.24084        BASE0008                       0.26027    0.08651   3.009  0.00343 **     as.factor(groupe)NAV:BASE0008 -0.26902    0.09503  -2.831  0.00576 **      Well, the interaction is still relevant, but it looks like there is not effect of the group and the quantitative variable is more important to determinate the output of the experiment. I want to know if my interpretation is accurate: What can I say about the group?


title: How to use Box-Muller transform to generate n-dimensional normal random variables

text: I'm trying to generate random variables. I read about **Box-Muller transform** which is a way to generate a pair of normal variables, 2-d normal distrubution. But how do I expand that transform to generate 3-d, 4-d, etc. normal variables with this approach? Or is there a different approach to consider?


title: How to derive Bayes prior given large amounts of sample data from different population members

text: In Poker, it is important to judge the likelihood of an opponents action, but for unknown opponents, we have few samples available. One method proposed in http://www.husng.com/content/interpreting-small-sample-sizes-bayesian- estimators was to use bayesian estimators, but the author used a simple heuristic for constructing the prior. Assume we have a database with a large number of opponents and their frequencies of the actions (however, the number of samples per opponent might be very different, e.g. there might be opponents where we have just 5 hands and others we have 10000 hands). In my understanding, each opponent can be described by a binomial distribution and our prior could be determined by "combining" those into a Beta distribution. The questions I have are:   1. Is the thought process so far correct? Especially, does it make sense to split the observations into several opponents or can we see all different ones as a single "default" opponent and it leads to the same prior?   2. Is it a problem that we have a different sample size for each opponent? I'm thinking especially on cases where we have a tiny sample size for an opponent.   3. How can calculate the parameters of our prior beta distribution with parameters $(\alpha, \beta)$ based on the data?


title: Hypothesis test on the Euclidean length of an unknown vector

text: ## Question Suppose I observe a vector $\mathbf{x}=[X_1 \ldots X_n]$, where each $X_i=m_i+n_i$, with $n_i$ being an independent zero-mean Gaussian random variable with variance $\sigma^2$ (i.e. $n_i\sim\mathcal{N}(0,\sigma^2)$ i.i.d.) and $\mathbf{m}=[m_1 \ldots m_n]$ is an unknown vector with a known Euclidean length to the origin $l(\mathbf{m})=\|\mathbf{m}\|_2=\sqrt{\sum_{i=1}^n m_i^2}$. Suppose that the Euclidean length of the unknown vector $\mathbf{m}$ can equally likely be either $l_0$ or $l_1$, and, without loss of generality, $l_0<l_1$. Now, given the vector of noisy observations $\mathbf{x}$ I need to decide whether the length of $\mathbf{m}$ is $l_0$ or $l_1$. Note that I do _not_ know (nor do I need to know for my answer) the values in $\mathbf{m}$, just the two possibilities for length. The intuitive procedure is to compute the Euclidean length of the observed vector $S=\sqrt{\sum_{i=1}^nX_i^2}$, and then select whichever of $l_0$ or $l_1$ is closer to $S$ (based on some threshold that is related to the probabilities of error that you are willing to tolerate). Since $S^2$ is a non-central chi-square random variable, one can (numerically) obtain (at least bounds) on the probabilities of error of this approach. However, is $S$ the "best" test statistic (in the sense of Neyman-Pearson optimality) for this problem? If it is not optimal, is a better test known? If $S$ is optimal, is there a proof of its optimality in the literature? * * * ## My prior effort This seems like a problem that should be well-studied. In fact, it is related to the problem of non-coherent detection in communications theory, where $S^2$ is used by the square-law detector. John Proakis has a proof on pages 304-306 of "Digital Communications" (4th edition) that $S$ as defined is the N-P test statistic for when $n=2$ (detection of a complex-valued symbol with arbitrary phase offset). He essentially projects the problem into polar coordinate system, and takes the expectation of the test statistic over the uniform distribution for the angle. I tried a naive approach of extending Proakis' proof to $n$-spherical coordinate system (since $\mathbf{m}$ can be thought of a coordinate on an $n$-sphere with radius $l_k$) and computing the likelihood function (from N-P) by taking the expectation over the uniform distribution on all the angles. However, while I can get a closed form expression $n=3$ which yields the N-P test statistic defined by $S$, the integration gets really nasty for $n>3$. Perhaps there is something more clever that one can do, maybe using the circular symmetry of joint distribution of i.i.d. Gaussians. This question is related to my previous question.


title: Intraclass correlation (ICC) per item instead of across items

text: I am using the intraclass correlation coefficiant to determine inter-rater reliability for a set of pictures that were rated by several raters on a scale from 1-9. I am using the two-way mixed model in SPSS and am interested in absolute agreement (rather than consistency). This works fine for all items (pictures) taken together, but I am also interested in the ICC for each picture individually, because I am planning to remove pictures with particularly poor agreement between raters. I have seen ICCs per item reported in several papers, but I don't know how to compute them in SPSS. When I use the ratings from only one picture for my analysis, SPSS aborts the process because the number of cases is too small. So my question is: How do I calculate the inter-rater reliability item-wise? Many thanks in advance!


title: Any test methods to show the linear behaviour of a data set?

text: I have a data set that has been collected by a set of empirical experiments. Having plotted them on a graph, it seems that the data behavior is linear. On the other hand, some colleagues claim that the data are not linear. How can I show them the confidence level of being linear for my data? In other words, is there any formal method for showing the linearity of data instead of plotting them?


title: ols regression on stationary series

text: I am trying to regress (OLS) some time series on stock returns. I am not interested in regressing the returns of those time series on my stock returns, but I want to include information about the level (the higher the better) of these time series. I read somewhere that the dependent variable has to be stationary in order to correctly regress it on my returns. Here's the approach I was thinking about: for series on the interval [-100;100] I use an indicator variable =1 for positive values and =-1 for negative values. Does it matter that I will have many 1's in a row, then many -1 in a row? It is by far not as erratic as for the returns. Is my approach ok? Or how could I modify while keeping the same idea?


title: Proportions with repeated measures

text: I wanted to ask what is the best way to analyze proportion data with repeated measures. I will give a simple example. Let's say you have a binary variable Y, which is a success/failure variable. You take a sample of n people. In normal circumferences, you would count how many success you had out of n, and get a point estimate along with a CI (Clopper-Pearson or other). But let's say, that some of the n people, gave 2 samples. For example let's say you are testing some cream that suppose to heal a wound. And some people come with two wounds and put the cream in 2 locations. You end up not with n samples, but N samples, some are correlated because they are coming from the same person. How would you now calculate a point estimate and CI for the success rate, taking into account the within-subject correlation ? Thank you !


title: Regression for a Rate variable in R

text: was tasked with developing a regression model looking at student enrollment in different programs. This is a very nice, clean data set where the enrollment counts follow a Poisson distribution well. I fit a model in R (using both GLM and Zero Inflated Poisson.) The resulting residuals seemed reasonable. However, I was then instructed to change the count of students to a "rate" which was calculated as students / school_population (Each school has its own population.)) This is now no longer a count variable, but a proportion between 0 and 1. This is considered the "proportion of enrollment" in a program. This "rate" (students/population) is no longer Poisson, but is certainly not normal either. So, I'm a bit lost as to the appropriate distribution, and subsequent model to represent it. A log normal distribution seems to fit this rate parameter well, however I have many 0 values, so it won't actually fit. Any suggestions on the best form of distribution for this new parameter, and how to model it in R? Thanks!


title: How does one get the gamma distribution parameters from a word problem description?

text: I am specifically confused as to the meanings of the shaping and scaling parameters of the gamma distribution and how they are used in a real context. Once I find the correct parameters, I have the equations to find everything else I need. The word problem that I am looking at is: > Audrey, an astronomer is searching for extra-solar planets using the > technique of relativistic lensing. Though there are believed to be a very > large number of planets that can be found this way, actually finding one > takes time and luck; and finding one planet does not help at all with > finding planets of other stars in the same part of the sky. Audrey is good > at it, and finds one planet at a time, on average once every three months.   >   > Find the expected value and standard deviation of the number of planets she > will find in the next two years. Through logic alone, one can figure out that the expected value is eight, which should be equal to alpha times beta, if I am understanding the gamma distribution correctly. Our teacher taught us gamma distribution with replacement parts. For instance, a system with 4 spare parts (plus the one already in the system) where each part lasts on average 4 months would be represented by a gamma distribution with shaping parameter 5, and scale parameter 4. On the other hand, I may have misinterpreted the problem and am trying to use the wrong distribution. That's just the distribution we have been using, so I assumed that's the one to use. Instead, I am thinking that this may be a Poisson distribution. Then the second part of the problem is: > When she finds her sixth new planet, she will be eligible for a prize. Find > the expected value and standard deviation of the amount of time until she is > eligible for that prize. I would think this would use the gamma distribution with alpha = 6 and Beta = 3.


title: Grouping strings by common themes?

text: I am attempting to group, for example, strings about programming with other strings about programming, strings about physics with other strings about physics, etc., for a wide range of topics. Despite the glaring theoretical linguistic aspect of the problem, I am looking to actually do this using programming/software. **The rundown:** Given a large number of strings, how would I go about grouping them by semantic theme? **The particular application:** I have ~200k trivia questions that I would like to categorize into common groupings (cars, computers, politics, Canada, food, Barack Obama, etc.). **What I've looked into:** Wikipedia has a list of natural language processing toolkits (assuming that what I'm trying to do is actually called NLP) so I have looked at a few but none seem to do anything similar to my needs. **Notes:** It has been pointed out that doing this requires additional knowledge (e.g. a Porsche being a car, C++ being a programming language). I assume then that training data is needed, but if I have only the list of questions and answers, how can I generate training data? And then how do I use training data? **More notes:** If the current formatting of my Q&As help (although it looks like JSON, it's basically a raw text file):               // row 1: is metadata     // row 2: is a very specific kind of "category"     // row 3: is the question     // row 4: is the answer     {       15343       A MUSICAL PASTICHE       Of classical music's "three B's", he was the one born in Hamburg in 1833       Johannes Brahms     }      But before someone points out that there already exists a category, note that there are ~200k questions and answers like this, and basically as many "categories". I am trying to group these into broader groups like the ones listed above. Also, this formatting can be changed for all the questions very easily, I do it programmatically. **And more notes:** I don't actually know how many categories I'll need (at least 10-20), because I haven't read through _all_ of the questions myself. I was partially expecting to have the finite number determined somehow during categorizing. In any case, I can always manually create a number of categories.


title: Binomial random variable conditional on another one

text: On the Wikipedia page for the Binomial distribution, the following property is mentioned (under the related distribution section): (paraphrased) If $X\sim \text{Bin}(n,p)$ and $Y|X \sim \text{Bin}(X,q)$, then $Y \sim \text{Bin}(n,pq)$ I interpret this the following way. The probability mass function for $X$ is: $$P(X=x) = \binom{n}{x}p^x(1-p)^{n-x}$$ The conditional mass function for $Y$ given $X=x$ is: $$P(Y=y|X=x) = \binom{x}{y}q^{y} (1-q)^{x-y}$$ The mass function of $Y$ is: $$P(Y=y) = \binom{n}{y} (pq)^y (1-pq)^{n-y}$$ There is no citation for this particular property. I have tried to prove it, but to no avail. I wrote the following R code to get a sense of the veracity of the claim.               # Observations of X & Y to be generated     obs <- 10000          n <- 10     p <- 0.6     q <- 0.4          X <- rbinom(obs, n, p)     Y <- X          for( i in 1:obs)     {       Y[i] <- rbinom(1, X[i], q)     }          # Simulated pmf of Y      hist(Y, breaks=obs)          # Theoretical/claimed pmf     Y_theoretical <- rbinom(obs, n, p*q)     hist(Y_theoretical, breaks=obs)      The two histograms generated are shown below: (The simulated pmf) ![Simulated distribution](http://i.stack.imgur.com/8wGSl.png) (Claimed pmf) ![Claimed distribution](http://i.stack.imgur.com/uSFrw.png) Both seem identical for the choice of $p$ and $q$. Can a proof of this claim be provided?


title: Alternatives to absolute error?

text: Let me explain my scenario in which I need to calculate absolute error. Lets say the X is the actual value. And X' is the value of X with some error 'e'. So X' = X + e'. Lets say i = 1 to 10000. I have X(i) as well as X'(i). Now I want to know (also want to plot on graph) that How much X' differs from X. One way to do is to calculate Absolute error and plot it on the graph. I want to know if there are other ways to see how much X differs from X' besides absolute error? and Which one will be more appropriate in my case as I have to plot it as a graph? And also easier to follow for those who will see that plotted graph? Thanks! :)


title: Visualize effects of a regression with categorical explanatory variables (3 levels) in R?

text: Using R, I want to run a linear regression to estimate the abnormal return on days with positive, negative and neutral news (CLASS). I'm a beginner in R, as well as in using regression models! First of all the data structure is as follows. CONTROLVAR just represents all the columns I use as control variables.               DATE <- c("1","2","3","4","5","6","7","1","2","3","4","5","6","7")     COMP <- c("A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B")     RET <- c(-2.0,1.1,3,1.4,-0.2, 0.6, 0.1, -0.21, -1.2, 0.9, 0.3, -0.1,0.3,-0.12)     CLASS <- c("positive", "negative", "aneutral", "positive", "positive", "negative", "aneutral", "positive", "negative", "negative", "positive", "aneutral", "aneutral", "aneutral")     LITIGATION <- c(1,0,0,0,0,0,0,0,0,0,0,0,0,0)     POLLUTION  <- c(1,0,1,1,0,0,0,0,1,0,0,0,0,0)     LAYOFF     <- c(0,1,0,0,0,0,0,0,0,0,0,0,0,0)     CONTROLVAR <- c("11","13","13","14","13","14","12","11","13","13","14","13","14","12")          mydf <- data.frame(DATE, COMP, RET, CLASS, LITIGATION, POLLUTION, LAYOFF, CONTROLVAR, stringsAsFactors=F)          mydf           #    DATE COMP   RET    CLASS LITIGATION POLLUTION LAYOFF CONTROLVAR     # 1     1    A -2.00 positive          1         1      0         11     # 2     2    A  1.10 negative          0         0      1         13     # 3     3    A  3.00 aneutral          0         1      0         13     # 4     4    A  1.40 positive          0         1      0         14     # 5     5    A -0.20 positive          0         0      0         13     # 6     6    A  0.60 negative          0         0      0         14     # 7     7    A  0.10 aneutral          0         0      0         12     # 8     1    B -0.21 positive          0         0      0         11     # 9     2    B -1.20 negative          0         1      0         13     # 10    3    B  0.90 negative          0         0      0         13     # 11    4    B  0.30 positive          0         0      0         14     # 12    5    B -0.10 aneutral          0         0      0         13     # 13    6    B  0.30 aneutral          0         0      0         14     # 14    7    B -0.12 aneutral          0         0      0         12      aneutral (neutral) will be the reference category. I would also like to see the effect of certain subjects of the article on the abnormal return. I'd like to include interaction, so my model looks like this:               mymodel <- lm(RET ~ CLASS * (LITIGATION + POLLUTION + LAYOFF)   # Interaction Variables                    + CONTROLVAR,     # Control Variables                    data=mydf)      Now I don't really understand how to interpret the coefficients I get and I would like to plot a regression line or anything that visualizes these effects to better understand the results. What's a good way of doing this for a three- class-problem like this? `abline()` doesn't seem to work, because there are too many variables. Thank you!


title: Does a high value for covariance suggest a strong relationship?

text: I have a data set consisting of about 100,000 rows and two columns. The values in one column range from 0 to 10^9 the values in the other column range from 0 to 245. I've calculated the covariance between these two values and the result is 5x10^8. I realise this suggests two things:               1) There is a correlation     2) It's a positive correlation      However, I'm now trying to read into this a little further. Does the magnitude of the value have any significance? Would a higher or lower (but still >0) value suggest a closer correlation or can the covariance not be used for this type of analysis? If it can't can anything other than linear regression?


title: How to promote a regression tree over a GLM?

text: Does anybody have any suggestions about promoting the use of a regression tree over a GLM when the two models fit the data almost exactly the same? My team's current arguments are a) a tree is easier for non-technical people to understand and b) when the two models' predictions are implemented -simplified as needed- as tables (necessary for use in actuarial software) the GLM results will be distorted more than the tree's predictions. Any other ideas? Background: Office politics. We want to use one of our own models instead of having our parent company's GLM forced on us without any comparison of the models, so to get our work a chance to be fairly evaluated we need to be promote our tree's advantages over a GLM.


title: Exception for sum of deviations from mean being 0

text: I was reading here > the sum of the deviations about the mean will be 0, except for possible > rounding. Could anyone explain me the what does it mean? I know about sum of deviations from mean being zero but what about this **_except for possible rounding_**?


title: Quicker way to optimize alpha parameter in croston model using R

text: I have the following R code:               library("forecast")     library("rbenchmark")          v <- c(95.3, 96.8, 97.2, 97.9, 98.2, 98.5, 99.3, 99.9, 102.7, 104.2, 107.2, 109.3, 109.8, 110.5, 111.2, 111.6, 112.6, 113.5, 114.3, 115, 117.4, 119, 119.3, 120.1, 120.7, 121.1, 121.9, 122.2, 122.9, 124.6, 129.1, 128.8, 129.3, 130.3, 131.1, 131.8, 132.9, 133.5, 134.5, 135.6, 139.4, 138.9, 140.9, 141.8, 142.8, 142.4, 142.8, 143.9, 144.8, 145.5, 149.6, 151, 152.1, 152.7, 154.5, 153.9, 153.8, 153.8, 154.5, 155.5, 158.9, 159, 158.7, 159.7, 159.7, 159.7, 160.1, 159.9, 161.3, 161.9, 164.4, 164.4, 164.7, 165.1, 165.2, 164.9, 166.9, 167.8, 169.4, 170.1, 171.6, 172.9, 173.7, 175.2, 175.8, 176.3, 177.1, 177.5, 178.8, 180.2, 183, 184, 184.7, 186.5, 187.3, 187.9, 187.9, 188.7, 190.2, 191.8, 199, 199.9, 205.4, 205.2, 206.4, 206.2, 208.2, 209.6, 212, 213.4, 218.9, 225, 225.8, 227.1, 227.3, 227, 227.1, 226.7, 229.2, 230.1, 230.2, 230.3, 231.3, 231.9, 232, 231.5, 231.2, 231.3, 234.6, 235.1, 241, 241.6, 242.7, 243.7, 243.1, 242.3, 241.9, 242.3, 244.5, 245.2, 245.1, 245.9, 246.8, 247.8, 248.3, 248.4, 248.4, 248.5, 250.7, 251, 251.3, 252.3, 253.3, 255, 255.3, 255.1, 254.8, 254.5, 256.2, 256.9, 255.6, 255.8, 257, 257.6, 257.3, 256.3, 255.7, 254.5, 256, 255.9, 254.6, 254.2, 255.2, 257, 257, 257.4, 257.3, 257.4, 259.8, 259.6, 256.9, 256.6, 257, 257.7, 258.1, 257.6, 257, 255.7, 256.8, 257.3, 256.2, 256.3, 257.3, 257.9, 258.3, 258.7, 257.6, 257.6, 259.4, 259.7, 257.5, 258.7, 259.9, 260, 261.3, 261.2, 260, 260.2, 262, 262.6, 261.7, 262.6, 264.6, 266.9, 268.7, 268.3, 266.9, 267.6, 269.9, 269.1, 268.8, 269.4, 271.8, 272.9, 273.6, 273.2, 272.3, 272.4, 274.5, 275.4, 276, 278.4, 279.8, 278.8, 278.5, 277.7, 276.8, 276.7, 278.7, 278.9, 278, 277.3, 279.4, 279.4, 280.1, 278.9, 278.5, 278.2, 280.2, 281, 277.9, 279.2, 279.8, 280.2, 280.3, 280.4, 279.4, 279.9, 281.9, 282.4)     t <- ts(v, frequency = 12)     o <- function(a){        accuracy(croston(t,h=12,alpha=a))["MAE"]      }     startTime <- proc.time()     optimize(o, 0.1, lower=0.1, upper=0.9, tol=0.1)     proc.time() - startTime      This yields an acceptable result however the optimization process is taking around 35 seconds. I really need to get it down to something more like 3 or 4 seconds. I tried playing around with the tolerance, but even with `tol=0.5` the optimization still takes 15 seconds (and no longer yields very useful results). I imagine the `optimize` function is using something like a Nelder-Mead algorithm. Is there maybe a better/faster optimization algorithm that could be used specifically to optimize the alpha parameter of the `croston` function?