title: Coefficient output in R for full model

text: I am having trouble understanding what exactly R is outputting when I look at the coefficient summary for the full model. For the lowbwt data, there are 4 variables: low, age, lwt, race (3 levels), and ftv. Low is the binary response variable. For example, if I am testing H0: beta_lwt=0 vs H1: beta_lwt=/=0, why don't I get the same estimates for lwt under h1 and m1? What is exactly is the difference between testing H0: beta_lwt = 0 vs H1: beta_lwt =/= 0 and testing M0: beta_age=beta_race=beta_ftv=0 vs M1: beta_age=beta_lwt=beta_race=beta_ftv=0?               > h1 <- glm(low ~ lwt, family=binomial, data=lowbwt)     > summary(h1)$coef                  Estimate  Std. Error   z value   Pr(>|z|)     (Intercept)  0.99831432 0.785290874  1.271267 0.20363370     **lwt         -0.01405826 0.006169588 -2.278639 0.02268856**     > m1 <- glm(low ~ age + lwt + race + ftv, family=binomial, data=lowbwt)     > summary(m1)$coef                    Estimate  Std. Error    z value   Pr(>|z|)     (Intercept)  1.29536575 1.071443130  1.2089916 0.22666607     age         -0.02382297 0.033729620 -0.7062924 0.48000632     **lwt         -0.01424461 0.006540719 -2.1778358 0.02941826**     race2        1.00389782 0.497858958  2.0164302 0.04375501     race3        0.43310843 0.362240236  1.1956387 0.23183760     ftv         -0.04930832 0.167239083 -0.2948373 0.76811818      Thanks so much in advance!


title: Are these 2 quasi poisson glm identical?

text: Approach 1:               lm.fit <- glm(response ~ 1, offset=log(lam), family="quasipoisson")     summary(lm.fit)      Approach 2: Feed summary.glm with a pearson dispersion.               lm.fit <- glm(response ~ 1, offset=log(lam), family="poisson")     summary(lm.fit,             dispersion=sum(residuals(lm.fit, type="pearson")^2)/df.residual(lm.fit))


title: How to convert predicted probabilities into predicted response?

text: How can I convert predicted probabilities of a logit model into predicted binary response ? Can I consider 0.5 as cut point to convert probabilities to binary variable (0,1). Or should I use binomial distribution to generate binary variable where predicted probabilities are used as success probability. By the way I am doing it for cross validation of the fitted model.


title: Comparing groups with multiple related variables

text: I have four groups with six subjects in each group. For each subject I have measured their mean blood pressure every minute for 3 hours. This means that I have 180 variables. The first 60 minutes represent "baseline" blood pressure, the next 60 minutes represent an intervention, and the last 60 minutes represents recovery. I want to determine if there is a difference in my baseline, intervention and recovery measurements. What is the best solution to this? My thought was to take an average of the variables in baseline, intervention and recovery periods for each subject. This would reduce my variables for each subject to three averages. Then I would compare these three averages across the groups. This means I use a subject average and compute a group average which I would compare. I'm not sure if this is the right approach? If relevant, I'm using SPSS for my statistics.


title: Convergence in probability of a product of RVs

text: I have found the following very useful theorem and I would appreciate some help comprehending it fully. **Theorem** Let $\{X_n \} $ be a sequence of random variables bounded in probability and let $ \{Y_n \} $ be a sequence of RVs which converge to $0$ in probability. Then $$X_n Y_n \rightarrow 0 \quad \text{in probability} $$ And its **proof** : For a given $\epsilon$ since $X_n$ is bounded in probability we can choose $N$ and a constant $B$ such that $$n \geq N \Rightarrow P \left[ |X_n| \leq B \right] \geq 1-\epsilon$$ Then $$ \overline{\lim_{n \to \infty}}P \left[|X_n Y_n| \geq \epsilon \right] \leq \overline{\lim_{n \to \infty}} P \left[ |X_n Y_n| \geq \epsilon, |X_n| \leq B \right]+ \overline{\lim_{n \to \infty}}P \left[|X_n Y_n| \geq \epsilon, |X_n| > B \right] \leq \overline{\lim_{n \to \infty}} P\left[|Y_n| \geq \epsilon/B \right] +\epsilon=\epsilon $$ If someone could help me understand how the last inequality is derived, given the information we have, I would be grateful. I understand the second to last inequality comes from the subadditivity property of limitsuperior so I merely need to comprehend the last one. Thank you in advance. **EDIT** : I think it is because $$P(A) \times P(B|A) \leq P(A) $$ and the complement of the event that is implied by the bound in probability.


title: Time series prediction with non-constant sampling interval

text: I have some data which can be modelled as such: each data sample $S$ is a series of discrete signal values $S(t_n) \in \{-1, 1\}$ measured at times $(t_{n, S})_{1 \leq n \leq N_S}$. The number of signal measurements $N_S$ per sample and the times at which the measurements were made all vary from sample to sample: $N_{S} \neq N_{S'}$ and $t_{n, S} \neq t_{n, S'}$ in the general case. In other words, each data sample is a series of "yes or no" answers asked at various times. I have some training data. Now, given a data sample $S$, I would like to predict the answer to the next question, in other words: the value of $S(t_{N_S + 1})$. Even better would be to have the probability of a "no" answer: $P(S(t_{N_S + 1}) = -1)$. I have no idea how to do this. Any hint? Starters? EDIT: if the signals were measured always at the same times ($t_{n, S} = t_{n, S'}$), and there were a constant number of signal measurements for each data sample ($N_S = N$), then I could produce for each data sample $S$ a vector of $N-1$ dimensions that would contain the $N-1$ first signal measurements: $S = [1, -1, -1, 1, 1, \dots]$. Measurement $S(t_N) \in \{-1, 1\}$ would be the label associated to sample $S$. Then, the problem would boil down to a classification/regression problem that can be solved e.g: with linear SVM. Unfortunately, the time differences between each signal measurement are important for the experiment.


title: Time series for Intermittent (1 month per year) data in R

text: I am going to analyze some data for an intermittent operation using R. Let's say I operate a Xmas tree stand from Black Friday to Christmas Eve every year. Let's say I operate 150 different Xmas tree stands every year. Let's say I have data for revenue-per-day for these 150 Xmas tree stands for the last 6 years. I want to analyze this data to figure out how to maximize my profit. Trends that I expect to exist include   * Daily trends in volume (higher volume on weekends, etc)   * Trend of # of days before Xmas (start low, max somewhere then decrease probably)   * Which day of week Xmas is probably has an affect   * Year over year (should show Great Depression of 2008, then recovery etc) I want to know how many stands to operate, and if I would do better with half the stands but move them around to different locations on different days to capture max revenue with minimum costs. For instance, I might find that for year when Xmas is on Wednesday, Location A should not be manned until 3 weeks before Xmas, should be manned on Thursday, Friday and Saturday for that week, and only on Tuesday and Saturday for the next 2 weeks. Bad example, but hopefully shows what I'm looking for. I want to try this in R. My limited understanding of time series is stumped when I have ~335 days between data sets which require a resolution of 1 day. Should I instead add dummy variables of day of week, days before Xmas, year and just try to fit on them? Any insight would be appreciated.


title: Predicting continuous variables from text features

text: I want to predict a continuous variable from text features. Lets say I have some student essays and I want to predict their quality, as measured by a human grader, using text features (mostly words they use). Linear regression is an obvious candidate, but if I have substantially more features than graded essays, this probably won't do well. If I wanted to classify them into good/bad, I might try the Naive Bayes classifier. I don't, but maybe I can draw inspiration from that. As I understand, Naive Bayes draws its power from assuming feature independence. Is there such a thing as Naive Multivariate Linear Regression, where you assume feature independence? I think this is the same as using Univariate Linear Regression for each regression coefficient. I would expect that to run into problems quickly, though. Is there something halfway between these two models? Putting a prior distribution on feature covariance that mostly expects conditional independence? Other models I should consider? Bayesian models preferred.


title: Is there a name for the variables in a regression equation that cause it to be overfit or near-overfit?

text: I call these "tweaker" variables -- the ones that pass some statistical test for inclusion, but whose contribution is minimal -- but I'm looking for a better or established label.


title: How can I interpret coefficients of categorical predictors in the negative binomial regression model?

text: I used some categorical variables as predictors to a negative binomial model. The dependent variable is numerical. I used glm.nb in R and the results show relative coefficients of one category respective to another category. Then I tried using lm.beta to standardize the coefficients, but still the results are relative. How can I interpret a positive versus negative coefficient here?


title: Software for ML/NLP/DM laboratory

text: suppose you might buy some software to be used in a ML/NLP/DM laboratory. What software package would you ask for? Let's say: MATLAB (with some toolboxes), SPSS, and what else? I know that there is a lot of free software one can use like R, Weka, Rapidminer, python packages, and so on. However, the above question targets commercial software. Thanks.


title: Mean and Median properties

text: Can somebody explain me clear the mathematical logic that would link two statements (a) and (b) together? Let us have a set of values (some distribution). Now, a) Median does not depend on every value [it just depends on one or two middle values]; b) Median is the locus of minimal sum-of-absolute-deviations from it. And likewise, and in contrast, a) (Arithmetic) mean depends on every value; b) Mean is the locus of minimal sum-of-squared-deviations from it. My grasp of it is intuitive so far. Thank you very much.


title: Newton Raphson Algorithm: negative semi definiteness

text: I would like to minimise the function $l(\theta|Y)$. Given the Newton's method below $$\theta^{(t+1)} = \theta^{(t)} - \left[l''(\theta\;|\;Y)\right]^{-1} l'(\theta^{(t)}\; | Y)\quad t = 0,1,...$$ where $l''(\theta^{(t)} | Y)$ is the observed Fisher Information matrix. The observed fisher information matrix (including the minus sign in front of it) would ideally be negative definite, in order to find the minimum. Correct me if I am wrong. My question is: what if the observed Fisher Information I obtained is **negative semidefinite**. What may happen is this case ? Would my algorithm fail miserably ?


title: Power calculation for cluster-level analysis in cluster randomized trials

text: I would like to solve for $\pi$1 in equation 7.14 of Hayes and Moulton's "Cluster Randomized Trials". I can't for the life of me remember how to do so. Here is a link to the equation. Here is the R script I wrote to calculate c as shown in equation 7.14:               # inputs     arms <- 2     alpha <- 0.05     tails <- 1     power <- .80     m <- 50             # observations per cluster     p0 <- 0.40          # true proportion in absence of treatment     p1 <- 0.60          # true proportion in presence of treatment     km <- 0.25          # between-cluster coefficient of variation within strata          num <- ((p0*(1 - p0))/m) + ((p1*(1 - p1))/m) + ((km^2)*((p0^2)+(p1^2)))     denom <- (p0 - p1)^2          x <- num/denom     c <- 2 + ((((abs(qnorm(alpha/tails))) + qnorm(power))^2)*(x))      Could anyone help me to re-write this to calculate p1 given c and all of the other variables? I'd be happy with a written formula that I can code if someone wants to help but not in R.


title: What is the best statistical method to compare case fatality rates of 17 years for 2 different states?

text: I have Case Fatality Rates (deaths per 100 cases) for 2 different states receiving different treatments for 17 years. What is the best statistical method to compare Case Fatality Rates ? The data are like this               Year   St.1 Cases  St.1 Deaths St.1 CFR St.2 Cases St.2 Deaths  State2 CFR       1994    1836        383         20.86     583        121         20.75     1995    1246        257         20.63    1126        227         20.16     1996    1450        263         18.14     896        179         19.98     1997    2953        407         13.78     351         76         21.65     1998    1161        149         12.83     1061       195         18.55     1999    2924        434         14.84     1371       275         20.06     2000    1729        169          9.77     1170       253         21.62     2001    1888        275         14.57     1005      199          19.80       2002    919         178          19.37      604      133         22.02        2003    865        142         16.42      1124      237         21.09     2004    1543       131          8.49      1030      228         22.14     2005    2887       336         11.64      6061     1500         24.75     2006    1484       108         7.28       2320     528          22.76     2007    1592       75          4.71       3024     645          21.33     2008    1920       53          2.76       3012     537          17.83     2009    1477       40          2.71       3073     556          18.09     2010    1534      26           1.69       3540     494          13.95


title: Difficulties with a Bayesian formulation of a model for human timing data

text: The Wing-Kristofferson model is a simple model of the behavior of a human trying to drum out a steady beat (that is, trying to mimic a metronome). Let $y_i$ be the $i$th interval between two drum beats then the model is: $$y_i \sim \text{Normal}(\mu + m_i - m_{i-1}, \sigma_y) $$ $$m_i \sim \text{Normal}(0, \sigma_m) $$ where $\mu$ is the mean interval of a time keeper and $m_i$ is the error of the $i$th drum stroke. What one wants is to estimate the parameters $\sigma_y$ and $\sigma_m$ ( while $\mu$ is considered more of a nuisance parameter). I've been trying to implement this model in a Bayesian framework by treating the $m_i$ as missing data. I've got two problems that I can't get my head around and that I would really appreciate some input on: **1. What would be appropriate "vague/objective" priors for the SD parameters $\sigma_y$ and $\sigma_m$?** I've been trying to to find "objective" priors for $\sigma_y$ and $\sigma_m$ that let the data dominate and that are relatively unbiased. Using simulated data I've tried a number of priors, non which have worked well so far. Flat priors on $\sigma_y$ and $\sigma_m$ tend to overestimate $\sigma_m$ but flat priors on $\log(\sigma_y)$ and $\log(\sigma_m)$ tend to pull $\sigma_m$ close to zero. **2. Is there a way to reduce the correlation between the parameters?** As $m_i$ are treated as missing data, and thus are parameters of the model, the number of parameters increase with the number of datapoints $y_i$. This wouldn't be too bad if it wasn't for the fact that both the parameters $m_i$ are heavily correlated both with each other and with $\sigma_m$. Another issue is that that $\sigma_m$ and $\sigma_y$ are heavily correlated. All these correlations makes the models hard to fit using MCMC frameworks such as JAGS or STAN. I've been tinkering around with different reparametizations of the model, but so far I haven't found anything that works. Still when the number of data points are low ($n < 60$) this isn't too much of a problem.


title: Convex Hull in R

text: I've to draw a polygon which will join the four vertexes on the plot in _R_. For this I need   minimum value of X which has minimum Y value.               X <- c(-62,  -40,   9,  13,  26,  27,  27)     Y <- c( 7, -14,  10,   9,  -8, -16,  12)     plot(x = X, y = Y)     abline(h = 0, v = 0, lty = 2.5, col = "green", lwd = 2)      ![enter image description here](http://i.stack.imgur.com/G0sv8.png) I've tough time to find four vertexes   * minimum of X and minimum Y (which is (-40, -14))   * minimum of X and maximum Y (which is (-62, 7))   * maximum of X and minimum Y (which is ( 27, -16))   * maximum of X and maximum Y (which is ( 27, 12)) in R.


title: How do I calculate the number of combinations, where a minimum number of elements must be different?

text: I have 54 names. I would like to know how many unique combinations of 6 names I can create, where at least 2 of the names are different from those in every other combination? My maths is extremely rusty but I understand how to calculate the number of unique combinations. I also understand how to calculate the number of unique combinations, where 1 name (or x names) must be include in every combination. I'm not at all sure how to go about solving this problem, except by brute- force for small numbers of names. An example based on 7 choose 3 would be 123,145,167,246,257,347,356. No combination contains 2 of the same numbers.


title: Importance of multiple linear regression assumptions when building predictive regression models

text: As far as I know, one can differentiate between two main goals of the regression analysis:   1. The goal is understanding causal relations between variables. Here, one has to check several common regression assumptions (main being linearity, normality, residuals zero mean, homoscedasticity, independence of errors). Collinearity needs to be analysed with care and removed if possible (either by deleting collinear predictors detected by variance inflation factors or applying principal component regression etc.).   2. The goal of the analysis is creation of predictive models. In contrast to the previous case, here one does not need to worry about collinearity at all (there are many online sources regarding this).  However, I am not sure what should one do with other "standard" regression tests, such as testing assumptions and hypotheses (p-values). Can one ignore these as well, just like collinearity, when the goal is building predictive models? In the other words: Could some "dummy" forward or backward feature selection be applied to select the features that maximize prediction R-square (using cross-validation or test sets), without taking care of regression assumptions or predictor significance?


title: What descriptive statistics are not effect sizes?

text: Wikipedia says > effect size is a measure of the strength of a phenomenon or a sample-based > estimate of that quantity. An effect size calculated from data is a > descriptive statistic that conveys the estimated magnitude of a relationship > without making any statement about whether the apparent relationship in the > data reflects a true relationship in the population. To understand it better, I was wondering what descriptive statistics are not effect size, except graphs and plots.


title: GLM conditional distribution from R GLM

text: I want to obtain the full distribution of a Gamma (or Inverse Gaussian) distributed $y_i$ given a vector of $\bar x_i$ that have been used in the linear predictor of a coefficient. Suppose also for the Gamma GLM I have used the log-link instead of the canonical link. Since both distributions are bi-parametric I know I can get the estimation of the mean parameter by `predict(glmObj, ..., type="response")` whichever distribution whichever link I have used. I'm not sure about the conditional variance. I know that $var\left(y_i\right)=\phi*V\left(\mu_i\right)$. My questions are:   1. Is it correct to estimate $\phi$, the dispersion parameter, as the square root of `glmObj$deviance/glmObj$df.residual`, regardless of the distribution and canonical link?   2. Is $V\left(\mu_i\right)$ dependent on the canonical link?


title: What is the proper way of converting ordinal values to numbers?

text: I have two painkillers and I have given them to two groups and recorded how the pain level changed in three categories: "Helped a lot", "Slightly better", "Did not help" Now I want to do a t-test and I was planning to convert each label to a number (say 3,2,1) and run a t-test on those numbers to see if any of the painkillers is significantly better than the other one. Now I have two questions:   1. If I go ahead with this ordinal to number conversion, what would be the methodical way of doing that? I can think of 3,2,1 or 8,4,2 and so on   2. Is there any test that allows me to do a hypothesis testing without converting to numbers? Thank you


title: Measuring spatial overlap in 2 different ways - paired t-test?

text: I have collected spatial data (GPS points and foraging heights) for a bunch of territorial animals. Until recently, these types of data were only analyzed in two dimensions (2D). However, 3D analysis is now available and I'm working on a paper that compares 2D and 3D estimates of spatial overlap. I define spatial overlap as being simply the area (2D) or volume (3D) of space overlapped divided by the total area or volume of the territory. This produces a proportion which I convert into a percentage. So basically for every bird I have a 2D percent overlap and a 3D percent overlap. I would like to compare these two different measures to see if one overestimates overlap compared to the other. It seems to make sense to me to use a paired t-test or the non-parametric equivalent (Paired-Sample Wilcoxon Signed Rank Test) as opposed to a normal t-test or Mann-Whitney U. Anybody have thoughts?


title: Chow's test and serially correlated model errors

text: how can one handle a time series with the Chow's test (in order to find a structural break) so that the assumption of independent model errors holds? I'm using the R function chow.test {gap}


title: Time series - correlation and lag time

text: I am studying the correlation between a set of input variables and a response variable, price. These are all in time series. 1) Is it necessary that I smooth out the curve where the input variable is cyclical (autoregressive)? If so, how? 2) Once a correlation is established, I would like to quantify exactly how the input variable affects the response variable. Eg: "Once X increases >10% then there is an 2% increase in y 6 months later." How should I implement this - in particular to figure out the lag time between two correlated occurrences? Example:![enter image description here](http://i.stack.imgur.com/F3Wtw.png) I already looked at: statsmodels.tsa.ARMA but it seems to deal with predicting only one variable over time. In scipy the covariance matrix can tell me about the correlation, but does not help with figuring out the lag time.


title: Predictive model & standardized variables

text: In a predictive model, I have standardized variables as predictors. Say I have to rescore the model on fresh data at some point in the future: do I use the means/stds as they were when I built the model to center and scale the new data, or do I use the means/stds as they are with the data I'm scoring. My take is to use the means/stds of the data I'm scoring, since I want the standardized variables to reflect distributions as they are at the time of scoring. Pros & cons of original means/stds vs. current means/stds? Thanks.


title: Decision network example

text: I am reading this example, but could you explain a little more. I don't get the part where it says "then we Normalize"... I know               P(sun) * P(F=bad|sun)  = 0.7*0.2 = 0.14     P(rain)* P(F=bad|rain) = 0.3*0.9 = 0.27      But where do they get               W    P(W | F=bad)     -----------------     sun     0.34     rain    0.66      ![enter image description here](http://i.stack.imgur.com/9ZWgQ.png) ![enter image description here](http://i.stack.imgur.com/fkFsh.png) ![enter image description here](http://i.stack.imgur.com/74J9o.png) Example from


title: Summary statistic for an error measurement

text: I have $n$ error measures that can be positive or negative values. Negative values are more important for my experiment than positive values. Of course, taking the mean of the $n$ values doesn't make sense because the mean could be $0$ and that won't reflect in one number what is happening. What summary statistic can you recommend me for this case? I would really appreciate any suggestion. Thanks in advance.


title: Logistic Regression and Design of experiment

text: I have seven factor and one response variable with 32 run and one replication i.e. 64 runs and I want to use logistics regression. I found in th analysis that 4 out 7 factors are significant. I want to check curvature. If curvature is not significant then should i make some other lay out of DoE only considering only significant factor? and If curvature is significant then should I make other layout with considering all factors?


title: Is there a theoretical basis for the shrinkage used in Boosted Regression Trees?

text: In Gradient Boosted Regression Trees, a shrinkage $\nu$ is often applied as: $$ f_t(x) \leftarrow f_{t-1}(x) + \nu h(x)$$ where $h$ is the regression tree learned by fitting the tree to the gradient. I've tried implementing this and found that this shrinkage is indeed necessary to prevent overfitting. The shrinkage required may vary by application but I found that anything greater than $\nu=0.01$ led to overfitting. Is there a theoretical justification for this kind of shrinkage? Are there more theoretically sound ways of regularizing GBRTs?


title: How to make a forest plot for single experimental group using event rate and sample size?

text: I am trying to create a forest plot for a small sample of case series. As the studies are case series there is only an experimental group, no control group. I would like to use event rate and sample size for each of the studies to create the graph. I am having difficulty as most software requires you input data for a control and experimental group. I am looking for a way to create the plot based soley off of experimental group event rate and sample size.


title: Sufficient statistical multivariate Gaussian

text: I have heard of this term supossing I have a matrix like $\begin{bmatrix} 1 & 12 & 13 & 14 \\ 2 & 22 & 23 & 24 \\ 3 & 32 & 33 & 34 \\ 4 & 42 & 43 & 44\\ 5 & 52 & 53& 54 \end{bmatrix}$ I count the number of rows $n$ wich is $5$ I compute the linear and quadratic sum like linear: $\sum_{i=1}^{n} x_{i}$ \begin{align} (&1+2+3+4+5,   12+22+32+42+52, \\ &13+23+33+43+53,14+24+34+44+54)\\ &=\left(15,160,165,170 \right ) \end{align} also the quadratic sum must be computed each as the quadratic sum of each column (I ommited because is a little long) Quadratic = $\sum_{i=1}^{n} x_{i}x_{i}^{T}$ Could you explain why are they call sufficient statistical could you point to some tutorial with more details? I have heard that this vector (linear sum) and the quadratic matrix (the quadratic sum of points) summarize the original data $x$


title: Designing a stopping rule using a hypothesis test

text: I am working on a stopping rule for an optimization algorithm that produces an upper bound and lower bound for the objective value of an optimization problem. In my case, the lower bound is deterministic, but the upper bound is an estimate derived from $N$ data points $UB_1, UB_2... UB_N$ with mean $\widehat{UB}$ and standard deviation $\hat{\sigma}_{UB}$, where: $$\widehat{UB} = \frac{1}{N} \sum_{i=1}^N UB_i$$ $$\hat{\sigma}_{UB}^2 = \frac{1}{N} \sum_{i=1}^{N} (UB_i - \widehat{UB})^2$$ Theoretically speaking, the algorithm needs to stop when the upper bound is equal to the lower bound. The stopping rule that I have in mind is to stop when $\frac{\widehat{UB}}{LB} = 1$. More formally, I would like to test: $$H_0: \frac{\widehat{UB}}{LB} = 1$$ $$H_A: \frac{\widehat{UB}}{LB} > 1$$ And stop when I fail to reject $H_A$. Broadly speaking, I would like to control for both type I and type II errors. That is, I would like to specify the probability that I am stopping too early (i.e. $H_0$ is accepted given $H_A$ is true) and also specify that I am stopping too late (i.e. $H_0$ is rejected given that $H_0$ is true). I'm wondering what the exact formulation should look like in this case? Or what type of test to do. I'm open to any ideas you may have, including a different test altogether so long as we can specify the probabilities I described above in some way.


title: Calcualting R2 in mixed models useing Nakagawa & Schielzeth's (2013) R2glmm method.

text: I have been reading about calculating R2 values in mixed models and after reading the R-sig FAQ, other posts on this forum (I would link a few but I don't have enough reputation) and several other references I understand that using R2 values in the context of mixed models is complicated. However, I have recently came across these two papers below. While these methods do look promising (to me) I am not a statistician, and as such I was wondering if anyone else would have any insight about the methods they propose and how they would compare to other methods that have been proposed. > Nakagawa, Shinichi, and Holger Schielzeth. "A general and simple method for > obtaining R2 from generalized linear mixed‐effects models." Methods in > Ecology and Evolution 4.2 (2013): 133-142. > > Johnson, Paul CD. "Extension of Nakagawa & Schielzeth's R2GLMM to random > slopes models." Methods in Ecology and Evolution (2014). The is method can also be implemented using the r.squaredGLMM function in the MuMLn package which gives the following description of the method. > For mixed-effects models, R^2 can be categorized into two types. Marginal > R^2 represents the variance explained by fixed factors, and is defined as: > >> R_GLMM(m)^2 = (σ_f^2) / (σ_f^2 + sum(σ_l^2) + σ_e^2 + σ_d^2 > > Conditional R^2 is interpreted as variance explained by both fixed and > random factors (i.e. the entire model), and is calculated according to the > equation: > >> R_GLMM(c)^2= (σ_f^2 + sum(σ_l^2)) / (σ_f^2 + sum(σ_l^2) + σ_e^2 + σ_d^2 > > where σ_f^2 is the variance of the fixed effect components, and sum(σ_l^2) > is the sum of all variance components (group, individual, etc.), σ_l^2 is > the variance due to additive dispersion and σ_d^2 is the distribution- > specific variance. In my analysis I am looking at longitudinal data and I am primarily interested in variance explained by the fixed effects in the model               library(MuMIn)      library(lme4)          fm1 <- lmer(zglobcog ~ age_c + gender_R2 + ibphdtdep + iyeareducc + apoegeno + age_c*apoegeno + (age_c | pathid), data = dat, REML = FALSE, control = lmerControl(optimizer = "Nelder_Mead"))          #Jarret Byrnes (correlation between the fitted and the observed values)     r2.corr.mer <- function(m) {        lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))        summary(lmfit)$r.squared     }          r2.corr.mer(fm1)     [1] 0.8857005          #Xu 2003     1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))     [1] 0.8783479          #Nakagawa & Schielzeth's (2013)     r.squaredGLMM(fm1)           R2m       R2c      0.1778225 0.8099395       Thanks.


title: How are propensity scores different from adding covariates in a regression, and when are they preferred to the latter?

text: I admit I'm relatively new to propensity scores and causal analysis. One thing that's not obvious to me as a newcomer is how the "balancing" using propensity scores is mathematically different from what happens when we add covariates in a regression? What's different about the operation, and why is it (or is it) better than adding subpopulation covariates in a regression? I've seen some studies that do an empirical comparison of the methods, but I haven't seen a good discussion relating the mathematical properties of the two methods and why PSM lends itself to causal interpretations while including regression covariates does not. There also seems to be a lot of confusion and controversy in this field, which makes things even more difficult to pick up. Any thoughts on this or any pointers to good resources/papers to better understand the distinction? (I'm slowly making my way through Judea Pearl's causality book, so no need to point me to that)


title: Column-wise matrix normalization in R

text: I would like to perform column-wise normalization of a matrix in R. Given a matrix `m`, I want to normalize each column by dividing each element by the sum of the column. One (hackish) way to do this is as follows:               m / t(replicate(nrow(m), colSums(m)))      Is there a more succinct/elegant/efficient way to achieve the same task?


title: Probability Distribution - which to use Normal or Hypergeometric

text: Good evening everyone, I am attempting a self-study question on probability and am faced with one which really quizzed me. > A study found that a family typically spends between \$500 and \$4500 a > month on expenses. Suppose the money spent is uniformly distributed between > these amounts. If we randomly select 10 families, what is the probability > that at least two families spend more than \$3000 a month on expenses. My approach to this question is that I present > x = families spending > \$3000 a month I then proceed on to form a normal distribution as follows: ![enter image description here](http://i.stack.imgur.com/9nPhc.jpg) I was thinking of using the Hypergeomatric distribution for this, but also thought that the normal distribution made sense. I am just puzzled how to move on. From what I understand, the normal distribution never touches the X-axis, so that will mean that it will be hard for me assert that 500 - 4500 is the range, and thus I was skeptical to convert it to the standard normal.


title: 100mn rows of "events" with 30 or so attributes, want to understand what affects one particular attribute - how? which software?

text: I have about 50-100 million rows of data of interesting "events". Each of these rows has about 32 attributes. One of these attributes is how much money we made :-) What's the best way to make sense of this data? I've heard that Decision trees might be a great way to understand this. To be clear, I don't have prediction needs for the future, I just want to know which attributes are highly correlated with the money being made. I would also like to know what kind of software should I be using to do this? Should I create a SQL database? My current data is in CSV, I've converted it to JSON. Also, is this a machine learning problem?


title: Question on Probability & Distribution

text: Good evening everyone, I am currently working on a self-study question: > There is 15 marbles in a jar. Six are blue, while the rest are green. Five > random marbles are selected. What is the probability that the number of blue > marbles selected is more than the number of green marbles selected. This question had really quizzed me. The following is my best attempt, however still falling short of going anywhere near the right answer. > P(Blue) = 0.4, P(Green) = 0.6. **Find P(Blue > Green)**. Appreciate any pointer and guidiance on how I can proceed please.


title: Self Study question on Conditional Probability

text: I am currently doing a self-study on Conditional Probability. I was faced with a question where I was provided with $P(a)$, $P(b)$, $P(c)$, $P(a \mid d)$, $P(b \mid d)$ and $P(c \mid d)$. The full context of the question: > A departmental store reports that 30% of payment is collected in cash, 60% > in credit card and 10% in debit card. 20% of cash purchases, 90% of credit > card purchases and 80% of debit card purchases are for more than 200 dollars > in purchases. What the probability that Sue paid cash if she purchased a new > bag that costs $98. Thus, I took $a = \mbox{payment in cash}$, $b = \mbox{payment in credit}$ and $c = \mbox{payment as debit}$. $P$ represented the probability that payments went above 200 dollars. I was attempting to find out what is $P(d)$. My attempt (which did not seem to yield the right answer): $$ P(a)P(a \mid d) + P(b) P(b \mid d) + P(c)P(c \mid d) $$ Appreciate some guidance please. **Edit** Based on the latest guidance, I still failed to answer the question with the final answer being 1.43567, however still deferring from the original answer of 0.750. My answer seemed grossly wrong. My solution: > Let $a = \mbox{cash}$, $b = \mbox{credit}$, $c = \mbox{debit}$ and $d = > \mbox{payment amounts more than 200 dollars}$. > > Therefore $P(a) = 0.3$, $P(b) = 0.6$, $P(c) = 0.1$, $P(d \mid a) = 0.2$, > $P(d \mid b) = 0.9$ and $P(d \mid c) = 0.8$. > > I will first attempt to find \begin{align*} P(d) &= P(a)P(d \mid a) + > P(b)P(d \mid b) + P(c)P(d \mid c) \\ &= 0.3 \cdot 0.2 + 0.6 \cdot 0.9 + 0.1 > \cdot 0.8 = 0.68. \end{align*} > > Then, I will attempt to find the value of $P(a \mid d)$ using the bayes > theorem which results in $(0.2+0.3)/0.68 = 0.73529$. > > Next, to find $P(a \mid d^c)$, we use the formula $P(d) = P(a)P(a \mid d) + > P(d^c) P(a \mid d^c)$. $0.68 = 0.3 \cdot 0.73529 + 0.32 P(a \mid d^c) = > 1.43567$. Since the probability can never be more than $1$, my answer looks terrible wrong. Appreciate and guidance please.


title: Fractional factorial design. only main effects considered

text: I have 15 factors which I want to optimize. Instead of running a fractional factorial design consisting of 32 runs, if I run a design consisting of 16 runs, will I only get main effects or any second order interactions too? I am sure that there are no third order int. but maybe 1 second order interaction will be present. If an interaction is significant, will the terms involved in the interaction be automatically significant although the interaction will not show up in the 16 run fractional factorial design?


title: What is the best way to analyze additional attributes in a market basket analysis?

text: I am performing a market basket analysis of customers and the products they purchased. I did an a priori analysis on all of the transactions on the products they purchased to determine which items are likely bought together. I have additional attributes (industry, geography, etc.) about the customers (B2B), and I would like to find purchasing trends. For instance, do customers in the financial services industry buy a group of products more frequently than customers in the retail industry (or all customers)? I can think of three ways to do this, but I am not sure which one is correct:   1. Add the attributes as an additional product in the transaction (so to speak), so the a priori algorithm itself will determine if there is a relationship between the _industry_ product and the _purchased_ products   2. Run the a priori analysis on pre-filtered customer transactions based on an attribute (e.g., run it on only retail customer transactions)   3. Run some other statistical analysis (chi-square?) to determine if a customer group differs


title: SVM classifier - grid search on C and gamma no affect

text: I'm working to find the best C and gamma values for the RBF on my dataset using an SVM classifier. Specifically, I'm using LibSVM & MATLAB.               g=[2e-5,2e-4,2e-3,2e-2,2e-  1,2e0,2e1,2e2,2e3,2e4,2e5,2e6,2e7,2e8,2e9,2e10,2e11,2e12,2e13,2e14,2e15];     C=g';     matrix_gC_aust=zeros(21,21);     matrix_gC_germ=zeros(21,21);     %%crossvalidate SVM test     for i=1:21         for j=1:21             model_1=svmtrain(aust_train(:,15),aust_train(:,1:14),['-g',num2str(g(i)), '-c',num2str(C(j)),'-h 0 -t 2']);             [pre1,acc1,decv1] = svmpredict(aust_test(:,15),aust_test(:,1:14),model_1);             model_2=svmtrain(germ_numeric_train(:,25),germ_numeric_train(:,1:24),['-g',num2str(g(i)), '-c',num2str(C(j)),'-h 0 -t 2']);             [pre2,acc2,decv2] = svmpredict(germ_numeric_test(:,25),germ_numeric_test(:,1:24),model_2);             matrix_gC_aust(i,j)=acc1(1,1);             matrix_gC_germ(i,j)=acc2(1,1);         end     end      The results from this show that altering C and gamma do not change the accuracy of the SVM classifier. Is this typical/possible? Has anyone else experienced this? One of my datasets is 1000 samples with 24 features, and the other is 690 samples with 14 features. I'm not sure how these sample/feature ratios compare to others. Thanks.


title: Reliability in Elicitation Exercise

text: I am trying to calculate the reliability in an elicitation exercise by analysing some test-retest questions given to the experts. The experts elicited a series of probability distributions which were then compared with the true value (found at a later date) by computing the standardized quadratic scores. These scores are the values that I am using to calculate the reliability between the test-retest results. Which reliability method would be appropriate here? I was looking mostly at Pearson's correlation and Chronbach's alpha (and got some negative values using both methods) but I am not sure this is the right approach. * * * **UPDATE:** Background information The data were collected from a number of students who were asked to predict their own actual exam mark in four chosen modules by giving a probability distribution of the marks. One module was then repeated at a later date (hence the test-retest exercise). Once the exam was taken, and the real results were available, the standardized quadratic scores were computed. These scores are proper scoring rules used to compare assessed probability distributions with the observed data which might be known at a later stage. The probability score _Q_ is defined as: ![Quadratic score](http://img717.imageshack.us/img717/9424/chart2j.png) where _k_ is the total number of elicited probabilities and _j_ is the true outcome. My question is which reliability method would be more appropriate when it comes to assessing the reliability between the scores of the repeated modules? I calculated Pearson's correlation and Chronbach's alpha (and got some negative values using both methods) but there might be a better approach.


title: Is a variable significant in a linear regression model?

text: I've got a linear regression model with the sample and variable observations and I want to know:   1. Whether a specific variable is significant enough to remain included in the model.   2. Whether another variable (with observations) ought to be included in the model. Which statistics can help me out? How can get them most efficiently?


title: What would the calculated value of the standard deviation of a uniform distribution be?

text: A colleague wants to compare models that use either a Gaussian distribution or a uniform distribution and for other reasons needs the standard devation of these two distributions to be equal. In R I can do a simulation...               sd(runif(100000000))     sd(runif(100000000,min=0,max=2))      and see that the calculated standard deviation is likely to be ~.2887 * the range of the uniform distribution. However, I was wondering if there was an equation that could yield the exact value, and if so, what that formula was.


title: How do I perform this complicated ANOVA type analysis in R?

text: I have the following type of dataset:                       |          |             |            variable_r         |     subject |  gender  |  age_group  |     Cond_1    |     Cond_2    |     --------|----------|-------------|---------------|---------------|        1    |    m     |      1      | r (A) | r (B) | r (A) | r (B) |        2    |    f     |      2      | r (A) | r (B) | r (A) | r (B) |     ...        8    |    f     |      2      | r (A) | r (B) | r (A) | r (B) |      So two genders, two age groups, two conditions (Cond_1 and Cond_2) under which the experiment was done and two ways the subjects were prompted (A and B). r is the numerical result from each experiment. So two within-subject variables (prompt A/B and Cond 1/2) and two between-subject variables (age group 1/2 and gender m/f) (right?). I should calculate the statistically significant effects of each variable and their interactions. How can I do this in R (or Python)? My googling found a lot of information about different types of ANOVA analyses, but I wasn't able to apply that information to my case. Thanks! e: the subjects were tested 4 times


title: R statistics Mixed ANOVA but outcome Poisson distributed - which model in R?

text: I have a mixed ANOVA with one between and one within factor: between factor: control versus treament within factor: ad1, ad2 (measures click rates on 2 ads)               aov(repdat~(within*between)+Error(subjcts/(withincontrasts))+(between),data.rep)      However my outcome repdat has a Poisson distribution and 0s (since sometimes users just dont click). How do I model this in R? I considered: 1.) using Zelig package with model Poisson, but I dont know how to properly write the mixed model formula here. Is it the same as for the mixed ANOVA? 2.) using lme4 package with model Poisson. Can I use both a between and within factor or even treat "ads" as a random variable? It would be great if someone could quickly show me example code for both 1. and 2. thanks


title: How does R handle missing values in lm?

text: I'd like to regress a vector B against each of the columns in a matrix A. This is trivial if there are no missing data, but if matrix A contains missing values, then my regression against A is constrained to include only rows where all values are present (the default _na.omit_ behavior). This produces incorrect results for columns with no missing data. I can regress the column matrix B against individual columns of the matrix A, but I have thousands of regressions to do, and this is prohibitively slow and inelegant. The _na.exclude_ function seems to be designed for this case, but I can't make it work. What am I doing wrong here? Using R 2.13 on OSX, if it matters.               A = matrix(1:20, nrow=10, ncol=2)     B = matrix(1:10, nrow=10, ncol=1)     dim(lm(A~B)$residuals)     # [1] 10 2 (the expected 10 residual values)          # Missing value in first column; now we have 9 residuals     A[1,1] = NA       dim(lm(A~B)$residuals)     #[1]  9 2 (the expected 9 residuals, given na.omit() is the default)          # Call lm with na.exclude; still have 9 residuals     dim(lm(A~B, na.action=na.exclude)$residuals)     #[1]  9 2 (was hoping to get a 10x2 matrix with a missing value here)          A.ex = na.exclude(A)     dim(lm(A.ex~B)$residuals)     # Throws an error because dim(A.ex)==9,2     #Error in model.frame.default(formula = A.ex ~ B, drop.unused.levels = TRUE) :      #  variable lengths differ (found for 'B')


title: Test for biological replicates

text: I am performing this simple experiment: I have one variety of grass and 8 different fungi (say #1 to #8). I am going to put 10-20 grass plants in each one of 18 containers, and then I will put each one of the fungi in 2 containers and let 2 with no fungus ("negative control"). After some time I will count how many plants have been killed by each one of the fungi in each one of two replicates for that fungus. I want to compare the results for each one of the fungi with respect to the negative control. I think I can do that using the Fisher exact test, but I don't know if I have to combine the counts of the two replicates into a single count (like using a single container with all the plants together) or if I can keep those count separated. Also, I would like to know if I have to run the test doing pairwise comparisons (controls vs #n) or if I should run a single test considering all the results at the same time. May I do that also with the Fisher exact test? Do I have to make a multiple comparison adjustment? Thank you for any help. I am very confused about this.


title: Confidence Interval for $\eta^2$

text: Following my question here, I am also looking at the difference between males and females and I have conducted linear regression in a general linear model setup for this purpose. My effect size for this part of the project is $\eta^2$ (which is 0.25). $\eta^2$ is equivalent to $R^2$ (which is the appeal for me in this case). I am wondering how do I calculate the confidence interval of $\eta^2$. Does anyone know of any online calculators that can I can use?


title: Checking $R$ functions arguments

text: I was wondering how to get arguments for _R_ functions. Is there any _R_ function which can be used to get all arguments for a certain function? As _R_ function glm has the following arguments               glm(formula, family = gaussian, data, weights, subset,     na.action, start = NULL, etastart, mustart, offset,     control = list(...), model = TRUE, method = "glm.fit",     x = FALSE, y = TRUE, contrasts = NULL, ...)      How these arguments can be obtained on _R_ console? Thanks


title: How can I calculate the standard error of a proportion?

text: I made a comparison of hatch success between 2 populations of birds using R's `prop.test()` function:               prop.test(c(#hatched_site1, #hatched_site2),c(#laid_site1, #laid_site2))      It gave me the proportions of each site as part of the summary. How can I calculate the standard error for each proportion?


title: Binary Logistic Regression Methods

text: I have data sample size of almost 15,000 cases. The dependent variable is a dichotomous variable stating whether the patient has the disease or not, Yes=1, and No=0. I have 12 more independent variables which are continuous as well as dichotomous. My question here is, before I apply the logistic regression, I need to know which parameters have a great impact on the DV; how do I do that? Secondly, there are 7 different methods of using binary logistic regression like enter, forward conditional, forward LR, etc; what is the basic difference between these? I couldn't find any matter on this on the internet. Thanks in advance


title: Are there terms to distinguish between the two types of outliers?

text: I hate talking about "outliers," because I view that term as encompassing two entirely different concepts. The first is when it refers to data that was incorrectly recorded or measured. For instance, a 400000 bedroom house for 4 dollars or a person that is 180 inches tall. There are various techniques to handle these that are outside the scope of my question. The second type of outlier refers merely to extreme points on the distribution. These might include things like Babe Ruth's home run totals or the net worth of Bill Gates. Although these can have adverse effects on certain types of analyses, they are nevertheless legitimate datapoints. They also have their own techniques that I'd rather not get into here. Are there any terms that are used to distinguish between the two cases of outliers?


title: Expectation Maximization, dice example, always converging in second iteration

text: I am simulating two loaded dice and trying to estimate individual die prior probabilities and probability mass functions for each of them using the EM algorithm. Below is my Matlab code. Likelihood and pmf values remain constant second iteration onwards. Why is that?               % dice pmfs pmf1 and pmf2          A = [1 2 3 3 2 1];     B = [5 4 3 2 1 4];     U = [1 1 1 1 1 1];          pmf1 = A./sum(A);     pmf2 = B./sum(B);     pmfU = U./sum(U);          for i = 1:6         cdf1(i) = sum(pmf1(1:i));         cdf2(i) = sum(pmf2(1:i));         cdfU(i) = sum(pmfU(1:i));     end          % generate N random numbers from dices     N = 1000;     dicePriors = [0.7 0.3];          for i = 1:6         if i == 1             rv1 = pmfU(i)*rand(1,ceil(dicePriors(1)*N*pmf1(i)));             rv2 = pmfU(i)*rand(1,ceil(dicePriors(2)*N*pmf2(i)));         else             temp = cdfU(i-1) + pmfU(i)*rand(1,ceil(dicePriors(1)*N*pmf1(i)));             rv1 = [rv1 temp];             temp = cdfU(i-1) + pmfU(i)*rand(1,ceil(dicePriors(2)*N*pmf2(i)));             rv2 = [rv2 temp];         end             clear temp;     end          % plot individual histograms     % bar(hist(rv1,6) ./ sum(hist(rv1,6)))     % figure;     % bar(hist(rv2,6) ./ sum(hist(rv2,6)))          rv = [rv1 rv2];     N = length(rv);          for i = 1:N         rv(i) = ceil(6*rv(i));     end          % EM estimation of pmfs          %  Initialization     pmf11 = [3 3 1 1 1 1];     pmf22 = [8 1 1 1 1 8];          pmf11 = pmf11./sum(pmf11);     pmf22 = pmf22./sum(pmf22);          dicePriors1 = [0.5 0.5];          iterNo = 10;          for iter = 1:iterNo     iter     likelihood(iter) = 0;          for i = 1:N        num = rv(i);        prob_dice1_given_num(i,num) = pmf11(num)*dicePriors1(1) / (pmf11(num)*dicePriors1(1) + pmf22(num)*dicePriors1(2));        prob_dice2_given_num(i,num) = pmf22(num)*dicePriors1(2) / (pmf11(num)*dicePriors1(1) + pmf22(num)*dicePriors1(2));        likelihood(iter) = likelihood(iter) + log(pmf11(num)) + log(pmf22(num));     end          for j = 1:6         pmf11new(j) = sum(prob_dice1_given_num(:,j))/sum(sum(prob_dice1_given_num));         pmf22new(j) = sum(prob_dice2_given_num(:,j))/sum(sum(prob_dice2_given_num));     end     for j = 1:6         pmf11new(j) = pmf11new(j)/(sum(pmf11new));         pmf22new(j) = pmf22new(j)/(sum(pmf22new));     end          pmf11 = pmf11new;     pmf22 = pmf22new;          dicePriors1(1) = sum(sum(prob_dice1_given_num)) / (sum(sum(prob_dice1_given_num)) + sum(sum(prob_dice2_given_num)));     dicePriors1(2) = sum(sum(prob_dice2_given_num)) / (sum(sum(prob_dice1_given_num)) + sum(sum(prob_dice2_given_num)));               dicePriors1     pmf11     pmf22     likelihood          end


title: Where is the uniform distribution with one parameter ($U(\theta, k \theta)$) useful for modelling?

text: I recently came across the distribution $U(\theta, k \theta)$ (where k is known) in the context of statistical theory (as a nice toy example for finding MLE and the likes). However, I was wondering if there are real use cases where one might wish to use this distribution for modelling something. But I couldn't seem to find any paper using this (or similar) distribution for actual data. My question is if you happen to know of such a use for this distribution (or one similar to it). And more generally, how can one search for it effectively?


title: What types of statistical analysis should be done for testing ISP connections?

text: I'm in the middle of analyzing data taken from testing two different ISP connections (cable and DSL). I've been running randomized testing on the two connection using three different tests which all measure the same basic thing: connection time (ms) and Download Speed (Mbit/s). Some of the tests also record Upload Speed (Mbit/s) and latency (ms). Since the testing is randomized, and due to outside forces, I don't have same number of samples per connection. What are some statistical analysis that I can perform to prove both the quality (Max Down/Up speed) and reliability (failures or latencies) of the connection using different sample sizes?


title: Statistic for count process and measurements

text: I am facing some problems concerning which statistical approach to use for my measurements. I have a sensors which counts how many times an events occured. I had to characterize two lots of sensors, and I tested 10 devices per lot. Each time I tested a sensor I collected more then 6000 counts in order to reduce the statistical error on the counts to <1%. I made 3 measurements on each sensor to verify the repeatibility of the measurements (which is less then 10 %). For each one of the 10 sensors I made the mean of the 3 measurements and then calculate the mean on 10 and the standard deviation for the lot. I made the same for the second lot. 1) is this process correct in your opinion? 2) I am using always a normal distribution, should I use a t student for low samples (n=10) ? 3) With which interval of confidence I can state that all the sensors in one lot follow the normal distribution with mean and average calculated as before? Thank you to all!!


title: How to compare Arellano-Bond models?

text: I want to compare different Arellano-Bond models, but I don't know exactly how to do it and which statistics to consider. For instance, should I compare $\chi^2$ of Sargan tests? What other options exist? PS: I am using `pggm` in `R` to fit models.


title: Multiple random effects in a generalized linear mixed model

text: I'm looking for an example of a data set that can be fitted using a GLMM with at least 3 separate random effects. I've taken a look at a few books on GLMM's including the one by McCulloch and Searle, but their examples always stop at 2 random effects. I don't know why I'm having so much trouble wrapping my head around a practical example of 3 or more random effects but I am. In the national youth tobacco survey for example, counties were the PSU's. Within sampled counties, schools were then randomly selected along with students. So if I understand it correctly it would make sense to perhaps treat both counties and schools as random. I suppose if they sampled classrooms within each school we could then have 3 nested random effects - one for each of county, school and classroom. Is that correct? Any help would be greatly appreciated.


title: Performing a Wald's group sequential test in R (for complex hypothesis)

text: The settings: we have a sequence of numbers generated from a normal distribution with $\mu=1;\sigma^2 = 1 $. And with every new observation collected, we'd like to test for: $H_0:\mu_0=0$ (two-sided). From what I've learned in school, there is the Wald's group sequential for when comparing two simple hypothesis, and I would like to know what is available for complex hypothesis. What I am looking for is:   1. Some reference text (article/book) that describes the solution for this (or a similar) case.   2. If there is some implementation in R for performing such a test. Thanks.


title: Is this outlier merely unusual or does it signify manipulation?

text: These are actual reported costs (60) for a series of projects for two years, say USD (though actually not) totalling $57,975,403.24:               (274,680.90)    3,900.13    103,999.61  512,918.71  1,079,734.06      (97,570.45)     4,218.45    110,499.74  516,501.94  1,249,999.48      (43,261.81)     9,951.36    116,287.74  531,650.58  1,398,872.64      (41,314.45)     14,191.48   124,999.48  551,420.65  1,999,996.52      (34,843.48)     34,999.09   145,808.39  574,996.39  2,249,270.19      (1,305.16)      35,218.84   215,999.87  605,869.55  3,255,268.13      (634.19)        42,007.48   227,045.03  715,369.04  3,758,841.42      0.00            47,045.03   262,369.29  749,881.94  3,897,417.42      60,000.00       54,999.09   347,790.19  794,996.65  3,925,001.81      960,000.00      74,175.23   359,948.78  800,321.94  4,266,898.32      39,999.74   *   77,998.97   359,998.06  824,731.87  7,494,344.91      59,999.48   *               404,998.07  849,997.80  10,646,324.13                                              919,938.97       Some negatives, the occasional zero and a high incidence of round thousands is to be expected. The two values asterisked are probably a single project. The values are all plausibly accurate but I would like to know whether I have grounds for suspecting the values have been manipulated as they seem to end just short of a round thousand far more often than I would expect: ![Distribution chart](http://i.stack.imgur.com/v4Eay.gif)! In particular, 15 of the 60 are less than USD 4 short of the next round thousand. Excluding 7 negatives, 1 zero and 2 round thousands, 50 values spread across ten centuries averages out at about 5 per band. A very slight propensity to fall just short of a round number is to be expected – either as a result of something like a currency conversion approximation or the rare purchase of a low value consumer item (such as a book or shoes that mostly seem to be priced $x.99) but there may also be a slight propensity to just exceed a band.   1. In this context, is the fact that the top band contains 20 examples at all significant?   2. And if not is there a sample size at which such an anomaly would be statistically significant (say 95%)?   3. Is there something like Benford’s Law that applies to trailing rather than leading digits?  * * * Update: I now have most of the values that aggregated as above. Of the 20 values that seemed odd to me, I do not have details for 4 but the sums of the components for 9 of the 20 _are_ round thousands. In effect, based on actual underlying values, 9 of the 20 in the range `900>999.99` should be moved to the range `Roundk` and hence change the chart appearance to something I would deem plausible. I have also investigated the distribution of the last integer in the details that were represented by the remaining 7 in the range 900>999.99 which excluding zero (as likely to be favoured), is: ![CV51022 example](http://i.stack.imgur.com/GCDbQ.gif) There were 973 data points of which 157 were ‘0’ and, though not quite horizontal, a linear trend line shows only relatively minor departure from a simple average of just under 91 instances per final integer. So, my conclusion is that the reporting of aggregate values (which takes place at intervals) may for one or more reporting periods have been rounded*. This I think confirms that the results were not as should have been expected but, in this case, provides a wholly innocuous explanation. Therefore I’d suggest A1: Yes, A2: n/a, A3: No – other than random might be expected (uniform distribution). *Current Total reported as Previous total + Current Period change, rather than as Current Cumulative Total.


title: code for ordered probit model

text: I have a data set with 7 predictor variables and one dependent variable. The dependent variable has 4 categories so it's not binomial. I need to fit a probit model. I need codes for probit model in R, SPSS or Matlab.


title: Regression with Time series model

text: I wanted to see if the different kinds of supports to small enterprises have contributed to the development/growth of small enterprises using 15 years period data (6 years before implementing these programs & 9 years after it). I used comparative research design & applied linear time series regression. At 10% significance level, the result gave me the following result:                           -R squared= 0.35, Adjusted R squared= -0.03                     -some coefficients are positive & some are negative                   -p-values for all the four independent variables' coefficients are >                      significance level                 -t-stat= 0.5, 0.2, -0.2, & 1.4       Therefore, 1) How can I interpret these results? 2) Can I use < 90% confidence interval to see at what % C.I. the support programs become significant to have effect on the growth? 3) Am I expected to provide explanations as to why it is unusually negative or insignificant? 4) What type of tests should I conduct in this time series model? 5) And finally, which software may I use- EVIEWS, SPSS or Stata? You may answer all or any of these 5 questions. -Here is the data: -               No. of     credit      Training    marketing      Business             Organizing services                 (in $)     (in No of   services       development service    (in No. of firms)                              Trainees) (in No.of firms)   (in No.of firms)     employees               45,351    44,062,000    2,990    1,322          1,920                19,134                   107,283   109,560,000   40,435   30,662         5,353                74,941                 20,131    36,050,000    8,375    7,500          4,848                8,225            64,945    83,600,000    16,538   18,779         7426                 6406     223,242   157,100,000   10,894   25,876         2479                 9909     134,611   229,912,620   33,547   58,687         6,642                27,971                                                                            39092     90,529    271,924,783   70,157   18,543         2,300                                       63,759    368,691,717   12,148   12,443         3,954                67,238     110,619   491,586,029   76,461   25,414         4,168                70,600


title: automating procedures in a complex data set

text: I have a data set on hands with 5,220 cases and I'd like to automate some procedures (but I have almost no programming knowledge). The data has some continuous variables that are grouped by 2 others: the name of species and the locality where they were collected. So, the samples are defined as 'each species on each locality'. For every sample I'd like to do multiple imputation (when applicable), test for the presence of outliers, standardize the variables, correct some species abundances, save individual samples to tab delimited text file, and assemble each individual sample (now, without NAs and outliers, corrected abundances, and with the new standardized variables) into a single data set. That task is pretty complex to me, since my programming knowledge is poor (and my free time to learn R programming is sparse). Could someone help me with that (I could provide you the data set and the script I have written to do that, sample by sample [ouch!])? [Below is an example is the codes I've used to accomplish my goals, sample by sample, which can exemplify the complexity of the procedures:               #Subsetting the data (v1-v11 are continuous "predictors"): species 1 at locality 1 (all data [5520 cases] are on a vector called 'morfo')     sp1.loc1<-morfo[which(spps=="sp1" & taxoc=="loc1"),] #getting only the observations of sp1 (species 1) at loc1 (locality 1)     str(sp1.loc1) #abundance -> 19 cases and the abundance variable ('abund') says 18…     sp1.loc1$abund<-rep(19,19)     summary(sp1.loc1) #missing values present; abundance for sp1 at loc1 corrected     attach(sp1.loc1)          #Dealing with NAs:     install.packages("mice", dependencies = T) #ok (R at: home & work)     library(mice)     imp <- mice(sp1.loc1)     sp1.loc1 <- complete(imp)     summary(sp1.loc1) #jaust checking... No more Nas!     attach(sp1.loc1)               #Detecting univariate outliers     z.crit <- qnorm(0.9999)          subset(sp1.loc1, select = id, subset = abs(scale(v1)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v2)) > z.crit)     morfo[47,6]     sort(v2[taxoc=="loc1"]) #the nearest observation close to 32.00 is 25.10     sp1.loc1[,6][sp1.loc1[,6]==32.00]<-25.10     subset(sp1.loc1, select = id, subset = abs(scale(v2)) > z.crit) #Rechecking for outliers (now, it's ok)          subset(sp1.loc1, select = id, subset = abs(scale(v3)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v4)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v5)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v6)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v7)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v8)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v9)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v10)) > z.crit)          subset(sp1.loc1, select = id, subset = abs(scale(v11)) > z.crit)          #Standardizing variables     v1.std<-with(sp1.loc1,(scale(v1)))     v1.pad<-v1.std[,1]          v2.std<-with(sp1.loc1,(scale(v2)))     v2.pad<-v2.std[,1]          v3.std<-with(sp1.loc1,(scale(v3)))     v3.pad<-v3.std[,1]          v4.std<-with(sp1.loc1,(scale(v4)))     v4.pad<-v4.std[,1]          v5.std<-with(sp1.loc1,(scale(v5)))     v5.pad<-v5.std[,1]          v6.std<-with(sp1.loc1,(scale(v6)))     v6.pad<-v6.std[,1]          v7.std<-with(sp1.loc1,(scale(v7)))     v7.pad<-v7.std[,1]          v8.std<-with(sp1.loc1,(scale(v8)))     v8.pad<-v8.std[,1]          v9.std<-with(sp1.loc1,(scale(v9)))     v9.pad<-v9.std[,1]          v10.std<-with(sp1.loc1,(scale(v10)))     v10.pad<-v10.std[,1]          v11.std<-with(sp1.loc1,(scale(v11)))     v11.pad<-v1.std[,1]               #Joining the new standardized variables to the sp1.loc1 data set          sp1.loc1<-data.frame(sp1.loc1,v1.pad,v2.pad,v3.pad,v4.pad,v5.pad,v6.pad,v7.pad,v8.pad,v9.pad,v10.pad,v11.pad)          attach(sp1.loc1)          write.table(sp1.loc1,"sp1.at.loc1.txt",quote=F,row.names=F, col.names=T,sep="\t")          detach(sp1.loc1)          #Subsetting the data (v1-v11 are continuous "predictors"): species 2 at locality 1...]


title: Which are differences between the hypergeometric distribution and chi-square distribution

text: As the title suggest...I have a very basic question. I have a case with the following data:               Universe: 18840 balls total     red balls in the universe: 6680     Sample: 382 balls total     red balls in the sample: 160      I would like to estimate if the percentage of red balls in my sample is significantly different from the percentage of reds in universe. In your opinion is it more correct to utilize a chi-square test or an hypergeometric distribution?


title: What to do when I have expected count <5 warning for a chi squared test?

text: I applied a survey consisting of 12 questions to 120 people and each questions include 4 nominal categories; I want to make comparison of people's answers according to their socio-demographic characteristics such as their educational level or socioeconomic status. All of my comparison criteria consist of nominal categories and number of categories change between 3 to 6. I cannot drop or combine comparisons; categories in other words number of categories are fixed. My question is when I compare the frequency distribution of people's answers based on their education levels (for example) through chi-squared test, I got a warning that says (for instance) 7 cells have less than count five; minimum expected count is; 43. I got this warning for almost my all questions - demographic characteristics comparisons. Shall I underestimate this warning and use my test results or is there a different test I should use? If I should use a different test, which one?


title: Examine SVM result by plotting histogram of decision values of training samples

text: I'm working for object detection(computer vision) and have some problems in SVM training. My training configuration is as below.   * Balanced training set (positive 3998/ negative 3998)   * The dimension of my feature set is 2592   * Training using LIBSVM with default options, except kernel function    * Kernel function is intersection kernel, or min-kernel; $K(x,y)=\sum_i\min(x_i,y_i)$. For more information, see ikSVM After training, I plotted density of the scores(decision values) of training samples. I expected the density function is unbiased, i.e. the center of the density function is near 0. But as you see figure, the distribution is biased to right. If the distribution is shitfed to left, then the training accuracy would be 100%. I can't understand the result and want to get a proper result. Is there anyone who can comment on that?


title: Comparing AIC of a model and its log-transformed version

text: _The essence of my question is this:_ > Let $Y \in \mathbb{R}^n$ be a multivariate normal random variable with mean > $\mu$ and covariance matrix $\Sigma$. Let $Z := \log(Y)$, i.e. $Z_i = > \log(Y_i), i \in \{1,\ldots,n\}$. How do I compare the AIC of a model fit to > observed realizations of $Y$ versus a model fit to observed realizations of > $Z$? * * * * * * _My initial and slightly longer question:_ Let $Y \sim \mathcal{N}(\mu,\Sigma)$ be a multivariate normal random variable. If I want to compare a model fit to $Y$ versus a model fit to $\log(Y)$, I could look at their log-likelihoods. However, since these models aren't nested, I can't compare the log-likelihoods (and stuff like AIC, etc.) directly, but I have to transform them. I know that if $X_1,\ldots,X_n$ are random variables with joint pdf $g(x_1,\ldots,x_n)$ and if $Y_i = t_i(X_1,\ldots,X_n)$ for one-to-one transformations $t_i$ and $i \in \{1,\ldots,n\}$, then the pdf of $Y_1,\ldots,Y_n$ is given by $$f(y_1,\ldots,y_n)=g(t_1^{-1}(y),\ldots,t_n^{-1}(y))\det(J)$$ where $J$ is the Jacobian associated with the transformation. Do I simply have to use the transformation rule to compare $$l(Y) = \log(\prod_{i=1}^{n}\phi(y_i;\mu,\Sigma))$$ to $$l(\log(Y))=\log(\prod_{i=1}^{n}\phi(\log(y_i);\mu,\Sigma))$$ or is there something else I can do? * * * [edit] Forgot to place logarithms in the last two expressions.


title: Finding related words

text: I have several files, each of which contains unique terms which are related to each other(without sentence structure). So for finding the word relationships I created a dictionary of bi-grams for every two words occurring in the document, and for every match found I increment the count for the bi-gram occurrence making it more important. For example - if I have football, goal, messi. I have 6 bi-grams football-goal, goal-football, messi-football, football-messi, goal-messi, messi-goal. Doing so in a large set of document I might get several occurrences of football-goal which might tell me that these words are related. And I can further develop the co-relation with 'goal' related words. I have two questions here:   1. Are there better algorithms to achieve this?   2. If I have to use this is algorithm. What are the possible way of visualizing the co-related terms?


title: Log linear models advantages

text: What are the advantages of log linear representation in opposite of table representation? Is it simply computational issue ( avoid overflowing)? For example, in a markov network A-B we can represent the factor P(A,B) as a table:               A B P(A,B)          0 0 10          0 1 1          1 0 1          1 1 10      Alternatively, if we represent the factor P(A,B) as log linear model: $P(A,B) = exp(\sum\limits_{i=1}^4\theta_i f_i(D_i))$ Here f are indicator function. Then basically each theta is the log of entry in the table representation. What would be the advantages of log linear here?


title: Method for calculating percentiles for a set of time series

text: I have a set of time series for a given quantity (e.g. CPU). The measurements are roughly evenly spaced, but the data points aren't synchronised between sets and some sets have missing measurements. The measurements often exhibit a strong pattern over the course of an hour or day. I'd like to compare the most recent time series to a set of previous ones for similar periods (e.g. the same hour from previous days), to see if the recent behaviour is similar to past behaviour. I want this for a data visualisation, rather than any measure of correlation. At the moment I'm taking a time slice from each previous set, combining them and calculating percentiles on each slice. The lack of synchronisation of measurements means that each time slice can have a varying number of measurements. Is there a better method to achieve this? Here's some sample data to illustrate the problem, though I need to be able to do this on other sets. The sample includes 9 sets of hourly CPU data for each of two machines. The `mapped_date` column aligns all series to the same hour. https://dl.dropboxusercontent.com/u/133024/crossvalidated/sample_cpu_data.csv


title: Plot the probability mass function

text: I am trying to plot the probability mass function of a sample of a discrete metric. If it was continuous, I know that using pandas it would be as simple as calling:               sample.plot(kind="density")      But I'm afraid that this is not enough (or not right) for my sample. Is there a function within matplotlib, scipy, numpy, etc. that I could use for plotting it?


title: Logistic regression, "sum" confidence interval

text: I have an existing fitted logit regression model. Model: $\hat{p}(x)=\frac{1}{1+e^{-\hat{\beta}x}}$ With parameter estimates $\hat{\beta}$, and observation $x$ Given a new set of datapoints $x_1,\ldots,x_n$ and $o_1,\ldots,o_n$, I want to construct a confidence interval for $E=\sum_i \hat{p}(x_i)$. To goal is to check if $O=\sum_i o_i$ is (not) significantly different. For a single observation, you can construct a confidence interval for the probability quite easily. Given $\hat{\beta}$ is normally distributed, and the confidence interval for these estimates is ($\hat{\beta}_{0.025}, \hat{\beta}_{0.975}$), the confidence interval for the predictor is: ($\frac{1}{1 + e^{-\hat{\beta}_{0.025}x}}$, $\frac{1}{1 + e^{-\hat{\beta}_{0.975}x}}$). Is there a good way to get the confidence interval for the "sum" of the predicted probabilities?


title: Cross validation in SPSS

text: I need to conduct cross validation of my data to check for predictive validity. This is what I did so far. I divided the sample into two sub- samples(50-50). In SPSS, I then used the split variable to instruct SPSS to keep the data divided into two-sub samples while running regression. Post this I got two regression outputs. And now I'm stuck here. The supervisor wants me to test for predictive validity and using each sub-sample I need to test models and compute r for predictions versus actual values in the hold out sample. How do I do it in SPSS? Thank you very much for your time.


title: Analyze a football match: similar players with DBSCAN and similar trajectories with TRACLUS

text: I'm trying to analyze a dataset that originates from sensors located near players' shoes in a match (http://www.orgs.ttu.edu/debs2013/index.php?goto=cfchallengedetails). I decided to look at clustering to identify:   1. Similar trajectories of players in the match by using the TRACLUS clustering algorithm   2. Similar players by counting some characteristics such as unsuccessful passages, unsuccessful crosses, shots and tackles. I thought to use DBSCAN to cluster them.   3. Group players that pass the ball to each other more often. How can I cluster them? Can I exploit something else from this type of dataset? Is there any other characteristic which I can use in point 2?


title: Regression says $H_1$ is accepted, One-Way ANOVA says $H_1$ is rejected. Are different results possible?

text: I am testing my hypotheses with a 2x2x2 between subjects design. Thus I am testing the moderation effects on the dependent variable. When I used regression to conduct this test, $H_1$ was accepted. However when I used One- way- ANOVA to further test this, $H_1$ was rejected. I am just wondering whether I have done something wrong, or is this normal? Update: Just to clarify my research has three moderating variables : Price, duration and Frequency and I am testing whether or not this has an moderating effect on consumer's enjoyment and attitude towards in app purchases. All questions are 7-point Likert scales.


title: Unbalanced Panel: pooled OLS vs FE vs RE - which method yield unbiased and robust estimators?

text: I am conducting an empirical study (relation between earnings and returns). I have an unbalanced panel with $N=449$ firms and $T=36$ time periods. Regarding the estimation methods I am confused:   1. In most of the papers they account for firm clustering and time issues and therefore cluster for firms and introduce time-dummy variable using OLS. => Is that the simple pooled OLS?   1. If there are entity and time issues I thought that one should use the within (Fixed effects) or at least a weighted average of within and between (Random effects) estimator.   2. I have run some diagnostic test: a) the Hausman test favors the fixed effects b) the Breusch-Pagan Lagrange Multiplier (LM) favors the OLS over the Random Effects So the first question that appears is: How can I check pooled OLS vs Fixed Effects? Additionally, my panel is unbalanced: firms in my panel might get bankrupt or merge with other companies. therefore, the question arises: In case of correlation of these pattern with the idiosyncratic errors, which is best to avoid biased estimtators? The next issue regards Stata.I do not fully understand the difference of the following regressions: > 1) reg y x i.time, cluster(id) > > 2) xtreg y x i.time, fe cluster(id) dfadj > > 3) xtreg y x i.time, re cluster(id) Is the only difference the used estimator? In most of the studies version 1) is used although I do not really get the advantage of version 1) over 2) or 3). In which case is pooled OLS preferred over the other methods? I´ve noticed that the SE are largest when fixing for firm and time using the fe-estimator. And the last question arises: If there are unobserved effects - firm and time: how can I determine if these effects are permanent or temporary (die away over time)? Important for me is to get an inference robust result.


title: How to obtain Item Information Curves for the Partial Credit Model?

text: I have estimated a Partial Credit Model in Stata following Zheng, X. and Rabe- Hesketh, S. (2007) "Estimating parameters of dichotomous and ordinal item response models using gllamm". I'd like to compute Item Information Curves but I can't figure out how to do it from what I found in Google or in other StackExchange posts. a) In Ayala's "the theory and practice of IRT" book, pg 200, it is explained that the item information function in the PCM model is obtained as the weighted sum of category information functions. But in the formula for item information function I don't see where do I sum across item category options. b) Is the item information formula for the Partial Credit model valid for the 2-paramenter version (Generalized Partial Credit Model)? If not, what is the formula for this? Any guidance would be much appreciated.


title: Estimating the time for completing a sequence of actions

text: In short: suppose I have observations for times taken to do some action. I want to estimate, how long will it take to complete a sequence of actions. The estimate should minimize the mean absolute error (L1 norm). === Details: I have observations for positive, non-identically distributed RV $x_i$, it can be assumed that RV are independent. I am interested in $y = \sum_{i=1}^n x_i$. For example, $x_i$ can be the time for action $i$, and $y$ the time for completing a sequence of actions. I would like to find an estimate $\hat{y}$ that would minimize the mean absolute error (L1 norm). So I need to estimate $\mathit{median}(y)$. For $n = 1$ the estimate is straightforward $\hat{y} = \mathit{median}(x_1)$. For large $n$ I have observed that $\hat{y} = \sum_{i=1}^n \mathit{mean}(x_i)$ gives a reasonably good estimate. I suppose this is explained by the Central Limit Theorem: for large $n$, $y$ approaches the normal distribution, thus $\mathit{median}(y) \approx \mathit{mean}(y)$, and $\mathit{mean}(\sum_{i=1}^n x_i) = \sum_{i=1}^n mean(x_i)$. **Question: Is there any analytical solution for approximating median of the sum of independent non identically distributed RV for small $n$, i.e. for approximating $\mathit{median}(y)$ for small $n$ in this context?** Practically, in the data $n$ is distributed: $n=1$ ($10\%$), $n=2$ ($20\%$), $n=3\ldots 5$ ($20\%$), $n=6\ldots 10$ ($10\%$), and $n>10$ ($40\%$), maximum $n$ is $70$. === So far my solution is as follows. I assume that that $\hat{y} = \sum_{i=1}^n \big((1-w)*\mathit{median}(x_i) + w*\mathit{mean}(x_i)\big)$, where $w \in [0,1]$ is a weight. When $n=1$, $w=0$. When $n\rightarrow \infty$, $w \rightarrow 1$. Then I model $w$ as a function of $n$. I estimate an empirical distribution of all observations (ignoring that the variables are non identically distributed), generate synthetic data from that distribution. On that data I learn $w = f(n)$ using some machine learning approach, e.g. ANN. The data is generated as follows. I sample $n$ uniformly at random. I compute the true $\mathit{median}(y)$. I let $w$ run from $0$ to $1$ and pick $w$ that gives the minimum absolute deviation of $\hat{y} = \sum_{i=1}^n \big((1-w)*\mathit{median}(x_i) + w*\mathit{mean}(x_i)\big)$ from $\mathit{median}(y)$. This gives me one data point. I generate, say 10000 such data points and learn a function on this data. This gives a slight improvement over simply using $w=0$ or $w=1$ for all $n$, but not much, and not always. === And finally, to illustrate my intuition that there may exist a nice functional form of $w$, here is an example of $w = f(n)$ learned on synthetic data sampled from log-normal distribution $ln {\cal N}(0,s)$. Solid lines denote estimated functions, circles denote out-of-sample test data. ![enter image description here](http://i.stack.imgur.com/lCScR.png)


title: Predictive models with data that have elements inside them?

text: I'm getting more experience in building predictive models like trees and random forests, but most of my experience is using data that is basically single observations (rows) with many variables (columns). How do I deal with making predictions about things that could potentially have a variable number of objects inside of it. Here is an example: Say I want to predict if an insurance claim is fraudulent. That claim could contain several items inside of it. One claim could be one for theft of an Xbox and a television. Another claim could be for theft of 12 individual pieces of jewelry. The claims have claim level data like total amounts of losses to be paid, but there are also several individual level data elements inside with their own values and attributes. What is the best way to structure the data or model to make a prediction about the entire claim? Do I need a claim level variable called HasXbox? or HasJewelry? Do I need a variable like that for every possible item that could be in the claim?


title: Resources for Computational Algorithms

text: The journal _Applied Statistics_ used to have a section on algorithms, which are available online. Griffiths and Hill also published a selection of these algorithms in a book in 1985. But I can't find any modern equivalent of these resources anywhere. By modern, I mean something that has the updated versions of the algorithms, includes more recent algorithms, and hopefully, the implementation is in a more modern programming language. Does it exist?


title: What interactions to include in my GLM model?

text: I realize this might be a too general question, so I'll describe what I'm doing right now first. I'm working for a virtual insurance company and I have this dataset. It has severity (meaning payment/number of claims), type of car, gender, marital status, age of car, risk type (high and low), company name (ours and others), deductible range. All the variables are categorical. The question I'm trying to answer is: **Is our company paying out more on claims than other companies do?** First I suppose I should choose main effects. The way I do it is that, if I can come up with a story or theory, then I include it. For example, I choose to include marital status because probably married people drive more cautiously since they have families. I include company because...well, this is what I need to study. Now here is my question:   1. After reading many posts, it occurred to me that merely looking at the test statistics is not a good idea. Sometimes even if the p-value is greater than 0.05, we should still keep the variable. My question is, is there a rule-of-thumb or something? For example, if the p-value o test statistics of _marital_status_ is 0.8, then I guess it's obvious that we should not include it. But what about 0.08? 0.07?   2. After selecting main effects, I might should start considering interactions between these variables. My question is, how do I even start? My variable _type_of_car_ has around 20 different values (20 different types of car). I thought about running a regression with all possible interactions, and then drop those with insignificant p-values. That didn't really work out because everytime I drop one, the p-value of the variables in the new model would alter drastically. I apologize for such a long post. I don't have a solid background in statistics (I'm a math major) and I really have to teach myself all the stuff. Any suggestion on how to proceed would be appreciated!


title: Which Models to prefer

text: (Model1) . reg lnQ lnL lnK COOP               ------------------------------------------------------------------------------              lnQ |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------              lnL |   .7682515   .0662688    11.59   0.000     .6364455    .9000576              lnK |   .0537095   .0495861     1.08   0.282    -.0449153    .1523342             COOP |   -.158698   .0660776    -2.40   0.019    -.2901238   -.0272723            _cons |   4.946906   .2614977    18.92   0.000     4.426797    5.467014     ------------------------------------------------------------------------------      (model 2) . reg lnQ lnL lnK COOP COOP_LNL               ------------------------------------------------------------------------------              lnQ |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------              lnL |   .7339817   .0719376    10.20   0.000      .590875    .8770885              lnK |   .0536394   .0494509     1.08   0.281    -.0447341     .152013             COOP |   -.690055   .4454689    -1.55   0.125    -1.576235    .1961245         COOP_LNL |   .0888071   .0736333     1.21   0.231    -.0576729    .2352872            _cons |   5.147985   .3095237    16.63   0.000     4.532244    5.763726              firm            float   %9.0g                 firm code (or firm ID)         va              float   %9.0g                 value added in millions of lire         labor           float   %9.0g                 permanent employees at year end         capital         float   %9.0g                 fixed assets at historical cost in millions of lire     COOP = a dummy variable for a worker cooperative. Conventional firms have code (or ID) numbers (given by FIRM) less than 100, while the code (or ID) numbers for worker cooperatives exceed 99.     COOP_LNL=COOP*lnL         lnQ = the natural logarithm of output (value added)     lnL = the natural logarithm of labor     lnK = the natural logarithm of capital      How to determine which models to prefer? My thinking is do a F-test on COOP and COOP_LNL to see if this p-value is bigger than the p-value of COOP for model 1. If this is the case, I will pick model 1 because it has less p-value...Is it the right way to approach it? Is there any better way or different way?


title: Theil-Sen estimation

text: Is the theil-sen estimation in robust regression only limited to a two dimensional problem or can you use it for more than one indepedent variable as well?


title: Tests of within subjects contrasts in R

text: I have run a split splot model in SPSS (via repeated measures function) and I would like to reproduce my results using R. To do so I used ezANOVA function from ez package to obtain sphericity tests and correction and type III SS. I have read that aov function does not give type III SS. So far the results are identical with those from SPSS but I cannot tell so when I run tests of within subjects contrasts. I use polynomial contrasts and my factor has six levels. What I manage to get is t-values instead of f-values or f-values but not the same as in SPSS. Here is my data and code `id subj treat m1 m2 m3 m4 m5 m6 1 1 1 455 460 510 504 436 466 2 2 1 467 565 610 596 542 587 3 3 1 445 530 580 597 582 619 4 4 1 485 542 594 583 611 612 5 5 1 480 500 550 528 562 576 6 6 2 514 560 565 524 552 597 7 7 2 440 480 536 484 567 569 8 8 2 495 570 569 585 576 677 9 9 2 520 590 610 637 671 702 10 10 2 503 555 591 605 649 675 11 11 3 496 560 622 622 632 670 12 12 3 498 540 589 557 568 609 13 13 3 478 510 568 555 576 605 14 14 3 545 565 580 601 633 649 15 15 3 472 498 540 524 532 583` m1-m6 represent measurements in 6 different time points. My code for the contrasts. `long.df<\- melt(data, id=c('subj','treat')) long.df<\- long.df[order(long.df$subj)] names(long.df)<\- c('subj','treat','time','meas') mod<\- lm(meas~time + time^2 + time^3 + time^4 + time^5, long.df) summary(mod) # this is how I obtain t-values for the polynomial contrasts # mod1<\- aov(meas~time + time^2 + time^3 + time^4 + time^5, long.df) summary(mod1, split=list(time=list("Linear"=1, "Quad"=2,'q'=3,'f'=4,'fif'=5))) # this is how i obtain f-values #` I think that the F-values correspond to type 1 SS (as I use aov function). How am I wrong? Is there a way to conduct the Tests of within subjects contrasts and obtain type III SS so I will have identical results to SPSS? Here are the SPSS results. Many thanks in advance! Source time Type III SS/ df/ Mean Square/ F/ Sig.   time Linear 123662.881/ 1/ 123662.881/ 83.591/ .000   Quadratic 5928.007/ 1/ 5928.007/ 18.968/ .001   Cubic 10462.676/ 1/ 10462.676/ 28.075/ .000   Order 4 798.193 / 1/ 798.193/ 4.010/ .068   Order 5 1702.743/ 1/ 1702.743/ 4.878/ .047


title: Interpreting simple linear regression

text: Good day, I will am currently doing some self-study on Simple Linear Regression. I understand the formula: $$\hat{Y}_i=\beta_0+\beta_1 X_i$$ But what does it mean when the value of $\beta_1$ is 0? Does it mean that $X$ and $Y$ have no relationship? Appreciate some advice and explanations please


title: Why do ANOVA contrast tests use N - k degrees of freedom?

text: I would have thought that the degrees of freedom would be the same as a regular t-test, i.e. N - 1, since in a contrast we are either comparing two groups, or two sets of groups. Why do we instead use the within SS degrees of freedom (N - k)? Here are some example SPSS tables from http://wiki.uva.nl/spss/index.php/One- way_ANOVA ![enter image description here](http://i.stack.imgur.com/RORFz.png) ![enter image description here](http://i.stack.imgur.com/kC0sw.png)


title: How do I correct for gender bias for Bayesian credible intervals?

text: Forgive me if this question is answered elsewhere - it may be all I need to be told is the correct terminology. I have a survey in which I ask a number of Yes/No questions (actually Yes/No/Don't know - but for simplicity I want to assume Yes/No for now). Say for one question I have 583 say yes and 237 say no. Suppose I had really no idea what the proportion of yes and no would be in the population as a whole (this was reasonably accurate), so that my prior belief can be represented by a uniform distribution over the proportion who would say yes in the population, then, if I have understood Bayesian inference properly, I obtain a Beta posterior distribution. I can compute credible intervals from that (using `scipy.stats.beta.interval(0.95, 583, 237)` for instance). (For Yes/No/Don't know I can use a Dirichlet distribution). But in fact my survey includes men and women. Consistently women are less likely to say yes than men. In my sample the yes v no numbers are: 172 v 118 (for women) and 411 v 119 (for men). As you can see more men answered the survey than women. I'd like to correct for that but I don't know how. Eg, I could fractionally weight the men's answers by the ratio of men/women. The problem is that I don't know how to compute credible intervals for the weighted numbers that would result. When I try to approach this from first principles, I get stuck working out what a good prior is. If `P(M)` is my prior belief that the proportion of men who would say yes in the population is M (and `P(F)` has the corresponding meaning) with `m` and `f` representing the sample values then I am computing `P(M, F | m, f)` via Bayes' theorem and I need a prior for `P(M, F)`. Although I have no real knowledge about Yes/No proportions for men or women individually I am quite sure they are correlated in some way. Are there well know priors that can incorporate that weak prior knowledge? Or can I avoid working out the whole problem and just do some kind of weighting? Telling me where to look would be very useful.


title: Why adding an NA indicator column instead of value imputation (for randomForest)

text: I have downloaded the sample codes from kaggle for the randomForest benchmark _[URL?]_ and there's this part that I don't understand.                 appendNAs <- function(dataset, cols) {       append_these = data.frame( is.na(dataset[, cols] ))       names(append_these) = paste(names(append_these), "NA", sep = "_")       dataset = cbind(dataset, append_these)       dataset[is.na(dataset)] = -1       return(dataset)       }      I understand that randomForest does not accept NAs and I usually just impute them with median or rfImpute. Adding extra columns as NAs indicator is new to me. I have looked at the help file for randomForest but didn't anything useful. Would be great if someone can explain how and why the extra columns work.


title: Spatial patterns analysis for images

text: i have done some research about spatial patterns analysis, using some spatial statistics mesures, all the cours that i found talks about spatial statistic in the context of GIS, in the context of image analysis, i am trying to use different technics such as Quadrat Analysis, Nearest Neighbors , Joins Count, Moran's I, in order to classify a pattern image, as randomly , clustred or uniform distribution, i have done some implementations in matlab, but i didn't anderstand how to interpret the results. especialy for test segnificance for exemple chi-square, Z statistic ..., in images i dont see here applying test segnificance, because ther is no Statistical inference !! Moran's I, try to mesure autocoorelation of one variable, but in my case i have points (pixels) with x & y coordinates, i am wondering how this could be applicable for images. the same issues for getting ranks from pixels coordinates to use for Spearman's rank correlation and Kendall's rank correlation. Thanks a lot.


title: ECM and Cointegration

text: Can i apply a standar OLS ECM after proving cointegration using johansen approch? The results for VECM are uncertain and one of the equation gives a negative R squared. If i want to test only one dependent variable in a simple bivariate model.