title: What is a user-friendly and powerful software package for econometrics modeling?

text: I would be using such software to run multiple regressions using macroeconomics variables as independent variables to estimate other macroeconomic ones as single dependent variables. I need to be able to run robust regressions including Weighted Least Squares(to resolve heteroskedasticity), Feasible Generalized Least Squares (to resolve autocorrelation of residuals). If using a regular regression, bypassing such models, would need to be able to run heteroskedastic resistant and/or autocorrelation resistant Standard Errors. Probably could also use having access to ARCH and GARCH models. I am not a computer programmer, coder, Visual Basic expert. I am just a pretty good power user of regular Excel functions and its Add-ins programs. Thus, I am interested in something very user friendly with a menu or window interface that would not be too difficult to pick up. Because of my ease of use requirements, I have no interest in SAS. Preliminarily, talking to a few economists they seem to gravitate towards EViews. I also hear good stuff about Stata. I sense EViews is a bit easier to use. But, Stata may be a bit more powerful. I hear SPSS is pretty good. But, I sense for econometrics models it probably would come as a distant second vs either EViews or Stata. I also hear a bit about Microfit. Do you have any experience with those programs? Can you give me some insights regarding their respective trade offs? What would you recommend? Is there another program I should consider?


title: Why do we use GLM?

text: I'm trying to justify using of GLM model in my project instead of a simple linear regression. A lot of sources that I've seen contain the statement that "GLM allow us to build regression models for response variables that are not normally distributed". I can't understand what that is supposed to mean since we don't make any assumptions about variable distribution. Could you, please, explain why it's often more appropriate (theoretically) to use GLM if your response variable is, say, count.


title: Intuition behind Dvoretzky Kiefer Wolfowitz inequality

text: I've been reading about the Dvoretzky-Kiefer-Wolfowitz inequality, in the context of confidence bands on empirical distribution functions. I think I understand the inequality at face value: that the probability that the empirical distribution function differs from the real distribution function by more than some amount is bounded by a function of sample size. (And a surprisingly aggressive function of sample size at that, at least to my mind). What I'm completely lacking, though, is any sense of _why_ the inequality holds. Given we're not (I think) making any assumptions about the real distribution function, apart from that it is a distribution function, I don't really see why we can say anything useful about it at all! It probably doesn't help that I've yet to find a comprehensible proof of the theorem (my definition of comprehensible is that which might be understood by a competent physicist!) So, can anyone point me towards a readable proof? Or perhaps supply a hint as to the intuition behind this inequality?


title: "the leading minor of order 1 is not positive definite" error using 2l.norm in mice

text: I am having a problem using the `2l.norm` method of multilevel imputation in `mice`. Unfortunately I cannot post a reproducible example because of the size of my data - when I reduce the size, the problem vanishes. For a particular variable, `mice` produces the following errors and warnings:               Error in chol.default(inv.sigma2[class] * X.SS[[class]] + inv.psi) :        the leading minor of order 1 is not positive definite     In addition: Warning messages:     1: In rgamma(n.class, n.g/2 + 1/(2 * theta), scale = 2 * theta/(ss *  :       NAs produced     2: In rgamma(1, n.class/(2 * theta) + 1, scale = 2 * theta * H/n.class) :       NAs produced     3: In rgamma(1, n.class/2 - 1, scale = 2/(n.class * (sigma2.0/H - log(sigma2.0) +  :       NAs produced      If I use the `2l.pan`, `norm` or `pmm` methods, the problem does not occur. The variable has the following distribution: ![enter image description here](http://i.stack.imgur.com/brJBH.png)                  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's         50.0   117.0   136.0   136.7   155.0   249.0    3124       Also, the class sizes have the following distribution:                  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.         3.00   50.00   80.00   88.52  111.00  350.00


title: Hypothesis test on data with confounding spatial clustering

text: This is a bit of an elaboration on a question I posted earlier, since I feel like my approach to the problem as a whole is probably quite flawed. Suppose I have a set of treatment and control cells, each with a position in space and some response value. I would like to know whether the treatment cells are significantly more similar in their responses than controls. However, the analysis is complicated by the fact that:   * Across all cells there is some tendency for the responses of pairs of nearby cells to be more similar   * My treatment cells are also significantly clustered in space So far, I have tried to find an upper bound on the probability of the null hypothesis using a bootstrap test:   1. For each possible pair of treatment cells, find a group of matched control pairs whose spatial distance is less than or equal to that of the treatment pair. The number of pairs that satisfy this criterion will vary depending on which treatment pair I'm considering.   2. Find the rank of each treatment pair within its distribution of matched control pairs and normalize it to between 0 and 1   3. Take the mean of the normalized ranks across groups as my 'real' score.   4. Bootstrap a null distribution by drawing randomly from each set of possible normalized ranks and taking the mean.   5. The normalized rank of the 'real' score within this null distribution gives me my p-value. Since I only consider control pairs whose distance is <= that of the treatment pair, this ought to be a very conservative test. However, if I were to consider control pairs whose distance is the same as the treatment pair to within +/- some tolerance, I will end up choosing control pairs that are, on average, more distant than the treatment pair (because the treatment pairs are spatially clustered there will tend to be a greater number of control pairs that are more distant than the treatment pair). Is there a better way of doing this analysis? I'm sure I must be missing something very obvious!


title: What is an unbiased estimate of population R-square?

text: I am interested in getting an unbiased estimate of $R^2$ in a multiple linear regression. On reflection, I can think of two different values that an unbiased estimate of $R^2$ might be trying to match.   1. **Out of sample $R^2$:** the r-square that would be obtained if the regression equation obtained from the sample (i.e., $\hat{\beta}$) were applied to an infinite amount of data external to the sample but from the same data generating process.   2. **Population $R^2$:** The r-square that would be obtained if an infinite sample were obtained and the model fitted to that infinite sample (i.e., $\beta$) or alternatively just the R-square implied by the known data generating process. I understand that adjusted $R^2$ is designed to compensate for the overfitting observed in sample $R^2$. Nonetheless, it's not clear whether adjusted $R^2$ is actually an unbiased estimate of $R^2$, and if it is an unbiased estimate, which of the above two definitions of $R^2$ it is aiming to estimate. Thus, my questions:   * What is an unbiased estimate of what I call above out of sample $R^2$?   * What is an unbiased estimate of what I call above population $R^2$?   * Are there any references that provide simulation or other proof of the unbiasedness?


title: Ratio of correlated sample variances (gamma distributed)

text: for $N$ samples of two correlated random variables $X \sim N\left(0,\sigma_X^2\right)$ and $Y \sim N\left(0, \sigma_Y^2\right)$ with correlation $\rho$, I am analyzing the ratio of the sample variances, $r=\frac{s_Y^2}{s_X^2}$. Currently, I am foccusing on $E\left[r\right] = E\left[\frac{s_Y^2}{s_X^2}\right]$. I have found that $s_X^2 \sim \Gamma\left(\frac{N-1}{2}, \frac{2 \sigma_X^2}{N-1}\right)$, $s_Y^2 \sim \Gamma\left(\frac{N-1}{2}, \frac{2 \sigma_Y^2}{N-1}\right)$ since $\frac{(N-1)S_X^2}{\sigma_X^2}\sim \chi_{N-1}^2$ (and likewise for $s_Y^2$). The ratio of the sample variances is therefore the ratio of two dependent gamma distributions. I have found some papers on this topic, for instance [1, 2]. These papers do however always assume the gamma distributions having different shape and same scale parameters (scale 1 to be precise) - just the opposite of my problem. I have tried using $\frac{N-1}{2 \sigma_X^2} s_X^2 \sim \Gamma \left(\frac{N-1}{2}, 1\right)$, $\frac{N-1}{2 \sigma_Y^2} s_Y^2 \sim \Gamma \left(\frac{N-1}{2}, 1\right)$ - the ratio of these should be beta prime distributed.   Using the formula for the expected value of a beta prime distribution this leads to   $E\left[\frac{s_Y^2}{s_X^2}\right] =E\left[\frac{\frac{N-1}{2 \sigma_X^2}}{\frac{N-1}{2 \sigma_Y^2}} \frac{\frac{N-1}{2 \sigma_Y^2}}{\frac{N-1}{2 \sigma_X^2}} \frac{s_Y^2}{s_X^2}\right] =\frac{\sigma_Y^2}{\sigma_X^2} E\left[\frac{\frac{N-1}{2 \sigma_Y^2}}{\frac{N-1}{2 \sigma_X^2}} \frac{s_Y^2}{s_X^2}\right] =\frac{\sigma_Y^2}{\sigma_X^2} \frac{\alpha}{\beta - 1}$ with $\alpha = \beta = \frac{N-1}{2}$,   which doesn't seem to be the correct result (shouldn't the ratio in some way depend on the correlation?) and which differs from the results of numerical experiments. How can the expected value of the ratio be calculated? At which point am I mistaken?   I would very much appreciate any help - thank you in advance. Sources:   [1] Lee et al.: Distribution of a ratio of correlated gamma random variables. SIAM Journal on Applied Mathematics 36.2 (1979): 304-320.   [2] Tubbs: Moments for a ratio of correlated gamma variates. Communications in Statistics-Theory and Methods 15.1 (1986): 251-259.


title: Multiple Imputation: Before or after case exclusion?

text: I've done quite some reading on multiple imputation (MI), but can not seem to figure out the next question: I have a dataset with missing values, some rows have many missing values, others have less missing values. They are missing at random. I want to exclude certain rows because there are missing too many values (columns) , though I'm not sure if I should first create a dataset only with the rows I want to analyze in the end, or use the complete dataset to fill in the missing values (with MI) and delete the rows afterwards. Does the second method help me create thrust worthier imputations?


title: temperature prediction algorithm

text: I found an interesting problem in a contest on temperature prediction: https://www.hackerrank.com/contests/expansion- challenge/challenges/temperature-predictions It is not about forecasting the temperature in the future, but estimate the missing data. I guess 1D interpolation may work. Also one can collect all the data from Jan, Feb, ....Dec and do an regression. What is the best algorithm for this kind of time-series problem (estimate the missing value as precise as possible)?


title: Naive Bayes: Imbalanced Dataset in Real-time Scenario

text: I am using scikit-learn Multinomial Naive Bayes classifier for binary text classification (classifier tells me whether the document belongs to the category X or not). I use a balanced dataset to train my model and a balanced test set to test it and the results are very promising. This classifer needs to run in real time and constantly analyze documents thrown at it randomly. However, when I run my classifier in production, the number of false positives is very high and therefore I end up with a very low precision. The reason is simple: there are many more negative samples that the classifer encounters in the real-time scenario (around 90 % of the time) and this does not correspond to the ideal balanced dataset I used for testing and training. Is there a way I can simulate this real-time case during training or are there any tricks that I can use (including pre-processing on the documents to see if they are suitable for the classifer)? I was planning to train my classifier using an imbalanced dataset with the same proportions as I have in real-time case but I am afraid that might bias Naive Bayes towards the negative class and lose the recall I have on the positive class. Any advice is appreciated.


title: Outlier problems

text: Having built a regression model with an ordinal response variable and predictors comprised of categorical and continuous nature, I have some questions that pertain to one of the final goals, i.e. suggest another class for certain observations (which ones is unknown) in the training set. I approach this problem with a sort of outlier analysis methodology:   1. Is there any way to know if the labels of certain observations in the training set were misplaced to begin with?   2. Even if I do find some observations that are likely "outliers", based on model's studentized residuals OR Mahalanobis Distance etc, how do I know that it is not due to model bias? Since the model is trained on mislabeled data.   3. How do I go about choosing another class label for these records? Thank you for any thoughts !


title: Statistics published in academic papers

text: I read a lot of evolutionary/ecological academic papers, sometimes with the specific aim of seeing how statistics are being used 'in the real world' outside of the textbook. I normally take the statistics in papers as gospel and use the papers to help in my statistical learning. After all, if a paper has taken years to write and has gone through rigorous peer review, then surely the statistics are going to be rock solid? But in the past few days, I've questioned my assumption, and wondered how often the statistical analysis published in academic papers is suspect? In particular, it might be expected that those in fields such as ecology and evolution have spent less time learning statistics and more time learning their fields. How often do people find suspect statistics in academic papers?


title: Confusion related to data normalization

text: I am trying to learn a linear regression model. However, I have some confusion related to the normalization of the data. I have normalized the features/predictors to zero mean and unit variance. Do I need to do the same for the target. If so why?


title: Why is it called white noise?

text: I know, that the white noise is called "white noise", because it comes from physics and has something to do with the spectral decompoistion (is that right?) of the white light? I am not familiar with signal theory, so could anyone explain me why it is called white noise, especially the relation to the terms light, spectral and so?


title: How to find out the best approximation of a probability distribution?

text: I have a sample from an unknown distribution. I checked with KS test the hypotheses that the sample represents lognormal, Weibull or Rayleigh distributions and the test shows that all these hypotheses cannot be accepted. Are there any means to find the distribution that approximates my sample in the best way?


title: Looking for Intuition Regarding Fixed/Random Effect Conundrum

text: I understand that we use random effects (or mixed effects) models when we believe that some model parameter(s) vary randomly across some grouping factor. I have a desire to fit a model where the response has been normalized and centered (not perfectly, but pretty close) across a grouping factor, but an independent variable `x` has not been adjusted in any way. This led me to the following test (using _fabricated_ data) to ensure that I'd find the effect I was looking for if it was indeed there. I ran one _mixed_ effects model with a random intercept (across groups defined by `f`) and a second _fixed_ effect model with the factor f as a fixed effect predictor. I used the R package `lmer` for the mixed effect model, and the base function `lm()` for the fixed effect model. Following is the data and the results. Notice that `y`, regardless of group, varies around 0. And that `x` varies consistently with `y` within group, but varies much more across groups than `y`               > data           y   x f     1  -0.5   2 1     2   0.0   3 1     3   0.5   4 1     4  -0.6  -4 2     5   0.0  -3 2     6   0.6  -2 2     7  -0.2  13 3     8   0.1  14 3     9   0.4  15 3     10 -0.5 -15 4     11 -0.1 -14 4     12  0.4 -13 4      If you're interested in working with the data, here is `dput()` output:               data<-structure(list(y = c(-0.5, 0, 0.5, -0.6, 0, 0.6, -0.2, 0.1, 0.4,      -0.5, -0.1, 0.4), x = c(2, 3, 4, -4, -3, -2, 13, 14, 15, -15,      -14, -13), f = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L,      4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor")),      .Names = c("y","x","f"), row.names = c(NA, -12L), class = "data.frame")      Fitting the mixed effects model:               > summary(lmer(y~ x + (1|f),data=data))     Linear mixed model fit by REML      Formula: y ~ x + (1 | f)         Data: data         AIC   BIC logLik deviance REMLdev      28.59 30.53  -10.3       11   20.59     Random effects:      Groups   Name        Variance Std.Dev.      f        (Intercept) 0.00000  0.00000       Residual             0.17567  0.41913      Number of obs: 12, groups: f, 4          Fixed effects:                 Estimate Std. Error t value     (Intercept) 0.008333   0.120992   0.069     x           0.008643   0.011912   0.726          Correlation of Fixed Effects:       (Intr)     x 0.000       I note that the intercept variance component is estimated 0, and importantly to me, `x` is not a significant predictor of `y`. Next I fit the fixed effect model with `f` as a predictor instead of a grouping factor for a random intercept:               > summary(lm(y~ x + f,data=data))          Call:     lm(formula = y ~ x + f, data = data)          Residuals:          Min       1Q   Median       3Q      Max      -0.16250 -0.03438  0.00000  0.03125  0.16250           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) -1.38750    0.14099  -9.841 2.38e-05 ***     x            0.46250    0.04128  11.205 1.01e-05 ***     f2           2.77500    0.26538  10.457 1.59e-05 ***     f3          -4.98750    0.46396 -10.750 1.33e-05 ***     f4           7.79583    0.70817  11.008 1.13e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.1168 on 7 degrees of freedom     Multiple R-squared: 0.9484, Adjusted R-squared: 0.9189      F-statistic: 32.16 on 4 and 7 DF,  p-value: 0.0001348       Now I notice that, as expected, `x` is a significant predictor of `y`. **What I am looking for** is intuition regarding this difference. In what way is my thinking wrong here? Why do I incorrectly expect to find a significant parameter for `x` in both of these models but only actually see it in the fixed effect model?


title: confusion related to the dual of svm

text: I have a confusion related to the dual of svm In the main objective function I have ![enter image description here](http://i.stack.imgur.com/qNxtL.png) Now to solve the dual of this objective function, I will minimize with respect to the primal variables first to get ![enter image description here](http://i.stack.imgur.com/VAEDx.png) I didn't get how substituting the value of w in the lagrangian will give me $-\frac{1}{2}\sum_{i,j}\alpha_i\alpha_jy^{(i)}y^{(j)}{x^{(i)}}^Tx^{(j)}$ Where did the minus sign come from in the above? I got the other terms by substituting the remaining terms. But I didn't get the above term how it came to give me ![enter image description here](http://i.stack.imgur.com/BWOOV.png)


title: How was this summation performed?

text: Say I have the following simple Bayesian network involving 3 r.v.s A, B, and C: $$ A \rightarrow C \rightarrow B $$ I am trying to prove that A and B are conditionally independent given by marginalizing out C and seeing if I get out $p(a, b) = p(a)p(b)$ in $$ p(a, b) = p(a) \sum_{c} p(b|c)p(c|a) $$ However, I am not sure how to simplify $\sum_{c} p(b|c)p(c|a)$. How do the $c$s go away?


title: Creating "demo" data from real data: disguising without disfiguring

text: _(I have no real idea what to tag this with because I'm no statistician and I don't know what field this falls into. Feel free to add more suitable tags.)_ I work for a company that produces data analysis software, and we need a decent set of data to test and demo our latest product with. We can't just fill the database with the output of a random number generator because the program's outputs would become nonsensical. One of the simplest ways to get such data is from a client; we have a large body of data from a trial we ran. Now, obviously we can't publish a client's actual data, so we need to alter it a bit, but we still need it to behave like real data. The aim here is to take their set of data, and apply a "fuzz" to it so that it can't be recognised as specifically theirs. My memory of statistical theory is itself a little fuzzy, so I'd like to run this by you guys: Essentially, the data we have (from the client) is itself a sample of all the data that exists (in the country, or the world). What I'd like to know is what type of operations can be applied to make the sample no longer strongly representative of the client's sample population, while still keeping it roughly representative of the world's population. For reference, as far as we're aware the data we have generally follows rough normal (Gaussian) distributions. The original dataset isn't widely available, but could theoretically be recognised from some regionally-specific characteristics (we don't know what those characteristics are, and it's doubtful whether _anyone_ does to a sufficient level, but we know that variations exist from place to place). Anyways, I'm more interested in the theory of this than the practice - I want to know whether an operation makes it impossible (or at least difficult) to identify the source dataset by parameter X, whether or not anyone has or could work out parameter X in the first place. The approach I've come up with is to separate the readings into the various types, (without giving much away, let's say a group might be "length" or "time taken to do X".) For each of those, calculate the standard deviation. Then, to each value, add a random value between the positive and negative values of (n * stddev) where n is some fraction that I can use to tune the result until the data is sufficiently "fuzzed". I didn't want to simply apply a static range (say, random between 90% and 110% of the original value) because some values vary much more or less than others - in some measurements, being 10% over the mean is barely noticeable, but in others it makes you a serious outlier. Is this sufficient to mask the source of the original data? If not, by which statistical measures would the data still be identifiable, and how would I mask those while still keeping the resultant data vaguely realistic?


title: Which statistical test to use?

text: I have acquired data from 3 different donors in which 6 different drugs were tested over a period of time. For a period of 1 week I have measured 1 parameter each day from each condition from each donor. I would like to compare two drugs with one and other. But when I plot the parameter values over time of each condition there is a height difference between the lines of the different donors (the shape of the line is the same and this for every condition). When I combine the 3 donors in 1 line the standard deviations are very large. Which statistical test can I use to compare 2 conditions and correct for the donor variability? thanks a lot


title: Cross Validation using NON-independent Data

text: i got a problem with cross-validation! I got two studies, t(N=1200) and t+1(N=200) respectively. The 200 corporations form t+1 are also part of the t0-Data. In t+1, there are two significant relationships using kendall's tau, but they are unexpected relationships. In order to check for errors, i want to do a cross validation using the whole data from t. But there's a problem: I can remove only 165 (t+1)-corporations from t, because 35 corporations cannot be identified. Thus, i cannot postulate independence between results and cross-validation. Is there a way around this? Thanks in advance!


title: Why do (a lot) fewer than 95% of my confidence intervals intersect the true parameter value in this simulation?

text: Why are my confidence intervals performing so poorly? I have a data frame `may` of two variables $C$ and $F$. I perform a linear regression               lm.1 <- lm(may$C ~ may$F)     lm.1$coef[2]         #  may$F     #  0.0161821       confint(lm.1)[c(2,4)]     # [1] 0.0151  0.0173      So far so good. I then want to simulate running a sample and getting the confidence interval from each sample, similar to this illustration from wikipedia. However when I do -- expecting to get 95% of my intervals to intersect the true parameter (`0.0161`) -- I find that only 55% (!!) of my intervals intersect the correct parameter.![100 confidence intervals with true parameter](http://i.stack.imgur.com/aLiFh.png)               sample.size <- 15          plot(c(1,n), c(0, 0.03), type="n")       lines(x=c(-1,n+1), y=c(0.016,0.016), lty=2, col=2)       x <- cbind(rep(0,n), rep(0,n))     for (i in 1:n){       may.temp <- may       lm.temp <- NA       set.seed(i*3)       may.temp <- may.temp[sample(1:nrow(may.temp), sample.size), ]       lm.temp <- lm(may.temp$C ~ may.temp$F)            x[i, 1] <- confint(lm.temp)[2,1]       x[i, 2] <- confint(lm.temp)[2,2]            lines(x=c(i,i), y=c(x[i,1], x[i,2]), col="blue", lwd=2)     }     y <- rep(0,n)     for (i in 1:n){       if( x[i,1] < 0.017 & x[i,2] > 0.015){       y[i] <- 1     }     }     sum(y)     # 55      I had two thoughts:   1. My sample size (15) is too low. But I thought the confidence intervals would simple become larger to take this into account.   2. My data has outliers. My data does have outliers, but (from my non-expert eye) these don't seem to bias the model. Where have I gone wrong?


title: Calculating parameters of a binomial distribution

text: Say I have a set of data which I know comes from a Binomial distribution (generated using the `rbinom` function in `R`). Is there any way to use the data to work backwards and work out the `p` and `n` parameters of the Binomial distribution from which it came? EDIT: I realise that there are infinite combinations: `n = 100` and `p = 0.1` will produce approximately the same distribution as `n = 10000` and `p = 0.0001`. What I'm looking for is an approximation - any combination will do.


title: R Time Series Issues

text: I would like to know more details about R issues mentioned in Time Series Analysis and Its Applications: With R Examples. For e.g. the first problem still exists in R version 3.0.1               # generate an AR(1) with mean 50     set.seed(66)      # so you can reproduce these results     x = arima.sim(list(order=c(1,0,0), ar=.9), n=100) + 50        mean(x)         [1] 50.60668   # the sample mean is close     arima(x, order = c(1, 0, 0))         Coefficients:                ar1  intercept  <--  here's the problem             0.8971    50.6304  <--  or here, one of these has to change       s.e.  0.0409     0.8365      *direct copy from the webpage Theoretical model of above simulation is: $X_t = 5 + 0.9 X_{t-1} + Z_t$, where $\{ Z_t\} $ is the white noise Is R following particular model convention or is R wrong? If it is wrong, it is still in R core package and I wonder why it isn't changed. Also, it would be very helpful if I get some hint on other five issues.


title: Non-linear regression references

text: I'm looking for non-linear regression analysis references. I'll be grateful if somebody can give me some references.


title: Showing that an apparent continuous variable is actually a 2-level categorical variable

text: I have some survey data and employees are divided into grades (from 1 to 7). At first glance, this appears to be a continuous variable. However, based on my knowledge of the organization, I know that levels 5 and above are managers, while those in grades 1 to 4 are lower level employees. From this variable, I would like to create a dummy varaible to use as an independent variable in a regression model. (I have done a number of t tests on various variables and shown that there are statistically significant differences between the two groups.) My question is this: what statistical test (if any) can I use to show that there are actually two different groups rather than 7?


title: Dynamic factor analysis vs state space model

text: The MARSS package in R offers function for dynamic factor analysis. In this package, the dynamic factor model is written as a special form of state space model and they assume the common trends follow AR(1) process. As I am not very familiar with those two methods, I come with two questions: Is the Dynamic Factor Analysis a special form of State Space Model? What is the difference between those two methods? In addition, the Dynamic Factor Analysis does not necessary assume the common trends as AR(1) process. Is there any package that allows the the common trends as seasonal ARIMA (or some other) process?


title: Iterating through test/train datasets

text: I am building a predictive model. I split 2/3 of the data into a training set, and the rest into the test set. I build a model based on the training set. I use the the model to test on the training set, and the test set. What are the advantages and disadvantages to iterating through all 3 combinations of test/train datasets with my data to get more results on how well my model performs on the data? Is there any reason to just test the model on this one partitioning of the data into test/train sets?


title: A priori selection of SVM class weights

text: I remember seeing/reading somewhere that for multiclass SVMs with unbalanced data, there was a way to determine the class weights from the training data (rather than X validation). Does anyone know what the method is or what paper its from? Thanks


title: Forecasting the target variable vs building a causal model and forecasting causal variables

text: I want to know the approaches people use to forecast lets say unemployment rate .... By itself it might not fit a time series model (ARMA) very well as the trend is dependent on many external factors. My hypothesis is that if a time series for the target doesn't fit well, but some of the underlying causal variables can be reasonably forecasted into the future it is better to use the forecasted causal variables and predict the target variable. Is that a common approach? I don't know about ARIMA but does it take into account external factors?..


title: Sample size for clinical trials, when effect size is predicted

text: How should I calculate sample size for a clinical trial when I need a power of 80% and I predict an effect size to be at lest 0.5? do I need to estimate the mean of outcome as well?


title: Numerical or computational method to generate density estimate

text: Suppose $X \sim U(0,1)$, and $Y[x] = g(x)$ where $g(.)$ is some complicated function. I want to calculate/plot the density of $Y$. I can do this analytically for simple enough $g$. I can also generate some large number $N$ samples of $X$ from the stardard uniform, calculate $y_i = Y(x_i)$ for each sample $x_i$, and then compute the histogram of these $ y_i, \forall i \in (1, ... N)$. This will turn out to be inefficient if, say the details of $p_Y(y)$ come from very small regions of $X$. Is there a better way to calculate/estimate $p_Y(y)$?


title: Hand calculated confidence intervals after logit regression, symmetric or not?

text: I've used Stata 12 to estimate a hierarchical logit model (xtmelogit) with random intercepts. The outcome is incarceration (1=yes) for a series of convicted criminal offenders. There are a slew of interaction terms and I'm trying to understand these by looking at the predicted odds at different levels of X while holding other values as observed. The margins command won't work here since it cannot include the random effects portion of the model. So I'm doing it by hand. After setting my values where I want them I calculate the fixed estimate (the logodds for the fixed portion of the model) and add that to the estimated random intercept to get a point estimate The code looks something like this: `predict randeffect, reffects` to fit with a random intercept and `predict randse, reses` to get the standard error for the random intercept The following code is done within a `forvalues` loop where I vary the value of some covariates and then estimate the fixed effect at each level... what I've posted is a simplification that demonstrates the principle without the complexity of the `forvalues` loop `predict fixedeffect, xb` which gets the point estimate for the fixed portion of the model (log odds) and `predict fixedse, stdp` which gets the standard error for the fixed portion of the model (log odds). To get the point estimate at each level I add the fixed effect to the `randeffect` and take the exponent of it to get the predicted odds. That is, `gen odds = exp(randeffect + fixedeffect)`. **And now the tricky part...** To get the confidence interval my first instinct was to do all the addition with the log odds and then take the exponent to get the upper and lower bounds. This results in very asymmetrical bounds because the exponent is (obviously) not a linear transformation of the values. The code is as follows: `gen lowerbound = exp((randeffect + fixedeffect) - (1.96 * (fixedse + randse)))` The above makes sense but the confidence intervals are not symmetrical. However, I see in the output of the logit command that all odds ratios have symmetrical confidence intervals. Alternatively, I can calculate the confidence interval and then take the exponent and then add and subtract that to the calculated odds ratio. This gives me nice a nice symmetrical confidence band around my point estimate when I graph it. This code looks like `gen lowerbound = odds - (exp(1.96 * randeffect + fixedeffect))` So which way is right? Do I exponentiate first and then add and subtract? Or do I add and subtract and then exponentiate? Thanks for any insight you can provide. EDIT: I wanted to follow up with a caveat for anyone considering doing the same in Stata. First, Stata can make predictions based on the fixed part of the model which basically ignores the intercept. Since the fixed effects are the same for all clusters, it's like estimating a model where all level 1 units are in the same cluster. It's like an apples to apples comparison. I've chosen this approach to understand the fixed effects and then separately consider variation in the intercepts across my level 2 units. It's much simpler this way. Second, from what I can tell the method described in my original post does not correctly model interaction effects which were the focus of the analysis. The problem, I think, is that when you vary the value for X2 to see how the effect of X1 on Y changes, you must also change/reset the value of X2*X1 (the interaction term). This sounds simple enough but the problem quickly becomes a thorny mess of code when your X1 variable is involved in more than one interaction (i.e. you have X2*X1 as well as X3*X1 in the model). So this is just something to be watchful for if you, after reading this, get the bright idea to take this route. I personally gave up on it, preferring instead to use the margins command with the predict (mu fixedonly) option               margins x1, dydx(x2) predict (mu fixedonly)      But I sincerely thank those of you who contributed to the effort!


title: Inter-rater reliability of time measurements

text: I am having some trouble trying to decide what measure of inter-rater reliability to use in a study. Part of a larger study involves accurately determining when participants began (onset) and stopped (offset) writing. I have video footage which captures the writing, and can determine the time of onset/offset to the nearest 0.01 second. The time is taken from a stopwatch which was running continuously from the start of each experiment, with multiple onset/offsets in each experiment. The onset/offset have been defined, and I would like to determine the inter-rater reliability for 2 raters. The data looks like this:               Participant A  trial one - onset: 21.35 offset: 24.55                    trial two - onset: 54.50 offset 57.75     Participant B  trial one - onset 10.35 offset 21.15                    trial two - onset 35.65 offset 39.90      and so on and so forth. Two raters would have data like the above and I would like to calculate if they are reliably similar. I looked at Kappa statistic, Krippendorff’s alpha, intra-class correlation, etc., but none seems appropriate for continuous time measurements.


title: What are 'aliased coefficients'?

text: While building a regression model in R (`lm`), I am frequently getting this message               "there are aliased coefficients in the model"      What exactly does it mean? Also, due to this `predict()` is also giving a warning. Though it's just a warning, I want to know how can we detect/remove aliased coefficients before building a model. Also, what are the probable consequences of neglecting this warning?


title: Estimating ratios of quantile means for sorted data

text: In studying the relationship between income and consumption, it is common to sort the observations by income and observe that high-income households have lower levels of consumption per dollar of income than low-income households, i.e., they have higher savings rates. This conforms to one's naive intuitions about the relationship, whether you think the higher income is the cause or consequence of the higher savings rate. I recently took a commonly used consumption data set (the U.S. Survey of Consumer Expenditure), and sorted it into deciles based on aggregate consumption. After this sorting, the decile with the highest consumption had the lowest savings rate. In fact, they dissaved by a considerable margin. I am looking at some real reasons why this might be true in the period in question (e.g., drawing down fictitious wealth from the housing bubble), but it seems to me that I might get the same result if income and consumption were both observed with error, that is, if the errors were not too highly correlated. Suppose that this is true, and that the errors (or the percentage errors) of household income and consumption are uncorrelated with the true errors. How might I estimate that relationship correctly? And more specifically, can I estimate it correctly groupwise, for, e.g., income deciles or consumption deciles or quintiles? (Here I use "deciles" to refer to the groups of sample observations divided by the decile break-points, rather than by the breakpoints themselves).


title: Combining expert estimates

text: I am trying to work on a process to improve how how well my team estimates. I want to look at using some statistics to help out and embrace the uncertainty in how we estimate tasks. If I have a group of experts and I ask them for a best (b), average (a) and worst (w) case estimate for a task, creating a triangular distribution. Using a monte carlo method I can generate 10,000 trials using the triangular distribution for each of the estimators and combine these together into one sample set. I can then use a MLE to create a beta distribution from the sample statistics, then use the cdf(x) functions to determine the likelihood of an estimated number of days, or use the inverseF to determine the number of days to maintain a certain likelihood. This seems to make sense to me logically, however I am unsure about the fitting of a Beta distribution to the samples. I am wondering if there is a better distribution I should be using, especially given that the samples are just randomly generated from N triangular distributions. Does anybody have any suggestions or experience in doing this that could point me in the right direction? It's been a few years since I left uni, so my statistics is a little rusty. There are plenty of examples of this process for one expert, however I would like to be able to do the same thing for multiple experts. Thank you, Aidan


title: What GLM family and link function for "proportion of time"?

text: A simple question to which I don't seem to find the answer anywhere. I have a response variable `duration of time spent doing A` of individuals tested for $\text{max duration}=X$. Therefore the variance is likely to be much wider at $x=0.5$ compared to $X=0$ and $X=1$. However, I am confused what family and link function I need as time is continuous but can not be less than $0$ and more than $X$. ps: I am using R so any code examples would be appreciated **Update**   I have heard about the tobit() function from "AER" which " _is a convenience interface to 'survreg' setting different defaults and providing a more convenient interface for specification of the censoring information._ " Therefore, it seems for my data I could run the model `tobit(duration.A ~ factor1*factor2, left=0, right=X)` However, it is unclear how A) I could do this for a mixed model i.e. a model with random factors, B) what it's assumptions are (they are not clear from the R documentation http://cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf), and C) how I could get residuals and fitted values for plotting.


title: What is this kind of tabular relation diagram called?

text: A relation between two sets can be visualized with a grid-like table in which the rows are labelled along the left with the elements of one set and the columns are labelled along the top with the elements of the other, and in which each cell of the table is occupied by a symbol (usually a dot) if & only if the elements for that row & column are related. What is this kind of diagram called? An example: ![Example](http://i.stack.imgur.com/UU3vV.png) (I made this image myself, but I'm positive I've seen these charts elsewhere before.)


title: Why does this excerpt say that unbiased estimation of standard deviation usually isn't relevant?

text: I was reading on the computation of the unbiased estimation of standard deviation and the source I read stated > (...) except in some important situations, the task has little relevance to > applications of statistics since its need is avoided by standard procedures, > such as the use of significance tests and confidence intervals, or by using > Bayesian analysis. I was wondering if anyone could elucidate the reasoning behind this statement, for example doesn't the confidence interval use the standard deviation as part of the calculation? Therefore, wouldn't the confidence intervals be affected by a biased standard deviation? EDIT: Thanks for the answers so far, but I'm not quite sure I follow some of the reasoning for them so I'll add a very simple example. The point is that if the source is correct, then then something is wrong from my conclusion to the example and I would like someone to point how how the p-value _doesn't_ depend on the standard deviation. Suppose a researcher wished to test whether the mean score of fifth graders on a test in his or her city differed from the national mean of 76 with a significance level of 0.05. The researcher randomly sampled the scores of 20 students. The sample mean was 80.85 with a sample standard deviation of 8.87. This means: t = (80.85-76)/(8.87/sqrt(20)) = 2.44. A t-table is then used to calculate that the two-tailed probability value of a t of 2.44 with 19 df is 0.025. This is below our significance level of 0.05 so we reject the null hypothesis. So in this example, wouldn't the p-value (and maybe your conclusion) change depending on how you estimated your sample standard deviation?


title: What is a "concordance score" for regression coefficients?

text: I came across this "concordance score" in a set of slides called Penalized regression methods for ranking variables by effect size, with applications to genetic mapping studies, by Ji Zhu: $$ \mathrm{CS} = \frac{ \sum_{j \neq k} \left( \mathbb{I}{\left[\hat{\beta}_j > \hat{\beta}_k\right]} + 0.5 \cdot \mathbb{I}{\left[\hat{\beta}_j = \hat{\beta}_k\right]} \right) \cdot \mathbb{I}{\left[\beta_j > \beta_k\right]} }{ \sum_{j \neq k} \mathbb{I}{\left[\beta_j > \beta_k\right]} } $$ where $\mathbb{I}$ is the indicator function. This on slide 7. The previous slide just says something like "consider $\operatorname{P}{\left(\hat{\beta}_{j\lambda}-\hat{\beta}_{k\lambda}>0\right)}$. A bit later, the slides suggest tuning $\lambda$ in penalized regression by maximizing CS, and they call the process "oracle tuning." What is the index $k$ referring to here, and what exactly is this measuring the concordance of? What kind of model does oracle tuning select?


title: Cross-Validation vs. AICc for LASSO

text: I was working on a research project in which I try to estimate the the individual contribution of a group of regional political leaders to local economic growth. The major challenge is that there is relatively few locality- year observations to the number of regional leaders. Overall I have about 300 political leaders and about 500 locality-year observations. The data matrix is also quite sparse as most regional leaders serve in one to three localities for an average of 3-4 years. I decided to use LASSO, which would grant me more degrees of freedom by throwing away the insignificant variables. Although I was told that AICc and cross-validation (and choosing the model that gives the minimum mean squared error) are theoretically equivalent, but in my dataset they perform very differently. LASSO gives much fewer non-zero individual coefficient estimates under AICc than under cross-validation (200 fold for a dataset wth 500 observations). So my question is: Why do I see different results under AICc and cross-validation? And which criteria should I be using in selecting the model and reporting the results? Any help on this will be greatly appreciated. And I'm happy to provide more details about my research if needed. Thanks!


title: Comparing two models in R

text: I'm using R to develop regression models, and I need to compare two different models' performance. The question that arises is, "Is Model 1 statistically better than model 2?" and I don't seem to have a way to answer that question. Background: Model 1 consists of Variable A regressed on the endpoint. Model 2 consists of Variables B, C, and D regressed on the endpoint. Both models are developed using lm - ordinary least squares, nothing too fancy here. Given that these are not nested models, I cannot compare them using ANOVA. I can look at the R2 of actual vs predicted for each model, and I see that Model 2 is better, but how do I determine if it is statistically significantly better? I've also used the Concordance Correlation Coefficient, but again, I can't find a way to prove significance. The best I've come up with is that the rho for Model 2 is better than Model 1, but that rho value is within the 95% confidence limits of the rho for Model 1. I should throw in there that my assessment of predicted vs actual has been on a 60 observation hold-out set (240 observations in the training set).


title: Residuals normal enough for LMM (R package nlme)?

text: My larger question is what might be the best analytical approach given the structure of my data (see data summary below) but presently, I'm wondering if the normality of the model residuals are normal enough to justify a linear mixed model or should I pursue a glmm? I did build lmm and glmm models and couldn't get perfectly normalized residuals in either. I'm not able to post images of my residuals but the lmm residuals are above the line for low and high values (chi-square skewed?). Is there a test within the nlme package? Do residuals need to be normal under a negative binomial distribution in glmmPQL? In the lmm model, the parameters were significant (intervals function), models were better with random intercepts (liklihood ratio test), and the plot of residual versus fitted values appears to meet the homogeneity of variance assumption with points symmetrical around zero. I've been pouring over Zuur et al. 2009 Mixed Effects Models and Extensions in Ecology with R. I did pursue unmarked (R), but had significant problems with the model output and I feel as though I can justify not accounting for detection probability given my research objectives. LMM:               full=lme(ABUNDHA.LOG~FAMILY+DAMCAT+RIVKM.SD+R1WT.SD+R3D.SD+R1WT.SD:R3D.SD+fSEASON,              method="ML",random=~1|REP/SITE_CODE/YEAR,data=data.all8)      GLMM:               glmmPQLnb1 <- glmmPQL(ABUND ~ FAMILY+RIVKM.SD+R1WT.SD+DAMCAT+fSEASON+AREAHA+offset(AREAHA),                            random = list(~1 | REP, ~1 | SITE_CODE, ~1 | YEAR),                            data = data.all8,                            family = negative.binomial(theta = 0.9111, link = log),                            na.action = na.exclude)      Data summary: The data set consists of surveys of river-associated birds, conducted from established sites (n = 80) of varying size (min = 0.3, max = 103.3, mean = 22.3, 1 s.d. = 24.3 ha) along the banks of multiple rivers (n = 11) in Maine. The number and species of birds were counted during the 20 minute survey within the set boundaries of the site. Rivers differ in the number of sites (min = 1, max = 38, mean = 7, s.d. = 10) and the number of surveys conducted at each site. Surveys (n = 4784) were conducted year-round approximately every two weeks in the spring and summer and every three weeks in the fall and winter, for four years (2009 – 2012). The variable “REP” was used to denote each survey at each site, 1….21 (for example site 1, could have REP = 2 in year = 2009, 2010, …). The response variable could be defined in two possible ways:   1. as count data with abundance of river-associated birds and size of survey site as an offset or,    2. as density (abundance of birds per hectare). River-associated bird species, could be grouped by “FAMILY” (GENERALIST, INSECTIVORE, PISCIVORE, SHOREBIRD, WATERFOWL) or by migration status, or by a single species model (BNAME = BAEA). There are 43 species under the BNAME variable and 7 general species categories (DUCK, FLYCATCHER, GULL, MERGANSER, SHOREBIRD, SWALLOW, and YELLOWLEGS). When a survey was conducted and no birds were observed, a value of “NOBI” was entered in the “BNAME” variable, and “NA” in the “FAMILY” and “MIGRAT” variables. There is also a “SEASON” (0-winter, 1-spring, 2-summer, 3-fall), “SITE_CODE” (each site has a unique number/code), “JULIAN” (julian date).  Explanatory variables include the following (all are site-specific except for WATERCAT (survey-specific):   1. DAMCAT (the position of a site above, below, or not at a dam) or   2. DAM01 (site does not (0) or does (1) contain a dam)   3. WATERCAT (survey was conducted during low (L), medium (M), or high (H) water levels).   4. RIVKM (distance from mouth of river to survey site in km)   5. TIDAL01 (survey site is not influenced by tide (0) or is (1))   6. R1WT (the percent of area within 100 m along the river bank that is wetland habitat)   7. R3D (the percent of area within 100 m along the river bank that is developed)


title: PACF and ACF for AR and MA

text: I once heard the following statement: > The PACF (partial autocorrelation function) for MA processes behaves much > like the ACF for AR processes; the PACF for AR processes behaves much like > the PACF for AR processes. How to understand the logic underlying this statement? In addition to the strict mathematical proof, are there any approaches to understand this statement, the inherent relationship of AR, MA along with their PACF/ACF, from time series properties, statistics, or any other high- level thoughts?


title: Limiting pmf as $n \to \infty$

text: Consider the simple pmf: $$p_n (x)=\begin{cases} 1\quad x=2+1/n \\ 0\quad \text{elsewhere} \end{cases}$$ Then my book states that $\lim_{n\to \infty} p_n (x)=0$ **for all values of $x$**. Is that really the case though? Why can we not say that the probability of $x=2$ equals $1$ as $n \to \infty$? Thank you.


title: Simple Average Method, calculating moving average and moving variance, how can I say if stationary or not?

text: Lets say I have a sample data (here is just 10 numbers, in real I have about 10000 measurement results). Then, I want to check if the data is stationary or not using Simple Average Method. For example, my data set of size `N = 10`:               X = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};      I calculated averages (window = 3):               SAM = {2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}      with this formula:               (1+2+3)/3 = 2     (2+3+4)/3 = 3     (3+4+5)/3 = 4     ...      and put them into SAM table above. Then, calculated the differences beetwen my averages, `SAM[i+1] = SAM[i]`, and I got a differences table: `1 1 1 1 1 1 1` from which I see that the difference between means (averages) is constant (its always 1). Can I assume that with this simple test my data `X = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};` is stationary?


title: Whether to use factor analysis based on binary multiple response data?

text: I have a survey where I have asked people which type of computer games they enjoy and whether they consider themselves a hardcore gamer. I allowed people to select multiple genres, but now I am unsure what to do with my data. I initially thought factor analysis, with the idea being if there were genres that belonged to a particular type, they would separate out and I would be able to see a pattern. However, since I have the data on whether they consider themselves hardcore or not, it seems like I should use it. If I did use it, would it make sense to test each genre individually, e.g. RTS-hardcore vs. RTS-non-hardcore, or should I look into combinations of genres? * * * [EDIT] @Srikant: Yes, I was planning to make the answers binary. It didn't make sense to me to score different answers with values depending on say the number of genres somebody chose, because they might not necessarily play each game evenly, and I would have no way to determine what ratio of gameplay each genre had, so binary seemed the most fair to me. @mbq & @chl: The aim is just to see if certain genres tend to be more "hardcore" than other genres. I was expecting to find RTS, MOBA and FPS to cluster towards hardcore seeing as they tend to have a higher learning curve than say music/rhythm games.


title: How to test whether a distribution follows a power law?

text: I have the data of how many users post how many questions. For example, [UserCount, QuestionCount] [2, 100] [9, 10] [3, 80] ... ... it means each of the 2 users posts 100 questions, each of the 9 users posts 10 questions, and so on. So, how to validate whether the UserCount-QuestionCount distribution follows a power law? I find the poweRlaw packages, e.g., http://cran.r-project.org/web/packages/poweRlaw/vignettes/examples.pdf. However, i can only pass one group of numbers to do the evaluation.(The provided exmaple of this package is word frequency.) So how do i use this package? Or i got something wrong? I also have the data of each user's question count, i.e., [100, 100, 10, 10, 10 ... ]. If i pass this data to the package, what i will get?


title: Logistic Regression: Transforming Variables

text: When transforming variables, do you have to use all of the same transformation? For example, can I pick and choose differently transformed variables, as in: Let, $x_1,x_2,x_3$ be age, length of employment, length of residence, and income.               Y = B1*sqrt(x1) + B2*-1/(x2) + B3*log(x3)      Or, must you be consistent with your transforms and use all of the same? As in:               Y = B1*log(x1) + B2*log(x2) + B3*log(x3)       My understanding is that the goal of transformation is to address the problem of normality. Looking at histograms of each variable we can see that they present very different distributions, which would lead me to believe that the transformations required are different on a variable by variable basis.               ## R Code     df <- read.spss(file="http://www.bertelsen.ca/R/logistic-regression.sav",                      use.value.labels=T, to.data.frame=T)     hist(df[1:7])       ![alt text](http://i.stack.imgur.com/5qPcd.jpg) Lastly, how valid is it to transform variables using $\log(x_n + 1)$ where $x_n$ has $0$ values? Does this transform need to be consistent across all variables or is it used adhoc even for those variables which do not include $0$'s?               ## R Code      plot(df[1:7])      ![alt text](http://i.stack.imgur.com/ShjOG.jpg)


title: Better Classification of default in logistic regression

text: Full Disclosure: This is homework. I've included a link to the dataset ( http://www.bertelsen.ca/R/logistic-regression.sav ) My goal is to maximize the prediction of loan defaulters in this data set. Every model that I have come up with so far, predicts >90% of non-defaulters, but <40% of defaulters making the classification efficiency overall ~80%. So, I wonder if there are interaction effects between the variables? Within a logistic regression, other than testing each possible combination is there a way to identify potential interaction effects? Or alternatively a way to boost the efficiency of classification of defaulters. I'm stuck, any recommendations would be helpful in your choice of words, R-code or SPSS syntax. My primary variables are outlined in the following histogram and scatterplot (with the exception of the dichotomous variable) A description of the primary variables:               age: Age in years     employ: Years with current employer     address: Years at current address     income: Household income in thousands     debtinc: Debt to income ratio (x100)     creddebt: Credit card debt in thousands     othdebt: Other debt in thousands     default: Previously defaulted (dichotomous, yes/no, 0/1)     ed: Level of education (No HS, HS, Some College, College, Post-grad)      Additional variables are just transformations of the above. I also tried converting a few of the continuous variables into categorical variables and implementing them in the model, no luck there. If you'd like to pop it into R, quickly, here it is:               ## R Code     df <- read.spss(file="http://www.bertelsen.ca/R/logistic-regression.sav", use.value.labels=T, to.data.frame=T)      ![alt text](http://i.stack.imgur.com/aVqtZ.jpg) ![alt text](http://i.stack.imgur.com/VQJDg.jpg)


title: Logistic Regression: Classification Tables a la SPSS in R

text: In SPSS output there is a pretty little classification table available when you perform a logistic regression, is the same possible with R? If so, how? ![alt text](http://i.stack.imgur.com/yFS6P.jpg)


title: Alternative to the Wilcoxon test when the distribution isn't continuous?

text: One of the assumptions for using the Wilcoxon sign-rank test is that the underlying distribution is continuous (see here.) However, there are cases (for example, when analyzing Likert scale data) where this assumption might not necessarily hold. In such cases, what test can you use? And how would you do it with R? (My only bet here is to use a randomization test on the median - which I imagine can be easily done using the boot package.)


title: Nonlinear regression

text: I have some functions of $x$, in the form of $d\sqrt{x}$ or $d\log(x)$ where $d$ is known. I would like to rewrite (approximate is fine) them in the form $a/(1 + bx^c)$, where $a$, $b$ and $c$ are arbitrary. I don't think there are $a$, $b$ and $c$ such that the two curves will match exactly, so I think maybe try to fit a nonlinear regression and find the closest $a$, $b$ and $c$. I tried the following R code, and I get the singular gradient error.               x     = seq(1, 50000, by=1000)     y     = 50*sqrt(x)     model = nls(y ~ a/(1 + b*x^c), start=list(a=1, b=-0.01, c=0.01))


title: Interpreting the relInterIntra {irr} results in R

text: I am using this for a reliability analysis test for quantitative values between two raters. I got my results using the "relInterIntra {irr}" package in R, but the reported values gave me several "rhohat" values. I am unsure which co-efficient the results represent. Could someone help enlighten me here? I am new to using R. or would it be advisable using the "icc" package instead for my reliability analysis? An example of the results I got using the package for a repeated trial, 2 rater test on 82 subjects is shown below:                   >relInterIntra(x=LCDScore_3,nrater=2,raterLabels=c('Rater1','Rater2'),rho0inter=0.6,rho0intra=0.8,conf.level=0.95)          Inter/Intrarater reliability          Subjects = 82     Raters = 2           rhohat = 0.774, 0.915, 0.882, 0.903, 0.904, 0.926, 0.891, 0.916      rhohat = 6.59, 2.10, 1.84, 2.46     rhohat = 1.11e-15, 4.59e-04, 3.17e-03, 3.42e-05     rhohat = 0.296, 0.845, 0.858, 0.840, 0.877     rhohat = 0.958, 1.024, 0.887, 1.561, 1.058


title: regularized loss function

text: I'm fitting a classifier with cross-entropy loss (i.e. Bernoulli likelihood). Some examples are very clearly associated with one class or the other, and despite some attempts at regularization, the classifier sometimes assigns probabilities $<10^{-9}$. While the true probability might actually be this small in some cases, it's still causing some headaches. Among other things, the gradient is basically vertical in this region, and I'm running into numerical errors. My question is this: in addition to regularizing the model's complexity, as I'm currently doing, is there anything wrong with also penalizing the model for making extreme predictions? I'd also be interested in knowing if there are published papers discussing this approach, and what kinds of penalties would make the most sense. It seems like the most natural approach would be a beta-distributed prior on the outcome variable (equivalent to a prior belief that the classifier shouldn't return extreme values), but I'm very open to alternatives.


title: Nonparametric ANOVA with interaction?

text: Is there any nonparametric version of ANOVA able to deal with multiple factors AND evaluate their interaction(s)? If the interactions are really impossible to evaluate, is there any generalization of the Friedman test to more than two factors?


title: categorical variables in regression analysis and interaction terms

text: I am learning multiple regression with categorical variables and in a book I came across this problem. For yield of corn suppose there are two factors affecting, nitrogen level and depth of ploughing. Say there are three nitrogen categories(1,2,3) and two depth categories(1,2). There's a interaction term of nitrogen*depth as well. If the introduced dummy variables are as   $E_{i1}=1 $ if ith observation has nitrogen level 1 ,0 otherwise.   $E_{i2}=1 $ if ith observation has nitrogen level 2 ,0 otherwise. And $D=1 $ if depth category is 1,0 otherwise . Then the model would be   $Y=\beta_0+\beta_1E_1+\beta_2E_2+\beta_3D+\beta_4[E_1.D]+\beta_5[E_2.D]+\epsilon$.   Is this correct? When the t-statistics are found for the category $[nitrogen=1*depth=1]$ significance value was 0.029 and for$ [nitrogen=2*depth=1]$ it was 0.290.   I was asked to interpret if the interaction term significantly affects the yield of corn. Under 5% confidence level clearly the coefficient of the variable $[nitrogen=2*depth=1]$ is not significant. But since the coefficient of the variable $[nitrogen=1*depth=1]$ is significant can I say that, **the interaction term significantly affects the yield of corn**.   If both had insignificant coefficients then I could have said that there isn't a significant effect from the interaction term right? Since the interaction term $[nitrogen=2*depth=1]$ had a insignificant coefficient then my fitted model would be $\hat Y=\hat\beta_0+\hat\beta_1E_1+\hat\beta_2E_2+\hat\beta_3D+\hat\beta_4[E_1.D]$ leaving out E2D.IS this correct?   Can someone please help me to figure this out.


title: Feature selection : how to select the Information Gain threshold?

text: I am trying to use **Information Gain** to select features when classifying text with a Support Vector Machine. For each word in our training data, we computed its information gain. Then, we should keep only the words with an IG superior to some threshold. When reading the litterature about this I did not find a clear explanation about how we should select this thresold. > Given a training corpus, for each unique term we computed the information > gain, and removed from the feature space those terms whose information gain > **was less than some predetermined threshold.** > > Source: A comparative study on feature selection in Text Categorization . > For each dataset we selected the subset of features **with non-zero > information gain**. > > Source: Information Gain, Correlation and Support Vector Machine When training our SVM with words having IG > 0 as suggested in this paper, the results are worse than without feature-selection. When using IG > 0.01 the results become better. **Is there a recommended way for selecting this IG threshold ?** Should we perform a grid-search (with Training and Cross Validation data set) to determine the best IG threshold ? In which case this threshold becomes a meta-parameter such as the C of the SVM model. Should we stick to the simple IG > 0 even if our generalization results are worse ?


title: Multinomial logistic regression low classification rate

text: I am running a multinomial logistic regression with SPSS and I have encountered a problem (?) with my data. I have a dependent variable (DV) with three categories, five independent variables (IV) as factors and four IVs as covariates. There are close to 4800 respondents in the data. The distribution of my DV is: First category (40,6 %), second category (28,1 %), third category (31,3 %). As I run the regression (last category as the reference), the classification table shows that only 1,3 % percent are classified in the second category. See below: ![enter image description here](http://i.stack.imgur.com/KGfmX.png) I have tried to find some solution to this dilemma, and the only thing I have come up is that the model is just not a good one considering the classification. The results of the logistic model otherwise seem "reliable", though. Does anybody have any idea of what might be wrong or what could be done alternatively to this problem? Any help would be much appreciated! Best regards


title: How to implement a 2-D Gaussian Processes Regression through GPML (MATLAB)?

text: I just touched Gaussian processes two weeks ago. I am not very familiar with the selection of a model and its hyperparameters. Here is the demo code that I run for a 2-D Gaussian processes regression. Its output is not what I expected.               % produce the training set for regression.     % Here, the regression target Y is the sum of input.     X1train = linspace(-4.5,4.5,10);     X2train = linspace(-4.5,4.5,10);     X = [X1train' X2train'];     Y = sum(X,2);          % produce the test set for regression     Xtest = [1 2];          % set the hyperparameters     covfunc = {@covMaterniso, 3};     ell = 1/4; sf = 1;     hyp.cov = log([ell; sf]);     likfunc = @likGauss;     sn = 0.1;     hyp.lik = log(sn);          % implement regression     [ymu ys2 fmu fs2] = gp(hyp, @infExact, [], covfunc, likfunc, X, Y, Xtest);      The result is: `ymu` = 0.131695275851991, but what I expected is `ymu` = 1 + 2 = 3


title: Assumptions of cluster analysis

text: Apologies for the rudimentary question, I am new to this form of analysis and have a very limited understanding of the principles so far. I was just wondering if many of the parametric assumptions for multivariate/univariate tests apply for Cluster analysis? Many of the sources of information I have read regarding cluster analysis fail to specify any assumptions. I am particularly interested in the assumption of independence of observations. My understanding is that violation of this assumption (in ANOVA and MAVOVA for example) is serious because it influences estimates of error. From my reading so far, it seems that cluster analysis is largely a descriptive technique (that only involves statistical inference in certain specified cases). Accordingly, are assumptions such as independence and normally distributed data required? Any recommendations of texts that discuss this issue would be greatly appreciated. Many Thanks.


title: Standard deviation larger than mean

text: I am trying to test if there has been a significant change in the weight of two groups, one without exercise and one with exercise. The two samples have large spread, which creates a standard deviation that is larger than the mean for each group. My question is, is it still valid to use the unpaired t-test to compare means when the standard deviation is larger than the mean. If not, what statistical test can I perform?


title: Working out probability with alpha and beta

text: Let's say in an experiment, the Null hypothesis is patients having Condition $\alpha$ (type I error) is $0.05$ and $\beta$ (type II error) is $0.15$. If I have a sample of 50 patients, 35 with condition A and 15 with condition B, how would I work out the following:   1. Probability of properly predicting the condition of all patients with condition A   2. Probability of properly predicting the condition of all patients with condition B   3. Probability of predicting every condition properly   4. Probability of incorrectly predicting at least 5 condition Bs. I'm really confused on this, but here's what I think:   1. Since the error rate is $0.05$, properly predicting it would be a chance of 0.95. So 0.95^number of Condition A = $0.95^{35}$. Is that right?   2. Reverse of 1, so $0.85^{15}$.   3. Multiplying 1 & 2   4. $0.15^{10}$. Or am I doing it totally wrong and it should be worked out as a binomial problem?


title: Sum of Variance or Standard Error?

text: My end goal is to present the average travel time for a corridor and the standard deviation (or standard error) of that travel time. I'm not sure what to present to be statistically correct, here is the information I have to summarize: I have the average travel time for a segment of road for a specific window of time each day of the week (say, segment A to B from 5:00 to 6:00 pm, I also have this information from B to C and C to D, etc.). I have 12 weeks of data so I basically have 60 data points for each segment (weekdays only). My first thought was to present the average and standard deviation, but then I realized I am presenting an average of averages so my standard deviation should probably be a standard error. My problem is, in some areas I have multiple segments of road that I need to sum up and present the total average travel time (So I have 60 average travel times for A to B and 60 Average travel times for B to C, but in some areas of town I need to present the average travel time from A to C or A to D, etc.). Would it be wrong to present standard deviations of everything (including those areas of town where I am just presenting A to B), even though my data points are averages? If I am to present the standard error, and I am adding multiple segments, would my "n" just be the average of the n's in each segment? [assuming SE=SD/sqrt(n) is the correct equation to use for standard error] This post seemed very relevant: How to 'sum' a standard deviation?, but I was worried that my situation may be different because I need to present the "A to D" information in some areas. Hopefully I presented that clearly enough, it's a challenge to describe. Any insight will be greatly appreciated!


title: What is a better way to construct a confidence interval for the probability of success in binomial distributions?

text: For a variable $B\sim \textrm{Bin}(p,n)$. I observe $m$ successes. I know I can estimate $\hat{p}$ with $$\hat{p} = \frac{m}{n}$$ and I can approximate the CI by using a normal approximation by CLT where I assume $$p\sim\mathcal{N}\left(\hat{p}, \frac{\hat{p}(1-\hat{p})}{n}\right)$$ However if $m =0$ then I am in trouble (sort of) as the CI's lower bound would be negative. I remember there was a transform I can use involving a matrix that will eliminate that problem. What is the tranformation in question?


title: What is the benefit of using permutation tests?

text: When testing some null versus alternative hypotheses by a test statistic $U(X)$, where $X = \{ x_i, ..., x_n\}$, apply the permutation test with the set $G$ of permutations on $X$ and we have a new statistic $$ T(X) := \frac{\# \{\pi \in G: U(\pi X) \geq U(X)\}}{|G|}. $$   1. What is the benefit of using the permutation test over not using it? I.e. what is it like when the permutation test works?    2. What conditions to make that happen? Such as some conditions on the test statistic $U$ and/or on the null hypothesis? For example,   1. Should $T(X)$ be equal to the p-value based on $U(X)$, for sample $X$? If yes, why? (references is also appreciated) The p-value of $U(X)$ is defined as $$\inf_{c \in \mathbb R: U(x) \geq c} \sup_{F \in H} P(U(X) \geq c | X \sim F)$$. If the permutation test is to estimate the permutation distribution of $U(X) | X=x$, how is $T(X)$ equals the p-value of $U(X)$ at $X=x$? Especially, there may be more than one distributions in the null $H$, and $T(X)$ doesn't consider the null distributions one by one and then take $\sup_{F\in H}$ and $\inf_{c: U(x) \geq c}$.   2. Should the permutation test make $T(X)$ distribution-free over the null hypotheses? What conditions will make that happen?   3. Should $T(X)$ be uniformly distributed over $[0,1]$? What conditions will make that happen? Notice that when $U(\cdot)$ is a constant function, $T(\cdot)$ is also constant at $1$ and the distribution of $T(X)$ is far from being uniform over $[0,1]$. Thanks and regards!


title: Can mutual independence be defined in terms of conditional distributions?

text: Suppose there are more than two random variables $X_1, X_2, \dots, X_n$. Can mutual independence between them be defined in terms of conditional distributions? Thanks!


title: Dealing with multicollinearity of explanatory variables in panel regression when the usual remedies fail

text: I am regressing firm characteristics on some stock trading-related measures in a panel dataset. Firm size is a highly significant control variable, independent of the estimation method etc. My focus variables are related to firm size though, either by construction (e.g. $focus variable = x / firmsize$) or because of an economic relationship. As a consequence, I am finding myself in a classic multicollinearity situation: If firm size is put in as a control variable, my focus variables become insignificant. If firm size is left out, the focus variables are highly significant. Any of the usual advice (e.g. http://en.wikipedia.org/wiki/Multicollinearity) is not helping: I cannot obtain more data, I cannot run my regression on principal components because I need interpretable coefficients etc. I have little experience with this kind of problem but with some imagination, I came up with the following two ideas:   1. Running the regression with firm size as a control variables and additionally including interaction terms between each focus variable and firm size.   2. Trying to strip away the firm size effect from both the dependent variable and the focus variables, e.g. by first regressing firm size on the dependent/focus variable and then using the residuals as the dependent/focus variable in the actual regression. Would either or both idea make any sense? Any comment or alternative ideas would be very welcome!


title: How to deal with kurtotic interaction variable?

text: I'm going to do a multiple regression analysis with two continuous IVs and their interaction. I've centered the IVs and then multiplied one with another for interaction. My problem is that the interaction variable is severely kurtotic (kurtosis 7,731, standard error=0,127, z=60.87). How should I deal with this? Removal of the most extreme cases of course helps, but what would be the most recommended way to define outliers? N=1479 with all cases, the range of centered IVs is around 50 and the range of interaction variable around 700. The two IVs are normally distributed. Let me know if any additional info or exact numbers would help.


title: How do you tell whether good performances come in streaks?

text: I solve Rubik's cubes as a hobby. I record the time it took me to solve the cube using some software, and so now I have data from thousands of solves. The data is basically a long list of numbers representing the time each sequential solve took (e.g. 22.11, 20.66, 21.00, 18.74, ...) The time it takes me to solve the cube naturally varies somewhat from solve to solve, so there are good solves and bad solves. I want to know whether I "get hot" - whether the good solves come in streaks. For example, if I've just had a few consecutive good solves, is it more likely that my next solve will be good? What sort of analysis would be appropriate? I can think of a few specific things to do, for example treating the solves as a Markov process and seeing how well one solve predicts the next and comparing to random data, seeing how long the longest streaks of consecutive solves below the median for the last 100 are and comparing to what would be expected in random data, etc. I am not sure how insightful these tests would be, and wonder whether there are some well-developed approaches to this sort of problem.


title: Alternatives to box plots for comparing clustering results

text: In my research paper, I wish to compare the accuracy of 2 clustering algorithms. I have conducted clustering using these 2 algortihms on 300 proprietary datasets. I use adjusted rand index as the evaluation metric. Since, clustering has been conducted 300 times, I have 300 pairs of adjusted rand indexes to be compared against each other. I need a method to graphically show the performance of both algorithms and compare them. **Due to space constraint, I could not use barcharts**. Another alternative is to use box plots. I have some apprehension in using box plots to compare clustering algorithms' performances in academic research papers as I have not seen many papers using this approach. My questions:   1. Is it OK to use box plots in this scenario in academic papers?   2. If not, what alternative plots/ techniques could be used to graphically compare the 2 algorithms? ( **NOTE: I just want to compare the adj. rand indexes and NOT the time consumed or algorithmic complexity of the 2 algorithms** )   3. It will be helpful if one could refer papers where similar comparisons are done. Thanks.


title: Market Segmentation Test for SEM

text: I am working with SEM and would like to know how I could conduct a market segmentation test. I am trying the same model on two different samples and would like to test whether it's better to have them modeled separately or modeled as one bigger sample. Can anyone help me with this? Thank you, Deuterium


title: Should we conduct post-hoc tests when there is a significant interaction but no significant main effects?

text: I have data showing outcomes of some treatment on different people grouped by sex and age (grouped into: infant, child, adult). I need to investigate if there is any difference in outcome among age groups or between sexes. I used:               model <- aov(data$outcome ~ data$sex*data$age)     summary(model)      and found that neither sex nor age were significant on their own, but the interaction term was significant. I read about the Tukey test so I tried:               TukeyHSD(model)      and found that none of the differences were significant. Was this the right thing to do? Should my conclusion be that there is no significant difference or is there more I should do?


title: Sample covariance matrix

text: $\newcommand{\X}{\mathbf{X}} \renewcommand{\S}{\mathbf{S}} \newcommand{\I}{\mathbf{I}} \newcommand{\1}{\mathbf{1}} $ I found an article with an unusual (for me) covariance matrix. Let $\X$ denote an $N\times T$ matrix of $T$ observations on a system of $N$ random variables representing $T$ returns on a universe of $N$ stocks. Then the sample covariance matrix is defined (p.\ 606) as: $$ \S= \frac{1}{T}\X\left( \I - \frac{1}{T} \1\1' \right) \X' $$ I find difficult to understand the above estimator. I would normally use something like: $$ \frac{1}{T-1} \left(\frac{1}{T}\X\I - \X\right) \left(\frac{1}{T}\X\I - \X\right)' $$


title: Predictive Probabilities

text: Other than a calibration plot, is there a way to decide how good one models' predictive probabilities as compared to another model. I'm not interested in error rates as I find them ineffective for the level of precision I'm looking for. The only quantity of interest is the predictive probability distribution, as I am pricing contracts using them. EDIT: I have no faith in scoring rules based on the experience below with several different classifiers. I've simulated data from a known model. Trained the known model and a worse model using the simulated data, and the brier and log rules don't agree that the known model is superior. The class probabilities are materially different.


title: A proof involving properties of moment generating functions

text: Wackerly et al's text states this theorem "Let $m_x(t)$ and $m_y(t)$ denote the moment-generating functions of random variables X and Y, respectively. If both moment-generating functions exist and $m_x(t) = m_y(t)$ for all values of t, then X and Y have the same probability distribution." without a proof saying its beyond the scope of the text. Scheaffer Young also has the same theorem without a proof. I don't have a copy of Casella, but Google book search didn't seem to find the theorem in it. Gut's text seems to have an outline of a proof, but doesn't make reference to the "well-known results" and also requires knowing another result whose proof is also not provided. Does anyone know who originally proved this and if the proof is available online anywhere? Otherwise how would one fill in the details of this proof? In case I get asked no this is not a homework question, but I could imagine this possibly being someone's homework. I took a course sequence based on the Wackerly text and I have been left wondering about this proof for some time. So I figured it was just time to ask.


title: Signal processing techniques for unevenly spaced and repeated measures series

text: I am considering using signal processing techniques to find the minimum on a noisy 1D response line. More specifically I have a simulation that requires one parameter, but also includes randomness, and I want to adjust the parameter to fit some data. By fixing the parameter, I can have many different answers by running different simulations. The variance of the different outputs is not constant in the parameter space. I am planning to generate random samples of the parameter, and do some filtering on the resulting error from the simulator output to the real data. But: a) the points in the parameter space (line) may be uneven. I could do an evenly spaced grid search in the parameter values but if there is a solution where I can just sample according to the prior on the parameter it would be better. b) there may be more than one measure for the same parameter value. I can substitute these multiple values for the same parameter value by their mean, but then I have to weight that data point more than non-average data points. Can someone point me in some directions/literature on signal processing techniques for a) uneven spaced data b) multiple valued or weighted data Thanks jacques


title: Experimental evidence supporting Tufte-style visualizations?

text: Q: Does there exist experimental evidence supporting Tufte-style, minimalist, data-speak visualizations over the chart-junked visualizations of, say, Nigel Holmes? I asked how to add chart-junk to R plots here and responders threw a hefty amount of snark back at me. So, surely, there must be some experimental evidence, to which I'm not privy, supporting their anti-chart junk position--- more evidence than just "Tufte said so." Right? If such evidence exists it would contradict a lot of psychological research we have regarding humans, their memory recall, and pattern identification. So I'd certainly be excited to read about it. A little anecdote: at a conference I asked Edward Tufte how he regards experimental evidence finding that junk animations and videos improve humans' understanding and memory recall [see research cited in Brain Rules]. His response: "Don't believe them." So much for the scientific method! P.S. Of course, I'm needling people a little here. I own all of Tufte's books and think his work is incredible. I just think that his supporters have oversold some of his arguments. NOTE: This is a re-post of a question I asked on StackOverflow. Moderators closed it because it wasn't programming-specific. CrossValidated might be a better home. UPDATE: There are some useful links in the comments section of my original question post---namely, to the work of Chambers, Cleveland, and the datavis group at Stanford. UPDATE: This question deals with similar subject matter.


title: Pre-PCA pooling in R

text: I was just wondering if it is possible to pool the data into groups before doing prcomp() in R? I have i.e. 100 individuals (rows) and 50 measurements(cols) with individuals being grouped in separate populations (20 in the first, 30 in the second, 15 in the third...etc). Now PCA can give scores of individuals or measurements on PC axes, but I need to see what PC score is associated with populations. On the graph I can only plot individual scores but I need just let`s say 5 or 6 populations and their corresponding score. Any ideas?


title: How do you calculate the probability density function of the maximum of a sample of IID uniform random variables?

text: Given the random variable $$Y = \max(X_1, X_2, \ldots, X_n)$$ where $X_i$ are IID uniform variables, how do I calculate the PDF of $Y$?


title: What is a good source to learn about multidimensional scaling?

text: I don't have much experience with statistics (my education included 8 credits of pure statistics, another 8 of related subjects). I mostly learned about distributions and significance tests for distribution parameters, and a bit about linear regression. I also had some factor and component analysis, but have completely forgotten them and don't have the textbook any more. Now we have to evaluate some empirical data and the other people on our team have even less knowledge in statistics. After looking at different methods, I think that multidimensional scaling is the right technique to use - the examples I see it used for are similar to what we want to do. I found a book about it, but I find it very hard to learn from it. First, it is not at my level - the second paragraph of the "basics" chapter starts with the words "If σij is the affinity between object i and object j" and I am not sure what the affinity in my case is. Second, it is very old, from 1979, and I think there could have been improvements since then. Third, there is of course nothing about software support in this book. Fourth, it is in German, and googling about things I don't know is difficult, because there are better sources in English online, but I don't know the correct English translations of specific terms. What I want is a better book. Here are my criteria for it.   * It should be a textbook aimed at students, not a reference tome for specialists.    * I would appreciate if it is written from an applied point of view, but it is more important for me to grasp the theory than to see a script which is supposed to spew out the result I am supposed to want.    * It won't hurt if it is a more general statistics textbook which includes a sizable section on MDS (I might need to refresh other knowledge), but please not an exhaustive brick which covers all the material for a B.Sc. in statistics and costs its weight in Euro notes.    * I would prefer it to be in English, but if the ultimative source is in another language, please tell me about it.    * I will use R for the evaluation, so if there is a book which uses R for the examples, it would be great. A book which is dependent on any tool - even if it is R - is a complete no-go (like a book which says "to obtain this result, use function X of package Y" instead of explaining the calculations which have to be made to obtain the result).  If you know of a book which meets at least the majority of these criteria, I'd love to hear about it.


title: How is it possible to turn out with a highly accurate prediction when all records were classified the same way?

text: We ran a CHAID decision tree model using the set up and process described in my related question here. We used the propensity scores to come up with a prediction. We measured the prediction at the end of the year simply as actuals/prediction and got 95% accuracy. This was last year. This year we have a lower accuracy at 80% but not too bad. I am coming up to speed with what a previous analyst did. But when I inspected his model I noticed everyone was classifed the same: "stayer" (part of a binary target). The min and max of propensity scores ranged from .035 to .41. I don't understand why the model was so accurate as measured at the end of the year. It was great to have those results but it just doesn't seem logical to me. **QUESTIONS#1** Was this just a fluke (two years in a row!?!)? Or is there something I am missing with respect to intpreting classification results (confusion matrix) versus actual model performance? I can post a link to a SPSS Modeler .str file (and related text files, all anonimized) if anyone wants to see it. **QUESTIONS#2:** Do folks who use decision trees for prediction ever just use the classifications assigned (i.e. count number of target class) as opposed to propensity scores? I suppose if they were accurate you might do so (unlike our case). But I think propensity scores should be better to use in either case, right? EDIT: Everyone is classifed as 0 ("Stayer")   ![enter image description here](http://i.stack.imgur.com/qRokh.png) This is what I mean by accuracy and how it is being calculated ![enter image description here](http://i.stack.imgur.com/ItDjr.png) This is a sample of what my data looks like. ![enter image description here](http://i.stack.imgur.com/Yb7jO.png)


title: Computation of hypergeometric function in R

text: I'm having tremendous difficulty evaluating $_2F_1(a,b;c;z)$ with the `hypergeo` package in R. In my case, values of $a$, $b$, $c$ are always positive real numbers. Even so, the hypergeometric function is incredibly sensitive to their values. I am not looking for extreme precision; I can use Excel to get a rough estimation of the Guass hypergeometric that is fine for my purposes. Any suggestions for an implementation in R that will give a fast, error-free, if not super accurate Gaussian hypergeometric computation of positive real numbers with a wide range of values? Edit: it seems there is far more code for this in MATLAB than R. Any thoughts as to why?


title: HMMs in protein or NA sequence alignments

text: I would be very grateful if someone could explain to me how hidden markov models are used to generate a sequence alignment. I was trying to read on it, however ended up in attaining basic concepts and I'm interested in the mechanism that takes place in between the two points: I input 2 or more sequences (first) and get the alignment on my screen (second). I don't really understand the way an HMM is generated in the case of performing sequence alignments. I'll be thankful even if someone could provide me with links to sources to read on the subject. I will be very grateful for any help, Thank you.


title: Comparing and contrasting, p-values, significance levels and type I error

text: I was wondering if anybody could give a concise rundown as to the definitions and uses of p-values, significance level and type I error. I understand that p-values are defined as "the probability of obtaining a test statistic at least as extreme as the one we actually observed", while a significance level is just an arbitrary cutoff value to gauge if the p-value is significant or not. Type I error is the error of rejected a null hypothesis that was true. However, I am unsure regarding the difference between significance level and the type I error, are they not the same concept? For example, assume a very simple experiment where I flip a coin 1000 times and count the number of times it lands on 'heads'. My null hypothesis, H0, is that heads = 500 (unbiased coin). I then set my significance level at alpha = 0.05. I flip the coin 1000 times and then I calculate the p-value, if the p-value is > 0.05 then I fail to reject the null hypothesis and if the p-value is < 0.05 then I reject the null hypothesis. Now if I did this experiment repeatedly, each time calculating the p-value and either rejecting or failing to reject the null hypothesis and keeping a count of how many I rejected/failed to reject, then I would end up rejecting 5% of null hypotheses which were in actuality true, is that correct? This is the definition of type I error. Therefore, the significance level in Fisher significance testing is essentially the type I error from Neyman-Pearson hypothesis testing if you performed repeated experiments. Now as for p-values, if I had gotten a p-value of 0.06 from my last experiment and I did multiple experiments and counted all the ones that I got a p-value from 0 to 0.06, then would I also not have a 6% chance of rejecting a true null hypothesis?


title: How to get predictors in bagging in R?

text: I am using `RPART` for Classification and regression tree. I was tempted to use `bagging` function in `ipred` package in R, after reading about it. Now I am puzzled how to get the final set of rules similar to the ones I was getting in the output of `rpart`?               library(ipred)     data("BreastCancer", package = "mlbench")     myformula <- Class ~ Cl.thickness + Cell.size + Cell.shape + Marg.adhesion +                  Epith.c.size + Bare.nuclei + Bl.cromatin + Normal.nucleoli +                  Mitoses     #rpart      rpart(myformula, data=BreastCancer)     #Bagging     mod <- bagging(myformula, data=BreastCancer, coob=TRUE)     print(mod)     summary(mod)


title: Does a univariate random variable's mean always equal the integral of its quantile function?

text: I just noticed that integrating a univariate random variable's quantile function (inverse cdf) from p=0 to p=1 produces the variable's mean. I haven't heard of this relationship before now, so I'm wondering: Is this always the case? If so, is this relationship widely known? Here is an example in python:               from math import sqrt     from scipy.integrate import quad     from scipy.special import erfinv          def normalPdf(x, mu, sigma):         return 1.0 / sqrt(2.0 * pi * sigma**2.0) * exp(-(x - mu)**2.0 / (2.0 * sigma**2.0))          def normalQf(p, mu, sigma):         return mu + sigma * sqrt(2.0) * erfinv(2.0 * p - 1.0)          mu = 2.5     sigma = 1.3     quantileIntegral = quad(lambda p: quantile(p,mu,sigma), 0.0, 1.0)[0]     print quantileIntegral # Prints 2.5.


title: Previous Hourly Spend versus Lagged Hourly Spend Koyck and Transfer Functions

text: There are several posts I have read that addresses hourly electric load forecasts, but none address the question I have now. If I am building a Dynamic sort of regression model that attempts to forecast hourly demand for a product, based upon hourly advertising impressions or GRP's, how do I model lead/lag effects in the case where we know the previous hours GRP's are more influential than yesterdays GRP's at 8am. Ad spend at 6am as a lower influence on demand than adspend at 7 am, which is lower than ad spend at 8am on demand at 8am. So, there is an expected inherent decay, as if I were modeling lags (which would be previous day's value of ad spend). My question is, would we treat past hour series of a covariate the same way we would treat a lagged value series of a covariate in a Dynamic Regression model? If so, is there away to save on degrees of freedom by imposing some kinda of distributed lag structure on the covariates where instead of have lagged values of adspend at 8 am, lag 1 is replaced with adspend at 7am, lagged 2 is replaced with 6 am. Or do I just run a regular ARIMA tossing the earlier hourly spend values on the right handside?


title: What is the physical meaning of the probability density function and cumulative distribution function?

text: I have started research in Electronic Engineering, where PDF & CDF take a core part in most of the applications. I have studied books on probability where they have discussed the PDF & CDF formulas and basic theory. But my objective is different, like thinking about their applications: e.g. If I take any graph like a unit step function graph (which have 3 piece-wise defined sub-functions -ve infinity to zero, at zero and zero +ve infinity) I should know about this that should I take/evaluate its conditional PDF or conditional CDF? So generally I want to know the physical meaning of both conditional PDF & CDF in engineering applications? If anyone can explain this to me or recommend any literature on it that would be great.


title: Variance covariance matrix of regression coefficients with a probit link

text: Suppose we are performing ordinal regression using a probit link function. The data are doses (log transformed) and responses. The responses are ordinal and can be from 0 to 4. Suppose that some of the responses are not seen in the data (e.g. 1 and 3). The variance covariance matrix from R looks like:                                       0|2      2|4             log10(dataframe[, 1])     0|2                   3           4                     8     2|4                   57          59                    25.3434     log10(dataframe[, 1]) 37         36                     18      Note that there is no ` 0|1 `, ` 1|2 `, or ` 2|3 ` since there were no ones and threes in the data set. I am interested in the variance of the intercept, the variance of the probit slope, and the covariance between the intercept and the probit slope. In the above matrix, would these values be the $(1,1)$, $(3,3)$ and $(1,3)$ entries respectively?