title: Intercepts (reference) in linear mixed effect model, what to choose?

text: I am working on a study that examines word processing with native and non- native speakers. We have three independent variables: Groups (NSs and NNS), word types (five conditions), and word relatedness (related and unrelated words), and one dependent variable; Reaction Time. I have successfully ran linear mixed effects model in R using (lme4) package and I was able to understand the output. However, there is something that I do not understand, the intercepts (reference level). Could someone explain what does the intercept or reference means? Does it mean all levels of the variables are compared to it? If so, should I adjust the intercept level or let lme4 uses the default? I am struggling to make this decision and I need reasons to back up any choice I make. Note that that the Word Types are not related to each other, but they are compared with regard to Word Relatedness. Interestingly, when I change the intercept level the whole output changes. If you see below the Related:WordType also changes when the intercept is changed -- shouldn't this be the same in the two outputs?               > contrasts(NS$WordType)     # default intercept by R      2 3 4 5      1 0 0 0 0      2 1 0 0 0      3 0 1 0 0      4 0 0 1 0      5 0 0 0 1          Linear mixed model fit by REML ['lmerMod']     Formula: RT ~ Related * WordType * (1 | Item) + (1 | Subject)         Data: NS           REML criterion at convergence: 59570.13           Random effects:      Groups   Name        Variance Std.Dev.      Item     (Intercept)  1099     33.15        Subject  (Intercept)  4002     63.26        Residual             16153    127.10       Number of obs: 4737, groups: Item, 120; Subject, 43          Fixed effects:                         Estimate Std. Error t value     (Intercept)          681.373     13.179   51.70     Related2              13.655      8.379    1.63     WordType2             1.220     12.729    0.10     WordType3             5.700     12.709    0.45     WordType4             5.199     12.696    0.41     WordType5           -48.050     12.690   -3.79     Related2:WordType2  -20.925     11.862   -1.76     Related2:WordType3   -6.780     11.869   -0.57     Related2:WordType4  -19.746     11.854   -1.67     Related2:WordType5   28.870     11.859    2.43          NS$WordType <- relevel(NS$WordType,"4")     # Changed the intercept to word type 4     contrasts(NS$PrimeType)       1 2 3 5     4 0 0 0 0     1 1 0 0 0     2 0 1 0 0     3 0 0 1 0     5 0 0 0 1          Linear mixed model fit by REML ['lmerMod']     Formula: RT ~ Related * WordType * (1 | Item) + (1 | Subject)         Data: NS          REML criterion at convergence: 59570.13           Random effects:      Groups   Name        Variance Std.Dev.      Item     (Intercept)  1099     33.15        Subject  (Intercept)  4002     63.26        Residual             16153    127.10       Number of obs: 4737, groups: Item, 120; Subject, 43          Fixed effects:                         Estimate Std. Error t value     (Intercept)         686.5723    13.1779   52.10     Related2             -6.0911     8.3859   -0.73     WordType1           -5.1994    12.6957   -0.41     WordType2           -3.9794    12.7283   -0.31     WordType3            0.5006    12.7083    0.04     WordType5          -53.2497    12.6895   -4.20     Related2:WordType1  19.7457    11.8541    1.67     Related2:WordType2  -1.1791    11.8656   -0.10     Related2:WordType3  12.9659    11.8734    1.09     Related2:WordType5  48.6158    11.8662    4.10      I have also came across afex package -- first intended to calculate p-values. I noticed it controlled and chooses a different default when running the model.               > contrasts(NS$WordType)     # default intercept by afex       [,1] [,2] [,3] [,4]     1    1    0    0    0     2    0    1    0    0     3    0    0    1    0     4    0    0    0    1     5   -1   -1   -1   -1               Linear mixed model fit by REML ['lmerMod']     Formula: RT ~ Related * WordType * (1 | Item) + (1 | Subject)         Data: NS          REML criterion at convergence: 59583.5           Random effects:      Groups   Name        Variance Std.Dev.      Item     (Intercept)  1099     33.15        Subject  (Intercept)  4002     63.26        Residual             16153    127.10       Number of obs: 4737, groups: Item, 120; Subject, 43          Fixed effects:                         Estimate Std. Error t value     (Intercept)         679.1560    10.2793   66.07     Related1             -4.9693     1.8766   -2.65     WordType1            9.0442     7.0943    1.27     WordType2           -0.1982     7.0930   -0.03     WordType3           11.3543     7.1003    1.60     WordType4            4.3708     7.0972    0.62     Related1:WordType1  -1.8580     3.7485   -0.50     Related1:WordType2   8.6044     3.7544    2.29     Related1:WordType3   1.5319     3.7579    0.41     Related1:WordType4   8.0148     3.7509    2.14      When I saw the negative numbers in the contrasts, I wonder if I could produce the same output as afex when I adjust the intercepts.               > NS$WordType <- relevel(NS$WordType,"5")   # changed the intercept to word type 5     > NS$Related <- relevel(NS$Related, "2")          # changed word relatedness to 2 (unrelated)     > contrasts(NS$WordType)                               # both adjusted to mimic afex output       1 2 3 4     5 0 0 0 0     1 1 0 0 0     2 0 1 0 0     3 0 0 1 0     4 0 0 0 1          Linear mixed model fit by REML ['lmerMod']     Formula: RT ~ Related * WordType * (1 | Item) + (1 | Subject)         Data: NS          REML criterion at convergence: 59570.13           Random effects:      Groups   Name        Variance Std.Dev.      Item     (Intercept)  1099     33.15        Subject  (Intercept)  4002     63.26        Residual             16153    127.10       Number of obs: 4737, groups: Item, 120; Subject, 43          Fixed effects:                         Estimate Std. Error t value     (Intercept)         675.8473    13.1857   51.26     Related1            -42.5247     8.3944   -5.07     WordType1           19.1802    12.6920    1.51     WordType2           -0.5246    12.6622   -0.04     WordType3           18.1004    12.7000    1.43     WordType4            4.6339    12.7001    0.36     Related1:WordType1  28.8701    11.8595    2.43     Related1:WordType2  49.7949    11.8727    4.19     Related1:WordType3  35.6499    11.8805    3.00     Related1:WordType4  48.6158    11.8662    4.10      As you can see, the output is different from afex output. To summarize, my questions are as follows: What is the intercept in Linear mixed effects models? What do I need to know in order to make a wise decision in adjusting the intercept? Why afex package differ in choosing the intercept? Would it be the "best" approach in choosing the intercept for this type of design?


title: Computing confidence interval of clustered data

text: _Hi, this is the first time I am posting a question here so I apologize if my questions are not written in the proper manner and corrections will be welcomed_. I have a dataset in which each instance is defined by a cluster (C) and two more fields: trials and successes forming the following instance vector {C,trials,successes}. I am trying to estimate the confidence interval of $80\%$ for $\hat{p}=\frac{successes}{trials}$ given that different clusters have different homogeneity of variance and a certain unknown $\hat{p'}$ (the instances are not distributed normally, but more like Exponential distribution most of the time). For example, there can be two instances both with $\frac{successes}{trials}=\frac{5}{200}$ but as they come from different clusters their confidence interval can be very different. I would like to combine somehow the Wilson score interval with the cluster real measured values that form $80\%$ interval (maybe for all instances with enough trails - not sure what should be the threshold at the moment) and get the expected confidence interval. I would be happy to receive ideas or references about how it can be done (any other ways that can lead to integrated confidence interval will also be welcomed).


title: Partitioned sum of squares

text: I'm trying to calculate partitioned sum of squares in a linear regression. In the first model, there are two predictors. In the second model, one of these predictors in removed. In the model with two predictors versus the model with one predictor, I have calculated the difference in regression sum of squares to be 2.72 - is this correct? If this is correct, why is the difference in sum of squares 2.72 (quite small) when the difference in r-squared is ~30% (quite large).               # model two predictors     mod <- lm(drat ~ hp + wt, mtcars)          # regression sum of squares two predictors     regressionSumSquares <- sum((predict(mod)-mean(mtcars$drat))*(predict(mod)-mean(mtcars$drat)))          # residuals two predictors     residualSumSquares <- sum((predict(mod)-mtcars$drat)*(predict(mod)-mtcars$drat))          # r-squared two predictors     totalSumSquares <- regressionSumSquares + residualSumSquares     rSquared <- regressionSumSquares/totalSumSquares                    # model one predictor      modJustHp <- lm(drat ~ hp, mtcars)          # regression sum of squares one predictor     regressionSumSquaresJustHp <- sum((predict(modJustHp)-mean(mtcars$drat))*(predict(modJustHp)-mean(mtcars$drat)))          # residual sum of squares one predictor     residualSumSquaresJustHp <- sum((predict(modJustHp)-mtcars$drat)*(predict(modJustHp)-mtcars$drat))          # r-squared one predictor     totalSumSquaresJustHp <- regressionSumSquaresJustHp + residualSumSquaresJustHp     rSquaredJustHp <- regressionSumSquaresJustHp/totalSumSquaresJustHp          # difference in sum of squares one vs two predictors     regressionSumSquares - regressionSumSquaresJustHp


title: Correct setup for leave-one-subject-out cross-validation

text: I've got a question concerning leave-one-subject-out cross-validation of a classifier and correct outlier handling in this case. Let's suppose I've got 5 subjects. Within each subject the features are standardised (mean subtracted and divided by standard deviation) missing values or outliers are represented by NaNs. Now I'm doing an outer leave-one-subject-out cross-validation. In every outer fold a classifier is trained on the data of 4 subjects and tested on the fifth held out subject. NaNs are filled with non-informational values i. e. the mean of the training or test set respectively. First question: **Setting NaNs to the mean of the training or test set is ok?** Second question: **Is it ok not to standardise the features within the training set?** So I have then standardised for every subject but then do not standardise them over subjects in the training set? Now within each outer fold I plan to tune a classifier's parameter with help of another cross-validation. Third question: **Does the inner cross-validation necessarily need to be leave-one-subject-out as well?** I suppose, but can't really see the differences and implications of using another folding for the inner cross- validation. Fourth question: **For the inner fold the NaNs again need to be set to non- informational values (mean) of the corresponding training "subset"?** Thanks in advance for help, suggestions and remarks concerning the correct setup for this!


title: When is a fixed effect truly fixed?

text: Consider a linear unobserved effects model of the type: $$y_{it} = X_{it}\beta + c_{i} + e_{it}$$ where $c$ is an unobserved but time-invariant characteristic and $e$ is an error, $i$ and $t$ index individual observations and time, respectively. The typical approach in a fixed effects (FE) regression would be to remove $c_{i}$ via individual dummies (LSDV) / de- meaning or by first differencing. **What I have always wondered: when is $c_{i}$ truly "fixed"?** This might appear a trivial question but let me give you two examples for my reason behind it.   1. Suppose we interview a person today and ask for her income, weight, etc. so we get our $X$. For the next 10 days we go to that same person and interview her again every day anew, so we have panel data for her. Should we treat unobserved characteristics as fixed for this period of 10 days when surely they will change at some other point in the future? In 10 days her personal ability might not change but it will when she gets older. Or asked in a more extreme way: if I interview this person every hour for 10 hours in a day, her unobserved characteristics are likely to be fixed in this "sample" but how useful is this?   2. Now suppose we instead interview a person every month from the start to the end of her life for 85 years or so. What will remain fixed in this time? Place of birth, gender and eye color most likely but apart from that I can hardly think of anything else. But even more importantly: what if there is a characteristic which changes at one single point in her life but the change is infinitesimally small? Then it's not a fixed effect anymore because it changed when in practice this characteristic is quasi fixed. From a statistical point it is relatively clear what is a fixed effect but from an intuitive point this is something I find hard to make sense of. Maybe someone else had these thoughts before and came up with an argument about when a fixed effect is really a fixed effect. I would very much appreciate other thoughts on this topic.


title: Forecasting of density function

text: I am doing some research about forecasting time series of probability density functions. We are aiming to forecast a PDF given historically observed (usually, estimated) PDF. The forecasting method we are developing performs pretty well in simulation studies. However, I need an numerical example from real applications to illustrate our method further. So, are there any proper examples in applications (finance, economics, biology, engineering, etc.) where a time series of PDFs are collected and it is important and difficult to forecast such a time series?


title: When is Fisher's z-transform appropriate?

text: I want to test a sample correlation $r$ for significance, using p-values, that is $H_0: \rho = 0, \; H_1: \rho \neq 0.$ I have understood that I can use Fisher's z-transform to calculate this by $z_{obs}= \displaystyle\frac{\sqrt{n-3}}{2}\ln\left(\displaystyle\frac{1+r}{1-r}\right)$ and finding the p-value by $p = 2P\left(Z>z_{obs}\right)$ using the standard normal distribution. My question is: how large $n$ should be for this to be an appropriate transformation? Obviously, $n$ must be larger than 3. My textbook does not mention any restrictions, but on slide 29 of this presentation it says that $n$ must be larger than 10. For the data I will be considering, I will have something like $5 \leq n \leq 10$.


title: How does a uniform prior lead to the same estimates from maximum likelihood and mode of posterior?

text: I am studying different point estimate methods and read that when using MAP vs ML estimates, when we use a "uniform prior", the estimates are identical. Can somebody explain what a "uniform" prior is and give some (simple) examples of when the MAP and ML estimators would be the same?


title: Initial descriptive analysis of a database join

text: I'd like to do an analysis to identify possible patterns in users demands against users demographics. I won't go into the details regarding the analysis, my question is mainly about the data cleaning/preparation. My data is coming from 2 tables in MySQL:   * user table (user_id, list of demographics), with 200 rows   * demand table (user_id,demand_id,product information) - with 3000 rows Initially I joined the two tables on user_id and imported in a unique query using the `<-dbGetQuery` command in R. While I was doing some initial descriptive analysis of the demographics variables, I realised that the whole exercise was wrong because the user data were repeated many time, as there my data contained now more than one demand for each user. As a result I'd assume that the distribution of the users age, for example, would be incorrect as the age group wich has more demands will count more than it should. So I stopped my analysis and started with a different approach:   * import user table only   * do the descriptive analysis of the user demographics    * import the demand table in a separate data frame   * join the two data frames in R   * proceed with the analysis Can you confirm if this is a better approach? Thank you for your help


title: Specifying nested terms in glmer

text: I have collected reproductive data from 393 trees from 51 provenances across the distribution of my study species. At each of these sites I have also recorded the level of fragmentation for each of the collected trees (1-4, where 1 is single paddock tree and 4 is an closed canopy forest). I want to test the influence of fragmentation on the reproduction data I collected using a GLMM (i.e. `lme4::glmer`). In these analyses I have treated fragmentation as a fixed factor and fitted the individual fragmentation scores of each tree as a nested factor within provenances in the following model: `glmm.fit1 <- glmer(response ~ fragScore + (1|prov), data = data, family = gaussian (link="identity"))` [NB. I realise the above code is equivalent to a `lmer` model when specifying `family = gaussian`, I have just kept it as a `glmer` as I will be running count/proportion response variables as well]. I am a bit confused on the correct specification of the nested term as another blog (see jpkelley response here) suggests `(1|prov)` is equivalent to `(1|prov/fragScore)` if prov is individually coded (which here it is). Whereas, Ben Bolker in his response to a similar question suggests `(1|prov/fragScore)`. Is this the correct syntax to specify tree nested within provenance or should I use `(1|prov/fragScore)`? Thanks in advance! PS. this is my first posting so apologies if it is a bit vague.


title: Gumbel Copula generation using nonparametric correlations like Kendall's tau

text: I have 2 different variates W,X. I want to compute Gumbel copula for these variates. I followed following steps for the same: 1\. To compute kendall's tau I used R's package Kendall. From kendall's tau I computed gumbel parameter(alpha) as alpha=1/(1-tau).   1. Then I computed gumbel copula using R's package copula. function used during the computation was gumbelCopula(alpha,dim=2).   2. After that I used rCopula() function to simulate copula values. Above method worked fine for 2 variates W,X. Now when I wish to compute gumbel Copula for 4 variates W,X,Y,Z using the same procedure, it is failing in 1st step itself as I cannot compute kandall tau for 4 variates. Kendall tau works only for 2 variates. Can someone please illuminate me with the method of computing Gumbel copula for 4 variates? Thanks in Advance


title: How to forecast time-series bounded by [0,1], i.e., forecast relative frequencies?

text: I am working with time series values which are all in the closed interval [0, 1]; these values represent relative frequencies, i.e., empirical probabilities. I would like to create a model such that all forecasted values are within [0, 1], but it would also be fine if the model's output was strictly within the open interval (0, 1). This answered question tackles the lower bound aspect of my question, but not the upper bound aspect: How to achieve strictly positive forecasts? I'd like to use the R `forecast` package if possible to achieve this, but I am open to other suggestions.


title: Marginalized local odds ratios structure in GEE

text: The `ordLORgee` command from the `multgee` package in R has an argument LORstr. With this argument the user can specify the marginalized local odds ratio structure. What is a marginalized local odds ratio structure? How do I determine which marginalized local odds ratio structure to use? I have subjects nested in countries, with scale and categorical independent variables at subject level and scale independent variables at country level. The focus of the research is to find the effects of the country level variables on the dependent variable. The author of the `multgee` package references the article below, but I cannot access it. Touloumis, A., Agresti, A. and Kateri, M. (2013). GEE for multinomial responses using a local odds ratios parameterization. Biometrics, 69, 633-640.


title: How to speed up clustering analysis in R?

text: I am conducting a simulation study, which involves several clustering methods such as model-based clustering method (MBCM). However, there is a big problem I find that MBCM runs so slowly (might be due to the EM algorithm). I have tried my best to avoid for loop in my simulation setting but it still takes me forever. For example, I simulate 100 datasets and each data set contains 1000 subjects. This tiny simulation took me around **879** minutes to complete. I can't image how much time will take me to complete 1000 simulated datasets. Does anyone have a better idea to speed up the clustering analysis in R? many many thanks Tu


title: How is the number of particles decided in particle filtering?

text: Say the observations are $x_i$ and the states are $y_i$ in a sequential model. I understand that particle filtering works by generating "particles" from $p(y_i | x_1,\ldots,x_i)$ for approximating $p(y_{i+1} | x_1,\ldots,x_{i+1})$. How do we decide on how many "particles" to use as we go along on the chain? Do we choose a fixed number in the beginning, and stick to it (one that works well experimentally), or do we change the number of particles used as the particle filtering algorithm proceeds?


title: A mixed model with gamma distribution: handling zeros

text: I have a panel cost data with some zero values (not missing values but rounded to 0). How should I handle zeros when I use SAS Proc Glimmix with Gamma distribution? Maybe change zero to very small values? Or is there an option for it (which I couldn't find) in PROC GLIMMIX? Janne


title: Sum of squared difference and Gaussian noise model

text: I have been reading that when the underlying error is distributed normally, then minimising the sum of squared difference between the observed data and the model is the appropriate cost function to do the model fitting. However, I am having trouble understanding why that is the case. I was wondering if there is an intuitive way to understand why this is the case? What would happen if I use other cost functions say point wise mutual information or some local measure like cross correlation.


title: Is the variance of each group unknown in one-way ANOVA?

text: I am trying to understand ANOVA. In one-way ANOVA to test equality of mean of several population groups, each group has a normal distribution with the same variance. I wonder if the common variance of each group is known and used somewhere or unknown and estimated from samples in ANOVA? Thanks and regards!


title: Meaning of covariance matrix row sums

text: Say I have an $n \times n$ covariance matrix for a sample set of $n$ random variables. Is there any meaning if the sum of the rows of this matrix? Is it a meaningful measurement of the contribution of each random variable to the variance of sum of the random variables?


title: Why might the means estimated by GLM ANOVA differ from those calculated by hand?

text: I'm currently analysing some EEG data using a repeated measures general linear model, and I've run into a problem. I have three within-subject factors; Scalp area (3), condition (2) and word (2). When I run a three-way analysis I get mean values for Word (separated by condition) which don't match up to the means I've calculated for myself, for condition 1 (those for condition 2 seem fine). When I take the relevant columns of data and enter them into a two-way analysis for just one condition (so Scalp area (3) x Word (2)), I get the correct means, and those for condition 1 don't match the ones from the three- way analysis. I've tried this in SPSS and Statistica with identical results; I've checked the data, and I'm pretty sure it's not an input error. I've checked the order in which I'm entering the variables (slowest-changing first), and that seems fine. So either I'm doing something wrong in the analysis, or I've misunderstood the output, or... something else? I may well be staring straight past an obvious error, but I've spent all afternoon tearing my hair out, so I'm hoping someone might be able to offer some insight! Thanks in advance!


title: Passing a single item and getting a prediction

text: I have completed a demo of a naive Bayes classifier for predicting whether an SMS message is spam or ham. My question is, how can I use this in a practical manner? It appears I would have to completely retrain the classifier every time there is a new message. Ideally, I would like to simply pass in a new string to a method and the output would tell me the probability of it's being spam or not. Here is how I am currently predicting:               # a lot of text manipulation to get the string values into a      # suitable format precede these lines          sms_classifier <- naiveBayes(sms_train, sms_raw_train$type)          sms_test_pred <- predict(sms_classifier, sms_test)      The results:               Total Observations in Table:  1390                             | actual         predicted |       ham |      spam | Row Total |      -------------|-----------|-----------|-----------|              ham |      1204 |        31 |      1235 |                   |    16.148 |   106.504 |           |                   |     0.998 |     0.169 |           |      -------------|-----------|-----------|-----------|             spam |         3 |       152 |       155 |                   |   128.660 |   848.596 |           |                   |     0.002 |     0.831 |           |      -------------|-----------|-----------|-----------|     Column Total |      1207 |       183 |      1390 |                   |     0.868 |     0.132 |           |      -------------|-----------|-----------|-----------|      **Edit for clarification:** I initially have 5k rows like this:               type,text     ham,Hope you are having a good week. Just checking in     ham,K..give back my thanks.     ham,Am also doing in cbe only. But have to pay.      I then perform a number of manipulations using `tm` package to be able to run the naive Bayes classifier to predict the type of the unseen data based on the text. I would like to simply pass a string to a method to utilize my existing classifier without having to retrain the entire thing, like so:               classifyNewItem <- function(text) {         # run new item against classifier here     }      **Is this how this is generally done? What is the preferred method of utilizing a classifier in a real world scenario such as a spam filter?**


title: Comparing Counts Between Two Independent Groups

text: I have 2 groups of patients and I'm looking at what percentage were discharged from the hospital at different time points after surgery.                           Day 1     Day 2     Day 3     Day 4     Group 1     12%       10%       25%       12%     Group 2     50%       20%       20%       15%      What is the most appropriate statistical test for this? These are two separate groups of patients and the total number of patients in each group is different. As far as the raw numbers, some of the values are less than 5 so I don't think I can use chi-square. I would appreciate your thoughts.


title: How to know whether the six spatial features values (in 3-time period) significantly differ from each other?

text: I am currently bothered on how to statistically analyze my simple data as shown below.               Date             Spatial Features               1      2      3      4      5      6      ----------------------------------------------     1988 0.7784 0.8310 0.8726 0.9307 0.9189 0.8814     1998 0.8956 0.9036 0.9368 0.9575 0.9283 0.9045     2009 0.9255 0.9479 0.9703 0.9797 0.9434 0.9555      **Note** : These values are Relative Entropy values. For my purpose, Im using Relative Entropy to measure compactness or dispersion. Entropy value ranges from 0 to 1, where 0 means compact, while 1 means dispersed. The entropy values were calculated for every spatial feature in 3-time period (1988, 1998, & 2009). As shown in the table, there are six spatial features being considered numbered 1, 2, 3, 4, 5, & 6. My goal is to prove whether there is a significant difference among these values of compactness/dispersion based on spatial features. Or whether compactness/dispersion differs significantly with spatial features. Does the time (year) also need to be taken into account in the analysis? Thank you very much and I would really appreciate any help or advice you could offer me on how to attain my goal, particularly on the appropriate stat test or method to be used.


title: Testing significance of peaks in spectral density

text: We sometimes use spectral density plot to analyze periodicity in time series. Normally we analyze the plot by visual inspection and then try to draw a conclusion about the periodicity. But has the statisticians developed any test to check whether any spikes in the plot are statistically different from white noise? Has the R-experts developed any package for spectral density analysis and for doing that kind of test? Great if someone could help. Regards,   P.


title: Measuring predictive accuracy for multiple dependent variables

text: In machine learning and in statistics there exist plenty of measures which estimate the performance of a predictive model. For example, classification accuracy, area under ROC curve ... for classification or $R^2$ for regression models. Is there any generalisation of such measures for multiple dependent variables in order to evaluate predictive accuracy of Multivariate Multiple Regression, Partial Least Squares Regression or even more complex predictive models where the dependent variables can be of mixed type?


title: Multiple dependant variables in logistic regression

text: I have a binary dependent variable, say food security status, where 0 = _food insecure_ and 1 = _food secure_ , and this response variable is recorded for two categories of households, _small farmers_ and _land less rural house holds_. Both these household categories may have food secure or food insecure. These dependant variables are a function of some common independant variables. I want to check which independant variable has more impact on farmers or landless households' food security. Is it possible to include it in one logistic regression model? Or do I need to run separate models for farmers and landless households?


title: Neural network with PROC logistic in SAS

text: I am wondering if somebody has used or has any reference to point as to how a neural network in conjunction with PROC logistic can be used in sas. A code would be an ideal thing for me but if unavailable, if somebody can pass me a link as to how can I benefit myself from using neural network alongside logistic ?


title: Does Pearson correlation require removal of bivariate or univariate outliers?

text: Does the Pearson's correlation estimator require no bivariate outliers, or no outliers in each of two individual vectors of data? The answer will impact on how I winsorize outliers before calculating my Pearson estimate.


title: Can I use standard deviation in Likert scale (1-7) data with non normal distribution?

text: I have 20 questions, each divided in 6 Likert scale items (1-7), representing 6 different values like Success, Oportunity, etc. The sample is rather small, 42 cases, and the distribution is not normal. Basically, I want to compare responses to values between individuals in certain groups consisting of 4-5 individuals, and see how much do the responses differ for each value within these mini-groups. Is it ok to use standard deviation? Please suggest any nonparametric tests that might be more suitable for this analysis. Any reccomandation is welcome!


title: Is high multicollinearity always an issue in OLS?

text: $$Y_t = a + bX_{1,t} + cX_{2,t} + dX_{3,t} + e_t$$ A high $R^2$ in $X_{1,t} = \alpha + \beta X_{2,t} + \gamma X_{3,t} + \varepsilon_t$ will always result in a higher standard error of the $b$ estimator, _all else held equal_. This is an example of the multicollinearity problem. It does not result in a biased or inconsistent estimator. _**Question_** : Assuming correct specification and the satisfaction of the other assumptions, if my statistical test is powerful enough to detect significance in $\hat{b}$ at an acceptable test size, **do I need to continue to worry about multicollinearity**?


title: Random music note generation

text: I'm working on a little project, which is in early development stage right now. What I'm trying to figure out is how do I generate random music, especially notes. By googling I've only managed to find _ransom_ note generator and lots of random music, where word random has nothing to do with definition of random variable :) I don't know if anybody have done this before and I wish to know about presence of such work to give me some guidance. Here some of my naive ideas: Given the root note of the chord generate a random note with mean of that note and some deviation, which might be controlled by part of the bar which you are currently in. I would prefer to set 0 as a reference to root note in every bar and generate integers, which would mean how many scale intervals up or down melody should deviate from the root. But this approach might make joining melody bars together harder, so another idea is to use the key of whole part as a reference note and then change the mean of random variable according to which bar corresponds to which key. So now, that I have explained a little on how do I expect the mean of the random note to be, lets move on to random number generation in the bar or a piece of music. What I need is a sequence of integers with mean $0$. Natural candidates could be   1. Sum of two Poisson distributed variables, one of them should be multiplied by $-1$   2. The same thing with binomial random variables. Maybe you can suggest something better? For the rhythm generation I will probably set the fixed number of maximum notes per bar and then generate sequence of zeroes and ones, the one will indicate that beat is changing. Maybe I will will set probability of zero conditional on how many zeroes were before. If anybody has ideas I would love to hear them!


title: What is the computational complexity of the EM algorithm?

text: In general, and more specifically for Bernoulli mixture model (aka Latent Class Analysis).


title: Benchmark data for Random Forest evaluation

text: I would like to make/explore several alternative formulations of random forest (link) and am looking for the current state of the science for evaluating performance. In the reference file Leo provides "microarray", "DNA", "glass", and "spectral" data to exhibit the performance of his method. Are these available? How do I get them in order to first reproduce the results, and explore performance of changed parameter settings? How general and authoritative are they? Other references that I have thus far found include (this), (this), and (this). Can you suggest rigorous, modern, and relatively complete list of quality classification benchmarks to be used to compare the performance of very similar random forests?


title: Understanding the definition of omnibus tests

text: From Wikipedia > Omnibus tests ... test whether **the explained variance** in a set of data > is significantly greater than **the unexplained variance** , overall. > > ... > > Omnibus test as a statistical test is implemented on an overall hypothesis > that tends to find general significance between parameters' variance, while > examining parameters of the same type, such as: Hypotheses regarding > equality vs. inequality between k expectancies $µ_1=µ_2=…=µ_k$ vs. > $µ_1≠µ_2≠...≠ µ_k$ in Analysis Of Variance(ANOVA) ; I was wondering what it means by a variance being explained or unexplained? How shall I understand what an Omnibus test is? What is a set that is not omnibus? Thanks and regards! * * * Added: I was wondering if the following likelihood ratio test for the global null is an omnibus test, quoted from Casella and Berger's Statistical Inference? ![enter image description here](http://i.stack.imgur.com/wnMj0.png) ... ![enter image description here](http://i.stack.imgur.com/eeUOZ.png)


title: Python module for change point analysis

text: I'm looking for a Python module that performs a change-point analysis on a time-series. There are a number of different algorithms and I'd like to explore the efficacy of some of them without having to hand-roll each of the algorithms. Ideally I'd like some modules like the bcp (Bayesian Change Point) or strucchange packages in R. I expected to find some in Scipy but I haven't been able to turn up anything. I'm surprised that there aren't any facilities in:   * statsmodels.tsa: Time series statistical analysis tools   * scikits.timeseries: Time series analysis tools to extend scipy   * scipy.signal: signal processing tools in scipy Are there any modules with change point detection algorithms in Python?


title: How to get a Standard Error for a Poisson Variable with Exposure?

text: I'm looking over some work performed by a consultancy and I'm unsure whether the standard error formula they have used is correct, and the subsequent conclusions they have drawn are erroneous. They have provided a Standard Error between the difference in observed claim frequencies and modelled claim frequencies, and then checked whether the modelled frequency lies between the observed 95% Confidence Interval. The Standard Error Formula they have used is: $\sqrt{\frac{Observerd Frequency * (1 - Modelled Frequency)}{Exposure}}$. However I would have thought a more appropriate Standard Error would be along the lines of the usual Standard Error for a Proportion: $\sqrt{\frac{Observerd Frequency * (1 - Observed Frequency)}{Exposure}}$. My problem is I don't know how appropriate is it to use the Modelled Frequency in the Standard Error Calculation? ## The data looks like this:   1. Sum Insured Band: \$0 - \$5,000    2. Exposure (Years): 27,233    3. Observed Claims: 1271    4. Modelled Claims: 1433    5. Observed Frequency: 4.7    6. Modelled Frequency: 5.3  I've had a look around at disease incidence and claims frequency standard error literature on google, but haven't really found anything concrete. Thanks for any help.


title: Probability of reaching a threshold after several trials

text: I was wondering if there wasn't some better way to get to the bottom of this question I had regarding the chances of reaching a certain threshold of guessing numbers correctly over a series of trials. This isn't my specialty, so would appreciate answers that are like you were talking to a 5 year old or at least google-able :) So say I have a game where I have to pick 10 out of 80 numbers and then 20 are drawn. For each one of my Picks that is chosen as one of the 20 I get a black marble. I need to know the odds that if I can get AT LEAST 20 black marbles in 5 games, picking 10 numbers each game. Each game is discrete in that the picked numbers are replaced and the 20 numbers that are drawn are totally random I can tell what the average amount of black marbles I will win per game (2.5) and therefore the average amount of black marbles I'll win in 5 games (12.5). I can do this because for a single game, hyper geometric distribution tells me the odds of each outcome of catches: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. I can then multiply each of these by their probabilities and sum them to get the average amount I will win in that game   e.g. $(HYPGEOMDIST(0, 10, 20, 80) * 0) + (HYPGEOMDIST(1, 10, 20, 80) * 1)$, etc Here's how I would do this by hand if I had to ( **this is the part I'd like to simplify / correct my assumptions on!!** ): Find each permutation of black marbles (0-10) I could win for each of the 5 games without regard to order. Filter each permutation so that if the set of 5 games does not sum to 20 it is discarded   For the remaining sets find and multiply the HYPGEOMDIST for each of the 5 games in each set, then sum every set Please excuse me if this has been asked somewhere else, I am unfamiliar with the terminology in order to search well enough on this site.


title: Real-life examples of moving average processes

text: Can you give some real-life examples of time series for which a moving average process of order $q$, i.e. $$ y_t = \sum_{i=1}^q \theta_i \varepsilon_{t-i} + \varepsilon_t, \text{ where } \varepsilon_t \sim \mathcal{N}(0, \sigma^2) $$ has some _a priori_ reason for being a good model? At least for me, autoregressive processes seem to be quite easy to understand intuitively, while MA processes do not seem as natural at first glance. Note that I am _not_ interested in theoretical results here (such as Wold's Theorem or invertibility). As an example of what I am looking for, suppose that you have daily stock returns $r_t \sim \text{IID}(0, \sigma^2)$. Then, average weekly stock returns will have an MA(4) structure as a purely statistical artifact.


title: Compare two teams who competed in slightly different competitions

text: Two running teams take part in competitions on the same course, but their competitions differ because of the number of teams entered in each competition. Results are based solely on position - times are not taken. Competition 1 has 7 teams of 8 runners Each team enters 2 of its runners in each of 4 races Points are awarded on the basis of 1 for first place down to 14 for fourteenth place The points for all 4 races are added together. The winning team had a total of 40 points. Competition 2 has 5 teams of 8 runners Each team enters 3 of its runners in each of two races (15 competitors in each) and a further 2 runners in a third race (10 competitors in this one). Points are awarded in each race on a similar basis to competition 1 (2 races award from 1 to 15 and the other, 1 to 10). The winning team has a total of 46 points Statistically, which is the better team? I'm not a statistician, but my logic lead me initially to compare each team's performance with what I would (probably incorrectly) call the mean number of points in each competition (average of the most and least number of points the team could score). However this would seem disadvantage the team in competition 2, because their best possible points score is 15 versus 12 for team 1. And the "mean" numbers are in any case similar. I did think about comparing each team's performance against the best case and worst case scenarios in their respective competitions and averaging those scores? The other thought is that I should somehow compare each runner's performance against the best & worst case? In reality, I'm a bit stuck. Help please!


title: Expectation of a product of multiple random (Bernoulli) variables

text: I'm doing a research and using random variables to model a random process. I'm defining a Bernoulli random variable as a product of several other Bernoulli variables (three or more). So, I have the following as the expectation: $E[X_1]=E[X_2X_3X_4]$ I cannot assume the independence of the variables $X_i$ for any $i$. The issue is to find a correct formula for the expectation to consider the correlation of the variables $X_2$ to $X_4$. I know that for a product of two variables, I can have: $E[X_2X_3]=E[X_2]E[X_3] + Cov(X_2,X_3)$. But this covariance will be a matrix in case of three or more variables. Then, how would one compute a value for $E[X_1]$? I appreciate if somebody can help with this or point me to a good resource (e.g. book).


title: What to do against sampling/selection bias?

text: Let's say in machine learning the test cases are biased towards one or more options. For example when there are 100 (just an example) test cases in total and every test case has a color to be suggested: 60 => blue   30 => green   9 => red   1 => yellow Most machine learning techniques will suggest blue just for probability reasons. To avoid this some write you should use the same amount of test cases for every answer. But that way I only have 4 test cases. Also when principal component analysis does not bring any insights, what can be done to help machine learning algorithms like neural networks or genetic algorithms to avoid the sampling bias towards the given training data? Would it help to create/train one for each color? Or what would you do?


title: Appropriate use of difference adjusted CIs in error bars for mixed ANOVA's time * group interaction

text: I have performed a mixed design ANOVA with time (pre-post) and group (treatment vs control) as variables. I'm wondering if it's sensible to use _between subjects difference adjusted CIs_ (from Baguley, 2012) in the interaction plot below. What I want to emphasize is the group difference at post-test (and of course non-significant difference at pre-test). The error bars seem to convey the statistical significance of the time x treatment interaction, which was _p_ = 0.0451. The error bars are calculated separately for pre- and posttest. Also, I'm thinking it'd be _incorrect_ to interpret these errorbars as evidence that, say, treatment-pre is significantly different than control-post, is this correct? Lastly, I'm aware that these error bars are not appropriate in interpreting differences within subject (and the figure caption will tell what the CIs represent). But is it sensible to use them like I've explained here, or is there some better way to graphically convey the significant group x time interaction? Thanks in advance for any input on this matter! ![enter image description here](http://i.stack.imgur.com/cCJiA.png) Baguley, T. (2012). Calculating and graphing within-subject confidence intervals for ANOVA. Behavior Research Methods, 44, 158-175.


title: Incorporating averaging models from AIC and still using k-fold cross validation?

text: I've a county/district that I've divided into ~300 grids that are 15km^2 in size attributed with various habitat and economic variables that have been summarized and standardized. I then have 2 types of data 1) presence/absence of leopard caused human deaths and 2) counts of those human deaths. Of the 300 grids 80 had >=1 human death. I came up with 14 a-priori hypotheses to test. Initially, I simply ran AIC model selection and averaged the best models. The reviewers said this was incomplete as I still need to validate the models. It was suggested I do k-fold cross validation (KCV) to validate the model and get at model performance or predictive accuracy. However, I am having difficulties figuring out how to do that when multiple models (say 3 models for example) are the "best" according to AIC and how you would do KCV over an averaged model in R. I've also been seeing that some people have used CV methods for model selection itself. That is, they select the "top model" based on the one that predicts the best instead of using AIC, BIC, etc. But, model selection using AIC and CV are asymptotically equivalent. That is, as sample size increase the AIC top model will converge to the "best predicting" top model selected using CV methods. How would I know if my sample size is "enough" to use CV instead of AIC? Can I still use AIC - and then only run the KCV on the absolutely best top model? Or run the KCV over all three best models and then somehow compare between them? I found the R package cvTools which has a cv.Select command, but I have yet to see any discussion of interpretation for the results or to know what it is exactly I am doing outside of the PDF that has the short short example. Basically k-fold type techniques were developed for assessing predictive accuracy for a single model... I've got multiple models according to AIC that work -- so what to do? Is the reviewer correct in saying I need to validate the model? If CV and AIC model selection are asymptotically equivalent. That is, as sample size increase the AIC top model will converge to the "best predicting" top model selected using CV methods -- would I even have to do this as CV would essentially give the same result?


title: Interpreting  negative binomial regression output in R

text: I've got a question concerning a negbin regression: suppose that you have the following commands               require(MASS)     attach(cars)     mod.NB<-glm.nb(dist~speed)     summary(mod.NB)     detach(cars)      (note that cars is a dataset which is available in R) and don't care if this model makes sense. What I'd like to know is, how can I interpret the variable theta (at the bottom of the summary). Is this the shape parameter of the negbin distribution and is it possible to interpret it as a measure of skewness? I appreciate every thought!


title: Finding similarity between the objects ( Actually how to justify the similarity of objects in mathematics form )

text: I have developed a system to trace the outlines of (images of) objects. Now I want to test whether two independent traces represent a common feature. Imagine two people (or machines) tracing the outline of a feature in an image, recording it as sequence of vertices. Inaccuracies in recognizing the feature boundaries and in specifying the vertices can be viewed as random errors in vertex positions. The problem is the two traces might use completely different (Cartesian) coordinate systems (set up on two digitizing tablets, for instance). The null hypothesis to test is that they represent a common feature. This is illustrated below. I drew a figure and recorded the $x$ and $y$ co- ordinates of its vertices in its coordinate system. Let's call this figure $M$. It is represented as a sequence $(x_i,y_i), i=1, 2, \ldots, m$. Then I drew the same figure in a bigger size in another coordinate system (with no known relationship to the first coordinates system) and recorded the $x$ and $y$ co-ordinates. Let's call this $N$, represented as a sequence $(x_i^\prime, y_i^\prime), i=1, 2, \ldots, n$. **Question:** Having these data, how can I test whether the figures $M$ and $N$ represent the same image features even though they have different sizes and co-ordinates? If the conclusion is that yes, they do represent a common object, then how can I estimate a similarity transformation between $M$ and $N$ so that I can work on formulas or equations to check the results with different figures? ![enter image description here](http://i.stack.imgur.com/VNR37.png)


title: Feature Selection

text: I've a facebook users dataset in which each user has a "huge" set of attribute, i.e about 220 attributes like age, hometown, religion, and a set of facebook liked pages to store the users tastes. Now I need an automatic tool (i've no time to develop an entire framework) to select only the relevant features in the dataset, in order to make agile the cluster algorithms that i've to use. My only constraint is that I need a **unsupervised** features selection technique, since I've no further information than my initial dataset. Any help?


title: How to get Cook's distance and carry out residual analysis for non-lm() and non-glm() models in R?

text: I usually use the `plot(lm())` or `plot(glm())` (combined with `par(mfrow=c(2,2)`) to analyze residuals. Unfortunately this is not possible whenn estimating a `tobit()` (package AER) or `zeroinfl()` (package pscl) model. Furthermore any Cooks distance command does not work. It's either not possible to make the residual analysis in a convienent way, or I'm just too dumb (this is probably the reason). Any help is appreciated! Here is a short example of what I am talking about:               require(AER)     require(stats)     attach(cars)     tob<-tobit(dist~speed)     summary(tob)     par(mfrow=c(2,2))     plot(tob)           #doesn't work     cooks.distance(tob) #doesn't work     detach(cars)


title: Inference on a probabilistic graphical model with observed continuous variable

text: Suppose we have the following probabilistic graphical model: $L \rightarrow X, L \rightarrow Y, F \rightarrow Y, X \rightarrow D, Y \rightarrow D$. That is, the joint is: $P(D,X,Y,L,F) = P(D|X,Y)P(X|L)P(Y|L,F)P(L)P(F)$. Here $L,F,D$ are discrete, $X,Y$ are continuous random variables (having proper pdf's). I want to compute $P(D=d,L=l | X=x, F=f)$. Note that we have already observed ($x$) for one of the continuous variables. Here is my solution attempt so far: $P(d,l | x, f) = \frac{P(d,l,x, f)}{P(x,f)}$ Let's ignore the denominator for now. The numerator: $P(d,l,x,f) = \int_y P(d|x,y) P(x|l) P(y|l,f)P(l)P(f)$ Now, $P(x|l)$ and $P(y|l,f)$ should be zero because they are continuous variables. But we have already observed the value of X. So, I am not sure how to proceed from now on. My guess is that we should use the pdf's instead of probabilities: $P(d,l,x,f) = \int_y P(d|x,y)\;f_{X|L}(x|l)\;f_{Y|L,F}(y|l,f)\;P(l)\;P(f)dy$ Is this correct? If this is correct, then another question is how to compute this integral. Would the following work? Sample a $y$ from $f_{Y|L,F}(Y|l,f)$, compute $P(d|x,y)$ and then compute $g(y) = P(d|x,y)\;f_{X|L}(x|l)\;P(l)\;P(f)$. Do this a large number of times and take the average of $g(y)$'s. ??? Any ideas, hints, directions would be highly appreciated.


title: Are there any good movies involving mathematics or probability?

text: Can you suggest some good movies which involve math, probabilities etc? One example is 21. I would also be interested in movies that involve algorithms (e.g. text decryption). In general "geeky" movies with famous scientific theories but no science fiction or documentaries. Thanks in advance!


title: Why we need single class classification?

text: I have started learning classification in machine learning. I face two terminologies, one is "single class classification" and a the other is "binary class classification". I am confused about when to use single class classification and when to use binary class classification. My question mainly concentrates on single class classification. **Why do we need single class classification?**


title: Correct scaling factor going from monthly to yearly data

text: I'm trying to compute the probability of an extreme event, using extreme value theory on my data. One of the typical problems with this kind of computations is the lack of sufficient data, which also occurs in my case. More specifically I would like to become a yearly occurance rate of a certain event, but I only have a few years of data. I do however have data for each of the months where I compute the rate from the number of occurances in the total population. However, as stated earlier, in the end I would like to obtain a yearly rate and as such I've opted to multiply my monthly rates by a factor 12 to obtain the yearly occurance rate in each month. However, this will result in an increased volatility which I would like to avoid to not make an overestimation of my final result of the yearly rate. Therefore, my question is the following: which scaling factor is optimal to obtain valid yearly rates each month, but without making an overestimate of the final exceedance rate due to the increased volatility I've introduced. Is this simply $sqrt(12)$ or is there a more appropriate value? ...and what about the fact that not all months have the same number of days?


title: How can I show that a random walk is not covariance stationary?

text: How can I show that a random walk ($y$ follows a random walk) is not covariance stationary? I tried to work on the formula below (with no results) Could you give me just a hint on how to proceed please? $$Cov(y_{t+h},y_t)=E(y_{t+h}\times y_t)-E(y_t)E(y_{t+h})$$ Is this approach right? Important: $\epsilon_t$, the shock, is an iid sequence with mean $0$ and variance $\sigma^2_\epsilon$. * * * If $y$ follows a RW we have $$y_t=y_{t-1}+\epsilon_t$$ then, $$Var(y_t)=Var(y_{t-1}+\epsilon_t)=Var(y_{t-1})+\sigma_\epsilon+2Cov(\epsilon_t,y_{t-1})$$ Now I see that the variance of $y_t$ depends on the variance of $y_{t-1}$. This should suggest me that we lack covariance stationarity.


title: 2 treatment levels and data from multiple time points - PERMANOVA or repeated measures MANOVA?

text: I have a dataset with continuous data (i.e. several time points) from an experiment with 2 treatments (i.e. independent variables) and 2 replicate measurements per time point and treatment. What is the right statistical analysis to apply? I tried MANOVA and PERMANOVA for distinct time points, but using this approach I miss dynamics from the entire time course. Is it possible to add "time" as another independent variable? I have the feeling that this would be wrong, as each time point is dependent on the previous time point. I read about repeated measures approaches, but not entirely sure how this can be combined with a multi-factorial analysis... I'd really appreciate help... Thanks, John


title: K-means clustering and R - what to do next?

text: I used `kmeans` command on my data-frame (as suggested in "R and Data Mining: Examples and Case studies"). Now my data is clustered into x number of cluster.   What they don't tell you in this book is what to do next.   How can I get `characteristics` specific for `cluster 1`.   Or how does `cluster 2` differ from `cluster 3`? Example of my data-frame after clustering.               Cluster   Char1   Char2   Char3   Char4   Char4     1         0.00    0.02    1.23    3.21    2.34     1         2.12    12.1    1.42    1.31    2.04     2         5.35    59.2    0.01    9.32    9.33     3         5.23    10.3    8.13    0.72    0.91        ...          Is there more informative guide on the internet? How can I investigate these things?


title: Clustering inertia formula in scikit learn

text: I would like to code a kmeans clustering in python using pandas and scikit learn. In order to select the good k, I would like to code the Gap Statistic from Tibshirani and al 2001 (pdf). I would like to know if I could use inertia_ result from scikit and adapt the gap statistic formula without having to recode all the distances calculation. Does anyone know the inertia formula used in scikit / know an easy way to recode the gap statistic using high level distance functions?


title: Is this a valid approach to testing a hypothesis about the relationship between two variables?

text: I am trying to test a hypothesis I have about consolidation in the real-world market for a certain machine. (My apologies in advance for obfuscating a bit here, but some of the data is proprietary and may not be redistributed.) 1) The market has two components: new machines sold and existing machines in use. Companies act as agents who sell new machines and service the old machines. Machines can be used for 10-15 years before becoming unfit for use. The pool of machines in use has risen steadily over the years. 2) The data show that over the decades the number of agents has fallen steadily, although over the same period the number of machines in use has risen steadily. If you divide the number of machines by the number of agents to get a "agents per machines in use" value, this has risen steadily over time. In effect, the industry has become more efficient: it sells and services more machines with fewer (but presumably larger) agents. 3) My hypothesis is that the driver of this increased efficiency has been the rise in the ratio of machines in use to new machines sold. My theory is that the increase in machines in use offered an opportunity to some agents in the form of increased service revenues. Those opportunities required a broader set of skills and more organisation than simply selling new machines, so some agents out-competed others. Over time - decades, in fact - the stronger agents drove the weaker ones out of the industry or consolidated them. This led to an increase in overall efficiency. I want to find out whether there is indeed a relationship between market maturity and market efficiency. My approach to testing this hypothesis has been as follows. (a) Calculate the ratio of machines in use to new machine sales as a proxy for market maturity; the higher the number, the more mature the market. (b) Calculate the ratio of machines to the number of dealers operating in the market as a proxy for market efficiency. This gave me 45 yearly data points for each of the two data series. I then made a scatter plot in Excel, with the machines per agent in use on the x axis and the ratio of machines in use to new machines sold on the y axis. (I planned to post the chart here but I don't yet have enough reputation to do include images in posts.) The linear trend line equation of the plot is y = 0.7277x + 6.528 with an R2 of 0.7488. This looks like a fairly strong relationship but not suspiciously high. However, I have no formal training in statistics, so I'm concerned that I may be making a gross methodological mistake. Q. My question is... is this approach valid or irredeemably flawed? I have tried to normalise the data somewhat but I am still comparing a ratio and an absolute figure. Is there something else I should do? (Another approach I considered was simply comparing the machines in use per dealer to the number of machines in the market, but the scales are different by a couple of orders of magnitude. I got a 45-degree line and an R2 of 0.99.) I use Stack Overflow (mostly for R) but this is my first time on CrossValidated. I realise this is a question that may have a discussion element to it rather than the preferred "answerable question". Happy to edit the question as required.


title: What is a difference between random effects-, fixed effects- and marginal model?

text: I am trying to expand my knowledge of statistics. I come from a physical sciences background with a "recipe based" approach to statistical testing, where we say _is it continuous, is it normally distributed -- OLS regression_. In my reading I have come across the terms: random effects model, fixed effects model, marginal model. My questions are:   * In very simple terms, what are they?    * What are the differences between them?    * Are any of them synonyms?    * Where do the traditional tests like OLS regression, ANOVA and ANCOVA fall in this classification?  Just trying to decide where to go next with the self study.


title: Should pruning be avoided for bagging (with decision trees)?

text: I came by several posts and papers claiming that pruning trees in a "bagging" ensemble of trees is not needed (see 1). However, is it necessarily (or at least in some known cases) damaging to perform pruning (say, with the OOB sample) on the individual trees in an ensemble? Thanks!


title: Can a factor be changed to binomial levels to achieve model validation and extract insignificant variables?

text: I have just run a NBGLM and want to know something. If I am aiming to drop the least significant explanatory variables until all explanatory variables are significantly correlated wih the response variable (i.e. model selection finding the strongest variables). What do I do when there is a factor within the model that has a level within it insignificant?                                    Estimate Std. Error z value Pr(>|z|)              (Intercept)             6.3081079  1.6497143   3.824 0.000131 ***               Height_            0.0247727  0.0081857   3.026 0.002475 **                Width             -0.0017481  0.0006614  -2.643 0.008215 **           Upper_Field_Layer      -0.0257683  0.0119503  -2.156 0.031061 *            MeanMin                -0.2361341  0.1420208  -1.663 0.096378 .            as.factor(Site_Treat)2 -0.7361044  0.1866798  -3.943 8.04e-05 ***          as.factor(Change)2     -0.5002666  0.1810585  -2.763 0.005727 **           as.factor(Change)3     -0.1910862  0.1821139  -1.049 0.294055              ---      0.29 in change 3 is the most obvious to drop but as it is a factor it means dropping all other levels also. Could I reformat the excel spreadsheet to look like this? Where C=Change. A row for change 1, change 2 & change 3. Would these then be okay to put in as discrete nominal variables whereby on insignificant levels can be removed?               C   C1  C2  C3     1   1   0   0     2   0   1   0     3   0   0   1     1   1   0   0     2   0   1   0     3   0   0   1     1   1   0   0     2   0   1   0     3   0   0   1     1   1   0   0     2   0   1   0     3   0   0   1      Best Wishes


title: Appropriate regression model when dependent variable is between 0 and 1?

text: I am performing a regression where my dependent variable is the value of a group's Simpson Diversity Index. This index value is constrained by $1/k$ and $1$ (where $k$ is the number of classes), though none of my values approach $1$. I know OLS regression is not suited for regression with a 'bounded' dependent variable, and my research on the appropriate method has pointed me in several directions, to include a logit transformation and a beta regression. Beta regression is well over my head, so I am considering the logit transformation, but still am looking for some advice on interpreting the resulting coefficients, and if this method is truly sufficient. Additionally, some other questions: Do I just transform the dependent variable and leave the independent variables alone? Do I transform both? (By the way, my dependent variables include percentages, integers, and dummy variables.) With the transformation, I have read that OLS would then be appropriate, but I have also seen suggestions for GLM.


title: Appropriate test for ordered response variable

text: As part of my research project, I'm studying the plumage of a species of bird called Brent Goose. The species can be divided into two subspecies, a dark- bellied form and a light-bellied form. The two subspecies are present side by side at a location in Northern England, whereas in Southern England, only dark-bellied is found. Ecological character displacement predicts that populations that winter side by side should be more different than populations that winter further from each other. Therefore, I want to make two comparisons:   1. Northern England dark-bellied vs Northern England light-bellied   2. Southern England dark-bellied vs Northern England light-bellied If ecological character displacement theory holds, the difference between comparison 1 should be larger than comparison 2. To test for this, I have scored 4 characters on a 3 point scale. As an example, I have scored the upper parts as follows: 1=blackish, 2=dark greyish-brown, 3=pale greyish-brown. I have two questions:   1. Using my ordered character scores, what would be an appropriate test to determine if the difference in comparison 1 is larger than comparison 2?   2. Do I need to analyse each of the 4 characters separately or can they be combined into a single test?


title: Distance metric for all categorical data

text: I have categorical data in the form of Y/N of 23 attributes and containing 200+ records. I need to compute the distance between the attributes with their class label. Which metric can be used?


title: Multiple regression problem S Shape Residuals

text: I have posted this question before but wanted to give a more thorough explanation of my problem so as to hopefully garner some assistance. I am also posting it on Talk Stats under the same title. I have a mutiple regression model with over 10,000 records and though the model has an R^2 of 91%, the lack of fit tests and slightly S-shaped normal probability plot of the residuals show suggest curvature in the data. My objective here is to figure out and correct/ transform the curvature so as to fit the fit appropriately. I have tried every method I have researched that could possibly help, though it is possible that I still missed something. Here is what I have done: 1) Added quadratic and cubic terms of each x. This only seemed to exacerbate the curavture. I also tried only added these terms for one x at a time and the results were the same. 2) Used the natural log of y and all x's. This again exacerbated the S shaped residuals plot. 3) I have 3 categorical variables and thinking that variance between them was an issue, I added weights for each category by 1/each category's variance. This also yielded no benefit to the residual plot or fit of the model. UPDATED PER MICHELLE'S REQUEST: \--I am trying to predict sales volume by the amount to which our price is more or less than our main compeitor. I am in a high volume commoditized industry so I expect to be able to get some good results. \--My response is sales volume (y) and my predictors are the 4 week average of the sales volume for that same day and half-hour period (x1) and the amount we are above or below (in cents) our main competitor (x2). Also the categorical variable is store and I have modeled both using this as a variable and without. \-- Here are my model results: The regression equation is               Gallons = 4.19 - 293 MUSAPriceMinusLowKey + 0.983 MUSA4wkgal          Predictor                 Coef   SE Coef       T      P     Constant                4.1932    0.7771    5.40  0.000     MUSAPriceMinusLowKey   -292.75     13.07  -22.41  0.000     MUSA4wkgal            0.982941  0.002580  380.92  0.000          S = 47.1163   R-Sq = 91.0%   R-Sq(adj) = 91.0%          Analysis of Variance     Source             DF         SS         MS         F      P     Regression          2  323496065  161748033  72861.37  0.000     Residual Error  14352   31860609       2220     Total           14354  355356674          Source                DF     Seq SS     MUSAPriceMinusLowKey   1    1380837     MUSA4wkgal             1  322115228          Source                DF   Seq SS      PriceMinusLowKey      1   177122      4wkgal                1  8611828          Lack of fit test          Possible curvature in variable MUSAPric  (P-Value = 0.000 )     Possible interaction in variable MUSA4wkg  (P-Value = 0.000 )          Overall lack of fit test is significant at P = 0.000      So I am at my wit's end now and am not sure what to do next. Is it possible that the S shape residual plot is common for a dataset with so many observations and it is not a concern? Thanks in advance for any help you can give.


title: What does the hidden layer in a neural network compute?

text: I'm sure many people will respond with links to 'let me google that for you', so I want to say that I've tried to figure this out so please forgive my lack of understanding here, but I cannot figure out how the practical implementation of a neural network actually works. I understand the input layer and how to normalize the data, I also understand the bias unit, but when it comes to the hidden layer, what the actual computation is in that layer, and how it maps to the output is just a little foggy. I've seen diagrams with question marks in the hidden layer, boolean functions like AND/OR/XOR, activation functions, and input nodes that map to all of the hidden units and input nodes that map to only a few hidden units each and so I just have a few questions on the practical aspect. Of course, a simple explanation of the entire neural network process like you would explain to a child, would be awesome. What computations are done in the hidden layer? How are those computations mapped to the output layer? How does the ouput layer work? De-normalizing the data from the hidden layer? Why are some layers in the input layer connected to the hidden layer and some are not?


title: Data collection and storage for time series analysis

text: I am designing a data capture method for a client for inplay sporting events and he wants to record the odds movements for later analysis in Excel once every half second. I want to get this right so that it's easy to use the data down the line for analysis in other packages. A bit more background and assumptions.   * Each event can have between 4 - 40 contenders (c)   * Each event has 10 variables that apply equally to all contenders (e)   * Each contender has 20 variables of same heading/type with values unique to contender (i) In essence I need to choose between   * **1\. Having 1 timeframe on 1 row, so each timeframe capture has** Columns required = e+max(c)i = 810 Rows required = 1 Good: Easy to manipulate, data on one row, 1 row describes all contenders in event per row. Bad: Huge number of columns, lots of blank column data if c is less than max(c), hard to search names across multiple columns or   * **2\. Having 1 timeframe on multiple rows, so each timeframe has** Columns required = e+i =30 Rows required = c Good: Less columns, easy to search/filter as all names in the same column Bad: Timeframes in different rows for different contenders Does it matter? Is it easy for packages to handle data in both forms? My client doesn't know the answer but wants the best solution! I'm tending towards 2. as it's much easier to manage and search in database terms but not sure about preparation for time series analysis? Can anyone one with experience offer some advice? Thanks Os


title: What is missing from this taxonomy  of methods used in pattern recognition?

text: I'm working up a taxonomy showing different methods used in pattern recognition and I'd be curious to hear about how it could be improved. The Mind Map groups different methods based on the discipline which influenced their development. ![taxonomy](http://bentham.k2.t.u-tokyo.ac.jp/media/zoo.png)


title: Confidence intervals which contain 1 but the z test and p value are significant

text: I calculated the confidence intervals for the regression coefficients using the formula: SE X 1.96. I then took this value away from the regression coefficient to get the lower confidence level and added it to the regression coefficient to get the upper confidence level. However, for those z tests and associated p values which are significant I am getting confidence intervals which include 1. For instance:               Predictor variable 1:          Coefficient   = 1.08     SE            = 0.47     Z test        = 2.3     p value       = 2.2e-02          Confidence interval calculation:      SE X 1.96   = 0.47 X 1.96 = 0.92      1.08 - 0.92 = 0.16     1.08 + 0.92 = 2          Lower and upper confidence levels = (0.16, 2)      I'm just wondering why this may be the case? Thanks so much!


title: How to calculate Newey West t-statistic in SAS 9.2?

text: I need to run Newey West t statistics in SAS 9.2. I already run regression, White's test, Breusch Godfrey test and Jarque-Bera normality test. Regression is simple. Number of observations = 522 Depended variable name: rtest 1 Independent variable name: rtest 2 Data are time series data. I found somewhere that i should start with:               PROC model;     PARMS B0 B1;     rtest1=B0+B1*rtest2;      I do not understand this part:               FIT rtest1/GMM Kernel=(BART,1,0.605);     RUN;      How to get numbers in parenthesis. What is meaning of those numbers? What does command `FIT` do to data? Do I need to add more commands?


title: How do I use the GPML package for multi dimensional input?

text: I have downloaded the Gaussian Processes for Machine Learning (GPML) package (gpml-matlab-v3.1-2010-09-27.zip) from the website, and I can run the regression example (demoRegression) in Octave. It works just fine. Now I have my own data for regression where the x (input) matrix is a 54x10 matrix (54 samples, 10 input vars), and the y (target) vector is 54x1. The problem is that I do not understand hov to calculate the meanfunc and the covfunc, the code provided for the regression example does not work for multiple input datasets. I do not understand enough Octave to decode the specific code used to calculate this. Are there anybody who has tried this, and maybe can show an example?


title: Grouping samples by clustering or PCA

text: If I have 5 binary variables with values for 100 observations to give me a 5x100 matrix.                   y1 y2 .... y100     x1  0  1  .... 1     x2  0  0  .... 0     x3  1  0  .... 0     x4  1  1  .... 1     x5  0  0  .... 0      I want to put the observations into fairly homogeneous groups based on their values from x1 to x5. I can try PCA. However, failing any obvious groupings from PCA, what's the best way to group my samples? Should I just have a go at the various means of clustering? How do I assess the homogeneity of the resultant groups? Sorry if my questions are a bit vague, so another question may be: What do I need to find out about my data in order to answer the questions above?


title: Can I use paired sample t-test for the following problem?

text: Suppose there is an experiment that in each trial returns a two vector of probabilities of size $n$. In each trial I select uniformly a number from ${1,\cdots, n}$ and compare that element in this two vectors. I want to see if in a large number of trials, that randomly selected element of first vector is greater than the second one. For this I am planning to compare their mean values according to paired sample $t$-test. I know that the vector elements have specified variance and expected value. With this information is it justified to use this hypothesis test? Thank you.


title: t-test equivalent for samples, given a distribution of unknown form?

text: If there are two samples from a distribution whose form is not known, and which is certainly not normal, how can one ascertain the probability that, by random chance, their means/medians are dissimilar? In other words, is there a t-test equivalent for samples coming from a distribution of unknown form?


title: Comparing cosine similarities for tf-idf vectors for documents with different length

text: I'm computing cosine similarities between 2 vectors. These vectors are information retrieval query and document representations respectively. They have been computed using tf-idf weights. Since my documents have different length, tf-idf weights are theoretically unbounded. The question is: is cosine similarity still a valid measure? Can I compare several cosine similarities for each doc?


title: Detecting a consistent pattern in a dataset via Decision Trees and cross-validation

text: Assume a classification problem where there are two classes and the aim is to detect a consistent pattern which successfully separates the input dataset regardless of how we divide it into training / testing. The broader question is does such a consistent pattern exist across different partitions of the dataset at hand ? Let's assume that the classification is done by decision trees (e.g. J48). Also, suppose that a k-fold cross-validation is done to evaluate how consistent the pattern is across different combinations of training / testing. In each fold the testing sets contains completely unrelated and never before seen samples by the training phase. Although some folds may yield the same tree, it's very likely that the majority of folds will yield a different decision tree, in other words a different pattern. I would like to ask the following questions :   1. If each fold reports a different decision tree but gives good classification accuracy, does this mean there is no consistent pattern ? Off course there may be a lot of ways to divide a dataset, not all of them covered by a k-fold validation.   2. Cross-validation aims to assess the performance of a classifier. But when each fold yields a completely different decision tree, does it make sense to assess the performance only by looking at the classification accuracy AUC, TP/FP rates ? Is the diversity/variability of models (usually pruning trees is responsible for that behavior) a sign of a poor classifier or an indication of a complex dataset ?


title: How to project video viewcount based on historicals?

text: My goal is to create a formula that can give an indication of how a YouTube channel's video will perform in the first 30 days of its lifespan and eliminate viral video / "lightening in a bottle" outliers that may be on the channel. The goal is to use the resulting number to price a video from a specific YouTube channel. As an example, a hypothetical YouTube channel uploads approx 10 videos a month. Variables:   1. Some videos get shared more and are more "viral"   2. Videos have a "fat head" and "long tail." Fat head refers to the largest chunk of viewership which in the case of established YouTube channels happens upfront, and the long tail refers to accumulated views over succeeding months. These viewcounts belong to videos that the same channel uploaded in the last 30 days (from most recent in descending order):                 351,170         770,783       1,183,166         154,645       1,568,569       2,564,857       1,023,498       1,409,113       1,006,203       1,244,092       So my questions:   1. Is there a formula I could plug in to my spreadsheet given this data that could accurately come up with a conservative estimation of how the video will perform in the first 30 days?    2. If not, how can I create one?    3. Because some of these videos are still generating a "fat head" (like the most recently published video with 351,170 views) would it make sense to instead gather and average videos uploaded in the last 30-60 days instead? (fat head has time to impact viewcount and settle)


title: Estimation of school effects using xtmixed (in STATA 12)

text: I am trying to estimate the school effects on student scores. I ran a model with student effects as a random effect and school effects as fixed effects as follows: xtreg studentscore lagstudentscore limitedEng poor `gradeyeardum'`schoolyeardum', re where student score is regressed on lag of student score, student's language proficiency and poverty status, and a set of grade-year dummy variables and a set of school-year dummy variables. The coefficients on the school-year dummy variables are the key fixed effects of interest. I want to estimate an equivalent model using xtmixed. I ran the code below but the MLE doesn't seem to converge - I got the "not concave" message fofor several iterations. Here's the code I ran: xtmixed studentscore lagstudentscore limitedEng poor `gradeyeardum' schoolyeardum', /// || schnamyearid : || studentid: I am not sure if I should estimate, both the individual level and school level standard deviations for the variance component. I am interested at the school level variation, so perhaps the following model makes more sense? xtmixed studentscore lagstudentscore limitedEng poor `gradeyeardum' schoolyeardum', /// || schnamyearid : I would appreciate any help I can get. Thanks, Robbie


title: Poisson regression

text: im analyzing an article for my studies with the hypothesis if a change in work motivation ist related with a change in mental well being (http://www.sciencedirect.com/science...01879113001541). Sadly i dont know much about poisson regression. The follow up measurement was 18 month later. Do you consider always the time, when you do a poisson regression?Im not quite sure if they did so in this study.. If i imagine a graph of this regression, what can i see on the x and what on the y axis? Thanks for your help


title: How much smaller can $p$ values from ANOVA's $F$-test be vs. those from multiple $t$-tests on the same data?

text: **Intro:** Having noted the attention received today by this question, "Can ANOVA be significant when none of the pairwise t-tests is?," I thought I might be able to reframe it in an interesting way that would deserve its own set of answers. A variety of incongruous results (at face value) can occur when statistical significance is understood as a simple dichotomy and judged on the mere basis of which is higher, the $p$ or the $\alpha$. @Glen_b's answer to the above question presents a useful example of a case where:   * An ANOVA $F$-test produces a $p_F<.05$ for one independent variable (IV) with four levels, but    * $p_t>.08$ for all two-sample $t$-tests that compare differences in the same dependent variable (DV) among observations corresponding to each pair of the IV's four levels.  A similar case arose despite Bonferroni corrections for post-hoc pairwise comparisons via this question: Anova repeated measures is significant, but all the multiple comparisons with Bonferroni correction are not? Previously mentioned cases with a slightly different test in multiple regression also exist:   * F and t statistics in a regression: $p_F<.001,p_{\beta t}>.09$   * How can a regression be significant yet all predictors be non-significant?     * In @whuber's answer, $p_F=.0003,p_{\beta t}>.09$ **I bet** that in cases like these, **some** (but not all) pairwise comparisons' (or regression coefficients' significance tests') **$p$ values must be fairly close to $\alpha$ if a corresponding omnibus test can achieve a $p <\alpha$**. I see this is the case in @Glen_b's first example, where $F_{(3,20)}=3.19$, $p_F=.046$, and the largest pairwise difference gives the smallest $p_t=.054$. Must this be the case in general? More **specifically** : * * * **Question:** If an ANOVA $F$-test produces a $p_F=.05$ for one polytomous IV's effect on a continuous DV, **how high could the lowest $p$ value be** among all two-sample $t$-tests that compare each pair of the IV's levels? Could the minimum pairwise significance be as high as $p_t=.50$? * * * **I welcome answers that address only this specific question**. However, to further motivate this question, I'll elaborate and throw in some potentially rhetorical questions. Feel welcome to address these concerns as well, and even to ignore the specific question if you like, especially if the specific question gets a definitive answer. **Significance:** Consider how much less important the difference between a $p_F=.04$ and a $p_t=.06$ would be if statistical significance were judged in continuous terms of the strength of evidence against the null hypothesis (Ron Fisher's approach, I think?), rather than in dichotomous terms as above or below an $\alpha=.05$ threshold for acceptable probability of error in choosing whether to reject the null wholesale. "$p$-hacking" is a known problem that partly owes its notoriety to an unnecessary vulnerability introduced by interpretation of $p$ values according to the common practice of dichotomizing significance into the equivalents of "good enough" and "not good enough." If one were to dispose this practice and focus instead on interpreting $p$ values as strength of evidence against the null on a continuous interval, might omnibus testing be somewhat less important when one really cares about multiple pairwise comparisons? Not useless necessarily, as any reasonably efficient improvement in statistical accuracy is of course desirable, but...if, for instance, the lowest pairwise comparison's $p$ value is necessarily within $.10$ of the ANOVA (or other omnibus test) $p$ value, doesn't this make the omnibus test somewhat more trivial, less compulsory, and even more misleading (in conjunction with preexisting misunderstandings), especially if one doesn't particularly want to control $\alpha$ across multiple tests? Conversely, if data may exist such that an omnibus $p=.05$, but all pairwise $p>.50$, shouldn't this further motivate omnibus and contrast testing throughout practice and pedagogy? It seems to me that this issue should also inform the relative merits of judging statistical significance according to a dichotomy vs. a continuum, in that the dichotomous interpretive system should be more sensitive to small adjustments when differences are "marginally significant", whereas neither system is safe from a failure to perform an omnibus test or adjust for multiple comparisons if this difference / adjustment can be very large (e.g., $p_t-p_F>.40)$ in theory. Other optional complexities to **consider or ignore—whatever makes answering easier and more worthwhile** :   * How high $p$s for $t$s could be if, for $F$, $p<.05$ instead (e.g., $p=.01, .001,\dots$)   * Sensitivity to number of levels in a polytomous IV   * Sensitivity to unevenness in significance of pairwise differences (while all $p_t>p_F$)     * whuber's answer indicates that including small differences can mask big differences.   * Differences among various omnibus tests' corrections for multiple comparisons     * See also: Correcting for multiple comparisons in a within subjects / repeated measures ANOVA; excessively conservative?     * With multiple IVs, it seems multicollinearity can exacerbate this issue.   * Restricted cases where data meet all assumptions of classic parametric tests optimally     * This restriction may be important to prevent this question from being somewhat moot.


title: Residuals for logistic regression and Cook's distance

text: 1. Are there any particular assumptions regarding the errors for logistic regression such as the constant variance of the error terms and the normality of the residuals?   2. Also typically when you have points that have a Cook's distance larger than 4/n, do you remove them? If you do remove them, how can you tell if the model with the removed points is better?


title: Why does -2*LL differ when using binary logistic regression vs GLM binary logistic in SPSS?

text: I'm comparing plausible models selected a-priori to predict a binary response variable. I used binary logistic regression in SPSS20 and obtained AIC=-2*LogLikelihood+2k where k is the number of parameters in the model. Using the GLM procedure (binary logistic) I get provided with AIC, AICc, BIC etc. straight away, however, the output differs from the binary logistic regression output (and results in different ranking of the models). Why do the outputs using binary logistic regression vs GLM (binary logistic) differ? And which approach should I use? Thanks heaps!!


title: If $\sigma$ follows Inverse-Gamma, can I just square the result to get $\sigma^2$ if I do simulation?

text: Suppose I have $\sigma|Y,Y*,\psi$~$Inverse-Gamma(A, B)$. How should I get $\sigma^2$. Can I just square the result from Inverse-Gamma? To me, this seem to work, if you think of it this way: $\sigma$ is the standard deivation and $\sigma^2$ is the variance, and this relationship is one-to-one. You are just squaring the number. But, sometimes the distribution changes, say $X\tilde{\ } N(0,1)$, but $X^2\tilde{\ } \chi^2(k)$, so, I don't know, and I am very confused. Thanks for your help!


title: Intuitive understanding covariance, cross-covariance, auto-/cross-correliation and power spectrum density

text: I'm currently studying for my finals in basic statistics for my ECE bachelor. While I think I have the math mostly down, I lack intuitive understanding what the numbers actually mean.(Preamble: I'll use rather sloppy language). I know the E[X] is the "weighted average" of all results of X weighted by their probability. Var[X] gives gives the expected variance from that E[X] squared, so tells us something about the "blurriness" of the distribution. The other properties I know the formulae, but lack any kind of intuition. Anyone have good explanations/resources to help with that?


title: Ways to measure distance from multivariate Gaussian (Mahalanobis distance)

text: I have a cluster of p-dimensional points and given a new p-dimensional point $x$ I want to determine whether or not it is likely to belong to this cluster. The cluster is made up of $n$ p-dimensional points, I am making the assumption that these points are drawn from a multivariate Gaussian distribution with sample mean ${\hat{\mu_X}}$ and sample covariance matrix $\hat{\sigma_X}$. Given a new point $x$ I am trying to decide if it is likely to belong to this cluster using the following theshold on the Mahalanobis distance: $$\frac{n}{(n-1)^2}\left(X_i-\hat{\mu}_X\right)'\hat{\sigma}_X^{-1}\left(X_i-\hat{\mu}_X\right)>B_{0.95}\left(\frac{p}{2},\frac{n-p-1}{2}\right)$$ However, some clusters have very few sample points $n$ in which case calculating the inverse covariance matrix $\hat{\sigma_X}^{-1}$ becomes impossible. Are there any other equivalent or more appropriate measure I can use in this case? Thanks


title: What insights can be found by using leverage plots?

text: I'm trying to figure out whether leverage plots can provide valuable information. See example here.


title: CAP and ROC measures

text: Why is the area under the curve measure better than the accuracy ratio? I know that the AR lies between 0,5 and 1 and the relationship $AUC=\frac{1}{2}(AR+1)$ holds, so the AUC lies between 0,75 and 1, but this does not really seem to be an advantage.


title: What is the expected range out of regression analyses?

text: I am computing multiple regression coefficients for these variables (for several cases):                        y,x1,x2,x3          log(y)= a(log(x1))+b(log(x2))+c(x3)+d          res=lm.fit(log(y)~log(x1)+log(x2)+(x3))      All values of `y,x1,x2,x3` are between `0 and 1`. The values of the coefficients range between `-10 to 10` for `a` The values of the coefficients range between `-5 to 5` for `b` The values of the coefficients range between `-5` to `5` for `c` The values of the coefficients range between `-3` to `3` for `d` I wonder if we should expect a certain range for each coefficient.or in other words,how to know that my regression analyses did a good job.


title: How to forecast time series from clusters?

text: I have a set of 1200 time series for which I need to build a forecasting model. Since it is too time consuming to build a model for each series, I separated the series into clusters and what I want to do now is to build a forecasting model for each cluster (using time series techniques). However, I don't know which is the best way to build that model for each cluster. Should I extract a series prototype and build the model on it? How do I do this? In a paper by E. Maharaj (1999) (Forecasting Time Series From Clusters) the author forecasts time series from clusters by "pooling" the series together, but I can't understand how this pooling is done.


title: How to estimate baseline hazard function in Cox Model with R

text: I need to estimate baseline hazard function $\lambda_0(t)$ in a time dependent Cox model $\lambda(t) = \lambda_0(t) \exp(Z(t)'\beta)$ While I took Survival course, I remember that the direct derivative of cumulative hazard function ($\lambda_0(t) dt = d\Lambda_0(t)$) would not be a good estimator because Breslow estimator gives a step function. So, is there any function in R that I could use directly ? Or any reference on this topic ? I am not sure if it is worth to open another question, so I just add some background why baseline hazard function is important for me. The formula below estimates the probability that the survival time for one subject is larger than another,. Under a Cox model setting, baseline hazard function $\lambda_0(t)$ is required. $P(T_1 > T_2 ) = - \int_0^\infty S_1(t) dS_2(t) = - \int_0^\infty S_1(t)S_2(t)\lambda_2(t)dt $


title: Choice of distance metric when data is combination text/numeric/categorical

text: I have a large table of attributes of different real-world movie theaters. I have classified them by the "true" physical entity to which they belong, so that there may be multiple records for a given movie theater entity. In this table, I have information such as names, how many screens they have, etc. Given some identifying information (part of the name, one line of the address and number of screens, say) I would like to classify the given information to the entity with which it is associated, and add it to the database. I was thinking of using an algorithm such as nearest neighbour but the choice of distance metric seems limiting. The only implementations I have seen use all numeric or all text information to calculate distance. How would I calculate a distance metric for data that may be numeric, text and categorical in nature?


title: Understanding Simpson's paradox

text: Andrew Gelman in one of his recent blog posts says: >   1. I do not think counterfactuals or potential outcomes are necessary for > Simpson’s paradox. I say this because one can set up Simpson’s paradox with > variables that cannot be manipulated, or for which manipulations are not > directly of interest. > >   2. Simpson’s paradox is part of a more general issue that regression coefs > change if you add more predictors, the flipping of sign is not really > necessary. > > > > Here’s an example that I use in my teaching that illustrates both points: > > I can run a regression predicting income from sex and height. I find that > the coef of sex is \$10,000 (i.e., comparing a man and woman of the same > height, on average the man will make \$10,000 more) and the coefficient of > height is \$500 (i.e., comparing two men or two women of different heights, > on average the taller person will make \$500 more per inch of height). > > How can I interpret these coefs? I feel that the coef of height is easy to > interpret (it’s easy to imagine comparing two people of the same sex with > different heights), indeed it would seem somehow “wrong” to regress on > height _without_ controlling for sex, as much of the raw difference between > short and tall people can be “explained” by being differences between men > and women. But the coef of sex in the above model seems very difficult to > interpret: why compare a man and a woman who are both 66 inches tall, for > example? That would be a comparison of a short man with a tall woman. All > this reasoning seems vaguely causal but I don’t think it makes sense to > think about it using potential outcomes. I pondered over it (and even commented on the post) and think there's something that begs to be understood with greater clarity here. Until the part on interpretation of gender it is so okay. But I do not see what's the problem behind comparing a short man and a tall woman. Here's my point: In fact it makes even greater sense (given the assumption that men are taller on average). You cannot compare a 'short man' and a 'short' woman for exactly the same reason, that the difference in income is explained in some part by the difference in heights. Same goes for tall men and tall women and even more so for short women and tall men (which is further out of the question, so to speak). So basically the effect of height is eliminated only in the case when short men and tall women are compared (and this helps in interpreting the coefficient on gender). Doesn't it ring a bell on similar underlying concepts behind the popular matching models? The idea behind Simpson's paradox is that the population effect might be different from the sub-group wise effect(s). This is in some sense related to his point 2 and the fact that he acknowledges that height should not be controlled for alone (what we say omitted variable bias). But I could not relate this to the controversy on the coefficient on gender. Maybe you might be able to express it more clearly? Or comment on my understanding?