title: Predicted Responses as Weighted Averages: Classical Regression Vs. Mixed Effect Models

text: Let's imagine that I have collected data in a longitudinal design that would lend itself to a mixed effects analysis. So suppose I perform a research study with human subjects and observe responses ($Y_{ij}$), where $i$ designates the subject and $j$ designates the trial, where the $j$th trial involves a possibly repeatable combination of covariates. In mixed effects models, the "predicted value" for the true response of the $i$th subject on the $j$th trial ($\widehat{Y_{ij}}$) is a weighted average of the response that was actually observed, $Y_{ij}$, and the population averaged mean-response profile, $X_i \widehat{\beta}$. How much each term is weighted depends, roughly speaking, upon the relatively "size" of the between-subject vs. within-subject variability. Is there an analogue for OLS linear regression? That is, is it similarly possible to express the "predicted value" for the response, $\widehat{Y_i}$ as a weighted average of the observed response, $Y_i$ and something else?


title: Ways to compare ordered binary datasets?

text: This question regards the best way to compare ordered binary data. My situation is as follows: I'm interested in evaluating how well a model conforms to human performance data on a set of validation problems. For a test set of 30 problems, I have a record of whether my model produced a correct response (denoted as a "1") or no response at all** (a "0") within a time limit. Because the model operates probabilistically, I have data for multiple runs of the model on the same test set (organized in a binary matrix of run # vs. problem #). I'd like to compare this accuracy data against the accuracy data produced by a group of human participants working on the same validation set. My question is: what test(s) can I use to measure the degree of coherence between the two datasets that takes into account the different questions in the validation set (i.e., that recognizes that data in the first column of the model accuracy matrix [corresponding to performance on the first question in the validation set] should only be compared against data in the first column of the human accuracy matrix, etc.)? Apologies if this has been answered elsewhere! ** Because the model is structured such that it never generates a response unless it is the correct one (i.e. it simply keeps running until it either (1) finds the correct answer or (2) runs out of time, whichever comes first), the response output doesn't include false positives or false negatives (so far as I understand those terms).


title: How to Code the Variance-Stabilized Bootstrap-t Method for Confidence Intervals

text: I'm interested in coding up the variance-stabilized bootstrap-$t$ approach for calculation of confidence intervals. However, I'm a little stuck and could use some help. I'm coding this in R and this is what I got so far:                   data = dataset              for(i in 1:B){                   boot = sample(nrow(data),replace=TRUE)              boot_data = data[boot,]              boot_stat [i]= f(boot_data)                   for(j in 1:M){                  boot_boot = sample(nrow(boot_data),replace=TRUE)                  boot_boot_data = boot_data[boot_boot,]                  boot_boot_stat = f(boot_data)              }                   se[i] = sum_over_M((boot_boot_stat-est_stat)^2)/(M-1)         }      Now I have B replicates of the se, and need to stabilize them as they're not independent of boot_stat.               temp = lowess(se~boot_stat)$y     inv = 1/temp     f = approxfun(boot_stat,inv,yleft=0,yright=0)     int = integrate(f,min(boot_stat),max(boot_stat))      Now I'm not sure how to proceed. I know I need to conduct another bootstrap, but I'm not sure what I'm bootstrapping. Something to do with the inverse transformation, but I'm still stuck. Any ideas?


title: Testing for non-random overlap of polygons

text: I have two polygon shapefiles and I want to see to what extent the observed amount of overlap is due to chance. I'm thinking of some kind of permutation test, but not sure of the best way to proceed. One idea is randomly relocating the polygons in one shapefile to new locations, measuring the total area of intersection, then repeating this a thousand times or so to get a distribution of areas of intersection. Then I can measure how likely the observed area is compared to the distribution of random areas. I have a recent GIS question relates to this approach, it was a comment on that question motivated this question. Another idea is to combine the two polgyon shapefiles into one shapefile, then randomly reassign the labels of the two types of polygon, measure the total area of intersection under the new labelling scheme, then repeat, etc. as above. The specific use-case is something like this: one polygon shapefile contains polygons of human activity areas (a few larger polygons). The other polygon shapefile contains polygons of rock outlines (many smaller polygons). The question is whether the observed arrangement of rocks is determined by the activity areas or not: did people put rocks on their activity areas or are we just looking at a natural distribution of rocks? Are either of the two ideas to test this on the right track, or should I be thinking of something completely different? Any insights would be most welcome!


title: Sum of Chi-squared Variables

text: What is the distribution of the variable $Y$ which is given by $$Y = \sum^N_{i=1} Z_i X_i $$ where : $Z_1,\dotsc,Z_N$ are independent and identically distributed and the density function $f_Z(z)$ is known (a shifted chi-square distribution), $X_i \sim\chi^2_L(\lambda)$, where $\lambda$ is the noncentrality paramter and $X_1,\dotsc,X_N$ are i.i.d. as well. The question is: (a) what is the distribution of $Y$ when $\lambda=0$. (b) what is the distribution of $Y$ when $\lambda\ne 0$.


title: Is it at all defensible to stratify a data set by the size of the residual and do a two-sample comparison?

text: This is something I'm seeing done as sort of an ad-hoc method and it seems very fishy to me but perhaps I am missing something. I've seen this done in multiple regression but let's just keep it simple: $$ y_{i} = \beta_{0} + \beta_{1} x_{i} + \varepsilon_{i} $$ Now take the residuals from the fitted model $$ e_{i} = y_{i} - \left( \hat{\beta}_{0} + \hat{\beta}_{1} x_{i} \right) $$ and stratify the sample based on the size of the residuals. For example, say the first sample is the bottom 90% of the residuals and the second sample is the top 10%, then proceed to do two sample comparisons - I've seen this done both on the predictor in the model, $x$, and on variables not in the model. The informal logic used is that perhaps points that have values far above what you'd expect under the model (i.e. a large residual) are different in some way, and that difference is investigated in this way. My thoughts on the subject are:   * If you see a 2-sample difference on a predictor in the model, then there are effects of the predictor not accounted for by the model in its current state (i.e. non-linear effects).    * If you see a 2-sample difference on a variable not in the model, then maybe it should've been in the model in the first place.  One thing I've found empirically (through simulations) is that, if you're comparing the mean of a predictor in the model $x$ and stratify in this way to produce the two sample means, $\overline{x}_{1}$ and $\overline{x}_{2}$, they are positively correlated with each other. This makes sense since both samples depend on the $\overline{y}, \overline{x}, \hat{\sigma}_{x}, \hat{\sigma}_{y}$ and $\hat{\rho}_{xy}$. That correlation increases as you move the cutoff down (i.e. the % you use to divide the sample). So at the very least, if you're going to do a two-sample comparison the standard error in the denominator of the $t$-statistic needs to be adjusted to account for the correlation (although I haven't derived an explicit formula for the covariance). Anyhow, my basic question is: Is there any rationale for doing this? If so, in what situations could this be a useful thing to do? Clearly I don't think there is but there may be something I'm not thinking of in the right way.


title: Lorenz curve and Gini coefficient for measuring classifier performance

text: I often use a ROC curve and the area under that curve as a measure of classifier accuracy in 2-class problems, e.g:               #Load a dataset     library(mlbench)     data(Sonar)          #Build a model     library(caret)     model <- train(Class~., data=Sonar, method='gbm', tuneLength=1, trControl=trainControl(method='cv'))     model          #ROC curve and AUC     library(pROC)     pMal <- predict(model, newdata=Sonar, type='prob')[,2]     roc(Sonar$Class, pMal, plot=TRUE)     >Area under the curve: 0.9705          #Lorez curve and gini?      ![enter image description here](http://i.stack.imgur.com/D6u2o.png) In a similar manner, I would like to be able to plot the lorenz curve and calculate the gini coefficient for my classifier. I know `Gini = 2*AUC-1`, but I'm not actually sure how to calculate it on it's own. Furthermore, every application of a lorenz curve I've seen looks at univariate data (e.g. income distribution). **How do I calculate a lorenz curve when I have 2 parameters: the predicted probability of the positive class, and the positive class itself?**


title: Pseudo-random orthogonal matrix generation

text: A little while back, J.M. suggested using the Stewart algorithm for generating $n$ by $n$ pseudo random orthogonal matrices in $O(n^2)$ time. He further noted that this methodology is implemented in Nick Higham's matrix computation toolbox (matlab). Now this package contains a bunch of .m files. Could anyone direct me to which one implements the pseudorandom orthogonal matrix generator? Alternatively, would anyone know of a R implementation of this algorithm?


title: What is the definition of "best" as used in the term "best fit" and cross validation?

text: If you fit a non linear function to a set of points (assuming there is only one ordinate for each abscissa) the result can either be:   1. a very complex function with small residuals   2. a very simple function with large residuals Cross validation is commonly used to find the "best" compromise between these two extremes. But what does "best" mean? Is it "most likely"? How would you even start to prove what the most likely solution is? My inner voice is telling me that CV is finding some sort of minimum energy solution. This makes me think of entropy, which I vaguely know occurs in both stats and physics. It seems to me that the "best" fit is generated by minimising the sum of functions of complexity and error ie               minimising m where m = c(Complexity) + e(Error)      Does this make any sense? What would the functions c and e be? Please can you explain using non mathematical language, because I will not understand much maths.


title: Clustering variables based on correlations between them

text: Questions:   1. I have a large correlation matrix. Instead of clustering individual correlations, I want to cluster variables based on their correlations to each other, ie if variable A and variable B have similar correlations to variables C to Z, then A and B should be part of the same cluster. A good real life example of this is different asset classes - intra asset-class correlations are higher than inter-asset class correlations.   2. I am also considering clustering variables in terms stregth relationship between them, eg when the correlation between variables A and B is close to 0, they act more or less independently. If suddenly some underlying conditions change and a strong correlation arises (positive or negative), we can think of these two variables as belonging to the same cluster. So instead of looking for positive correlation, one would look for relationship versus no relationship. I guess an analogy could be a cluster of positively and negatively charged particles. If the charge falls to 0, the particle drifts away from the cluster. However, both positive and negative charges attract particles to revelant clusters. I apologise if some of this isn't very clear. Please let me know, I will clarify specific details.


title: Testing homoscedasticity for multiple regression in SPSS

text: I am conducting a multiple regression with 1 DV and 6 IVs. I am trying to test Homoscedasticity on SPSS using a scatterplot since all my variables are scales. I conducted a the residual vs predictor value scatterplot and I think it might be a little heteroscadestic. How do I know which variable is the one causing the problem? And what should the next step be to try to make my data homoscedastic?


title: Normalizing or detrending groups of samples

text: How do I detrend or normalize multiple series of data so that I can inter- compare between the series? * * * Specifics below may not be appropriate for this forum. Please let me know and I can remove or re-phrase, but I think it might be helpful to fully understand the generic question above. I have a data set that I would like help analyzing. I think this question belongs here and not in at http://gis.stackexchange.com/ Specifically, I have the following situation: Each series is collected from an airplane flying a path, and has a variable number of (value, lat,lon,time) tuples. I have several of these flightpaths, each at a different time, and flying different paths (sometimes crossing, sometimes not). Flights happened months apart, at different times of day, and due to natural phenomena, the data (thermal in this case) varies. Part of the region flown over by multiple flights may or may not have an anomalous temperature signature. This is what I want to investigate. I am seeking an algorithm to detrend or normalize all of the flights so that I can increase my SNR and determine if there is a temperature anomaly over a sub- region.


title: Classification rate

text: I have used ordinal logistic regression where the response variable is NPS and the explanatory variables are different questions. I have found a model where all the questions in it are significant. Then I did a classification table to see how well the model can predict the observations to the right NPS group. The classification rate I got was 70% and I'm wondering if that is a good classification rate. What is considered a good and a bad classification rate?


title: Demonstrate difference in growth over time

text: **Short version** Is there a difference **per treatment** given time and this dataset? **Or** if the difference we're trying to demonstrate is important, what's the best method we have for teasing this out? **Long version** Ok, sorry if a bit _biology 101_ but this appears to be an edge case where the data and the model need to line up in the right way in order to draw some conclusions. Seems like a common issue... Would be nice to demonstrate an intuition rather than repeating this experiment with larger sample sizes. Let's say I have this graph, showing mean +- std. error: ![p1](http://i.stack.imgur.com/eIKeF.png) Now, it looks like there's a difference here. Can this be justified (avoiding Bayesian approaches)? The simpleminded man's approach would be to take Day 4 and apply a _t-test_ (as usual: 2-sided, unpaired, unequal variance), but this doesn't work in this case. It appears the variance is too high as we only had 3x measurements per time-point (err.. mostly my design, p = 0.22). **Edit** On reflection the next obvious approach would be ANOVA on a linear regression. Overlooked this on first draft. This also doesn't seem like the right approach as the usual linear model is impaired from heteroskedasticity ( _exaggerated variance over time_ ). **End Edit** I'm guessing there's a way to include **all** the data which would fit a simple (1-2 parameter) model of growth over time per predictor variable then compare these models using some formal test. This method should be justifiable yet accessible to a relatively unsophisticated audience. I have looked at `compareGrowthCurves` in statmod, read about grofit and tried a linear mixed-effects model adapted from this question on SE. This latter is closest to the bill, although in my case the measurements are not from the **same subject** over time so I'm not sure mixed-effects/multilevel models are appropriate. One sensible approach would be to model the rate of growth per time as linear and fixed and have the random effect be **Tx** then test it's significance, although I gather there's some debate about the merits of such an approach. (Also this method specifies a linear model which would not appear to be the best way to model a comparison of growth which in the case of one predictor has not yet hit an upper boundary and in the other appears basically static. I'm guessing there's a generalized mixed-effects model approach to this difficulty which would be more appropriate.) Now the code:               df1 <- data.frame(Day = rep(rep(0:4, each=3), 2),                   Tx = rep(c("Control", "BHB"), each=15),                   y = c(rep(16e3, 3),                   32e3, 56e3, 6e3,                   36e3, 14e3, 24e3,                   90e3, 22e3, 18e3,                   246e3, 38e3, 82e3,                   rep(16e3, 3),                   16e3, 34e3, 16e3,                   20e3, 20e3, 24e3,                   4e3, 12e3, 16e3,                   20e3, 5e3, 12e3))     ### standard error     stdErr <- function(x) sqrt(var(x)) / sqrt(length(x))     library(plyr)     ### summarise as mean and standard error to allow for plotting     df2 <- ddply(df1, c("Day", "Tx"), summarise,                  m1 = mean(y),                  se = stdErr(y) )     library(ggplot2)     ### plot with position dodge     pd <- position_dodge(.1)     ggplot(df2, aes(x=Day, y=m1, color=Tx)) +      geom_errorbar(aes(ymin=m1-se, ymax=m1+se), width=.1, position=pd) +      geom_line(position=pd) +      geom_point(position=pd, size=3) +      ylab("No. cells / ml")      Some formal tests:               ### t-test day 4     with(df1[df1$Day==4, ], t.test(y ~ Tx))     ### anova     anova(lm(y ~ Tx + Day, df1))     ### mixed effects model     library(nlme)     f1 <- lme(y ~ Day, random = ~1|Tx, data=df1[df1$Day!=0, ])     library(RLRsim)     exactRLRT(f1)      this last giving                   simulated finite sample distribution of RLRT.  (p-value based on 10000         simulated values)          data:       RLRT = 1.6722, p-value = 0.0465      By which I conclude that the probability of this data (or something more extreme), _given the null hypothesis that there is no influence of **treatment** on **change over time_** is close to the elusive 0.05. Again, sorry if this appears a bit basic but I feel a case like this could be used to illustrate the importance of modelling in avoiding further needless experimental repetition.


title: Combining evidence from different sources

text: I am computing similarities between words using similarity measures based on wordnet. I have obtained similarity measures for two words using three different methods - JNC (Jiang and Conrath), BNP (Banerjee and Pederson) and LIN (Lin). These similarity measures give values between 0 and 1 for two words (1 for exactly similar words). I need to combine the similarity measures from the above three methods, so that I can get a better measure of similarity. I want to combine the results because each of them has a drawback, and combining them would give a better measure. Any ideas on how I can do that? I read about Dempster-Shafer, but I want to use some other approach, because I'm finding it difficult to understand and apply it according to my requirement.


title: Probability Questions

text: I've got a 2 advanced probability questions that I'm having trouble with so just asking to confirm an answer/method   1. A computer company provides an insurance policy for one of its systems. If the system fails during the first year the policy pays £3000. The benefit decreases by £1000 each year until it reaches zero. If the system has not failed at the beginning of a year, the probability that it fails during the year is 0.1. How much should the company charge for the insurance policy so that, on average, its net gain per policy be £100?   2. A man aiming at a target receives 10 points if his shot is within 1cm of the center of the target, five points if it is between 1cm and 3cm, and three points if it is between 3cm and 5cm. What is the expected number of points scored if the man's shot is uniformly distributed in a circle of radius 8cm centers on the target?


title: Recommend classification algorithms to try

text: I am working on a binary classification problem that is reasonably-sized (100k observations). I extracted 60 numerical features; the classes in the training set are well balanced. There are some significant linear patterns, but after that the patterns seem very nonrandom and so I need classifier models that can deal with this. I am really looking to squeeze out the best possible (estimated) accuracy, at the sacrifice of computational effort, so I am considering creating an ensemble classifier. So far, I have received pretty good results with:   * a Random Forest classifier (90% CV accuracy)   * a radial basis SVM classifier (87% CV accuracy, still busy tuning it on a finer grid). I am now wondering if there are any other potentially interesting algorithms that I could add to the mix (three would be nice for majority voting, for example). I hope diverse models will help me shed some potentially remaining bias and improve accuracy a little bit. Preferably I'd use algorithms available through R's caret package. I am looking at Gaussian processes right now. My background in machine learning is not very theoretical; I really only have intimate experience with SVM, decision trees and random forests, so the list of algorithms in caret is pretty daunting and I am having a hard time finding applied studies that compare them. I know it's hard to predict relative performance on particular data, but I'm willing to burn through a few of them!


title: Time Series Similarity : Differing Lengths with R

text: I am experimenting with creating a distance matrix between time series for clustering and similarity searching. The main reference I am using is for the Similarity procedure in SAS (Paper). I would like to conduct the analysis in R using the `dtw` package. What I am confused about it the application of DTW to series of different lengths. 1) Is this possible? 2) Any hints on how to do it with the R package? Regarding question #2, trying to calculate a distance matrix on a matrix of series that differ in length immediately fails:               dist_mat<-dist(appliances_t,method="DTW")          Error in dtw(distance.only = TRUE, ...) :        No warping paths exists that is allowed by costraints      > dput(appliances_t)               structure(c(1L, 14L, 1L, 1L, 2L, 1L, 1L, 7L, 1L, 33L, 20L, 1L,      1L, 8L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 14L,      0L, 1L, 0L, 1L, 1L, 6L, 1L, 32L, 20L, 1L, 2L, 8L, 0L, 0L, 2L,      1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 19L, 0L, 3L, 6L, 1L, 1L,      7L, 1L, 42L, 27L, 1L, 3L, 10L, 0L, 1L, 3L, 1L, 3L, 1L, 1L, 1L,      0L, 0L, 1L, 1L, 22L, 1L, 7L, 4L, 1L, 5L, 7L, 1L, 51L, 32L, 5L,      4L, 12L, 1L, 1L, 9L, 5L, 7L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 33L,      1L, 6L, 4L, 3L, 5L, 5L, 1L, 80L, 49L, 5L, 5L, 19L, 1L, 1L, 9L,      5L, 7L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 28L, 1L, 7L, 8L, 3L, 5L,      6L, 3L, 63L, 41L, 5L, 6L, 15L, 1L, 1L, 9L, 5L, 8L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 30L, 1L, 10L, 22L, 10L, 8L, 7L, 5L, 70L, 44L,      8L, 7L, 16L, 1L, 1L, 12L, 8L, 11L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,      32L, 2L, 13L, 30L, 10L, 10L, 11L, 12L, 74L, 47L, 10L, 8L, 18L,      2L, 1L, 16L, 10L, 14L, 1L, 1L, 1L, 1L, 2L, 1L, 5L, 23L, 4L, 6L,      30L, 10L, 5L, 9L, 12L, 53L, 33L, 5L, 9L, 12L, 5L, 1L, 9L, 5L,      8L, 5L, 5L, 5L, 3L, 5L, 5L, 7L, 27L, 7L, 9L, 22L, 14L, 7L, 11L,      12L, 61L, 36L, 7L, 10L, 14L, 7L, 1L, 11L, 7L, 11L, 8L, 8L, 8L,      6L, 7L, 8L, 14L, 27L, 14L, 13L, 38L, 19L, 10L, 7L, 15L, 61L,      37L, 10L, 11L, 14L, 14L, 1L, 14L, 10L, 14L, 12L, 12L, 12L, 9L,      14L, 12L, 14L, 38L, 14L, 14L, 44L, 10L, 11L, 9L, 20L, 86L, 54L,      11L, 12L, 20L, 14L, 1L, 16L, 11L, 14L, 22L, 22L, 22L, 16L, 14L,      22L, 14L, 27L, 14L, 16L, 42L, 13L, 13L, 12L, 11L, 61L, 39L, 13L,      13L, 14L, 14L, 1L, 19L, 13L, 18L, 22L, 22L, 22L, 16L, 14L, 22L,      20L, 28L, 19L, 10L, 40L, 16L, 9L, 12L, 14L, 66L, 41L, 9L, 14L,      16L, 20L, 1L, 13L, 9L, 12L, 23L, 23L, 23L, 18L, 20L, 23L, 27L,      16L, 27L, 10L, 60L, 19L, 8L, 11L, 20L, 39L, 23L, 8L, 15L, 9L,      27L, 4L, 12L, 8L, 11L, 30L, 30L, 30L, 22L, 27L, 30L, 14L, 24L,      14L, 14L, 136L, 21L, 11L, 4L, 21L, 56L, 33L, 11L, 16L, 12L, 14L,      3L, 16L, 11L, 14L, 42L, 42L, 42L, 31L, 14L, 42L, 18L, 23L, 18L,      16L, 206L, 14L, 14L, 3L, 25L, 54L, 33L, 14L, 17L, 12L, 18L, 2L,      20L, 14L, 18L, 22L, 22L, 22L, 16L, 18L, 22L, 24L, 28L, 23L, 25L,      398L, 14L, 20L, 6L, 16L, 65L, 40L, 20L, 18L, 16L, 24L, 4L, 29L,      20L, 28L, 28L, 28L, 28L, 21L, 24L, 28L, 27L, 18L, 27L, 20L, 380L,      19L, 16L, 5L, 16L, 40L, 25L, 16L, 19L, 10L, 27L, 3L, 25L, 16L,      21L, 39L, 39L, 39L, 28L, 27L, 39L, 32L, 19L, 30L, 21L, 406L,      22L, 18L, 3L, 20L, 42L, 27L, 18L, 20L, 10L, 32L, 3L, 27L, 18L,      24L, 43L, 43L, 43L, 32L, 32L, 43L, 20L, 21L, 20L, 22L, 504L,      33L, 19L, 7L, 27L, 49L, 31L, 19L, 21L, 12L, 20L, 2L, 28L, 19L,      27L, 49L, 49L, 49L, 37L, 20L, 49L, 20L, 20L, 20L, 16L, 682L,      28L, 14L, 7L, 41L, 48L, 30L, 14L, 22L, 12L, 20L, 2L, 20L, 14L,      18L, 33L, 33L, 33L, 25L, 20L, 33L, 27L, 13L, 27L, 19L, 374L,      30L, 14L, 8L, 32L, 29L, 19L, 14L, 23L, 6L, 27L, 3L, 21L, 14L,      20L, 32L, 32L, 32L, 23L, 27L, 32L, 32L, 20L, 32L, 19L, 489L,      32L, 14L, 7L, 33L, 47L, 28L, 14L, 24L, 11L, 32L, 3L, 21L, 14L,      20L, 42L, 42L, 42L, 32L, 32L, 42L, 49L, 16L, 49L, 27L, 628L,      23L, 21L, 13L, 36L, 35L, 21L, 21L, 25L, 8L, 49L, 4L, 32L, 21L,      30L, 51L, 51L, 51L, 40L, 49L, 51L, 41L, 27L, 41L, 19L, 791L,      27L, 15L, 14L, 27L, 62L, 39L, 15L, 26L, 14L, 41L, 1L, 22L, 15L,      20L, 80L, 80L, 80L, 61L, 41L, 80L, 44L, 33L, 44L, 20L, 898L,      27L, 16L, 16L, 29L, 77L, 48L, 16L, 27L, 18L, 44L, 1L, 25L, 16L,      21L, 63L, 63L, 63L, 48L, 44L, 63L, 47L, 21L, 46L, 13L, 439L,      38L, 10L, 21L, 30L, 51L, 32L, 10L, 28L, 12L, 47L, 2L, 14L, 10L,      14L, 70L, 70L, 70L, 51L, 47L, 70L, 33L, 32L, 33L, 18L, 515L,      27L, 14L, 22L, 43L, 75L, 47L, 14L, 29L, 18L, 33L, 2L, 20L, 14L,      20L, 74L, 74L, 74L, 54L, 33L, 74L, 36L, 33L, 36L, 16L, 450L,      28L, 14L, 17L, 30L, 76L, 48L, 14L, 30L, 18L, 36L, 3L, 20L, 14L,      18L, 53L, 53L, 53L, 40L, 36L, 53L, 37L, 33L, 37L, 20L, 726L,      16L, 16L, 14L, 32L, 78L, 48L, 16L, 31L, 19L, 37L, 3L, 24L, 16L,      21L, 61L, 61L, 61L, 45L, 37L, 61L, 54L, 48L, 53L, 13L, 1069L,      24L, 10L, 17L, 19L, 109L, 68L, 10L, 32L, 27L, 54L, 3L, 14L, 10L,      14L, 61L, 61L, 61L, 45L, 54L, 61L, 39L, 51L, 37L, 13L, 889L,      23L, 11L, 19L, 28L, 120L, 76L, 11L, 33L, 28L, 39L, 3L, 16L, 11L,      14L, 86L, 86L, 86L, 63L, 39L, 86L, 41L, 43L, 40L, 16L, 1083L,      28L, 13L, 22L, 28L, 97L, 62L, 13L, 34L, 22L, 41L, 3L, 19L, 13L,      18L, 61L, 61L, 61L, 46L, 41L, 61L, 23L, 50L, 23L, 15L, 1110L,      18L, 13L, 25L, 32L, 118L, 73L, 12L, 35L, 28L, 23L, 2L, 19L, 13L,      17L, 66L, 66L, 66L, 49L, 23L, 66L, 33L, 80L, 33L, 10L, 803L,      19L, 7L, 32L, 20L, 185L, 116L, 7L, 36L, 43L, 33L, 7L, 11L, 7L,      11L, 39L, 39L, 39L, 28L, 33L, 39L, 33L, 63L, 33L, 14L, 828L,      21L, 12L, 33L, 20L, 146L, 91L, 12L, 37L, 34L, 33L, 9L, 18L, 12L,      16L, 56L, 56L, 56L, 42L, 33L, 56L, 40L, 68L, 40L, 11L, 1001L,      20L, 9L, 29L, 25L, 157L, 97L, 9L, 38L, 35L, 40L, 11L, 14L, 9L,      13L, 54L, 54L, 54L, 41L, 40L, 54L, 25L, 80L, 24L, 20L, 1281L,      13L, 16L, 32L, 23L, 181L, 113L, 16L, 39L, 43L, 25L, 3L, 22L,      16L, 21L, 65L, 65L, 65L, 48L, 25L, 65L, 27L, 68L, 27L, 24L, 918L,      20L, 20L, 32L, 14L, 155L, 96L, 20L, 40L, 35L, 27L, 1L, 29L, 20L,      27L, 40L, 40L, 40L, 29L, 27L, 40L, 31L, 85L, 30L, 16L, 974L,      16L, 14L, 45L, 23L, 194L, 120L, 13L, 41L, 47L, 31L, 4L, 20L,      14L, 18L, 42L, 42L, 42L, 31L, 31L, 42L, 30L, 76L, 29L, 23L, 524L,      27L, 19L, 51L, 18L, 172L, 109L, 19L, 42L, 40L, 30L, 7L, 28L,      19L, 27L, 49L, 49L, 49L, 36L, 30L, 49L, 19L, 74L, 19L, 23L, 486L,      33L, 20L, 40L, 30L, 168L, 105L, 19L, 43L, 39L, 19L, 3L, 28L,      20L, 27L, 48L, 48L, 48L, 36L, 19L, 48L, 28L, 80L, 28L, 24L, 380L,      21L, 20L, 51L, 39L, 185L, 115L, 20L, 44L, 43L, 28L, 3L, 29L,      20L, 28L, 29L, 29L, 29L, 21L, 28L, 29L, 21L, 131L, 21L, 33L,      456L, 32L, 27L, 54L, 27L, 300L, 188L, 27L, 45L, 70L, 21L, 3L,      40L, 27L, 38L, 47L, 47L, 47L, 34L, 21L, 47L, 39L, 139L, 38L,      37L, 360L, 33L, 30L, 63L, 37L, 315L, 198L, 30L, 46L, 74L, 39L,      3L, 47L, 30L, 41L, 35L, 35L, 35L, 27L, 39L, 35L, 48L, 197L, 48L,      30L, 488L, 33L, 24L, 90L, 38L, 449L, 280L, 24L, 47L, 104L, 48L,      8L, 36L, 24L, 33L, 62L, 62L, 62L, 47L, 48L, 62L, 32L, 175L, 32L,      36L, 622L, 48L, 29L, 86L, 39L, 401L, 249L, 29L, 48L, 95L, 32L,      3L, 44L, 29L, 40L, 77L, 77L, 77L, 58L, 32L, 77L, 47L, 147L, 47L,      58L, 522L, 51L, 47L, 63L, 53L, 336L, 211L, 47L, 49L, 79L, 47L,      3L, 69L, 47L, 64L, 51L, 51L, 51L, 40L, 47L, 51L, 48L, 212L, 48L,      47L, 306L, 43L, 36L, 76L, 60L, 483L, 304L, 35L, 50L, 114L, 48L,      3L, 54L, 36L, 50L, 75L, 75L, 75L, 56L, 48L, 75L, 48L, 175L, 48L,      48L, 628L, 50L, 39L, 106L, 48L, 401L, 249L, 39L, 51L, 95L, 48L,      0L, 58L, 39L, 52L, 76L, 76L, 76L, 57L, 48L, 76L, 68L, 186L, 67L,      56L, 320L, 80L, 47L, 85L, 59L, 426L, 266L, 47L, 52L, 100L, 68L,      1L, 68L, 47L, 63L, 78L, 78L, 78L, 59L, 68L, 78L, 76L, 144L, 75L,      48L, 646L, 63L, 39L, 69L, 92L, 332L, 207L, 39L, 53L, 77L, 76L,      2L, 58L, 39L, 52L, 109L, 109L, 109L, 80L, 76L, 109L, 62L, 140L,      61L, 61L, 790L, 68L, 48L, 83L, 72L, 323L, 200L, 48L, 54L, 76L,      62L, 3L, 72L, 48L, 66L, 120L, 120L, 120L, 91L, 62L, 120L, 73L,      170L, 72L, 54L, 556L, 80L, 43L, 96L, 78L, 390L, 244L, 43L, 55L,      91L, 73L, 6L, 63L, 43L, 60L, 97L, 97L, 97L, 73L, 73L, 97L, 116L,      157L, 114L, 51L, 552L, 68L, 43L, 113L, 91L, 357L, 224L, 43L,      56L, 84L, 116L, 7L, 63L, 43L, 58L, 118L, 118L, 118L, 88L, 116L,      118L, 91L, 225L, 91L, 58L, 806L, 85L, 47L, 96L, 78L, 512L, 321L,      47L, 57L, 119L, 91L, 5L, 68L, 47L, 65L, 185L, 185L, 185L, 138L,      91L, 185L, 97L, 103L, 96L, 95L, 855L, 76L, 75L, 108L, 96L, 236L,      148L, 75L, 58L, 56L, 97L, NA, 112L, 76L, 104L, 146L, 146L, 146L,      109L, 97L, 146L, 113L, 128L, 112L, 97L, 1050L, 74L, 78L, 116L,      86L, 294L, 184L, 78L, 59L, 69L, 113L, NA, 118L, 78L, 110L, 157L,      157L, 157L, 117L, 113L, 157L, 96L, 140L, 96L, 140L, 1120L, 80L,      112L, 123L, 83L, 321L, 199L, 112L, 60L, 75L, 96L, NA, 167L, 112L,      155L, 181L, 181L, 181L, 137L, 96L, 181L, 120L, 146L, 119L, 126L,      908L, 131L, 99L, 192L, 92L, 334L, 209L, 99L, 61L, 78L, 120L,      NA, 150L, 99L, 138L, 155L, 155L, 155L, 117L, 120L, 155L, 109L,      150L, 107L, 105L, 1068L, 139L, 84L, 373L, 150L, 343L, 214L, 83L,      62L, 80L, 109L, NA, 126L, 84L, 117L, 194L, 194L, 194L, 146L,      109L, 194L, 105L, 158L, 104L, 151L, 1696L, 197L, 120L, 195L,      158L, 363L, 227L, 120L, 63L, 84L, 105L, NA, 181L, 121L, 168L,      172L, 172L, 172L, 129L, 105L, 172L, 115L, 160L, 115L, 125L, 1658L,      175L, 99L, 213L, 224L, 366L, 228L, 99L, 64L, 86L, 115L, NA, 150L,      100L, 138L, 168L, 168L, 168L, 126L, 115L, 168L, 188L, 179L, 186L,      132L, 1872L, 147L, 106L, 215L, 199L, 408L, 256L, 106L, 65L, 96L,      188L, NA, 160L, 106L, 149L, 185L, 185L, 185L, 138L, 188L, 185L,      198L, 148L, 195L, 103L, 2262L, 212L, 81L, 232L, 167L, 338L, 212L,      81L, 66L, 80L, 198L, NA, 125L, 81L, 115L, 300L, 300L, 300L, 226L,      198L, 300L, 280L, 145L, 278L, 100L, 2120L, 175L, 80L, 318L, 242L,      332L, 209L, 80L, 67L, 138L, 280L, NA, 120L, 80L, 112L, 315L,      315L, 315L, 236L, 280L, 315L, 249L, 151L, 248L, 122L, 2215L,      186L, 97L, 184L, 199L, 348L, 226L, 97L, 68L, 137L, 249L, NA,      146L, 97L, 137L, 449L, 449L, 449L, 337L, 249L, 449L, 211L, 166L,      208L, 111L, 1756L, 144L, 91L, 195L, 213L, 402L, 434L, 91L, 69L,      190L, 211L, NA, 134L, 91L, 124L, 401L, 401L, 401L, 299L, 211L,      401L, 304L, 159L, 299L, 160L, 2010L, 140L, 128L, 243L, 165L,      422L, 380L, 128L, 70L, 149L, 304L, NA, 192L, 128L, 177L, 336L,      336L, 336L, 253L, 304L, 336L, 249L, 188L, 248L, 74L, 2278L, 170L,      60L, 242L, 160L, 492L, 501L, 59L, 71L, 127L, 249L, NA, 89L, 60L,      81L, 483L, 483L, 483L, 364L, 249L, 483L, 266L, 119L, 263L, 92L,      3064L, 157L, 74L, 264L, 195L, 403L, 359L, 74L, 72L, 128L, 266L,      NA, 110L, 74L, 102L, 401L, 401L, 401L, 300L, 266L, 401L, 207L,      237L, 206L, 99L, 3001L, 225L, 80L, 199L, 177L, 892L, 719L, 80L,      73L, 291L, 207L, NA, 119L, 80L, 112L, 426L, 426L, 426L, 319L,      207L, 426L, 200L, 383L, 199L, 104L, 4429L, 103L, 83L, 173L, 256L,      1418L, 710L, 82L, 74L, 342L, 200L, NA, 126L, 83L, 117L, 332L,      332L, 332L, 247L, 200L, 332L, 244L, 163L, 242L, 109L, 4118L,      128L, 86L, 207L, 118L, 1011L, 541L, 86L, 75L, 223L, 244L, NA,      128L, 86L, 118L, 323L, 323L, 323L, 241L, 244L, 323L, 224L, 463L,      221L, 113L, 3112L, 140L, 91L, 152L, 145L, 1268L, 436L, 91L, 76L,      276L, 224L, NA, 137L, 91L, 125L, 390L, 390L, 390L, 293L, 224L,      390L, 321L, 325L, 318L, 115L, 3968L, 146L, 92L, 173L, 160L, 1112L,      286L, 92L, 77L, 278L, 321L, NA, 138L, 92L, 125L, 357L, 357L,      357L, 268L, 321L, 357L, 148L, 579L, 145L, 127L, 3395L, 150L,      102L, 146L, 166L, 1606L, 382L, 102L, 78L, 253L, 148L, NA, 151L,      102L, 141L, 512L, 512L, 512L, 385L, 148L, 512L, 184L, 431L, 181L,      106L, 3693L, 158L, 84L, 158L, 171L, 1221L, 432L, 83L, 79L, 324L,      184L, NA, 126L, 85L, 118L, 236L, 236L, 236L, 178L, 184L, 236L,      199L, 833L, 198L, 103L, 3038L, 160L, 82L, 145L, 179L, 814L, 210L,      81L, 80L, 243L, 199L, NA, 125L, 82L, 115L, 294L, 294L, 294L,      220L, 199L, 294L, 209L, 580L, 207L, 109L, 2303L, 179L, 87L, 156L,      182L, 952L, 280L, 87L, 81L, 244L, 209L, NA, 131L, 88L, 120L,      321L, 321L, 321L, 240L, 209L, 321L, 1063L, 654L, 212L, 118L,      2959L, 148L, 95L, 162L, 202L, 920L, 332L, 95L, 82L, 202L, 214L,      NA, 142L, 95L, 130L, 334L, 334L, 334L, 249L, 214L, 334L, 1535L,      950L, 225L, 113L, 2509L, 306L, 91L, 297L, 167L, 907L, 393L, 91L,      83L, 201L, 227L, NA, 137L, 91L, 125L, 343L, 343L, 343L, 258L,      227L, 343L, 1705L, 1252L, 227L, 135L, 2718L, 397L, 109L, 332L,      166L, 804L, 320L, 108L, 84L, 121L, 228L, NA, 162L, 109L, 150L,      363L, 363L, 363L, 271L, 228L, 363L, 1203L, 1192L, 253L, 95L,      1813L, 481L, 76L, 451L, 174L, 869L, 308L, 75L, 85L, 91L, 256L,      NA, 113L, 76L, 104L, 366L, 366L, 366L, 275L, 256L, 366L, 1176L,      1127L, 210L, 140L, 2103L, 441L, 111L, 425L, 188L, 872L, 447L,      111L, 86L, 160L, 212L, NA, 167L, 111L, 154L, 408L, 408L, 408L,      306L, 212L, 408L, 943L, 1168L, 206L, 176L, 2156L, 600L, 141L,      346L, 179L, 980L, 478L, 141L, 87L, 169L, 209L, NA, 213L, 141L,      197L, 338L, 338L, 338L, 256L, 209L, 338L, 685L, 1469L, 215L,      129L, 2212L, 414L, 104L, 373L, 214L, 1195L, 764L, 104L, 88L,      234L, 217L, NA, 157L, 104L, 143L, 332L, 332L, 332L, 248L, 217L,      332L, 951L, 690L, 234L, 164L, 1995L, 358L, 130L, 522L, 150L,      849L, 441L, 129L, 89L, 145L, 236L, NA, 196L, 130L, 182L, 348L,      348L, 348L, 261L, 236L, 348L, 779L, NA, 227L, 147L, 1995L, 648L,      118L, 398L, 222L, NA, NA, 118L, 90L, NA, 228L, NA, 178L, 118L,      165L, 378L, 378L, 378L, 283L, 228L, 378L, 853L, NA, 266L, 167L,      2252L, 352L, 133L, 437L, 282L, NA, NA, 135L, 91L, NA, 270L, NA,      202L, 136L, 188L, 364L, 364L, 364L, 274L, 270L, 364L, 593L, NA,      186L, 213L, 2142L, 415L, 160L, 346L, 208L, NA, NA, 170L, 92L,      NA, 188L, NA, 257L, 170L, 236L, 430L, 430L, 430L, 324L, 188L,      430L, 844L, NA, 276L, 129L, 1933L, 584L, 135L, 423L, 261L, NA,      NA, 104L, 93L, NA, 280L, NA, 157L, 104L, 143L, 301L, 301L, 301L,      227L, 280L, 301L, 1055L, NA, 352L, 150L, 1973L, 606L, 177L, 507L,      236L, NA, NA, 119L, 94L, NA, 354L, NA, 180L, 119L, 167L, 447L,      447L, 447L, 334L, 354L, 447L, 773L, NA, 258L, 159L, 1656L, 1012L,      281L, 399L, 269L, NA, NA, 126L, 95L, NA, 261L, NA, 191L, 127L,      175L, 567L, 567L, 567L, 424L, 261L, 567L, 722L, NA, 324L, 167L,      1883L, 617L, 267L, 430L, 340L, NA, NA, 134L, 96L, NA, 326L, NA,      202L, 134L, 187L, 417L, 417L, 417L, 311L, 326L, 417L, 687L, NA,      293L, 169L, 1788L, 530L, 275L, 444L, 206L, NA, NA, 136L, 97L,      NA, 295L, NA, 205L, 136L, 189L, 521L, 521L, 521L, 392L, 295L,      521L, 711L, NA, 334L, 186L, 1763L, 597L, 234L, 460L, 240L, NA,      NA, 147L, 98L, NA, 338L, NA, 222L, 147L, 205L, 473L, 473L, 473L,      354L, 338L, 473L, 889L, NA, 422L, 198L, 1337L, 446L, 303L, 584L,      253L, NA, NA, 166L, 99L, NA, 427L, NA, 237L, 158L, 220L, 540L,      540L, 540L, 405L, 427L, 540L, 635L, NA, 258L, 310L, 1749L, 558L,      396L, 600L, 269L, NA, NA, 239L, 100L, NA, 261L, NA, 373L, 247L,      345L, 683L, 683L, 683L, 512L, 261L, 683L, 669L, NA, 298L, 359L,      1960L, 590L, 608L, 588L, 271L, NA, NA, 372L, 101L, NA, 303L,      NA, 431L, 288L, 398L, 416L, 416L, 416L, 311L, 303L, 416L, 712L,      NA, 315L, 214L, 1383L, 576L, 355L, 893L, 295L, NA, NA, 256L,      102L, NA, 319L, NA, 257L, 170L, 236L, 482L, 482L, 482L, 361L,      319L, 482L, 727L, NA, 334L, NA, 1631L, 730L, NA, 864L, 316L,      NA, NA, NA, 103L, NA, 336L, NA, NA, NA, NA, 509L, 509L, 509L,      382L, 336L, 509L, 887L, NA, 339L, NA, 1513L, 736L, NA, 580L,      496L, NA, NA, NA, 104L, NA, 340L, NA, NA, NA, NA, 539L, 539L,      539L, 403L, 340L, 539L, 794L, NA, 369L, NA, 1687L, 647L, NA,      632L, 572L, NA, NA, NA, 105L, NA, 264L, NA, NA, NA, NA, 546L,      546L, 546L, 408L, 325L, 546L, 838L, NA, 392L, NA, 2145L, NA,      NA, 712L, 340L, NA, NA, NA, 106L, NA, 357L, NA, NA, NA, NA, 594L,      594L, 382L, 443L, 462L, 594L, 880L, NA, 616L, NA, 1522L, NA,      NA, 700L, NA, NA, NA, NA, 107L, NA, 456L, NA, NA, NA, NA, 634L,      634L, 424L, 476L, 447L, 634L, 1030L, NA, 712L, NA, 1578L, NA,      NA, 703L, NA, NA, NA, NA, 108L, NA, 873L, NA, NA, NA, NA, 996L,      996L, 497L, 746L, 518L, 996L, 842L, NA, 256L, NA, 1546L, NA,      NA, 753L, NA, NA, NA, NA, 109L, NA, 337L, NA, NA, NA, NA, 1151L,      1151L, 515L, 862L, 465L, 1151L, NA, NA, NA, NA, 1685L, NA, NA,      699L, NA, NA, NA, NA, 110L, NA, NA, NA, NA, NA, NA, 685L, 685L,      310L, 512L, NA, 685L, NA, NA, NA, NA, 1726L, NA, NA, 857L, NA,      NA, NA, NA, 111L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, 1561L, NA, NA, 1099L, NA, NA, NA, NA, 112L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      1411L, NA, NA, 1156L, NA, NA, NA, NA, 113L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1714L, NA, NA, 3876L,      NA, NA, NA, NA, 114L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, 2258L, NA, NA, 1507L, NA, NA, NA, NA,      115L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, 1442L, NA, NA, 1524L, NA, NA, NA, NA, 116L, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      2073L, NA, NA, NA, NA, 117L, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2725L, NA, NA, NA,      NA, 118L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, 2255L, NA, NA, NA, NA, 119L, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, 1634L, NA, NA, NA, NA, 120L, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1360L, NA, NA,      NA, NA, 121L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, 1270L, NA, NA, NA, NA, 122L, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, 1328L, NA, NA, NA, NA, 123L, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1470L, NA,      NA, NA, NA, 124L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, 1040L, NA, NA, NA, NA, 125L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1095L, NA, NA, NA, NA, 126L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1823L,      NA, NA, NA, NA, 127L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 3221L, NA, NA, NA, NA, 128L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 3986L, NA, NA, NA, NA, 129L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3632L,      NA, NA, NA, NA, 130L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 3760L, NA, NA, NA, NA, 131L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 3225L, NA, NA, NA, NA, 132L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3473L,      NA, NA, NA, NA, 133L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 3359L, NA, NA, NA, NA, 134L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 2032L, NA, NA, NA, NA, 135L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 3337L,      NA, NA, NA, NA, 136L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 2474L, NA, NA, NA, NA, 137L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1748L, NA, NA, NA, NA, 138L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1581L,      NA, NA, NA, NA, 139L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 2590L, NA, NA, NA, NA, 140L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 2769L, NA, NA, NA, NA, 141L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2134L,      NA, NA, NA, NA, 142L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 2312L, NA, NA, NA, NA, 143L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1822L, NA, NA, NA, NA, 144L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1890L,      NA, NA, NA, NA, 145L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 1770L, NA, NA, NA, NA, 146L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1664L, NA, NA, NA, NA, 147L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1533L,      NA, NA, NA, NA, 148L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 1628L, NA, NA, NA, NA, 149L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1682L, NA, NA, NA, NA, 150L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1416L,      NA, NA, NA, NA, 151L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 1721L, NA, NA, NA, NA, 152L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1732L, NA, NA, NA, NA, 153L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1748L,      NA, NA, NA, NA, 154L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA, NA, NA, 1312L, NA, NA, NA, NA, 155L,      NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,      NA, NA, NA, 1659L, NA, NA, NA, NA, 156L, NA, NA, NA, NA, NA,      NA, NA, NA, NA, NA, NA, NA), .Dim = c(25L, 156L), .Dimnames = list(         c("units_1", "units_2", "units_3", "units_4", "units_5",          "units_6", "units_7", "units_8", "units_9", "units_10", "units_11",          "units_12", "units_13", "units_14", "units_15", "units_16",          "units_17", "units_18", "units_19", "units_20", "units_21",          "units_22", "units_23", "units_24", "units_25"), NULL))


title: Do I need heteroscedasticity consistent Standard Errors in LMER when applying Heckman's two stage procedure?

text: When including an Inverse Mill's ratio to account for selection in a mixed model using LMER in R (following Heckman's [1979] two-stage procedure), do I still need to estimate robust SEs or does the mixed model result in heteroscedasticity consistent SE?


title: How to handle Regression data thats not linear

text: I'm new to stats and am using Python 2.7 to fit a regression model (Random Forest). When I plot the percentile plot of the prices before and after a `log` transformation, I get **Percentile plot** _(Left: No transformation. Right: Log transformed)_ ![enter image description here](http://i.stack.imgur.com/qz8G2.png) I notice that the plot became more linear, except at the beginning where it's a little strange. What will you do under such circumstances, and what may have caused that abnormality? When I plot only the prices that have higher values (> 400000), the plot becomes ![enter image description here](http://i.stack.imgur.com/RccZy.png) Is it sufficiently linear now? What does a linear plot imply? How should I treat values under 400000 which I have excluded from being used to fit the regression model? How do you decide whether to log transform the variables too? By doing a similar percentile plot? * * * **Quantile-normal plot** _(Log transformation of y variable)_ _Left: Training Set. Right: Test Set_ ![enter image description here](http://i.stack.imgur.com/DyKt1.png) **Quantile-normal plot** _(_ **No** _transformation of y variable)_ _Left: Training Set. Right: Test Set_ ![enter image description here](http://i.stack.imgur.com/CxscG.png) Why does log transformation not do anything?


I use Lloyd's algorithm for clustering. Since it relies on a random initialization and Lloyd's algorithm can get stuck in local optima of the k-means objective function, I have to run it several times. How many different random initializations should I perform with Lloyd's algorithm to obtain the optimal clustering with X% of confidence? Or at least is there any rule of thumb advising for a decent number of random initializations?


title: How to determine if % selecting Option A is statistically different from selecting Option B in a survey?

text: I ran a survey on pricing and asked respondents to choose between option A or Option B. 61% chose option A and 39% chose Option B. The sample size is 90. How do I determine if the % selection Option A is statistically significant? What test do I run and how? The cross tab tables show stat tests across the rows (e.g., the % of Group 1 choosing option A is statistically different from group 2 choosing option A). However, I want to compare within the column (i.e., Is the % choosing Option A statistically different from the % choosing Option B at a confidence interval of 95% and 90%) I would like to do in Excel. Not sure how to use Binom.Dist formula. I have a feeling that would be the right test but I am not certain.


title: Time series forecasts of appointments with pre-registration

text: Looking for some tips and ideas. I get a list every day of the number of appointments for each day for the next two weeks for a clinic. I have quite good history of these list, and the actual number of visits that were seen each day. So for each day I can look at how many visits were actually seen, the total number of appointments scheduled for that day one day earlier, two days earlier, ..., back 14 days earlier. My desire is to develop forecasts for the actual number of appointments they can expect for the next 14 days. The simplest answer I can come up with is to develop seasonal time series models based on the actual visits and use the future appointments as a floor (e.g. if 10 appointments are scheduled for a 7 days from now, I should not predict that they will only have 8 appointments). But, this is not ideal, as it does not take into account that if there are 10 appointments 1 or more may not show. I would somehow like to incorporate the history data I have into the forecasts. Using the future appointments as factors does not seem like a good ideas as they are not independent. I have not found any examples of anyone performing forecasts with this type of data. Is there a way to potentially use this data in a better way? FYI: I am not necessarily bound to any technology, I can use either R or SAS.


title: Comparing linear regression models created from different data sets

text: I have one linear regression model [Mold] created from 12 points where I can calculate a single value of RMSE between the predicted values and the actual observed values. This model is then used to predict another dataset [Box, 42 points]. I have a second linear regression model [Box] in which all combinations of 12 points out of a 42 point sample are selected and RMSE between the predicted values and actual observed values is calculated for all combinations. The linear regression model in this scenario predicts the other 30 points not selected to construct the regression model. I am trying to figure out a way to compare the [Mold model] and the [Box model] with a statistical technique. My idea at the moment is to compare the single RMSE value from the [Mold Model] to a CDF of the RMSE from the multiple combinations of [Box Model]. Therefore I could say **XX%** of the time, the [Box Model] has an RMSE of **XX** better/worse than the [Mold Model]. If someone can please comment on the best technique you think I should be using it would be greatly appreciated. My master's thesis depends on it! Thanks!


title: Confidence intervals on a matrix of data in SAS

text: I have the following matrix of data, which I am reading into SAS:               1         5        12        19        13     6         3         1         3        14     2         7        12        19        21     22        24        21        29        18     17        15        22         9        18      It represents 5 different species of animal (the rows) in 5 different areas of an environment (the columns). I want to get a Shannon diversity index for the whole environment, so I sum the rows to get:               48        54        68        79        84      Then calculate the Shannon index from this, to get:               1.5873488      What I need to do, however, is calculate a confidence interval for this Shannon index. So I want to perform a nonparametric bootstrap on the initial matrix. Can anyone advise how this is possible in SAS?


title: Moderated regression for categorical, within-subjects variable

text: I am planning to run a moderated regression on my dataset with two IVs.   * One of them is a straightforward interval data-type, between subjects IV (e.g. anxiety level).    * The other is a categorical data-type, within subjects IV (e.g. type of social situation).    * DV is outcome scores on a variable (e.g. psychological adjustment).  I am predicting that anxiety level is a moderator of the effect of social situation on psychological adjustment. My question is:   * How can I run moderated regression since all participants will have the same scores on the categorical IV, since it is a within-subjects variable. ie. all participants will have 0 for social sitution A, and 1 for social situation B and so on?   * Am I better off running a two-way mixed design ANOVA?


title: How to reduce dimension for text document dataset?

text: I have a basic doubt in dimension reduction for text dataset eg. 20Newsgroup, rcv1 etc. Initially I extract the number of word occurrence in each document, i.e word x document matrix would be $n \times d$ where $n$ is the number of documents and $d$ is the dimension. I would like to reduce the dimension, say $d_1 << d$. What is the standard technique of reducing the dimension?   1. Choose the top $d_1$ feature from the original word occurrence matrix$( n \times d)$ and then calculate TF-IDF for the reduced matrix $(n \times d_1)$, or   2. Calculate the TF-IDF matrix for $n \times d$ matrix and then select top $d_1$ features.  Also, it is mentioned in many literature that top feature are selected. I wanted to know what is selecting the top features means? How do they define it?


title: How to calculate fitted values when using AdaBoost via gbm package in R?

text: First of all I'm not good at English. So please be patient to interpret my words. I really need your help. I used the R package GBM as probably my predictive modeling. I made a model with adaboost distribution. In these cases you would use `predict.gbm` code to predict some new data. But I have to predict them without R (maybe with SAS IML etc). So I want to know how to calcultate fitted values based on the model that I made. Data is as below (tab delimited). Variable 'sex' and 'houstype'(APT, DANDOK, ETC, VILLA) are categorical variables and y is a target variable that has 0,1 categories.               sex age income  area    housprice   houstype    y     M   30  72  31  99  APT 1     M   40  73  1   97  DANDOK  0     M   40  34  1   90  DANDOK  0     M   29  45  4   86  DANDOK  0     F   29  37  39  86  DANDOK  0     F   29  53  4   98  DANDOK  0     M   35  63  4   61  DANDOK  0     F   27  47  4   74  DANDOK  0     M   28  38  3   79  DANDOK  0     F   30  45  31  99  APT 0     M   31  32  4   82  DANDOK  0     M   30  47  64  98  APT 1     M   35  34  4   92  DANDOK  0     M   26  56  3   83  DANDOK  0     M   32  43  15  85  VILLA   1     M   24  32  4   69  DANDOK  0     M   31  35  3   68  DANDOK  0     M   40  47  8   52  DANDOK  0     F   26  40  56  14  DANDOK  0     F   29  45  61  63  VILLA   0     M   32  43  25  24  DANDOK  0     M   37  36  58  47  APT 1     M   36  49  49  53  APT 1     M   27  42  31  71  APT 1     M   33  53  7   48  DANDOK  0     M   29  60  8   77  VILLA   0     M   31  43  4   88  VILLA   0     M   28  35  3   83  DANDOK  0     M   38  53  12  58  VILLA   1     M   39  30  3   61  DANDOK  0     F   35  32  24  63  APT 0     F   26  36  31  97  APT 0     M   31  30  31  10  DANDOK  1     F   24  39  31  18  VILLA   0     M   38  96  4   23  APT 1     M   37  69  39  56  APT 1     M   39  37  58  17  APT 1     M   31  42  26  69  APT 1     M   32  64  64  47  APT 1     M   31  54  61  8   VILLA   1     F   29  35  18  88  APT 0     M   38  79  7   40  DANDOK  0     M   26  39  7   48  DANDOK  0     M   32  57  76  24  DANDOK  1     M   32  89  18  74  DANDOK  0     F   25  31  3   38  DANDOK  0     M   38  74  93  38  APT 1     F   34  37  76  80  APT 0     M   32  49  39  67  APT 1     M   33  60  18  69  APT 1     F   29  61  18  69  APT 0     M   37  30  83  73  APT 1     M   36  38  4   43  DANDOK  0     M   28  45  27  63  APT 1     M   39  89  39  76  APT 1     M   39  66  17  83  APT 1     M   29  67  1   61  DANDOK  0     M   30  41  1   69  DANDOK  0     M   32  44  1   83  DANDOK  0     M   32  57  49  41  VILLA   1     M   32  49  25  73  DANDOK  1     M   30  39  12  38  DANDOK  0     M   36  37  17  40  DANDOK  0     M   28  48  1   97  DANDOK  0     M   29  51  8   69  VILLA   0     M   29  50  8   69  VILLA   0     M   30  48  31  93  APT 1     M   40  86  39  98  APT 1     M   29  56  18  93  APT 1     F   29  79  8   92  APT 0     F   30  88  8   92  APT 0     M   29  34  2   69  DANDOK  0     M   34  41  25  54  DANDOK  0     M   27  34  2   92  DANDOK  0     F   34  33  49  39  VILLA   0     F   30  56  49  6   DANDOK  0     F   34  54  8   62  APT 1     F   23  37  31  38  VILLA   0     M   35  56  31  65  VILLA   1     F   26  42  39  60  APT 0     F   27  60  12  23  APT 0     F   29  52  8   37  APT 0     M   26  73  61  16  VILLA   1     F   33  45  49  39  APT 0     M   35  34  39  69  APT 1     F   25  30  8   55  APT 0     M   27  70  46  43  APT 1     M   30  43  18  81  APT 1     F   31  33  15  48  APT 0     M   32  80  24  84  APT 1     M   33  59  76  80  APT 1     M   23  32  27  32  VILLA   0     M   34  87  76  66  APT 1     F   30  42  31  46  VILLA   0     M   26  64  12  11  DANDOK  0     F   34  36  56  19  APT 0     M   35  50  4   7   DANDOK  0     M   32  75  4   7   DANDOK  0     F   29  58  69  72  APT 0     M   32  33  8   7   DANDOK  0     F   31  71  89  69  DANDOK  0     M   33  45  12  23  APT 1     M   26  56  4   76  ETC 0     M   38  48  3   44  DANDOK  0     M   37  62  39  69  APT 1     M   32  34  48  41  APT 1     M   32  54  89  20  APT 1     M   36  42  4   3   DANDOK  0     M   34  38  30  30  APT 1     M   37  51  18  38  DANDOK  0     F   27  32  24  23  APT 0     F   34  38  18  27  DANDOK  0     M   29  41  4   4   ETC 0     M   30  56  39  17  APT 1     M   32  68  8   77  APT 1     M   37  89  64  1   APT 1     M   31  39  3   57  DANDOK  0     M   38  36  83  73  VILLA   1     M   35  41  8   17  APT 1     M   33  30  49  1   DANDOK  1     M   32  32  4   4   DANDOK  0      Let the dataname as 'dataflc1'. And code is as below. Because I just want to know the method, I simplified the code.               dataflc1<-read.table("data.txt",header=TRUE)          as.character(dataflc1$sex)     as.character(dataflc1$houstype)          library(gbm)          dataflc1.output<-gbm(y~sex+age+income+area+housprice+houstype,     data=dataflc1, distribution="adaboost", n.trees=2,     keep.data=TRUE,     interaction.depth = 3)      Decision trees are as below.               pretty.gbm.tree(dataflc1.output,i.tree = 1)     pretty.gbm.tree(dataflc1.output,i.tree = 2)      which yields                 SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight    Prediction     0        3 10.0000000000        1         2           9     17.7855736     60 -0.0002277533     1       -1 -0.0008718615       -1        -1          -1      0.0000000     24 -0.0008718615     2        0  0.0000000000        3         7           8     29.1469266     36  0.0002016521     3        2 48.5000000000        4         5           6      0.9517364     23  0.0008808467     4       -1  0.0007259475       -1        -1          -1      0.0000000     10  0.0007259475     5       -1  0.0010000000       -1        -1          -1      0.0000000     13  0.0010000000     6       -1  0.0008808467       -1        -1          -1      0.0000000     23  0.0008808467     7       -1 -0.0010000000       -1        -1          -1      0.0000000     13 -0.0010000000     8       -1  0.0002016521       -1        -1          -1      0.0000000     36  0.0002016521     9       -1 -0.0002277533       -1        -1          -1      0.0000000     60 -0.0002277533            SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight    Prediction     0        5  1.0000000000        1         5           9      20.248252     60 -0.0001990321     1        0  2.0000000000        2         3           4      28.015752     30  0.0003333333     2       -1  0.0010000000       -1        -1          -1       0.000000     20  0.0010000000     3       -1 -0.0010000000       -1        -1          -1       0.000000     10 -0.0010000000     4       -1  0.0003333333       -1        -1          -1       0.000000     30  0.0003333333     5        3 28.0000000000        6         7           8       2.523166     30 -0.0007313975     6       -1 -0.0010000000       -1        -1          -1       0.000000     20 -0.0010000000     7       -1 -0.0001941924       -1        -1          -1       0.000000     10 -0.0001941924     8       -1 -0.0007313975       -1        -1          -1       0.000000     30 -0.0007313975     9       -1 -0.0001990321       -1        -1          -1       0.000000     60 -0.0001990321      Fitted values are as below.               dataflc1.output$fit      which yields                 [1] -0.2249587 -0.2288306 -0.2288306 -0.2288306 -0.2281529 -0.2288306 -0.2288306 -0.2288306 -0.2288306      [10] -0.2289587 -0.2288306 -0.2252328 -0.2288306 -0.2288306 -0.2272328 -0.2288306 -0.2288306 -0.2288306      [19] -0.2281529 -0.2281529 -0.2272328 -0.2252328 -0.2249587 -0.2252328 -0.2288306 -0.2288306 -0.2288306      [28] -0.2288306 -0.2269587 -0.2288306 -0.2289587 -0.2289587 -0.2264270 -0.2281529 -0.2268306 -0.2249587      [37] -0.2252328 -0.2252328 -0.2249587 -0.2261529 -0.2289587 -0.2288306 -0.2288306 -0.2261529 -0.2269587      [46] -0.2288306 -0.2249587 -0.2289587 -0.2249587 -0.2249587 -0.2289587 -0.2252328 -0.2288306 -0.2252328      [55] -0.2249587 -0.2249587 -0.2288306 -0.2288306 -0.2288306 -0.2261529 -0.2269587 -0.2272328 -0.2272328      [64] -0.2288306 -0.2288306 -0.2288306 -0.2252328 -0.2249587 -0.2249587 -0.2288306 -0.2288306 -0.2288306      [73] -0.2272328 -0.2288306 -0.2281529 -0.2281529 -0.2288306 -0.2281529 -0.2261529 -0.2289587 -0.2289587      [82] -0.2288306 -0.2261529 -0.2289587 -0.2252328 -0.2288306 -0.2249587 -0.2252328 -0.2289587 -0.2249587      [91] -0.2249587 -0.2272328 -0.2249587 -0.2281529 -0.2269587 -0.2289587 -0.2288306 -0.2288306 -0.2289587     [100] -0.2288306 -0.2281529 -0.2252328 -0.2288306 -0.2288306 -0.2249587 -0.2252328 -0.2249587 -0.2288306     [109] -0.2252328 -0.2269587 -0.2289587 -0.2289587 -0.2288306 -0.2249587 -0.2268306 -0.2249587 -0.2288306     [118] -0.2264270 -0.2268306 -0.2264270 -0.2288306      I have some questions. First, the SplitCodePred value of variable 'houstype' that has index 5 is 1.0000000000 in the second decision tree. Which category does 1.0000000000 belong to? APT? DANDOK? ETC? VILLA? How should I interpret the result? Second I want to know how to calculate 'ErrorReduction', 'Weight', 'Prediction' in each of the trees and fitted values based on the two trees. Third, What is 'MissingNode' in the results about decision trees? By the way, all variables don't have missing values. I sincerely hope your advices. Thank you.


title: Value of Information for a simple investment problem

text: Assume the following problem: You're deciding whether to invest into an opportunity with uncertain cost $c$ and value $v$. The cost has been estimated to be normally distributed with 90% CI between 1 and 5 million. Similarly, the value has been estimated to be normally distributed with 90% CI between 1 and 20 million. The default decision is to invest. However, you're interested in finding the value of reducing the uncertainty about both the cost and value in order to be certain about your decision. I intuitively understand that here we're after the Expected Value of Perfect Information or $EVPI$ which is equal to the potential loss in case we take the wrong decision. **Questions:**   1. How do we represent the problem above mathematically? Intuitively, I would write something like: $EVPI=\int_c\int_v(v-c)p(c)p(v)$   2. How do we solve the problem for a generic case? Is MCMC or a simple Monte Carlo a good choice here? I would prefer a layman-level answer as I'm not very strong with stats.


title: The meaning of probability when large number law fails

text: When the large number law fails, the frequency of an event does not converge to the probability, so what is the meaning of probability in these cases? PS: Please adds references to the answers (if it possible).


title: Converting monthly log data into annual growth rates

text: I have a time series of monthly log data pertaining to gross wages, going from Jan 2002 to Aug 2008. I need to use the Box-Jenkins modeling strategy to specify the correct ARIMA model, but for _annual wage growth rates_. How do I convert the data into this form? Thanks!


title: Multinomial Logit with mlogit and Yogurt Data

text: I am using the mlogit package in R to estimate a choice model, I started with the famous Yogurt dataset.               library(mlogit)          data(Yogurt)     names(Yogurt)          Yogdata <- mlogit.data(Yogurt, choice = "choice", shape = "wide", varying = c(2:9),                             sep=".")          Yog1 <- mlogit(choice ~ feat+price , data = Yogdata)          summary(Yog1)     phat      <- Yog1$probabilities     Yname     <-colnames(phat)     Y1hat     <- max.col(phat)     Y1hatname <- Yname[Y1hat]          preds = predict(Yog1,newdata=Yogdata)     preds[100:120,]     Yogurt[100:120,c('choice')]      The output `preds` includes probabilities that sum up to 1 for each row.                        dannon     hiland    weight   yoplait      [1,] 0.3881972 0.03053851 0.2200040 0.3612603      [2,] 0.6834572 0.01580043 0.1138286 0.1869137      [3,] 0.4136804 0.03137184 0.2344462 0.3205015      [4,] 0.4136804 0.03137184 0.2344462 0.3205015      [5,] 0.4136804 0.03137184 0.2344462 0.3205015      [6,] 0.4136804 0.03137184 0.2344462 0.3205015      [7,] 0.4136804 0.03137184 0.2344462 0.3205015      [8,] 0.4136804 0.03137184 0.2344462 0.3205015           [1] dannon  dannon  yoplait yoplait dannon  dannon  dannon  yoplait dannon      [10] dannon  yoplait yoplait dannon  dannon  dannon  dannon  dannon  dannon      [19] dannon  dannon  dannon       The choice would be the highest probability in each row. Would this be correct? When I compare this choice against real choices coming from Yogurt, for the same rows, I see matches are not too great. Is there a way to improve this model? Is this the right way to judge the success of this model?


title: Online resources for philosophy of causation for causal inference

text: Can you recommend any books, articles, essays, online tutorials/courses, etc that would be interesting and useful for an epidemiologist/biostatistician to learn about the philosophy of causation/causal inference? I know quite a bit about actually doing causal inference from an epi and biostats framework, but I would like to learn something about the philosophy which underlies and motivates this work. For example, it's my understanding that Hume first talked about ideas that could be interpreted as counterfactuals. I have basically no training or experience with philosophy, so I need something relatively introductory to start of with, but I would also be interested in recommendations for more complex but important/foundational texts/authors (but please indicate that they are not introductory). I hope this isn't too off-topic for cross-validated, but I'm hoping that some of you will have been in the same boat as me before and able to share your favorite resources.


title: Bayesian meta-analysis

text: A colleague is conducting a Bayesian meta-analysis of some large-scale drug studies. Although she has experience with 'conventional' meta-analysis, she in new at the Bayesian approach. We have read several papers with related studies and approaches, and some authors report mean drug concentrations, while others report medians. Is one more appropriate than the other in the Bayesian context?


title: Wrong prediction for data with weekly and annually cycles

text: I have one year daily data which has weekly and annually seasonality. There is a problem with seasonality definition. If my weekly seasonality has 7 cycles (one for each day) and annual seasonality has 52*7 cycles (52 weeks in a year each having 7 days), I will have 364 annual cycles but my data includes 366 days in the sample year. the prediction is not right because I'm including 364 seasonality but having 366 data. Hence predictions might be off by 2 days (Honestly I don't know if my explanation makes sense as I myself am new to forecasting and trying to understand how seasonality work).


title: Results with and without interaction

text: I'm working on an analysis with another person. First we did a logistic regression with study group and variable X. They were both significant. Then we added the interaction between study group and X only the study group and the interaction were significant. Actually, there are 3 groups and only one of the groups and one of the group*X interactions is significant. I am having a hard time to convince this person that we should not present the odds ratios for X from the first analysis that do not account for this interaction. This person thinks that we can show the results of the first analysis AND then the second analysis. I think this is misleading and will confuse people. I am not the best at explaining things (as you can tell), so any help with how to argue my point would be appreciated.


title: Empirical Bayes vs "non-informative" priors

text: I am familiar with the mechanics with both methods, but don't know what factors I should consider when choosing between these two approaches for adjusting a prior. I would imagine that, on a case by case basis, one could run posterior predictive checks, compare model complexity, etc. to see which model gives the best simplicity-fitting tradeoff. But in the general case, how does one make the decision of whether to use non- informative priors vs Empirical Bayes?


title: How to compare AR and ARIMA models?

text: Relatively new to stats. I use linear regression and get R^2, which is quite low. **MODEL 1**                   lmoutar=lm(formula = ts_y ~ ts_y_lag + ts_x)      So switched to arima with external regressor. Using "auto.arima", I formulate arimax model **MODEL 2**                   fitarima <- auto.arima(ts_y, xreg=ts_x)         arimaout<-arima(ts_y,order=c(2,0,5),xreg=ts_x)      How can I compare the explanability of AR model with arima model. From the thread How can I calculate the R-squared of a regression with arima errors using R?, I understand R^2 is not an option for ARIMA. From the thread Model comparison between an ARIMA model and a regression model, AIC/BIC is not the right criteria and MSE from forcast/predict can be possible criteria for comparison across AR and ARIMA model. Is MSE the best option for model comparison, if so how would I generate MSE for AR and ARIMA? I tried to compare the above ar and arima model using anova, but I get following error message               anova.lm(lmoutar,arimaout)        Warning message:         In anova.lmlist(object, ...) :                 models with response ‘"NULL"’ removed because response differs from model 1      What does this error message mean? **_EDIT_** Thanks for the response so far and insight that AR is nested within ARIMA. How would one answer this question, if I rephrase as "How to compare AR, ARIMA and General Linear Models?". The first model I listed has AR(1) and independent variable; it is a general linear model. So how would I compare a GLM versus ARIMAX model? Any thing else besides MSE that I could use to judge between GLM and ARIMAX


title: Posterior Predictive Checks

text: I understand what the posterior predictive distribution is, and I have been reading about posterior predictive checks, although it isn't clear to me what it does yet.   1. What exactly is the posterior predictive check?   2. Why do some authors say that running posterior predictive checks is "using the data twice" and should not be abused ? (or even that it is not Bayesian)? (e.g. see this or this)   3. What is this check exactly useful for? Can it really be used for model selection? (e.g. does it factor in both, fitness and model complexity?)


title: How is lasso an M-Estimator?

text: The definition of an M-estimator is an estimator (from Casella and Berger) of the form $$\hat{\theta}=\min \sum_{i=1}^n \rho(X_i-\theta),$$ where $X_1,X_2, \cdots, X_n$ is the data for some function $\rho$. The LASSO estimate is defined as $$\hat{\theta}=\min ||AX-\theta||_2^2 + \lambda ||\theta||_1,$$ for some matrix A, and a vector parameter $\theta$. Many papers refer to LASSO as an M-estimator, which does not appear to be true. Am I missing something here? The definition seems to imply there can not be any M-estimator for a vector parameter.


title: princomp loadings

text: My main goal is a model of the form `Y ~ D + N` . But lets say N and D correlate strongly so I decide to do a PCA on them I can do..               DN <- structure(list(D = c(0.7, 0.7, 0.8, 0.6, 0.5, 0.6, 0.6, 0.43,      0.43, 0.32, 0.6, 0.6), N = c(84, 62, 83, 67, 72, 69, 68, 59,      81, 52, 63, 83)), .Names = c("D", "N"), row.names = c(NA, -12L     ), class = "data.frame")          PCA <- princomp(DN, cor = TRUE)     summary(PCA)     #Importance of components:     #                         Comp.1    Comp.2     #Standard deviation     1.233169 0.6923106     #Proportion of Variance 0.760353 0.2396470     #Cumulative Proportion  0.760353 1.0000000       i feel happy with these results. And I can get my new vectors `Comp.1` and `~Comp.1` from `PCA$scores` which will allow me to do `Y ~ Comp.1 + Comp.2` . However if I want to know which comp explains mostly `D` or `N` I think I use `PCA$loadings`.               #Loadings:     #  Comp.1 Comp.2     #D -0.707  0.707     #N -0.707 -0.707          #               Comp.1 Comp.2     #SS loadings       1.0    1.0     #Proportion Var    0.5    0.5     #Cumulative Var    0.5    1.0      I dont understand why the loadings matrix are all equal...Is comp2 mostly D ? Questions:   1. Is this the correct approach for what I am trying to achieve?   2. How can I show which component is mostly explained by D and which by N?


title: How should I determine if a person differs from typical performance?

text: I have a set of data that looks like this:               Person 1 [48 total records]     2, 2, 2, 1, 2, 2, 1, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 2, 2, 1, 2, 1, 1, 1, 1,         1, 1, 2, 2, 1, 2, 1, 1, 2, 2, 2, 2, 1, 1, 2, 2, 1, 2, 1, 1, 1,     --     Person 2 [56 total records]     1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1, 1,         2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1, 1, 1,         1, 2, 1,     --     Person 3 [18 total records]     1, 2, 1, 2, 1, 1, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1,      A '1' indicates a false answer, and a '2' indicates a correct one. I want to be able to compare records to see if any people in the data set are performing significantly above or below average. I've heard of using z-scores and standard deviations, but I'm not sure if that's the correct approach, or even how I would go about doing the calculation. I also need to find out the minimum number of records I would need in order to have sufficient confidence in the results (if that's the correct term). Basically, I want to make sure I'm only performing this analysis on people with enough records to minimize the role chance would play. My math skills are pretty limited, so a simple explanation would be greatly appreciated.


title: How do you calculate the expected value of mixed lognormal distribution?

text: Suppose $X=\log(Y)$ can be modeled by a mixture of two normal distributions with proportion $p$ of $X_1$ and proportion $1-p$ of $X_2$, where $X_1\sim\mathcal N(U_1, \sigma^2_1)$ and $X_2\sim\mathcal N(U_2, \sigma^2_2)$. How do you calculate $E(Y)$; i.e., $E(\exp(X))$ where $X$ is a mixture of two normals?


title: Case-mix adjustment versus risk adjustment, what are their differences in practice and objective?

text: I have encountered in swathes of medical literature the use of the terms "case-mix" and "risk" adjustment without any citations or explanations of their exact usage and motivation from a modeling perspective. I understand the principles of covariate adjustment in multiple regression modeling to address confounding (bias) and stratification (efficiency). However, I can't seem to find a reference with discussion on the definition of these terms, their impact on analyses, and the objectives of their use. Can someone explain generally how and why case-mix and risk adjustment would be used in practice?


title: What is so cool about de Finetti's representation theorem?

text: From Theory of Statistics by Mark J. Schervish (page 12): > Although DeFinetti's representation theorem 1.49 is central to motivating > parametric models, it is not actually used in their implementation. How is the theorem central to parametric models?


title: Desirability approach to multiple optimization: drawbacks and alternatives

text: When fitting some relations $y_j=f_j(x_1, \ldots, x_p)+\epsilon_j$ between a multivariate response $y=(y_1, ..., y_k)$ and covariates $x_i$, a common problem is to determine for which values of the covariates $x_i$ the responses $y_j$, or the model functions $f_j$, attains an optimum. This problem is classical in the literature about "response surface" models. Of course it is not possible in general to find a simultaneous optimum, that is, a solution $(x_1, \ldots, x_p)$ for which _every_ function $f_j$ attains its optimum. Based on experimental data the desirability approach to this problem consists in:   * define some kind of "average" response $z=d_1(y_1)\cdots d_k(y_k)$, fit an univariate model $z=g(x_1, \ldots, x_p)+\epsilon$ and optimize the univariate fitted function $\hat g$;   * or firstly fit $p$ univariate regression models $y_j=\hat f_j(x_1, \ldots, x_p) + \epsilon_j$, $j=1,\ldots,p$, and then seek values of the covariates $x_i$ for which the functions $\hat f_j$ attain some kind of "simultaneous" optimum by defining a univariate "average" fitted model $\prod d_j(\hat f_j)$ and finally find its optimum. I am not attracted by this approach for the following reasons:   * The choice of the "desirability functions" $d_j$ seems too delicate and arbitrary to me. Even with strong prior knowledge and clear objectives (a rare practical situation). I would not be suprised that many users try numerous choices _a posteriori_ and finally choose the one that suits them best (hence their choice is data-dependent)   * I don't know whether there exists an appealing way to quantify the uncertainty about the solution provided by the desirability approach (even when assuming the desirability functions $d_j$ are the good ones I don't know whether it is possible to obtain some kind of confidence interval about the solution)   * This approach doesn't take into account the correlation between the univariate responses $y_j$ I would be interested in discussing these points and any other comments. Moreover I am interested in the alternatives to the desirability approach. Peterson and other authors (see for instance this powerpoint presentation and references given therein) proposed a more appealing approach: they fit a Bayesian multivariate linear model and then they optimize the posterior predictive probability that the multivariate response $y$ belongs to some prespecified set $A$. There's no need to be Bayesian for this approach: this probability can also be considered as an estimate of $\Pr(y \in A)$ and we can also provide a confidence interval about the optimized probability $\Pr(y \in A)$. The user must choose $A$ but with clear objectives this choice is not as arbitrary as the choice of the desirability functions $d_j$ (for instance we want a temperature to be higher than $20$, a concentration between $2$ and $4$...), and it well defined _a priori_. Do you know other alternatives to this multiple optimization problem ? **EDIT 06/09/2012:** The introduction of this paper by Peterson & al is an excellent discussion about the present topic.


title: What is the maximum likelihood estimate of the covariance of bivariate normal data when mean and variance are known?

text: Suppose we have a random sample from a bivariate normal distribution which has zeroes as means and ones as variances, so the only unknown parameter is the covariance. What is the MLE of the covariance? I know it should be something like $\frac{1}{n} \sum_{j=1}^{n}x_j y_j$ but how do we know this?


title: R obtaining a VAR forecast when some variable outcomes are already known

text: I have four variables and would like to construct a VAR model I would then like to make a VAR forecast on one of the variables using my own data for the forecasts of the other three variables. **Is there an R package or process from the model fit that allows me to do this?**


title: What is the estimation bias of the top estimate in a list sorted by value?

text: Let's make the problem as simple as possible. Assume two related random variables, $X_1$ and $X_2$. On the basis of some data we estimate their true means $\mu_{X_1}$ and $\mu_{X_2}$ by sample means $\hat\mu_{X_1}$ and $\hat\mu_{X_2}$. These estimates are unbiased. But now let's sort our two random variables by their sample means and look at the variable with the highest sample mean. Now for this top-of-the-list random variable the sample mean is now a **biased** estimator of its true mean (under some reasonable assumptions, e.g. that the means of these random variables are themselves distributed in a certain way and that distribution has a mean) -- that's easy to verify by Monte-Carlo. For obviousness, take not two but a thousand random variables and make their true means similar. The question is what is this bias and how do I analytically calculate it? I'd also appreciate some conceptual discussion on how does estimation bias arise out of sorting by estimated values.


title: Calculating standard error for funnel plots

text: I'm a programmer, with a decent but not-expert knowledge of stats, and I'm working through these instructions for how to create funnel plots. I know that there are better ready-made tools (e.g. R) for generating funnel plots, but these instructions are helping to understand the basics of funnel plots. So, my question: **Step 4 inthe instructions refers to a 'special kind' of standard error. Why is it special?** And, can I validly copy the way standard error is calculated in these instructions to build a funnel plot for any dataset that includes raw population count and raw incident count?


title: Is it valid to use a difference score as an independent variable in a regression analysis

text: I would like to see if the difference in the number of BPD symptoms from baseline to follow-up two years later can predict psychosocial functioning at the 2 year follow-up. So I wanted to do a linear regression using a difference score (BPD score T1- BPD score T2) as the IV and different measures of functioning at T2 as the DV's (e.g. SOFAS score as global functioning; peer rating scale as interpersnal functioning; and some binary measures such as work or no work at T2; so will need to complete a series of regressions). Is it valid to use the difference score. I have read a lot about using a difference score as a dependent variable but am unsure if the same information would correspond to using a difference score as a IV? I guess an additional question would be is it more approriate to use T1 measures of the outcome variables as covariates in the regression equation in order to measure change in that variable also? Or is that a seperate question? I appreciate any comments as I am quiet unsure of how to proceed and am going around in circles in how best to get some clarity.


title: RandomForest factor with too many levels

text: I have a model with about 200,000 training observations, where I am regressing, with 4 factors and 2 continuous variables. One of my features has 927 levels, which is causing the R implementation of randomForest to fail (it has a limit of 32 levels for any feature). Unfortunately, I don't see a simple way to avoid using this factor, or to decompose it into a series of continuous variables. Since my predictors are a mix of categorical and continuous, I thought of trees. Can anyone suggest a different implementation (package or language), ML approach, or a better way to pre-process or massage my inputs?


title: When should dimensional-reduction be used?

text: Yesterday I asked this question in which I had 180 subjects with 500 features each. While I was sure that dimensional-reduction is a must in this case (500 features), most of the answers I got said that 500 are not too many. **So, My question is:** Is there any rule of thumb when one should use dimensional-reduction before the classifier? How many features is too many? (I guess it is depends upon the ratio between the number of subjects and features. Isn't it?)


title: Understanding factor potentials in PyMC

text: I'm trying to understand factor potentials from the PyMC documentation, but need some help on the implementation piece--or it may turn out that I am misunderstanding how potentials work altogether. Imagine that we are building a Poisson switchpoint model as specified in PyMC's documentation tutorial. Now we want to introduce a factor potential exactly as described in the documentation, such that the difference between early and late Poisson means is less than 1. When I tried to implement the potential as in the code below, I see that values outside the allowed constraint appear in the posterior distribution. Why does the posterior distribution of potentialCheck contain these values? Obviously, I am doing something wrong...               from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform, potential     import numpy as np          disasters_array =   \      np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,                    3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,                    2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,                    1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,                    0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,                    3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,                    0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])          switchpoint = DiscreteUniform('switchpoint', lower=0, upper=110, doc='Switchpoint[year]')          early_mean = Exponential('early_mean', beta=1.)     late_mean = Exponential('late_mean', beta=1.)          @deterministic(plot=False)     def rate(s=switchpoint, e=early_mean, l=late_mean):         ''' Concatenate Poisson means '''         out = empty(len(disasters_array))         out[:s] = e         out[s:] = l         return out          @potential(plot=True)     def examplePotential(em=early_mean, lm=late_mean ):         if abs(em-lm) < 1:              return e         else:              return 1.          @deterministic(plot=True)     def potentialCheck(e=early_mean, l=late_mean):         '''Replicate the Potential to plot and check if the constraint has held.'''         return abs(e-l)          disasters = Poisson('disasters', mu=rate, value=disasters_array, observed=True)          disaster_model = [switchpoint, early_mean, late_mean, examplePotential, potentialCheck]          from pymc import MCMC     M = MCMC(disaster_model)     M.sample(iter=10000, burn=1000, thin=10)          from pymc.Matplot import plot as mcplot          mcplot(M)


title: Functional principal component analysis (FPCA). what is all about?

text: FPCA is something I have stumbled upon and never got to understand what it all about ? please see this paper, and I'm citing: > PCA runs into serious diﬃculties in analyzing functional data because of the > “curse of dimensionality” (Bellman 1961). The “curse of dimensionality” > originates from data sparsity in high-dimensional space. Even if the > geometric properties of PCA remain valid, and even if numerical techniques > deliver stable results, the sample covariance matrix is sometimes a poor > estimate of the population covariance matrix. To overcome this diﬃculty, > FPCA provides a much more informative way of examining the sample covariance > structure than PCA, I just don't get it: The drawback which is this paper mentioned above: wasn't PCA is ultimate method for handled the situation like the “curse of dimensionality” etc' ?


title: How to verify that the model is real?

text: Here is a block diagram which I'm using when I want to verify that my model is real. ![enter image description here](http://i.stack.imgur.com/22LrK.jpg) In each round a fold of 11/12 percent of the data is used to bulid the model (e.g. eigenvectors of the PCA) After 12 rounds I check that the models (e.g the eigenvectors of the PCA) from each round **are not** statistically different and if they indeed don't I declared the model as stable. Where the idea is that if a model is stable in a k-fold sense it is some sort of indication that the model is real I more or less thought about it myself (mainly on the last verification step) thus, I would like to know what do you think about it ? are you aware of other ways to "verify that my model is real" ?


title: Neural networks for simplistic image classification

text: I want to train a neural network to classify a few simple, cartoony images like the ones below (for the moment I only have the classes _house_ , _tree_ , and _sword_ ). ![house](http://i.stack.imgur.com/F9Pra.png) ![tree](http://i.stack.imgur.com/Tc7hD.jpg) ![sword](http://i.stack.imgur.com/Jt8QL.png) The images I am (currently) using are downsized to 32x32 pixels, and the feed- forward network architecture I use is 1024-512-256-3. This means that I end up with a total number of weights (excluding biases) of               1024*512 + 512*256 + 256*3 = 656128      That is a huge number and some function optimization algorithms end up depleting the available memory because of it. Obviously I'm doing something wrong. Should I use a different architecture or a different type of neural network (not MPL)? Should I reduce the image sizes further? When people train neural nets for complex object recognition tasks, how to they avoid using up all the memory?


title: How can I discern whether a classifier's outcome is significantly different?

text: This is a question I have posted here some while ago. and I would like to know if you think of more solutions for it from the ML perspective? Unfortunately, I can't use McNemar's test due to the small sample size (hence the values are not normally distributed) * * * How can I compare whether a classifier's output is significantly different? I have a relatively small data set which I'm trying to classify. My training set consists of 24 items which come from two different groups, 12 items from each group. Each item has two properties. I have 12 new items which I would like to classify based on this training set. (not sure whether the numbers are so important to my question... ) I have tested 5 different classifiers, and their classification result is as follows:               1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0     0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1     0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1     0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1     0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1      How can I test whether those results are significantly different from each other?


title: Forming an unbiased estimator of the maximum of several parameters, given independent estimators of each parameter?

text: Say I have $K$ independent normals, $X_i \sim \mathcal{N}(\mu_i, \sigma_i)$ for $i = 1,...,K$. How can I form an unbiased estimator of $\max_i \mu_i$ using $X_i$'s?


title: Higher order tensors in Statistics

text: I am looking for problems in Statistics or application areas in Statistics which deal with three-way data or higher order data or 3-Tensors. Is there a regression setting where the covariates are expressed as a higher order tensor tensor instead of a matrix (2-Tensor)?


title: Tail bounds on Euclidean norm for uniform distribution on $\{-n,-(n-1),...,n-1,n\}^d$

text: What are known upper bounds on how often the Euclidean norm of a uniformly chosen element of $\:\\{-n,~-(n-1),~...,~n-1,~n\\}^d\:$ will be larger than a given threshold? I'm mainly interested in bounds that converge exponentially to zero when $n$ is much less than $d$.


title: Regressing price on volume

text: What are the basic pitfalls of regressing stock price on volume of shares traded? This is a time series dataset. The model I'm using is: $$\ln(price)=b_0+b_1t+b_2\ln(volume)+\epsilon$$ As you can see I include a time trend $t$. The R2 and p-value appears to come out good, but I'm more accustomed to cross- sectional analysis so don't really know what I'm doing. I'd appreciate it if someone with more time series knowledge could let me know if this is totally naive. Thanks


title: T-test shows no differences, but the experiment group shows tendency more benefit in all variables measured than control group

text: I've just finished an animal experiment. I compared 1 control group and 1 experimental group, the only difference between the two is type of diet. For statistical analysis I used the independent groups t-test, and the result showed no significant differences between the two groups. However, the data shows the tendency that the experimental group has more benefit in all variables measured. So, what should I say about my data? All data are normally distributed. My supervisor said that maybe because I used very small sample (each group n=8) that I could not find any significant differences. He suggested me to do some "probability test" or something to extrapolate my data (unfortunately, I don't have any clue what he was talking about). So, is there any statistical analysis that I can use like what my supervisor told me to do?


title: Fixed effects models vs. models where years serve as instruments

text: I've been looking for a good explanation of the difference between running a regression with year dummies as fixed effects vs. year dummies as instrumental variables. A prominent example of this comes from Josh Angrist's 1991 paper, "Grouped- data estimation and testing in simple labor-supply models" http://www.sciencedirect.com/science/article/pii/030440769190101I


title: If you have a really big N, is using a regular vs. an m-n bootstrap ever a good idea?

text: I know that the regular bootstrap (assuming that the necessary assumptions are valid) is more efficient than the m-n bootstrap. However, supposing that you have access to really large data (such as a 1% sample of Google's daily searches), does it ever make sense to use the regular bootstrap whose assumptions may be violated?


title: In structural econometrics, what is meant by "simulation variance"

text: When people speak of "simulation variance", what does this mean? Does simulation variance disappear as $N \rightarrow \infty $? Are models that have less simulation variance for a fixed N considered more "robust"?


title: Example of estimation vs. calibration

text: What's a good example to demonstrate the difference between estimation and calibration? Edit: I'm looking for something that has the same definition of "calibration" and "estimation" that Krugman uses here (i.e. that they are in some way substitutes): > “calibration” — which basically means tweaking the parameters of your model > until it fits some aspects of the data, rather than flat-out estimating the > model


title: Derivation of Regularized Linear Regression Cost Function per Coursera Machine Learning Course

text: I took Andrew Ng's course "Machine Learning" via Coursera a few months back, not paying attention to most of the math/derivations and instead focusing on implementation and practicality. Since then I have started going back to study some of the underlying theory, and have revisited some of Prof. Ng's lectures. I was reading through his lecture on "Regularized Linear Regression", and saw that he gave the following cost function: $$J(\theta) = \frac{1}{2m}[\sum_{i=1}^m(h_\theta (x^{(i)}) - y^{(i)})^2 + \lambda\sum_{j=1}^n\theta^2_j]$$ Then, he gives the following gradient for this cost function: $$\frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m}[\sum_{i=1}^m(h_\theta (x^{(i)}) - y^{(i)})x^{(i)}_j - \lambda\theta_j]$$ I am a little confused about how he gets from one to the other. When I tried to do my own derivation, I had the following result: $$\frac{\partial}{\partial \theta_j}J(\theta) = \frac{1}{m}[\sum_{i=1}^m(h_\theta (x^{(i)}) + y^{(i)})x^{(i)}_j + \lambda\theta_j]$$ The difference is the 'plus' sign between the original cost function and the regularization parameter in Prof. Ng's formula changing into a 'minus' sign in his gradient function, whereas that is not happening in my result. Intuitively I understand why it's negative: we are reducing the theta parameter by the gradient figure, and we want the regularization parameter to reduce the amount that we are changing the parameter to avoid overfitting. I am just a little stuck on the calculus that backs this intuition. FYI, you can find the deck here, on slides 15 and 16.


title: Number of observations in groups - linear mixed effects model

text: I would like to fit linear mixed effects model to my dataset, but I was wondering if quantity of observations in groups matter? I have some groups with about 60 observations in each, but there are also some with only 1, so Im curious if there will be any influence on my linear mixed effects model because of that. Thanks in advance