title: Is it possible to combine the F-test and t-test to test for differences in both the variance and the mean?

text: I would like to test $H_{0} :$ $\mu_{1} = \mu_{2}, \ \sigma_{1}^{2} = \sigma_{2}^{2}$ vs $H_{0} :$ $\mu_{1} \neq \mu_{2}, \ \sigma_{1}^{2} \neq \sigma_{2}^{2}$      Is there a test statistic that would allow me to detect both differences in the mean and differences in the variance? Is there a way that I could possibly combine the $F$ - $test$ and the $t$ - $test$ to derive this test statistic? I am assuming that the population is normally distributed but what if it wasn't?      I am very stuck so any help would be greatly appreciated.


title: Estimating the mean prediction

text: Suppose I have realized n predictions $p_1, \dots, p_n$ using a classification method (for instance a logistic regression or a linear discriminant analysis). Each prediction is associated with a standard error $\sigma_1, \dots, \sigma_n$. I want to estimate the mean prediction and a confidence interval for it at a certain level. I'm wondering what the best approach is to tackle this problem. Should I compute a classic sample mean and a sample standard deviation from the predictions? Or would it be more accurate to compute a kind of pooled estimate using the individual standard errors? Another possibility to get the bounds of the confidence interval would be to compute the empirical quantiles but I'm not sure it would be statistically right. What do you think? * * * That's right, my n predictions refer to the same target which is a binary variable. But I may have forgotten to precise that it was the same model (in this case a fitted logistic regression) that produced all the predictions, that's why I wasn't sure whether a pooling method was appropriate. Thanks for the information. I didn't think of considering these predictions realizations of a mixture model. I rather looked into pooling methods and the fact is I found many different formulas to combine the predictions into a "consensus" so I was a bit confused. I also had difficulty with building the confidence interval for $p_c$ since I didn't know the quantiles of the distribution of the predictions. Of course, I thought of using the quantiles of the standard normal distribution, but in some cases, the predictions reduced to a pointwise estimation (for a specific sub-sample of the population, all individuals were similar and the prediction was the same for all of them), so this approximation should be very bad. Can the normal quantiles be used in this very situation? Another thing bugged me: on some sub-samples, the discretized predictions returned only 0's or only 1's so I changed them back to probabilities in order to get lightly different values and to be able to estimate a mean prediction. But I'm not sure this was the right approach to follow...


title: Risk table for a Kaplan Meier plot in R

text: I need to make a Kaplan Meier plot with a risk table beneath it. Otherwise stated, I need a table of the number of subjects at risk at different time points aligned below the figure. I found a website that explains how to do this for a plot that contains multiple subgroups. This is the `ggkm` function, the code for which is available here.               # Example of a plot like this     library(survival)     data(colon)     fit <- survfit(Surv(time,status)~rx, data=colon)     ggkm(fit, timeby=500, ystratalabs=c("Obs","Lev","Lev+5FU"))      ![Example of a plot like this](http://i.stack.imgur.com/9iaSp.png) Looking at the data that was used to make the above plot, we see:               Call: survfit(formula = Surv(time, status) ~ rx, data = colon)                     records n.max n.start events median 0.95LCL 0.95UCL     rx=Obs         630   630     630    345   1723    1323    2213     rx=Lev         620   620     620    333   1709    1219    2593     rx=Lev+5FU     608   608     608    242     NA      NA      NA      **My question:** My plot only has one group. The call looks like this:               Call: survfit(formula = Surv(Recur_day/365.242, Recur) ~ 1, data = study_data)          records   n.max n.start  events  median 0.95LCL 0.95UCL          440     440     440      92      NA      NA      NA       ![My plot](http://i.stack.imgur.com/QqMz4.png) When I try to use the `ggkm` function I get an error like so:                ggkm(survfit(formula = Surv(Recur_day/365.242, Recur) ~ 1, data = study_data), timeby = 2)     Error in data.frame(time = sfit$time[subs2], n.risk = sfit$n.risk[subs2],  :        arguments imply differing number of rows: 1, 0     In addition: Warning message:     In max(nchar(ystratalabs)) :       no non-missing arguments to max; returning -Inf      Does anyone have an idea of what I'm doing wrong? Or is there an easier way for me to do this by hand? I am open to just putting the data in a text box in PowerPoint myself, I'm just not certain of how to get the number of subjects at risk at each timepoint (I'd like to do at two years, four years, six years, etc.)


title: PCA components changing over time

text: If I do the PCA on the whole dataset I get 7 components that can explain 90% of the variance, if I split the dataset into 2 (sorted by time), the number of significant components in the first half goes to 5 (with 15 variables present in one or more components) and in the second half goes to 8 (with 21 variables present in one or more components), can we infer that some of these variables become more significant in the latter half compared to first half?


title: Impact of regression normality-assumption on model comparison & prediction?

text: This question is a continuation of the discussion here: How to test the statistical significance for categorical variable in linear regression? Following Macro's suggestion, I started a new thread. The new question is not limited to only study the inclusion/exclusion of categorical variable. It's about general model comparison and prediction. I found that my data is highly non-normal. The QQ plot is as follows: the curve is all below the straight 45 degree line. The curve is tangent to that straight line. And the curve looks like the curve of f(x)=-x^2 ( shape-wise). It's not the entire set of points that's under the 45 degree line. It's the curve with the shape of f(x)=-x^2 is "tangent" to the 45 degree line. By "tangent" I should have meant that those points around the "tangent" point are actually above the 45 degree line, very slightly though. Therefore, visually speaking, most of the data (~98%) are below the 45 degree line... These are the residual QQ plot coming out from "plot(lmModel)"... Using `qqnorm(lmModel$res); qqline(lmModel$res)` I got exactly the same curves and lines. My questions are:   1. If my end-goal is to use yhat to do prediction onto a wide data-set, does the non-normality of data matter?   2. For Model Comparison, the approaches pointed out by Macro probably won't work; any other alternative approaches that don't assume Gaussian distribution?   3. What shall I do to fix the non-normality problem? (data-size: 10 variables, 1700 observations).  Thank you!


title: Dealing with a large number of predictors in Logistic Regression

text: Let's say I have a logistic regression model which predicts whether a consumer will buy an item based on about 10 consumer characteristics. $$\begin{array}{rcl}Buy &=& B_0 + B_1\times Gender + B_2\times CreditType + B_3\times Education + B_4\times OwnsHome \\\\\phantom{Buy} && \+ B_5\times CarMake + B_6\times CarYear + B_7\times State + B_8\times Income + B_9\times Insurance \\\ \phantom{Buy} &&\+ B_{10}\times CarAccidents\end{array} $$   1. Is there ever an issue with including too many predictors in a logistic regression model? I'm not talking about insignificant variables or ones that may be related, but just the sheer number of variables included in a model.    2. With a larger number of predictors, how should one present the regression results in a meaningful manner? Is it just a matter of plotting the probability curve for $Y=1$, or are there "better" ways of doing this. I'd be doing this in R, so any help on that end would be appreciated.


title: Gradient descent and elastic-net logistic regression

text: I'm currently in the process of trying to understand the paper Regularization Paths for Generalized Linear Models via Coordinate Descent by Friedman et al. with regard to the regularization of logistic regression. Unfortunately, I was not able to figure out the exact algorithm that is used for optimization. I was looking for other papers that use coordinate descent in this context and came across Sparse logistic regression for text categorization by Genkin et al. They go into more detail and in section 3.1 they mention > Let $f(\beta_j)$ be the objective expressed as the function of only one > component $\beta_j$ with all the rest being fixed, and $Q(\beta_j, > \Delta_j)$ be an upper bound on the second derivative of $f (\cdot)$ in the > $\Delta_j$-vicinity of $\beta_j$ , $\Delta_j > 0$:   >  $Q(\beta_j , \Delta_j ) \geq f (\beta_j + \delta)$ for all $\delta \in > [−\Delta_j , \Delta_j ]$. Later on they use the following formula as upper bound for L2-penalized logistic regression: > $Q(\beta_j , \Delta_j ) = \sum_i x_{ij}^2 F(\beta^T x_i y_i, \Delta_j > x_{ij}) + 2\lambda$   >  $F(r, \delta) = 0.25$ if $|r| \leq |\delta|$   >  $F(r, \delta) = 1 / (2 + \exp(|r| - |\delta|) + \exp(|\delta| - |r|))$ > otherwise. I was surprised to see that this apparently necessary upper bound was not mentioned in the paper by Friedman et al. In addition, the bound differs between L1 and L2 penalty. What would the bound look in the generic framework proposed by Friedman et al. which is applicable for a wide range of penalties? My second question is related to the upper bound proposed by Genkin et al. They did not provide or reference a prove that validates there choice. Am I missing something obvious here?


title: Scatterplot with contour/heat overlay

text: I saw this plot in the supplement of a recent paper and I'd love to be able to reproduce it using R. It's a scatterplot, but to fix the overplotting there are contour lines that are "heat" colored blue to red corresponding to the overplotting density. How would I do this? ![enter image description here](http://i.stack.imgur.com/fAzpe.png)


title: Proof for principal component score

text: Can we proof that two principal component scores are not equal or one is better than the other in terms of variation? For example, the 1st principal component score is defined as $z_1 = a_{11}(y_1 - \mu_1) +...+ a_{1p}(y_p - \mu_p)$ with the eigenvalue $\lambda_1$. We all know that the proportion of total population variance due to $j^{th}$principal component is determined by its eigenvalue, such as $w_j = \frac{\lambda_j}{\lambda_1 + ... +\lambda_p}$. If we calculate a new score such as $s = z_1\times w_1 + ... + z_p\times w_p$. In addition, we also know that, Var$(z_1) = \lambda_1$, Cov$(z_i, z_j) = 0$, total population variance = $\lambda_1 + ... + \lambda_p$, and $\Sigma{w_p} =1$. My question is can we proof that Var$(s) = \lambda_1\times {w_1}^2 + ... + \lambda_p\times {w_p}^2 \geq$ or $\leq \lambda_1$? Many thanks in advance.


title: Two slopes between 0 and 1, neither different from 0, could they be different from each other?

text: I'm reading a paper that does not report the coefficients from two OLS regressions. In both cases there is 1 response variable and 1 predictor variable. The predictor variable is the same in both cases. I know the subject matter of the paper well, which leads me to believe that the means of the slopes are almost certainly between 0 and 1. Although the slopes for these two regressions are not reported, the author does report that neither slope is significantly different from 0 (p ≥ 0.05). If neither slope is different from 0, but both slopes are between 0 and 1, could the slopes be different from each other? To try to figure this out, I did a quick test in R. I used two slopes that were very different (0.99 and 0.01), but chose s.e.'s for each that would make them barely "insignificant". To compare the slopes, I used the formula from the answer to THIS question.               pnorm(        (0.99 - 0.01)/ #difference between means       sqrt(0.61^2 + 0.0062^2), #sqrt of sum squares of s.e.'s       lower.tail=FALSE     )      OK, so this quick-and-dirty test suggests that the two slopes in the author's analysis can't be different. Is it necessarily true that if two slopes are between 0 and 1, and neither different from 0, that they cannot be significantly different from each other?


title: How does ties.method argument of R's rank function work?

text: I am using `rank(a, ties.method="max")` to rank a. But I am not quite sure what does `ties.method="max"` do. Can you please help?


title: Significance testing of A/B test with poisson distribution

text: I am evaluating results of an A/B-style experiment where we have different groups of people that can perform a target action. They can do the target action as many times as they want, but most people don't do this very often (poisson distribution). I'm interested in how I can best calculate the statistical significance of the impact on the number of times they do the target action, similar how a normal A/B test works. How can I do this? Can I use a normal t-test?


title: Testing for differences in very small proportions

text: This is a somewhat typical a/b test setup. However the proportions that I'd like to test for differences are very small (< 1%). Consequently a 20% difference, which is meaningful, is 0.001. That's a small effect size! A power test of the chi-squared test for differences in proportions (same sample size) using the `pwr` package in `R` yields               difference of proportion power calculation for binomial distribution        (arcsine transformation)                     h = 0.001              n1 = 1e+05              n2 = 1e+05       sig.level = 0.05           power = 0.05574725     alternative = two.sided      So power grows very slowly. How can I get around this? The group `a` proportion is lower than group `b` (which contains the treatment we're interested in). As expected the group `b` proportion is higher than group `a`. Say group `b`s proportion is 0.007. Can I use an exact binomial test on group `a` with the alternative based on being greater than group `b`s proportion? When I do this I get a p-value of 0.97, indicating that it's very unlikely that group `a` will increase in conversion to trump group `b`s current position. I guess the question becomes one of power with this test? Any comments on the validity of this?


title: Statistical test for a gradually self-improving suggestion system

text: I built a suggestion system, which provides a list of candidate items to users. As users choose some items from the suggested list, the system remembers those choices, and the next time the system reorder the candidate list so the items that were heavily selected before will appear at the top of the list. So the more user use and select items, the better the ranking in the future. I asked some users to try out the system and collect the time each user spent in selecting items. From my raw calculation, I know that on average, the first half of the users spent 20% more time than the second half of the users. This may suggest that the self-improvement of the system is effective. But I don't know which statistical test to use in order to check this hypothesis. Some of my problems include:   * Which distribution to use   * Since the system is changing as more and more users try it, how to deal with this?   * I only have 12 users now, it is not that many, how to deal with this?


title: Does anyone know how this visualization was made?

text: I've been tasked with implementing a visualization something like this one from the New York Times: http://www.nytimes.com/interactive/2012/08/05/sports/olympics/the-100-meter- dash-one-race-every-medalist-ever.html?hp Mine doesn't have to be nearly this amazing and sophisticated, but I'd like to capture at least some of the presentation elegance. I can't find any sort of behind-the-scenes discussion of this graphic. Does anyone have any idea how this might have been done?


title: Looking for metric to measure linear and constrained relationship between 2 variables

text: Data are 2D x-y pairs, both of which have [0-1] possible range. I need to calculate how well these data fit a linear relationship, where the intercept is 0 and the slope is 1. In other words, a "perfect" metric score should be achieved if the data form a straight line between (0,0) and (1,1), and always worse than this otherwise. The data points are skewed such that most (~98%) are in the range (0.95,0.95)-(1,1). However, I am much more interested in the shape and slope of the scatterplot, rather than how well the points are correlated overall. Gaps may be present in the data on either axis, and there are 100s to thousands of data points. Variance occurs in y-variable, but does not follow any standard distribution shape. What single-valued metric could I use to measure how well the shape of the scatterplot fits the ideal linear relationship, while remaining unaffected by the huge skew in values?


title: Binomial GLMM: Model validation & ceiling effect

text: My data has a binary response `acc`(correct/incorrect), one continuous predictor `score`, three categorical predictors (`race`, `sex`, `emotion`) and a random factor `subj`. All predictors are within-subject. By selecting the random effects first and then the fixed effects, I ended up with this model: `M<-glmer(acc ~ race + sex + emotion + sex:emotion + race:emotion + score +(1+sex|subj), family=binomial, data=subset)` I need help on interpreting validation plots, figuring out if they show a "ceiling effect" in `acc`, and fix any problems that need to be fixed. * * * To validate the model I get the residuals and fitted values                fitted<-predict(M,type="response")      resid<-resid(M,type="pearson")      And plot the residuals against the categorical predictors                plot(subset$race,resid)          plot(subset$sex,resid)      plot(subset$emotion,resid)      Those three plots show a slight pattern of more negative and dispersed residuals in "easy" conditions. The pattern looks slight to me (i may be wrong). I plot the residuals against the continuous predictor                plot(subset$score,resid)      ![enter image description here](http://i.stack.imgur.com/WZbMp.png) This plot of residuals against the continuous predictor is worrying and shows a clear pattern of more negative and dispersed residuals when `score` increases (the task becomes easier).                plot(fitted,resid)       ![enter image description here](http://i.stack.imgur.com/LSTXm.png) This plot is also worrying showing a clear pattern of more negative and dispersed residuals when the probability of a correct answer increases (either for y=0 or y=1, not sure which one). Apparently these patterns may simply be coming from the log() in the link function. I further tried to plot a regression line as shown in here: link. ![enter image description here](http://i.stack.imgur.com/YX5EU.png) Supposedly it should be straight. Are these patterns strong enough to abandon the model? I would think that they are not, since the plots look very much like the ones from the links, except there is a general tendency to predict more "y=1" i gather. I know there is a ceiling effect in my data, with some easy conditions having almost only correct responses (y=1). This is why I am being maybe overly skeptical about my model. Are these patterns a symptom of this?


title: How can I apply LR with non-linear features?

text: My training data looks like this: ` win gender age country ... 1 F 20 US ... 1 M 21 CA ... 0 M 22 CA ... 0 M 20 CN ... ` Now I want to apply `Logistic Regression` to predict the win rate of a specified person, I understand the feature `age`, but for `country` and `gender`, I don't think they fit for my need (I think only linear features should be used), is there anyway to use them in LR(or other simple algorithms) ?


title: How to compare multiple (17) frequencies (as a "combined trait") of different measured features between two groups

text: I have two groups of samples ($A$, $n=100$ and $B$, $n=30$). For each groups we measured whether a trait is present or absent in each sample (0 or 1), resulting in frequency of "ones" for each trait. This we have done for 17 traits. Each individual trait can be easily compared using Chi square and we found that 3 out of 17 traits are significantly different between $A$ and $B$. Now, is there a way to test all 17 frequencies together and to show that populations $A$ and $B$ are different as a population when considering all 17 traits, although individual samples from $A$ and $B$ will look very similar when considering 1,2,... or 17 traits. Or is it enough to say that since 3 out of 17 traits are significantly different between $A$ and $B$, these populations ARE indeed significantly different?


title: Alternatives to ANCOVA when assumptions are violated

text: I am doing research on organisational change and looking at how certain types of change (e.g. closures of offices, layoffs, investment on activities/functions) have affected seven psychological variables (commitment, trust in leaders, procedural fairness, perception of organisational support, etc.). Now, I want to know how the changes (related to the types of organisational change described above) in the content of the psychological contract (measured through four of the psychological variables) impact organisational commitment (if at all). My original idea was to run an analysis of covariance in SPSS, but in checking for the basic assumptions I found that the regression slopes are not homogeneous, and that is the end of ANCOVA according to several authors. I am left without a clear idea of how to circumvent this problem. What kind of analysis can be performed when:   1. The covariates and the factor (type of change) are related, and   2. The regression slopes are not homogeneous? Please, bear in mind that mine is NOT an experimental design, rather observational, change happens without my intervention and I just look at an internal survey that studies the psychological and demographic variables.


title: Pairwise comparisons

text: I have a situation like this: > >               Sample1   Sample2   Sample3    .....   >     Gene1      0.38        0        0.13     .....   >     Gene2      0.11       0.12       0       .....   >     Gene3      0.55        0        0.9      .....   >     .....     .....      .....     .....     .....   >   I need a statistical test able to compare Gene by Gene and Sample by Sample the numerical values reported. In other words I would like to perform a statistical test able to tell me how statistically significant are the measurements: `0.38` (`Gene1`) and `0.11` (`Gene2`) for `Sample1`, as well as `0.55` (`Gene3`) and `0.38` (`Gene1`) for `Sample1`, and `0.55` (`Gene3`) and `0.11` (`Gene2`) for `Sample1` and so on. 0 means "not_reported" so I suppose I will ignore comparisons with 0. Can anyone help me please? I'm not an expert statistician.


title: Validating questionnaires

text: I am designing a questionnaire for my dissertation. I am in the process of validating the questionnaire I have applied a Cronbach's alpha test to the initial sample group. The responses to the questionnaire are on a Likert scale; can anyone suggest any further tests to apply to help test its validity. I am not an expert on statistics so any help would be appreciated. I have been doing some research and it appears I can do a Rasch analysis has anyone got any free software sites to apply this test and advice?


title: Sources of within-cluster correlation other than "random shocks"

text: We are often interested in estimating the limiting distribution of a parameter in situations where the data exhibit dependence within clusters. For example, a study of the effects of a household-level treatment on household-level outcomes must contend with the possibility that households within villages will have correlated outcomes. Thus, when computing standard errors for the treatment effects, we typically take this clustering into account by using a "cluster robust" covariance estimator or perhaps a "random effects" model. (See http://bit.ly/bAah5L for an example.) The properties of these covariance estimators are typically studied by assuming the source of the clustering are common "shocks" that occur within a group. That is for groups indexed by $g$ and units within groups indexed by $i$, we typically write down a model of the form, $y_{ig} = \alpha_g + \epsilon_{ig}$, where the $\alpha_g$'s denote the group level shocks. We then typically assume that the $\alpha_g$ are independent across groups with group level variance $\sigma^2_{\alpha}$ and the $\epsilon_{ig}$ are independent draws with variance $\sigma^2_{\epsilon}$, in which case the intra-cluster correlation is the familiar $\frac{\sigma^2_{\alpha}}{\sigma^2_{\alpha}+\sigma^2_{\epsilon}}$. Here's what I am wondering: does this linear "group shocks" model implicitly rule out important forms of within-cluster correlation? (I am limiting the question to correlation for now, and not other types of non-correlation dependence.) That is, are there ways that within-cluster correlation can arise that cannot, through algebraic manipulation, be expressed in terms of a group- level shock? Or is the group level shock characterization general with respect to within group correlation? I am inclined to try working out some of the examples, but I thought I'd put this up to see if someone else has already done the heavy (or maybe it's not so heavy) lifting. If this is not a well-posed question (i.e. unanswerable), I'd also welcome comments about why. (Note I am not putting "clustering" on the tag list here because it's not about cluster methods as the term is usually applied.)


title: Significance of $r^2$ value

text: I know about $r^2$ tells you about the amount of variation that can be explained by the predictor variables. I have run a model in which the rsquare has value 0.3010 but has false positive rate of around 15.60%. So this model which is logit in nature predicts 84% cases right. I want to know two things:   1) Is this false positive rate significant, i.e. is this prediction good enough?   2) Does it matter if I $r^2$ is as low as 0.3 if my model is good in prediction as we know that $r^2$ is a value used for explaning and not predicting?


title: Training a model

text: I have heard the terms training and validating a model. I know that we select variables which are most statistically significant and we look for other things like multicollinearity. My question is: what does traning a model involves more than this ?


title: Cross-validation and weights with glmnet

text: Are any weights supplied to glmnet`::cv.glmnet` taken into account for the calculation of the figure of merit (AUC for logistic regression in my case) on the "held-out" dataset, or only used directly in the fitting?


title: What is the difference between loadings and correlation loadings in pca and pls?

text: One common thing to do when doing Principal Component Analysis (pca) is to plot two loadings against each other to investigate the relationships between the variables. In the paper accompanying the pls R package for doing Principal Component Regression and PLS regression there is a different plot, called the 'correlation loadings plot' (see figure 7 and page 15 in the paper). The correlation loadings, as it is explained, is the correlation between the scores (from the pca or pls) and the actual observed data. It seems to me that the loading and correlation loadings are pretty similar, except that they are scaled a bit differently. A reproducible example in R, with the built in data set mtcars is as follows:               data(mtcars)     pca <- prcomp(mtcars, center=TRUE, scale=TRUE)          #loading plot     plot(pca$rotation[,1], pca$rotation[,2],          xlim=c(-1,1), ylim=c(-1,1),          main='Loadings for PC1 vs. PC2')          #correlation loading plot     correlationloadings <- cor(mtcars, pca$x)     plot(correlationloadings[,1], correlationloadings[,2],          xlim=c(-1,1), ylim=c(-1,1),          main='Correlation Loadings for PC1 vs. PC2')      ![loadingplot](http://i.stack.imgur.com/GJA4c.png) ![correlationloadinsplot](http://i.stack.imgur.com/y32sb.png) What is the difference in interpretation of these plots? And which plot (if any) is best to use in practice?


title: Practical issues with dynamic panel data modeling

text: Unfortunately for me, I've got a situation where I need to control for the lag of a dependent variable as a robustness check against an alternative interpretation of my main regression. The baseline specification is $$ y_{it} = \alpha_i + \delta y_{i,t-1} + X'\beta + \epsilon_{it} $$ The lagged dependent variable renders OLS (via the within (de-meaning) transformation) inconsistent, so it seems that the standard approach is the Arellano-Bond estimator, described in Chapter 8 of [this book].1 Basically an instrument matrix is constructed from all available lags of the DV starting with the second lag (i.e. there are more instruments in later than earlier observations), and applied as a 2-step GMM estimator to the first-differenced model of interest. I'm still digesting these methods, and it seems clear that this is an active area of econometrics research. I've got a few practical questions. Many of them might not have answers yet.   1. What do you do when the AR2 test fails? By construction, these models will have AR1 serial correlation. But they shouldn't exibit AR2. What if they do? Can I include more than one lag of the DV? How then should I construct the instrument matrix?    2. When one has a huge dataset, will the AR2 test be too-powerful -- detecting small degrees of correlation with no practical significance?    3. Is it possible to include more than one variable among the instruments? Say I'm interested in an interaction of $y_{t-1,i}$ with some variable $C$. This will be endogenous and require an instrument -- probably of the same structure as the instrument for the lagged DV itself -- but it isn't clear to me how I'd generalize the AB instrument matrix to incorporate this.   4. Finally, what are some readable books and papers on these methods? I get the sense that this is a fresh enough research area that little effort is yet expended on cleaning up the research for the purpose of exposition.


title: What are appropriate validation methods for a Bayesian network model with low sample size?

text: I am currently using a Bayesian network model with 20 variables and 210 data points, with 15 locations measured at 14 different time points each. There are also some restrictions on what types of connections are allowed. I have looked at **leave-one-out cross-validation** methods (and arguably leave-two-out would be computationally feasible with this sample size), and I am currently using that method. However, I want to know which of the following, if any, are appropriate means of going about this:   1. **k-fold cross-validation** , leaving out a single location (14 points). My main concern is that if the model were to be applied to locations not sampled, it is quite possible that correlations within a single locations's set of measurements may be unnecessarily influential on a model. Additionally, it's been suggested that k-fold cross-validation produces better results for simpler models such as linear regression. [1] Since my sample size is small, however, I'm not sure how much I can afford to cut out.   2. **Bagging** (bootstrap aggregating)[3]. I've mostly seen this method applied to decision trees, but it's used to avoid overfitting.    3. **Jacknifing** to determine the properties of the estimators being derived.  Additionally, one interesting property (at least for the datasets used in Zuk et al.2) of the Bayesian Networks is that underfitting decays exponentially fast with sample size, whereas overfitting decays as a power of N. Keeping that in mind, I am trying to decide what methods are most appropriate for my scenario. Sources used: 1 FAQs.org: What are cross-validation and bootstrapping? 2 Zuk et al. "On the Number of Samples Needed to Learn the Correct Structure of a Bayesian Network". 3 Elidan, Gel. "Bagged Learning Structure of Bayesian Networks". Edit: In addition to the above methods used, it has been suggested to perform **confirmatory factor analysis** as an additional method to determine if groups of connections make sense. It doesn't tell you anything about the signs/magnitudes of relationships, though.


title: R GLM verifying performance statistics for validation

text: I do not want to say my AUC is 0.77 and find out I am overlooking a lot. Below is my code and a two question at the bottom:               setwd("g:/docs/Model")     mydata <- read.csv("TModel_JUL16.csv")     summary(mydata)     mydata$ZIP <- factor(mydata$ZIP)     mydata$RACE <-factor(mydata$RACE)     mydata$GENDER <-factor(mydata$GENDER)     mydata$ADM_SOURCE <-factor(mydata$ADM_SOURCE)     mydata$TYPE <-factor(mydata$TYPE)     mydata$SVC <-factor(mydata$SVC)     mydata$WARD <-factor(mydata$WARD)     mydata$DISPO <- factor(mydata$DISPO)     mydata$STATUS <- factor(mydata$STATUS)     myglm <- glm(RA~ZIP+WARD+SVC+AGE+RACE+GENDER+ADM_SOURCE+TYPE+DISPO+STATUS+MEDICARE+LOS+IN+ED+PCI+CABG+DIALYSIS+AMI+CHF+PN+COPD+PVD+CVD+DEM+CPD+RhD+PUD+DIA+HePa+RD+SLVD+AIDS+DRUG+PHY+DEP+ALCH+PC+SYNC+SEP+ANEM+MALNEO+METCNCR+COAG+STK+NUTDEF,family=binomial,data=mydata)     summary(myglm)     prob <- predict(myglm, type="response")     mydata$prob=prob     library(pROC)     g <- roc(RA~prob,data=mydata)     plot(g)      These are my results:               Call:     glm(formula = RA ~ ZIP + WARD + SVC + AGE + RACE + GENDER + ADM_SOURCE +          TYPE + DISPO + STATUS + MEDICARE + LOS + IN + ED + PCI +          CABG + DIALYSIS + AMI + CHF + PN + COPD + PVD + CVD + DEM +          CPD + RhD + PUD + DIA + HePa + RD + SLVD + AIDS + DRUG +          PHY + DEP + ALCH + PC + SYNC + SEP + ANEM + MALNEO + METCNCR +          COAG + STK + NUTDEF, family = binomial, data = mydata)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.1163  -0.5366  -0.4152  -0.1993   3.3030            Coefficients:                           Estimate Std. Error z value Pr(>|z|)         (Intercept)         -1.689e+00  3.420e-01  -4.938 7.89e-07 ***     ZIP27546            -2.739e-01  3.302e-01  -0.830 0.406772         ZIP28301            -2.902e-01  2.303e-01  -1.260 0.207606         ZIP28302            -7.106e-01  3.673e-01  -1.934 0.053053 .       ZIP28303            -2.099e-01  2.302e-01  -0.912 0.361849         ZIP28304            -2.671e-01  2.296e-01  -1.164 0.244602         ZIP28305            -2.592e-01  2.400e-01  -1.080 0.280171         ZIP28306            -2.382e-01  2.292e-01  -1.040 0.298556         ZIP28307            -3.045e-01  3.266e-01  -0.932 0.351269         ZIP28309            -2.157e-01  3.672e-01  -0.587 0.556900         ZIP28311            -3.222e-01  2.307e-01  -1.396 0.162574         ZIP28312            -3.215e-01  2.318e-01  -1.387 0.165421         ZIP28314            -3.420e-01  2.300e-01  -1.487 0.137042         ZIP28318            -5.677e-01  2.739e-01  -2.073 0.038160 *       ZIP28320            -4.234e-01  3.033e-01  -1.396 0.162745         ZIP28323            -2.376e-01  3.488e-01  -0.681 0.495835         ZIP28326            -4.029e-01  3.134e-01  -1.285 0.198625         ZIP28328            -2.879e-01  2.830e-01  -1.017 0.309107         ZIP28331            -3.942e-01  3.676e-01  -1.073 0.283493         ZIP28334            -5.230e-01  3.004e-01  -1.741 0.081674 .       ZIP28337            -7.301e-01  2.773e-01  -2.633 0.008465 **      ZIP28339            -7.499e-01  5.351e-01  -1.401 0.161082         ZIP28344            -5.818e-01  3.578e-01  -1.626 0.103930         ZIP28348            -2.847e-01  2.305e-01  -1.235 0.216773         ZIP28356            -2.548e-01  2.606e-01  -0.978 0.328270         ZIP28357            -4.823e-01  2.845e-01  -1.696 0.089980 .       ZIP28358            -5.426e-01  3.033e-01  -1.789 0.073592 .       ZIP28360            -7.954e-01  4.333e-01  -1.836 0.066423 .       ZIP28371            -3.188e-01  2.474e-01  -1.289 0.197564         ZIP28372            -8.082e-01  4.966e-01  -1.627 0.103637         ZIP28376            -3.329e-01  2.355e-01  -1.414 0.157397         ZIP28377            -5.883e-01  2.924e-01  -2.012 0.044251 *       ZIP28382            -5.989e-01  2.901e-01  -2.065 0.038942 *       ZIP28384            -3.009e-01  2.497e-01  -1.205 0.228286         ZIP28385            -5.568e-01  3.771e-01  -1.476 0.139834         ZIP28386            -3.303e-01  3.212e-01  -1.028 0.303847         ZIP28390            -3.056e-01  2.364e-01  -1.293 0.196018         ZIP28391            -2.482e-01  2.453e-01  -1.012 0.311606         ZIP28392            -2.709e-01  3.952e-01  -0.686 0.493006         ZIP28395            -4.320e-01  2.730e-01  -1.582 0.113550         ZIP28399             1.235e-01  3.079e-01   0.401 0.688349         ZIP28433            -6.764e-01  3.964e-01  -1.706 0.087937 .       ZIP28441            -6.724e-01  3.890e-01  -1.729 0.083870 .       ZIP99999            -6.239e-01  2.378e-01  -2.623 0.008709 **      WARD1YSU             1.606e-01  2.082e-01   0.772 0.440365         WARD2SOU            -7.323e-02  2.365e+00  -0.031 0.975300         WARD3SOU            -5.522e-04  2.365e+00   0.000 0.999814         WARD3VPN            -4.266e-03  2.366e+00  -0.002 0.998561         WARD3VPS             1.675e-03  2.366e+00   0.001 0.999435         WARD4SOU            -1.826e-01  2.365e+00  -0.077 0.938458         WARD4VPN            -9.387e-02  2.367e+00  -0.040 0.968371         WARD4VPS            -1.046e-01  2.370e+00  -0.044 0.964807         WARD5NOR            -1.089e+00  2.747e+00  -0.397 0.691615         WARD5SOU            -2.720e-02  2.365e+00  -0.012 0.990821         WARD6SOU             1.019e-01  2.365e+00   0.043 0.965635         WARD7SOU            -1.450e-01  2.366e+00  -0.061 0.951123         WARD8SOU             5.411e-02  2.365e+00   0.023 0.981749         WARDAMUN             2.164e-01  2.389e+00   0.091 0.927827         WARDAUNA            -8.279e-02  2.366e+00  -0.035 0.972087         WARDAUNB            -2.720e-01  2.365e+00  -0.115 0.908445         WARDCICU            -1.689e+00  2.397e+00  -0.705 0.481115         WARDCPCU            -2.284e-01  2.365e+00  -0.097 0.923078         WARDCSIC            -6.011e-01  2.380e+00  -0.253 0.800643         WARDCSOH             1.874e+00  2.462e+00   0.761 0.446472         WARDCVOU             1.647e+00  2.547e+00   0.647 0.517838         WARDCVRU            -2.009e-01  2.590e+00  -0.078 0.938169         WARDGATE             1.330e-01  2.370e+00   0.056 0.955249         WARDLDRU             7.628e-01  2.371e+00   0.322 0.747704         WARDMICU            -2.455e+00  2.412e+00  -1.018 0.308850         WARDNEPH             1.173e-01  2.365e+00   0.050 0.960456         WARDOBRR             3.546e+00  2.731e+00   1.298 0.194237         WARDOBSB            -2.015e-01  2.367e+00  -0.085 0.932144         WARDOBSC            -1.468e-02  2.368e+00  -0.006 0.995052         WARDREBI            -1.046e+01  2.272e+02  -0.046 0.963262         WARDRGEN            -9.732e+00  2.276e+02  -0.043 0.965887         WARDSICU            -2.785e+00  2.421e+00  -1.150 0.249997         SVCCAN              -1.211e+01  3.248e+02  -0.037 0.970243         SVCCAR               7.815e-02  2.367e+00   0.033 0.973658         SVCCTS               1.624e-01  2.368e+00   0.069 0.945298         SVCDEN              -9.862e+00  2.292e+02  -0.043 0.965679         SVCDET              -4.493e-01  1.940e-01  -2.316 0.020561 *       SVCEBS              -1.049e+01  1.080e+02  -0.097 0.922647         SVCEDC              -9.275e+00  1.762e+02  -0.053 0.958009         SVCEGN               4.178e-01  2.623e+00   0.159 0.873433         SVCEGS              -7.110e-01  2.605e+00  -0.273 0.784884         SVCEMS               1.271e+00  2.674e+00   0.475 0.634671         SVCENT              -1.088e+00  2.574e+00  -0.423 0.672522         SVCERS              -1.021e+01  1.622e+02  -0.063 0.949780         SVCEYS               5.866e-01  2.727e+00   0.215 0.829683         SVCFAM               1.585e-01  2.365e+00   0.067 0.946573         SVCGAS              -1.160e+00  2.472e+00  -0.469 0.638878         SVCGYN              -1.024e-01  2.368e+00  -0.043 0.965523         SVCHEM               1.286e+00  2.400e+00   0.536 0.592043         SVCLDO              -1.005e+01  3.248e+02  -0.031 0.975306         SVCMED               1.140e-01  2.365e+00   0.048 0.961554         SVCNEP               2.326e-01  2.374e+00   0.098 0.921941         SVCNEU               4.451e-02  2.402e+00   0.019 0.985213         SVCNRS               1.361e-01  2.370e+00   0.057 0.954209         SVCOBA              -9.934e+00  1.621e+02  -0.061 0.951132         SVCOBD              -1.295e+00  2.370e+00  -0.546 0.584748         SVCOBN               2.103e+00  2.368e+00   0.888 0.374650         SVCOBS               6.057e-01  2.404e+00   0.252 0.801028         SVCONC               9.306e-01  2.367e+00   0.393 0.694249         SVCOPH              -9.914e+00  3.248e+02  -0.031 0.975646         SVCORT              -4.223e-01  2.368e+00  -0.178 0.858437         SVCPED               1.826e-01  2.751e+00   0.066 0.947061         SVCPLS               1.412e+00  2.634e+00   0.536 0.591885         SVCPOD              -1.153e-01  2.375e+00  -0.049 0.961262         SVCPUL               5.282e-01  2.380e+00   0.222 0.824376         SVCSUR               4.083e-02  2.366e+00   0.017 0.986228         SVCURO              -4.731e-01  2.378e+00  -0.199 0.842283         AGE                  8.543e-04  1.034e-03   0.827 0.408519         RACEAMERICAN INDIAN  7.895e-02  9.533e-02   0.828 0.407607         RACEASIAN           -1.006e-01  1.444e-01  -0.697 0.485921         RACECAUCASIAN       -3.327e-02  2.774e-02  -1.200 0.230327         RACEOTHER           -1.619e-01  6.037e-02  -2.681 0.007344 **      GENDERM              4.273e-02  2.687e-02   1.590 0.111745         ADM_SOURCEOTHER      1.652e-01  4.486e-01   0.368 0.712637         ADM_SOURCERA        -5.346e-02  4.963e-02  -1.077 0.281430         TYPEJ               -3.104e-02  2.346e-01  -0.132 0.894746         TYPEQ                4.522e-02  2.310e-01   0.196 0.844763         DISPOHOME           -4.442e-01  3.863e-02 -11.499  < 2e-16 ***     DISPOHOSPICE        -8.949e-01  9.628e-02  -9.295  < 2e-16 ***     DISPOOTHER          -9.028e-01  5.873e-02 -15.370  < 2e-16 ***     DISPOSNF            -4.500e-02  4.718e-02  -0.954 0.340221         STATUSMARRIED       -1.442e-01  3.970e-02  -3.631 0.000282 ***     STATUSOTHER          4.810e-02  6.469e-01   0.074 0.940723         STATUSSEPERATED      4.294e-02  6.882e-02   0.624 0.532630         STATUSSINGLE        -6.439e-02  4.105e-02  -1.568 0.116793         STATUSWIDOWED       -8.910e-02  4.675e-02  -1.906 0.056670 .       MEDICARE             1.305e-01  2.955e-02   4.417 1.00e-05 ***     LOS                  1.800e-02  1.815e-03   9.915  < 2e-16 ***     IN                   2.472e-01  1.237e-02  19.994  < 2e-16 ***     ED                   3.874e-02  5.256e-03   7.370 1.71e-13 ***     PCI                 -6.748e-03  6.773e-02  -0.100 0.920647         CABG                -8.911e-02  1.604e-01  -0.555 0.578575         DIALYSIS            -2.061e-01  1.090e-01  -1.891 0.058589 .       AMI                  4.675e-02  6.219e-02   0.752 0.452175         CHF                  1.293e-01  3.723e-02   3.473 0.000515 ***     PN                  -1.917e-02  4.010e-02  -0.478 0.632604         COPD                 7.144e-02  5.151e-02   1.387 0.165506         PVD                 -1.145e-02  7.589e-02  -0.151 0.880032         CVD                 -3.169e-02  7.305e-02  -0.434 0.664416         DEM                 -1.982e-01  1.141e-01  -1.737 0.082310 .       CPD                  7.009e-02  4.825e-02   1.453 0.146283         RhD                  1.712e-01  7.599e-02   2.253 0.024249 *       PUD                 -4.804e-02  7.202e-02  -0.667 0.504767         DIA                  1.380e-01  3.193e-02   4.321 1.55e-05 ***     HePa                 3.654e-02  1.240e-01   0.295 0.768183         RD                   2.044e-01  8.159e-02   2.505 0.012246 *       SLVD                 1.630e-03  5.426e-02   0.030 0.976036         AIDS                 1.104e-02  1.582e-01   0.070 0.944372         DRUG                 2.828e-01  4.687e-02   6.033 1.61e-09 ***     PHY                  1.747e-02  4.766e-02   0.366 0.714022         DEP                  8.486e-02  3.800e-02   2.233 0.025548 *       ALCH                 2.984e-02  1.049e-01   0.284 0.776109         PC                  -8.581e-02  1.433e-01  -0.599 0.549288         SYNC                -4.742e-02  8.126e-02  -0.584 0.559487         SEP                 -5.585e-02  4.355e-02  -1.283 0.199657         ANEM                 1.775e-01  3.449e-02   5.147 2.65e-07 ***     MALNEO               9.727e-02  6.513e-02   1.493 0.135309         METCNCR              9.645e-02  9.398e-02   1.026 0.304741         COAG                 1.308e-01  4.110e-02   3.184 0.001454 **      STK                  3.670e-02  8.925e-02   0.411 0.680945         NUTDEF               4.451e-02  4.605e-02   0.967 0.333744         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 54757  on 67295  degrees of freedom     Residual deviance: 46202  on 67131  degrees of freedom     AIC: 46532          Number of Fisher Scoring iterations: 11      ROC:               Call:     roc.formula(formula = RA ~ prob, data = mydata)          Data: prob in 57805 controls (RA 0) < 9491 cases (RA 1).     Area under the curve: 0.7731      Using lrm I get:                myglm <- lrm(RA~ZIP+WARD+SVC+AGE+RACE+GENDER+ADM_SOURCE+TYPE+DISPO+STATUS+MEDICARE+LOS+IN+ED+PCI+CABG+DIALYSIS+AMI+CHF+PN+COPD+PVD+CVD+DEM+CPD+RhD+PUD+DIA+HePa+RD+SLVD+AIDS+DRUG+PHY+DEP+ALCH+PC+SYNC+SEP+ANEM+MALNEO+METCNCR+COAG+STK+NUTDEF,x=TRUE,y=TRUE,data=mydata)      I get the following output:                                    Model Likelihood     Discrimination    Rank Discrim.                                 Ratio Test            Indexes          Indexes            Obs        67296    LR chi2    8554.47    R2       0.214    C       0.773          0         57805    d.f.           164    g        1.246    Dxy     0.546          1          9491    Pr(> chi2) <0.0001    gr       3.477    gamma   0.549         max |deriv| 0.07                          gp       0.129    tau-a   0.132                                                   Brier    0.103                         when I try to run:               library(rms)     myglm <- lrm(RA~ZIP+WARD+SVC+AGE+RACE+GENDER+ADM_SOURCE+TYPE+DISPO+STATUS+MEDICARE+LOS+IN+ED+PCI+CABG+DIALYSIS+AMI+CHF+PN+COPD+PVD+CVD+DEM+CPD+RhD+PUD+DIA+HePa+RD+SLVD+AIDS+DRUG+PHY+DEP+ALCH+PC+SYNC+SEP+ANEM+MALNEO+METCNCR+COAG+STK+NUTDEF,x=TRUE,y=TRUE,data=mydata)     my.valid <-validate(myglm,method="boot",B=1000)     my.calib <-calibrate(myglm,method="boot",B=1000)     par(bg="white", las=1)     plot(my.calib, las=1)      I get a singular information matrix error on SVC. Is my only option to remove the SVC variable? I do not receive any such error when using glm. Thanks!


title: How and why would MLPs for classification differ from MLPs for regression? Different backpropagation and transfer functions?

text: I'm using two 3-layer feedforward multi-layer perceptrons (MLPs). With the same input data (14 input neurons), I do one classification (true/false), and one regression (if true, "how much")¹. Until now, I've lazily used Matlabs patternnet and fitnet, respectively. Lazily, because I haven't taken the time to really understand what's going on — and I should. Moreover, I need to make the transition to an OSS library (probably FANN), that will likely require more manual setup than the Matlab NN Toolbox. Therefore, I'm trying to understand more precisely what's going on. The networks created by `patternnet` and `fitnet` are nearly identical: 14 inputs neurons, 11 hidden neurons, 1 target neuron (2 for the `fitnet`, but only 1 piece of informatio). But, they're not completely identical. The differences by default are:   * Matlab uses a Scaled conjugate gradient backpropagation for the classification network (`patternnet`), and Levenberg-Marquardt backpropagation for the regression network (`fitnet`).   * The classification network uses a hyperbolic tangent sigmoid transfer function between the input and the hidden layer, and between the hidden and the output layer. The regression network (`fitnet`) uses the Hyperbolic tangent sigmoid transfer function between the input and the hidden layer, and a purely linear transfer function between the hidden and the output layer. Should those differences be? What kind of backpropagation functions are optimal for classification, and what kind for regression, and why? What kind of transfer functions are optimal for classification, and what kind for regression, and why? * * * ¹The classification is for "cloudy" or "cloud-free" (2 complementary targets), the regression is for quantifying "how much cloud" (1 target).


title: Confusion related to the dual being concave

text: I have some confusion related to the dual function or the Lagrangian being concave for a general convex optimization problem with objective function $f_0(x)$ and equality and non equality constraints given by $h_j$ and $f_i$ $$ g(\lambda, \mu) = \text{inf}_{x \epsilon D} L(x,\lambda,\mu) = \text{inf}_{x \epsilon D}(f_o(x) + \sum_{i=1}^{m}\lambda_if_i(x) + \sum_{j=1}^{m'}\mu_ih_j(x)) $$ Why is it concave? Further, they have said that it is a pointwise infimum of a set of affine functions. I know what affine function is but what is pointwise infimum in this case? I thought we were trying to minimize w.r.t. $x$, so that is why we have $\text{inf}_x$. Can anyone please clarify? Well after some googling I found how they explained it. Infimum which means we are trying to minimize over the domain of x which means all possible xs in the domain D. So the infimum over the family of affine functions if concave. Let me elaborate if I have understood it correct. Lets jus say x=1,2 in the domain D. Then I have two family of affine functions of $\lambda, \mu$. So I am taking infimum over these two families one with x=1 and the other with x = 2. However, my question is are we taking infimum/ minimizing over $\lambda,\mu$ or x here? As far as I know we are minimizing over x not $\lambda, \mu$


title: Unstandardized effect size for two non-equal size samples for experimental design in R

text: I am trying to calculate the required (unstandardized) difference in between a test and control group of different sizes. I have historical data with transactions for   * the test group (n=272, mean transactions = 63.21, sum of transactions = 17194)    * the control group (n=831, mean transactions = 54.12, sum of transactions = 44971) The pooled standard deviation of transactions is 257.436. I am now trying to calculate the required difference in means for an experiment that exposes the test group with treatment that is supposed to increase the no. of transactions. In R, I am calculating:               power_test <- pwr.t2n.test(n1=272,              n2=831,              sig.level=0.05,              power=0.8              )      Which yields a standardized effect size of 0.196. As need the unstandardized effect size, I am multiplying with the pooled standard deviation:               d_unstandardized <- as.numeric(power_test["d"])*257.436      Which yields 50.42. How can this be interpreted? For a test with high power and significance, a total of 50 additional transaction in the test group would suffice? This is surprising, as the sums of the transactions are very high compared to this small estimated difference. Help is appreciated!


title: CAR spatial models in JAGS

text: WinBUGS comes with the GeoBUGS add-on, which contains a number of predefined model structures that are suitable for modelling spatial data structures e.g. geostatical structures (spatial.exp), conditional auto-regressive, CAR, structures (car.proper, car.normal) amongst others. Unfortunately there are no comparable structures for JAGS. This discussion: http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/280820cb on the mailing list indicates that I probably shouldn't hold my breath either. How then can I implement a CAR model in JAGS? Does anyone have examples that they could share? Cheers, Mark


title: Probability of dependent second event without knowing result of first event?

text: I have this question in my textbook I am having trouble answering. You have a fleet of 19 vehicles, 4 of which are vans and 15 of which are cars. Say you need to choose two vehicles at random. Regardless of which vehicle you choose first, what is the probability that second vehicle you choose is a van? I've tried several answers but none of them make sense. The one I thought was correct is to sum the two possible ways that the second vehicle be a van. $$(4/19)\times(3/18) + (15/19)\times(4/18)$$ But I'm told that's wrong. Could somebody please explain the way to do this? Thank you very much


title: Validating a model with 2 random effects and fixed effects of different levels

text: I have been working on this problem for the last little while and trying to write out the code to the below math/description for an ecological dataset in `lmer`, but I wanted to check with a few people to see if the code I've written out adequately reflects the math (which I am certain of!), as I could not find any equivalent examples in Gelman and Hill, 2007 or Zuur, 2009, and no one in my department has dealt with a a similar problem in `lmer` before. Thank you in advance for taking the time to look at this, I really do appreciate it very much! For individual i of species j in lake k,               Hgijk ~ b0j[i] + b1k[i] + b2j[i]*lengthi + ei     ei ~ N(0,s2e)     b0j ~ N(mb0,s2b0)     b1k ~ go + g1*chlak + eb1     eb1 ~ N(0,s2b1)     b2j ~ N(mb2,s2b2)      This model says that the Hg for fish i of species j from lake k is a function of the species, the lake, and the length of the fish. The species effects are like intercepts – the average Hg for a particular species at average length (if the length data are centered) in an average lake – and they come from a normal distribution. The lake effect is itself a linear function of chlorophyll concentration, so that lakes with higher chlorophyll can potentially have higher or lower average Hg concentrations. The length effects vary by species, on the assumption that Hg concentration may change with size differently in different species. These length effects come from a normal distribution. So in `lmer` I wrote it out like this:               M1<-lmer(THg~Chlo-a+CenteredLength+ (1|Lake) + (1+CenteredLength|Species), data=data)      Is this right? Does this reflect the math I wrote out above? Would this type of fully nested or partially nested structure be better?               M2<-lmer(THg~Chlo-a+CenteredLength+(1|Lake:Species)+ (1|Lake) + (1+CenteredLength|Species), data=data)     M3<-lmer(THg~Chlo-a+CenteredLength+(1+CenteredLength|Lake:Species)+ (1|Lake), data=data)      Also, if any of this is correct, when running these models, I come up with the problem that Length and Species have a correlation of 1.00. Based on this post, this makes be believe that the model without the varying slopes for length would be better. Is this an ok way to conclude this?


title: Heteroskedasticity in residuals vs. fitted plot

text: I am testing whether price per ounce of beer (continuous variable, range of values mostly between 0.1 and 0.5 dollars) and the presence of promotion, advertisement, and display (all binary) have effect on the total amount of ounces purchased (continuous variable). Here is my residual vs. fitted plot before the log transformation of y: ![before log transformation of y](http://i.stack.imgur.com/fUFfl.png) This is the residuals vs. fitted plot after the log transformation of y: ![after log transformation of y](http://i.stack.imgur.com/DdqM0.png) Heteroskedasticity is very high (White's general t statistics is nearly 800). This is the histogram of the transformed y: ![enter image description here](http://i.stack.imgur.com/C2bgi.jpg) Any ideas or suggestions on how to improve my model or where to look for errors in order to improve the problem of heteroskedasticity are greatly appreciated.


title: How to analyse my time series experiment

text: I have a question about how I can analyze the data from a time series experiment. I'm still a student and therefore am not really familiar with all the procedures and analyzes in SPSS, so i hope someone with the proper knowledge can help me with this problem. The goal is to show a relation between advertising en sales. To show this we have conducted a experiment with an experimental and control group. The experimental group is a province in which the people gets to see more advertising. In the control group the advertising stays the same as usual.   The dataset contains weekly sales data from several stores throughout the country from the control area and experimental area. And how much advertising is used in all the periods. The data goes back for 1.5 years so there are 70 measures before the experiment, and 10 measures in/after the experiment. So my question is, how can I prove a relation between those two factors. I know I can do a simple t-test to test the difference between the two groups. But this doesn’t take care of any trends of long/short term effects.


title: Should the mean of the bootstrapped distribution always be asymptotically equal to the sample estimate?

text: Suppose I bootstrap the distribution of the sample mean. Normally, one would use the mean of the bootstrapped distribution as point estimate of the parameter and the s.d. as its standard error. The mean of the bootstrapped distribution is asymptotically equal to the sample estimate (i.e. for a large number of iterated draws). Now suppose the mean, or more generally, parameter, has an asymetrical bootstrap distribution, so that the sample estimate of the parameter and the mean of the bootstrap distribution are likely to be unequal for a moderate amount of iterated draws. Should I still expect both to be asympotically equivalent, so as I increase the number of draws, both will be equal? And if this is so, is it customary practice to increase the number of iterations until they are equal before I report statistics? In my practical case, both deviate after 1000 iterated draws. So I am unsure whether I should report the sample estimate or the mean of the bootstrap distribution of the parameter.


title: R-squared of a vector error correction model

text: I hope somebody can help me with this. I estimated vector error correction model (and Johansen's cointegration test) and my R-squared is only 23%. Does R-squared tell me my model is bad in this case too?


title: Improving SVM classification

text: I have a classification problem (bioinformatics domain) where I have around 333 features. Currently, I am first selecting features (using importance feature of random forest) and then pushing the same through RBF kernel for SVM. On the results front, I have section of data which gives bad results on the classifier. There is feature rscindex, which separates the data into well classified data and incorrectly classified data,i.e data points with rscindex < 0.8 are performingly badly during classification stage. Please give me ideas on how to improve my results on this badly performing section of data. Will using Fisher kernel help (if so is there a library I can use) or is there a sampling technique that I should be looking into.? Please help. Thanks


title: Can I use optimally scaled variables for a factor analysis to account for rotation? If I can then how?

text: I have discussed this issue several times in this site, but I am asking it again for a final justification from the experts of our community. I wanted to extract four factors (I should call dimensions here I think) from a CATPCA along with the factor scores (note that, factor scores are not available if I use polychoric correlation matrix and factor analyze it) from a data set containing only categorical variables. I have 22 categorical variables overall. Now, as far I know, CATPCA doesn't use rotation which can be handy to extract factors in standard factor analysis or PCA. I am having a little problem to get clear factors although it seems like a rotation may improve understanding of the factors. Or, in simple terms, I am a little desperate to extract neat factors so that I want to employ different rotations. (I should not be though!) :D However, we can save the transformed variables in some out file and use them in a standard factor analysis as the variables are now quantified. Note that we need to mention the number of dimensions once while getting the optimally scaled (i.e. quantified) variables. I have reckoned that different number of dimensions mentioned produces different quantification. Using these quantified variables into further factor analysis will require to mention the number of factors again. Is this a problem? Because as I want 4 factors to be extracted, while asking for the optimally scaled variables I am mentioning 4 dimensions and then planning to use the quantified variables in a standard factor analysis (to get the facility of rotation) where I need to mention the number of factors once again. If this is a problem, then can we take the number of dimensions to be 22 for CATPCA to get optimally scaled variables (as there are 22 variables overall)? This is only to make sure we have no loss in percentage of variance explained. Then we can use the quantified variables from that for a standard factor analysis and get factor scores in SPSS with mentioning 4 factors this time. I totally don't know if it is a weird idea! So, waiting for your kind direction.


title: Why applying model selection using AIC gives me non-significant p-values for the variables

text: I have some questions about the AIC and hope you can help me. I applied model selection (backward, or forward) based on the AIC on my data. And some of the selected variables ended up with a p-values > 0.05. I know that people are saying we should select models based on the AIC instead of the p-value, so seems that the AIC and the p-value are two difference concepts. Could someone tell me what the difference is? What I understand so far is that:   1. For backward selection using the AIC, suppose we have 3 variables (var1, var2, var3) and the AIC of this model is AIC*. If excluding any one of these three variables would not end up with a AIC which is significantly lower than AIC* (in terms of ch-square distribution with df=1), then we would say these three variables are the final results.   2. A significant p-value for a variable (e.g. var1) in a three variable model means that the standardized effect size of that variable is significantly different from 0 (according to Wald, or t-test). What's the fundamental difference between these two methods? How do I interpret it if there are some variables having non-significant p-values in my best model (obtained via the AIC)?


title: Confirming the detailed calcs for NPS Margin of error

text: Based on this very good post How can I calculate margin of error in a NPS (Net Promoter Score) result? I've pulled together the detailed calculations to perform the test for the generalised NPS problem. I suspect that this question is lower order than generally posted here but my last true stats course was 20 years ago and I'm hoping that someone will validate the steps. To determine if there is a significant (95%) probability of a real difference between two sets of NPS scores: Where:               #P is the number of Promoters     #N is the number of Neutrals     #D is the number of Detractors.     #T = #P + #N + #D      First calculate NPS in the normal way:               NPS = #P/#T - #D/#T      Now determine the Variance of the sample NPS:               Var(NPS) = (1 – NPS) ^ 2) * #P/#T + (0-NPS) ^ 2 * #N/#T + (-1-NPS) ^ 2 * #D/#T)      Now calculate the Margin of Error (MoE) for the sample:               MoE = SQRT(Var(NPS)) / SQRT(#T)      If you do this for both of your samples you will end up with two MoE, one for each sample. To determine if your score probably (95%) changed:               If ABS(NPS1 – NPS2) > 2 * SQRT (MoE1 ^ 2 + MoE2 ^ 2) then the difference is real      Thanks.


title: Comparing data files

text: We have several weather files for a year's data sampled hourly. In each we have several variables (up to ten), temperature, wind speed, solar intensity etc. I would like to try and develop a system that could determine if these files are statistically different. If I was using a single variable, I could use a KS-test, but with a collection of variables how should I approach this problem?


title: Prediction with not significant covariate in logistic regression

text: I have a logistic regression model with several variables and one of those variables (called x3 in my example below) is not significant. However, x3 should remain in the model because it is scientifically important. Now, x3 is continuous and I want to create a plot of the predicted probability vs x3. Even though x3 is not statistically significant, it has an effect on my outcome and therefore it has an effect on the predicted probability. This means that I can see from the graph, that the probability changes with increasing x3. However, how should I interpret the graph and the change in the predicted probability, given that x3 is indeed not statistically significant? Below is a simulated data in R set to illustrate my question. The graph also contains a 95% confidence interval for the predicted probability (dashed lines):               > set.seed(314)     > n <- 300     > x1 <- rbinom(n,1,0.5)     > x2 <- rbinom(n,1,0.5)     > x3 <- rexp(n)     > logit <- 0.5+0.9*x1-0.5*x2     > prob <- exp(logit)/(1+exp(logit))     > y <- rbinom(n,1,prob)     >      > model <- glm(y~x1+x2+x3, family="binomial")     > summary(model)          Call:     glm(formula = y ~ x1 + x2 + x3, family = "binomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.0394  -1.1254   0.5604   0.8554   1.4457            Coefficients:                 Estimate Std. Error z value Pr(    >|z|)         (Intercept)   1.1402     0.2638   4.323 1.54e-05 ***     x1            0.8256     0.2653   3.112  0.00186 **      x2           -1.1338     0.2658  -4.266 1.99e-05 ***     x3           -0.1478     0.1249  -1.183  0.23681         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)          Null deviance: 373.05  on 299  degrees of freedom     Residual deviance: 341.21  on 296  degrees of freedom     AIC: 349.21          Number of Fisher Scoring iterations: 3          >      > dat <- data.frame(x1=1, x2=1, x3=seq(0,5,0.1))     > preds <- predict(model, dat,type = "link", se.fit = TRUE )     > critval <- 1.96     > upr <- preds$fit + (critval * preds$se.fit)     > lwr <- preds$fit - (critval * preds$se.fit)     > fit <- preds$fit         >          > fit2 <- mod$family$linkinv(fit)         > upr2 <- mod$family$linkinv(upr)         > lwr2 <- mod$family$linkinv(lwr)         >          > plot(dat$x3, fit2, lwd=2, type="l", main="Predicted Probability", ylab="Probability", xlab="x3", ylim=c(0,1.00))     > lines(dat$x3, upr2, lty=2)         > lines(dat$x3, lwr2, lty=2)      ![enter image description here](http://i.stack.imgur.com/ljW7W.png) Thanks! Emilia


title: Is it possible to build a more "controllable" decision tree like below?

text: I have 2 real and 1 discrete input variable whereas the output variable takes either of the 2 nominal values (i.e. 2 class problem). First I used Weka to train C 4.5 decision tree in a 10-fold cross validation setup. As a result I get a single decision tree whose True positive, False positive rates don't change much with the pruning factor. What I want is a single variable whose value I can threshold (in addition to class identified by decision tree) to do classification. At the same time, I want readability of decision trees. Is it possible to do this with regression trees ? Any starting points ? Basically I want to know if it's possible to build a decision tree like in below image - ![enter image description here](http://i.stack.imgur.com/HJNFC.jpg)


title: Analyzing multiply imupted data from Amelia in R: Why do results from zelig and mice differ?

text: Motivated by my answer to this question, I played around with analyzing mulitply imputed data from the `Amelia` package in R. As I have explained in my answer, the multiply imputed datasets can be analyzed using the combined `Zelig` and `mitools` packages or using a combination of `Zelig` and `mice`. Now, to me it seemed rather inconvenient to fit the linear model using `zelig()` as `mice`, too, provides the `with.mids()`-function to fit linear models to multiply imputed datasets. However, I found that the results differ depending on the function used for fitting. For the analysis using `with.mids`, I first had to circumvent a bug in the current `mice`-package by defining the following function, as has been explained in another question:               as.mids2 <- function(data2, .imp=1, .id=2){       ini <- mice(data2[data2[, .imp] == 0, -c(.imp, .id)], maxit=0)       names  <- names(ini$imp)           if (!is.null(.id)){             rownames(ini$data) <- data2[data2[, .imp] == 0, .id]           }           for (i in 1:length(names)){             for(m in 1:(max(as.numeric(data2[, .imp])) - 1)){               if(!is.null(ini$imp[[i]])){                 indic <- data2[, .imp] == m & is.na(data2[data2[, .imp]==0, names[i]])                 ini$imp[[names[i]]][m] <- data2[indic, names[i]]           }         }       }       return(ini)     }      Once I had done this, I used `Zelig`:               library("Amelia")     data(freetrade)     amelia.out <- amelia(freetrade, m = 15, ts = "year", cs = "country")          library("Zelig")     zelig.fit <- zelig(tariff ~ pop + gdp.pc + year + polity, data = amelia.out$imputations, model = "ls", cite = FALSE)     zelig.results <- lapply(zelig.fit, function(x) x$result)          library("mice")     zelig4mice <- as.mira(zelig.results)     zelig.mice.res <- summary(pool(zelig4mice, method = "rubin1987"))      Then I tried the same thing using only `mice`:               imp.data <- do.call("rbind", amelia.out$imputations)     imp.data <- rbind(freetrade, imp.data)     imp.data$.imp <- as.numeric(rep(c(0:15), each = nrow(freetrade)))     mice.data <- as.mids2(imp.data, .imp = ncol(imp.data), .id = NULL)          mice.fit <- with(mice.data, lm(tariff ~ polity + pop + gdp.pc + year))     mice.res <- summary(pool(mice.res2, method = "rubin1987"))      These are the results:               > zelig.mice.res                       est       se     t    df Pr(>|t|)     lo 95     hi 95 nmis   fmi lambda     (Intercept)  3.18e+03 7.22e+02  4.41  45.9 6.20e-05  1.73e+03  4.63e+03   NA 0.571  0.552     pop          3.13e-08 5.59e-09  5.59 392.1 4.21e-08  2.03e-08  4.23e-08   NA 0.193  0.189     gdp.pc      -2.11e-03 5.53e-04 -3.81 329.4 1.64e-04 -3.20e-03 -1.02e-03   NA 0.211  0.206     year        -1.58e+00 3.63e-01 -4.37  45.9 7.11e-05 -2.31e+00 -8.54e-01   NA 0.570  0.552     polity       5.52e-01 3.16e-01  1.75  90.8 8.41e-02 -7.58e-02  1.18e+00   NA 0.406  0.393          > mice.res                       est       se     t     df Pr(>|t|)     lo 95     hi 95 nmis    fmi lambda     (Intercept)  3.42e+03 8.87e+02  3.86   8.01 4.80e-03  1.38e+03  5.47e+03   NA 0.7599 0.7066     pop          3.20e-08 5.25e-09  6.10 504.30 2.10e-09  2.17e-08  4.24e-08    0 0.0927 0.0891     gdp.pc      -2.09e-03 5.31e-04 -3.93 189.23 1.19e-04 -3.13e-03 -1.04e-03    0 0.1543 0.1454     year        -1.70e+00 4.46e-01 -3.83   8.02 5.02e-03 -2.73e+00 -6.78e-01    0 0.7594 0.7061     polity       5.74e-01 3.60e-01  1.59  13.93 1.34e-01 -2.00e-01  1.35e+00    2 0.5907 0.5358      From these data it is apparent, that the linear models fit by the two methods differ and so do the determined degrees of freedom. **Why do these results differ? What is the correct analysis procedure?**


title: Extreme value theory for count data

text: I am aware of extreme value theory for continuous distributions. I need to fit an extreme value distribution to the maximum observation of number of events on a day, per month. This seems to be the block maxima problem, which is approximated by the GEV family of distributions for continuous distributions. How do I do this for count data? As a secondary question, let's assume the basic count process is ~ Poisson. Then does this lead to a different answer to the original question? Thanks!


title: Adjusting data for measurement errors

text: I have approximately 500,000 data points that were mechanically measured over about 2000 days. The equipment is manually set up each day, which likely induces some error. What I am attempting to do is use a linear mixed effects model to estimate the fixed effect of each day's potential measurement error while accounting for other known random effects and covariates, and then adjust that day's measurement based on the estimated effect. I am using lme4 in R. Given that lme4 cannot estimate fixed effects with so many levels, I have broken my data into 50 day chunks. I then compare each day to 50 other randomly selected days that I treat as one level. My first question is: is this a reasonable thing to do? The only other thread I could find dealing with such an issue regarding lm4 suggested that the author probably didn't actually want to estimate a model with so many levels. Given that each day could theoretically have a different fixed measurement error based on how the equipment is set up though, it seems to make sense to me in this case. My second question is: assuming I have reasonable estimates of each day's fixed effect, how would I go about adjusting my data? If I was just trying to determine which day's effects are likely to be significant I would perform some adjustment for multiple-comparisons. As I see it, I am trying to minimize the error over two sources - the measurement error from the equipment, and the error in estimating each fixed effect.


title: Question about a normal equation proof

text: How can you prove that the normal equations: $(X^TX)\beta = X^TY$ have one or more solutions without the assumption that X is invertible? My only guess is that it has something to do with generalized inverse, but I am totally lost.


title: biological replications with one extreme outlier

text: My problem is I have 3 biological replications (reps) withe each having 4-5 technical reps with two of the biological reps having comparable results for the treatment and control. These values represent the colonization of donor stem cells. (Treatment is about 50% of control) whereas in the third rep I still have ~50% drop from control but the values for both treatment and control are drastically reduced. Is there a statistical method to analyze this data because using a t-test obviously results in a non- significant p value? Or would it be best to just do an additional rep and throw out the last rep as an outlier as results from other experiments are comparable to the last 2 reps? Example - Biological reps below, average of technical reps   Rep 1 control 300 treatment 150   Rep 2 control 230 treatment 112   Rep 3 control 24 treatment 11


title: Correlating ranked lists

text: Let's assume I have 10 users who rank a list by preference.   1. A   2. B   3. C etc. My questions are:   1. What is the best statistical method to find which lists are the most correlated? If they are the same, 100% correlation, or 1 item switched with another in the list(say they swapped 3 and 4) are also very close, but what if 8/10 are exact but 1 list has the most important and the other list has the same as the least important(i.e, switched 1 and 10 but the rest match).   2. Beyond that, what if I also wanted to weight each category? So if each list entry was of specific importance, how could I weight each one as a percentage of the total? Say I think "A" is by far the most important so I want it to carry 30% of the total list weight, for example.


title: Testing whether error variances of two regression models are equal

text: Is it possible to test the following? Assume you have two linear regression models, one regressing $Y_1$ on $X_1$ the other $Y_2$ on $X_2$. This gives error variance $\sigma_1^2$ and $\sigma_2^2$. From a sample all variables are available. I would like to evaluate $$H_0: \sigma_1^2=\sigma_2^2,$$ or one- sided $$H_0: \sigma_1^2 \le \sigma_2^2.$$


title: Compare effects of paired explanatory variables within a (binomial) GLMM

text: I am using a GLMM to model the probability of visiting a 400m radius area. I have 135,000 observations (areas) and for each area I have the proportion of the habitat types within it. I want to examine the effect of habitat type on probability of visitation. However, I have split the habitat types into whether they are protected (sssi) or not, and I want to see how the influence of EACH habitat type on probability of visitation changes when protected. The approach I have taken is to split each habitat type into unprotected (not.sssi) and protected (sssi) and included all as explanatory variables in one model. Then I've plotted a clustered bar chart of coefficients to see how they change for each habitat (e.g. broadleaf.sssi vs broadleaf.not.sssi). But now I want to know if there is a 'significant' change in the influence of the habitat when protected vs not. A small sample of data is below, but a larger sample with more variables is available at http://pastebin.com/Dn3A4JNw The sample below doesn't include my region factors for the mixed effects model as there are not enough observations.               structure(list(visit = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), broadleaf.not.sssi = c(0,      0.004962779, 0.178082192, 0, 0, 0.019900498, 0, 0, 0.1105919,      0.00621118, 0.067247821, 0, 0.133333333, 0.350310559, 0.326276463,      0.117647059, 0.007453416, 0.217391304, 0, 0.546134663, 0.034825871,      0, 0.485714286, 0, 0.073474471, 0, 0.01120797, 0.00872818, 0,      0.12923463, 0, 0.080745342, 0, 0.04691358, 0, 0, 0, 0.114035088,      0, 0.232587065, 0, 0.162531017, 0, 0, 0.080946451, 0.047204969,      0, 0.260325407, 0, 0, 0, 0.199751861, 0, 0, 0, 0.376558603, 0.020075282,      0.108074534, 0, 0.009937888, 0.054794521, 0, 0, 0.054794521,      0.049689441, 0.049689441, 0.045996593, 0.054794521, 0.10931677,      0.127791563, 0.02484472, 0, 0.02484472, 0, 0.031094527, 0.008684864,      0.348692403, 0, 0.506815366, 0, 0.118453865, 0.074027604, 0.02864259,      0.013647643, 0.164383562, 0, 0.136645963, 0.067924528, 0.023602484,      0.067924528, 0.001243781, 0.018656716, 0.530510585, 0.001242236,      0.058457711, 0, 0, 0.022332506, 0.032338308, 0), conifer.not.sssi = c(0,      0, 0, 0, 0, 0, 0, 0, 0.00623053, 0, 0, 0.00248139, 0, 0, 0, 0,      0, 0, 0, 0.03117207, 0, 0, 0.514285714, 0, 0, 0, 0.01743462,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.080845771, 0.022332506,      0, 0, 0, 0, 0.01242236, 0, 0, 0, 0, 0, 0, 0, 0, 0.28358209, 0.163341646,      0.003764115, 0, 0, 0, 0, 0.018656716, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0.535491905, 0, 0.120198265, 0, 0.063591022,      0, 0, 0, 0.01743462, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00996264,      0, 0, 0, 0), broadleaf.sssi = c(0, 0, 0, 0, 0, 0, 0, 0, 0.102803738,      0, 0, 0, 0, 0, 0, 0.013767209, 0, 0, 0, 0.295511222, 0, 0, 0,      0, 0, 0, 0, 0.002493766, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.064676617, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.02739726,      0, 0, 0, 0.264339152, 0, 0.0124533, 0, 0, 0, 0, 0, 0, 0, 0.047263682,      0, 0, 0.047204969, 0, 0, 0, 0, 0, 0), conifer.sssi = c(0, 0,      0, 0, 0, 0, 0, 0, 0.020249221, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.211442786, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0.052369077, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0)), .Names = c("visit", "broadleaf.not.sssi", "conifer.not.sssi",      "broadleaf.sssi", "conifer.sssi"), row.names = c(1L, 2L, 3L,      4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L,      18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L,      31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L,      44L, 45L, 46L, 47L, 48L, 49L, 50L, 127894L, 127895L, 127896L,      127898L, 127900L, 127901L, 127902L, 127903L, 127905L, 127906L,      127907L, 127910L, 127913L, 127914L, 127915L, 127916L, 127917L,      127918L, 127920L, 127921L, 127922L, 127923L, 127924L, 127925L,      127927L, 127928L, 127929L, 127930L, 127932L, 127933L, 127935L,      127940L, 127941L, 127942L, 127943L, 127944L, 127945L, 127946L,      127947L, 127948L, 127949L, 127950L, 127951L, 127952L, 127953L,      127954L, 127956L, 127957L, 127961L, 127963L), class = "data.frame")      My code is as follows (with all variables found in the pastebin):               library(MASS)     sssi<-read.table("http://pastebin.com/Dn3A4JNw")     glmmPQL.sssi<-glmmPQL(visit ~ broadleaf.not.sssi+conifer.not.sssi+arable.not.sssi+impro.grass.not.sssi+broadleaf.sssi+conifer.sssi+arable.sssi+impro.grass.sssi+path.lgth+mean.elev+random=~pop.sqrd|region, family=binomial, data=sssi)      Can I use the coefficients and standard errors to test for a significant difference between each pair ("habitat.type".sssi vs "habitat.type".not.sssi). Or as someone else suggested should I restructure my model, entering habitat type without splitting by protection status and using an interaction with 'protected or not'? Our thought about this was that if we just include the proportion protected within the 400m radius area, this isn't indicating which habitat type is within the protected area, and this is our question, hence measuring this explicitly. I am using R version 3.0.2. on Windows 7


title: Classifier success rate and confidence intervals

text: Suppose we measure the classifier error on a test set and obtain a certain success rate - say, 75%. Now, of course, this is only one measurement - how to calculate the "true" success rate? Sure it will be close to 75% but how close? I understand it's related to confidence intervals but now I'm lost in confidence intervals. I think my example is similar to this one on wikipedia where they look at weight distribution of margarine cups. (Sorry, math is not rendering here so I created a screenshot - you might also want to flick through the corresponding section in the wikipedia article).        ![enter image description here](http://i.stack.imgur.com/YZhf3.png) I have the following questions:   1. Why they use the above standard error formula?   2. Where does this Ф^{-1}(0.975)=1.96 come from?   3. To solve my "true success rate" problem, should I repeat the estimation N times and then apply the same reasoning as they do with margarine cups?


title: Dynamic recommender systems

text: > **Possible Duplicate:**   >  Dynamic recommender systems A Recommender System would measure the correlation between ratings of different users and yield recommendations for a given user about the items which may be of interest to him. However, tastes change over time so _old ratings_ might not reflect _current preferences_ and vice versa. You may once have put "excellent" to a book you would now rate as "not too disgusting" and so on. Moreover, the interests themselves do change as well. **How should recommender systems work in a changing environment?**   1. One option is to cut off the "old" ratings, which may work just fine assuming you correctly define "old" (you can even say ratings never expire and pretend that the problem doesn't exist). But it's not the best possible option: of course tastes evolve, it's a normal life flow, and there's no reason why we cannot use the extra knowledge of once correct past ratings.    2. Another option is to somehow accommodate this extra knowledge. Thus we could not just find an "instant match" for your current interests but suggest you the things you may like _next_ (as opposed to the things you may like _now_ ). I'm not sure if I'm explaining this well enough. Basically I'm in favor of the second approach and am talking about a Recommender System which would measure the correlations of taste _trajectories_ and yield recommendations which would cater for.. well, let's call it personal growth - because they will be coming from people whose "tastes trajectory" (and not just "tastes snapshot") once have been similar to yours. **Now the question:** I wonder if something similar to the "option 2" already exists and, if it does, I wonder how it works. And if it doesn't exist, you're welcome to discuss how it should work! :) * * * _(asked this one atMachine Learning Beta but no answers yet, so reposting here)_


title: Cook's Distance

text: The formula of Cook's distance is $$D_i=\frac{(\hat Y-\hat Y(i))^{\prime}(\hat Y-\hat Y(i))}{p\times MSE}$$ where, $\hat Y$ is the prediction from the full regression model and $\hat Y$ is a $(n\times 1)$ matrix $\hat Y(i)$ is the prediction from a refitted regression model in which observation $i $ has been omitted. That is , $\hat Y(i)$ is a $((n-1)\times 1)$ matrix. ### So the order of $\hat Y$ and $\hat Y(i)$ is not same. Then _How can i subtract this two matrix since for subtracting we need equal order of two matrices._


title: What is the equation for random forest?

text: I need an equation for random forest so that I can score fresh data I receive every week, based on beta estimates I got after building model using this ensemble methodology. Every week I do not want to build random forest model again and then score it using following commands in Python. Because I have build 15 such models and I need to score them every week. So easiest way is to get the equation (betas I already have) and score.               # importing the ensemble libraries     from `numpy` import `genfromtxt`     from `sklearn.metrics` import `classification_report`          #importing the dataset          dataset = genfromtxt(open('~//win_5050_6oct.csv','r'))[1:]            target = [x[1] for x in dataset]     train = [x[2:] for x in dataset]     test = genfromtxt(open('~\\win_act_6oct.csv','r'))[1:]          val = [x[2:] for x in test]     y_act = [x[1] for x in test]     Gst_i=[x[0] for x in test]          #building random forest model          rf = RandomForestClassifier(n_estimators=100)     rf.fit(train, target)


title: What does initial level of GDP mean?

text: This might come as a very trivial thing and way below standards of this group but I am struggling to figure out what do the authors mean when they say that they have used initial level of income in growth regressions? assuming that the first obervation in the sample is 120, is it something like this:   Year GDP   1991 120   1992 120   1993 120 or something like this: Year GDP   1991 120   1992 .   1993 . or is it just the lagged value of levels? an example of this type of models is the following: Benhabib, J., & Spiegel,M.M. (1994). The role of human capital in economic development: Evidence from aggregate cross-country data. Journal of Monetary Economics, 34(2), 143–173. Thank you very much in advance.


title: Power analysis for binomial data when sample sizes are different

text: I'm running an A/B test where I want to compare two different layout of a button and want to see if one layout lead to a higher Click-Through Rate with a certain statistical power. I want to use a significance level of 0.05. Layout 1 has 25302 visualization and received 40 clicks. Layout 2 has 27924 visualization and received 50 clicks. The distribution is binomial, considering the nature of the event I'm studying. Every visualization can in fact lead to a click or not. I was trying to calculate it in R, using the `pwr.2p2n.test()` function included in the `pwr` library. The downside of this function is that it's asking the effect size as an input, which I don't know how to calculate.


title: How to rank and compare results of screening tests?

text: I have a question about how to do analysis of an experiment that has already been done, I hope you can help me with some advice! I will try to keep it as simple as possible, but will give some detail so you know what I'm talking about! What has been done is a "screening trial" to look at the activity of about 50 subjects (fungi) as antagonists (against a pest), the 50 individuals are members of groups (species), but some groups have many more members than others I have results of several types of screening tests for each of the 50 individuals, with reps of each. The screening tests look at different aspects, like growth rate, direct effects, and indirect effects. I can rank the isolates by their results in each screening test, and there looks like a lot of variability. I want to be able to report the findings of screens, for each screening test and also to see if some individuals are in top ranks in different screening tests (and also the opposite, if some are great at some tests but not at others). I think what I want is to know if the results of the tests correlate for each individual....? I am not sure how to say - this individual is the best - how can I tell if it is different than the next in the rank? If I list the top ten from each screening test, I would like to know that they are statistically different from those I excluded from the list. I would also like to compare them as groups, to be able to say, this species was the best, but with different numbers of representatives within the species, I dont think I could do this (please advise) This seems like it would be a common research experiment, for example, for testing drugs in medical experiments, so I am looking for examples of what others have done to present this type of result. I have seen a similar experiment to what I have to analyse but that had been done on a small scale, and the researchers used ANOVA to test differences among individuals and among groups, and some posthoc test to give each group little letters designating their means different than other groups. This seems to be unwieldy for 50 subjects, and I'm not sure about this.... I think I need some kind of mixed model regression to put all the test results in a model to test for correlation/covariance, but my understanding is weak! I have been learning R and would like to do this analysis using R. Can you give me advice/suggestions? I would appreciate any help in understanding and clarifying this problem and solutions!


title: How to plot Optics Clustering result in Matlab (reachability plot)

text: I modified the following script for Optics clustering ( http://chemometria.us.edu.pl/download/OPTICS.M ) in order to work with DTW distance instead than Euclidean's. I obtained the Order vector (order), Reachability distances (RD) vector and Core distances (CD) vector. As I read from Wikipedia (http://en.wikipedia.org/wiki/OPTICS_algorithm#Extracting_the_clusters) > Using a reachability-plot (a special kind of dendrogram), the hierarchical > structure of the clusters can be obtained easily. It is a 2D plot, with the > ordering of the points as processed by OPTICS on the x-axis and the > reachability distance on the y-axis. Since points belonging to a cluster > have a low reachability distance to their nearest neighbor, the clusters > show up as valleys in the reachability plot. The deeper the valley, the > denser the cluster. The image above illustrates this concept. In its upper > left area, a synthetic example data set is shown. The upper right part > visualizes the spanning tree produced by OPTICS, and the lower part shows > the reachability plot as computed by OPTICS. Colors in this plot are labels, > and not computed by the algorithm; but it is well visible how the valleys in > the plot correspond to the clusters in above data set. The yellow points in > this image are considered noise, and no valley is found in their > reachability plot. They will usually not be assigned to clusters except the > omnipresent "all data" cluster in a hierarchical result. I can represent the results on a scatter plot using the Order as X-axis, and RD as Y-axis. How can I obtain a hierarchical dendrogram from my data? How can I obtain a spanning tree from my data (and display it)?


title: In stochastic gradient descent, is there only one update to $\theta$ for each iteration?

text: I have read that the update equation for stochastic gradient descent is as shown below, for each iteration, k. Does one iteration correspond to one training example? So for each example is there only one update to $\theta$? $ \theta^{k+1} = \theta^k - \epsilon_k \frac{\partial L(\theta^k,z)}{\partial \theta^k} $ Update: Is it different for Online learning?


title: Where can I find the testing result of Netflix contest

text: I am using Netflix data as my dataset (although he competition is closed, I can still find the whole data set) But I could not find the result data. Does anyone have that testing result data ? Thanks


title: Pros and Cons of different types of regression?

text: I'm currently studying Linear Regression and am wondering if somebody would help me out with a high-level comparison of different methods. Just to list some: Lasso, Ridge, Elastic Net, Principal Component Regression, Least Angle Regression, Kernel Ridge Regression, Partial Least Squares, etc. I'd like to know are some data cases that would urge you to consider one over the other (relationship between sizes of p and n, multicolinearity, interpretability, etc.) Any help would be appreciated, I'm just trying to get a much better grasp on the material.


title: Bayesian model selection in PyMC

text: I'm trying to do model selection using PyMC (v2.2), but having difficulty assessing the models using various Information Criteria and/or Bayes Factor. My model is similar to a typical regression, with several parameters (~10) with priors modelled by uniform distributions, and a single observation modelled by a normal distribution with a deterministic mean and uniform standard deviation (through a precision deterministic). The mean (dynamic) response of the system is actually generated by an ordinary differential equation and typically yields around 1500 data points. I can get reasonably accurate results using either Adaptive or Non-Adaptive MCMC, with around 50K samples following 50K burn-in samples. Based on issues reported with BF for complex models, I started looking at the DIC values produced after performing the MCMC analysis, but for models consisting of various combinations of the true parameter set, they were quite large with very little difference between them. For example, -14623.9 and -14624.8. Are the DIC, and other similar criterion such as BPIC, normally so insensitive to different (sub-)models? Using some code from the sandbox, I also computed the sample likelihood for my models in an attempt to exploit Bayes Factor. However, the log-likelihoods produced were very large, ranging from -1e5 to 1e5. Regularization further decreased the minimum and exponentiating resulted in overflow (see function `weight`), so the array `exp(loglikes[m])` ended up comprising all zeros and one one! Why would my log-likelihoods, calculated by `logp` after a call to `draw_from_prior`, be so large? I'm fairly new to Bayesian estimation, so any help would be greatly appreciated!


title: Expectation under similar distribution to inverse Gaussian

text: I have the following integral that I wish to evaluate: $\int_0^\infty x\,\, p(x)dx$. Where, $p(x)\propto \exp\left(-ax+\frac{b}{x+k}\right)$ for $a,b\ge0$ Firstly what would the normalising constant be, and secondly what would the expectation be. In the event that a=0, this could be evaluated using the logarithmic integral function: http://en.wikipedia.org/wiki/Logarithmic_integral_function.


title: Understanding Differences of Mean between Two Groups

text: I have a variable `test score` that is measured for an organization consisting of 2050 people. The goal is to test if men outperform women (or vice versa). Unfortunately, only 50 of the 2050 people are women. I was told that in order to run inference when there is such a disproportion of one group within a sample that I should pursue some sort of matching strategy. First, is this true, that if I match on a few characteristics that I can run inference (e.g., ordinary least squares, linear probability models) and be less worried about the fact that there are so few females--since many males will drop out? Second, is there anything I can do to compare the groups before the match even as an anecdote? This is the data:                          Test Score    Test Score          N        Mean        Variance    Coef. Var.     Gender            2000      26.12         10.89         0.13         Male                  50      56.10         25.01         0.09        Female


title: Distribution of Product of Normal and Poisson?

text: Suppose that X is distributed Poisson with a known rate and Y is a normal distributed with a know mean and variance. My goal is to approximate the distribution Z where P(Z) = P(X) * P(Y), where Z is a non-negative integer. I could get a good approximation by sampling, but I'd really like to have a fast solution, ideally closed-form.


title: In finding the moment generating function why do we multiply by $e^{tx}$ for each pmf term?

text: The moment generating function that is associated with the discrete random variable $X$ and pmf $f(x)$ is defined as: $$M(t) = E\left[e^{tX}\right] = \sum_{x \in S} e^{tx} f(x).$$ Where does this $e^{tx}$ come from? This vaguely looks like an integrating factor from differential equations. Also, I find it strange that when we take the derivative of $M(t)$ we don't say "lets take the partial derivative with respect to t" but isn't this what we're doing, treating $x$ as a constant? Why do we use the derivative notation instead of the partial derivative operator?