title: What is the 95% confidence interval / ellipse in a PCA plot telling me?

text: I am new to statistics and within the last two days I tried to get my head around PCA plots. Now, I kind of understand what they are showing but I am still not sure about the 95% confidence ellipse that is very often shown in such plots. The 2 dimensional PCA plot displays the two biggest variances (whatever these are) in the data but I don't know what the ellipse is trying to tell me and what it means if a sample/dot (whatever is displayed) is lying outside that ellipse. Help is very much appreciated.


title: Use of Monte Carlo vs predictive modeling?

text: I am building a scoring model for rating an individual on trustworthiness. The parameters of the model are not fixed and I have no historical data to test. So to test the validity of the parameters and their weightings, I have built random samples using Monte Carlo simulation. I have now came across the term predictive modeling (logistic regression, etc). So here, does Monte Carlo or any other technique work?   Could anyone clarify when to use simulation and when to use predictive modeling?


title: definition of deep belief network

text: I was studying Deep Belief Network (DBN) and have questions. 1) According to the definition of DBN, DBN is formed by stacking RBM on top of each other such that the hidden layer in a lower layer becomes the input layer in the above layer. However, when I read the papers by Geoff Hinton (for example, "a fast learning algorithm for deep belief nets"), his DBN doesn't seem to have multiple RBMs. His DBN has only one RBM (i.e., undirected Markov Random Field with two layers) which sits on top while the other layers have only directional edges. I am confused by this difference. Are these two architectures essentially same? When the term, "DBN" is used in literatures, does it refer to Geoff Hinton's DBN? Or does it refer to a general class of multi-layer architectures with one or multiple RBMs. 2) Following up the above question, does anybody know why Geoff Hinton's DBN has undirected RBM (called associative memory) on top? What role does it play? What would happen if, instead of undirected edges, directed edges are used on top? Thanks,


title: Hypothesis test between ordered categorical independant var and continuous dep. var

text: Is there a standard hypothesis test for differences in the mean of a continuous dependent variable with respect to a single ordered categorical variable. By specifying that the independent variable is ordered, I also mean to impose the assumption that the dependent variable is monotonic in the independent variable.


title: How to find out if an online poker-site is fair?

text: Last week I had an interesting discussion with a good friend of mine. He had been playing some online poker and suggested that there is a relationship between new subscription/additional money transfer and the cards that you're dealt, i.e. you get good cards to get hooked. The sites would probably be risking a lot if this was true but the problem still fascinates me. My first approach to this was to ask my friend to define "good cards" and do a simple binomial test. My friend had though a hard time defining what exactly good cards are. If he gets really bad cards he knows to fold while if he gets good cards he knows to go all in - the bad cards are the ones in between. My other approach would be to calculate the exact probability of each given hand and then to see if it differs from the expected, perhaps using a Wilcoxon signed-rank test since this should detect a different distribution shape as well as a true shift. I guess the hard part is to calculate the exact probability. The data would consist of the first 0-100 dealt cards compared with 300-400 dealt cards a week later (or a friend that's been on the site for a while). **Question** : How would you suggest to approach the issue? **How Texas hold'em works** I'm no expert gamer (I've only played Texas hold'em 3-4 times) but it's fairly simple, you can find more details on the Wikipedia page here. The main difference from regular poker is that you only get 2 cards at start. You don't get to switch these cards. On the table are another 5 cards face down. By combining your two with the tables 5 you choose the best possible 5 card-poker hand. For instance if you get 2 aces you have a good start and you will probably go in strong, likewise if you have a 7 and a 2 your chances to win are very slim and you quickly fold. The hard part is perhaps a queen and a 9 where you might end up without anything although your cards are above the "average". You can find a list of the poker hands here.


title: Sum of variances from regression coefficients, larger then total variance. Why?

text: For a linear regression fit for a problem with p variables X_i ranging between 0 and 1, where p>20 (I don't know if that is relevant or not), and the number of samples is about 1000, I wanted to estimate the variance contribution for each of the variables using the regression coefficients. If I understood correctly var(A*X) = A^2*var(X), and therefore I thought that taking the square of the regression coefficients and multiplying that with the variance for each of the variables should give a vector containing all the variance contributions of the different variables. The problem is that I expected the sum of those variance contributions to be equal to the sum of the total regression model variance, but it isn't. The sum of the variances is some times up to 50% larger, then the total variance of the regression model. here some `Matlab` like pseudo code to explain the problem more in detail.               X    %sample matrix         Y    %output sample matrix            Linmodel=polyfitn(X,Y)   %fit a model          for ii=1:nr_colsX           VARCONT(ii)=var(Linmodel.COEF(ii)*X(:,ii))) %variance of the contributions               VARRC(ii)=(Linmodel.COEF(ii)).^2*var(X(:,ii))  %variance based on Reg. Coef.          end           SVC=sum(VARCONT(ii))  %sum of the variance contributions          SVSRC=sum(VARRC)          VY=var(Y)        %sum of the variance of the samples          VYmod=var(polyvaln(Linmodel,X))    %sum of the variance of the model on the samples          XR=rand(100000,nr_colsX)     %sum of the variance of the model with large number of samples          VYmodR=var(polyvaln(Linmodel,XR))      for one of the models that are supposed to be almost linear, VY is almost equal to VYmod. but SVC is about 50% largen then that. and VYmodR comes more in the direction of SVC. 1) Could some body please explain me why the sum of the variance contributions from the Regression coefficients, can be quite a bit larger then the variance of the regression model? 2) If this is so as it seems to be the case should there then not be some sort of upper bound for the sum of the square of the regression coefficients, such that the sum of their squares, multiplied by the variances of the input, should not be larger then the total variance of the output? Because it seems strange to me that the output of an interpolation model could result in a larger variance, then the variance of the data output points used for the interpolation. Any help is highly appreciated, but any help that is written in a way, such that also just a silly engineer as me can understand it is appreciated even more.


title: What does y argument do in chisq.test in R?

text: The `chisq.test` function in R includes a `y =` argument, which is to `NULL` by default. The help page doesn't explain what this argument does, and playing round with numbers doesn't give any clues, for example all these give exactly same results:               chisq.test(x=c(20, 10, 5, 3, 2, 1, 0), y  = 100:106)     chisq.test(x=c(20, 10, 5, 3, 2, 1, 0), y  = 0:6)     chisq.test(x=c(20, 10, 5, 3, 2, 1, 0), y  = 5:11)


title: Monotonically increasing transformation function

text: The dependent variable of my problem is highly concentrated around zero. Here is a stemplot                 The decimal point is 6 digit(s) to the right of the |            -2 | 511       -1 | 92221       -0 | 87777666666665555555555555555444444444444444444444444444444433333333+2428        0 | 00000000000000000000000000000000000000000000000000000000000000000000+2113        1 | 00000000000000000000000000000000000011111111111111111111111111111122+110        2 | 0000000000111111233333444444555666666677778889        3 | 00112457778889        4 | 11233456999        5 | 0000389        6 | 01477        7 | 259        8 | 033        9 | 002356       10 | 9       11 |        12 |        13 | 069       14 |        15 |        16 |        17 | 13      Normally when I have a dependent variable (DV) that looks like this I apply a logarithmic transformation for reasons both economical and mathematical. But obviously this will not work in this case as I have values below zero. Is there another monotonically increasing transformation function that will reduce the peakedness of this distribution?


title: How to rank film production companies based on weighted ranks?

text: IMDB uses a Bayesian estimate to compute its top 250 films. Let's say I know all the variables needed to compute this Bayesian estimate for these films. I also have the production company associated with each film. How could I fairly rank the production companies based on the weighted rank (and review counts) of their films? I'd be using all films with a more than m reviews. As a refresher, IMDB uses `weighted rating (WR) = (v / (v+m)) × R + (m / (v+m)) × C`. It probably doesn't matter, but I'm using R to do the analysis in case there are some packages that can help solve this.


title: For a research project I'm trying to estimate how many videos currently exist on Youtube

text: Official stats on how many videos (rather than how many hours or what volume of data) seem quite hard to come by. My current idea is something like this: A youtube URL is like: http://www.youtube.com/watch?v=w1sjRD7NSec where                id = w1sjRD7NSec      Assuming that each character in the id can be upper (26) or lower (26) case letter or number (10) the number of combinations should be                = (26 + 26 +10 ) ^ 11           = 62^11           = 5.2 x 10 ^19      So if I write a script to try random urls of 100/1000/whatever videos, and Z% are successful, then will                Z% x 5.2 x 10 ^19      give me the number of videos that exist (or at least are downloadable), albeit with a very low confidence?


title: lme: random effects for replicated growth curves

text: I am measuring the evolution of the brain response to a visual stimulation over time. The measures are done every seconds from 1 second to 14 seconds (each measure at time t gives a value summarizing the magnitude of the brain response from time 0 to time t). I have 8 subjects and 2 experimental conditions. For each subject and condition I replicate the measurement 12 times. I obtain therefore for each subject and condition 12 growth curves. My aim is to assess the influence of the experimental condition on the response profile. The responses appear to be fairly linear in time. Additionally, a noticeable aspect of the obtained results is the increase of the response variability over time (for each subject individually as well as between subjects). I am not sure about how to specify properly my random effect for the linear mixed model (using the nlme package). My initial idea was to start with the following model (time being coded as an integer and condition and subject as factors):`lm0 <\- lme( response ~ time*condition, random = ~time | subject/condition, data = psdData)`. Now due to the increase of variability of the response over time, I also get an increase of variability of the residuals over time. I initially wanted to try correcting for this using correlation structures, unfortunately, due to the replications, it seems that I cannot use the correlation option; when trying `lm1 <\- update(lm0, corr = corCAR1(form = ~time))` or `lm1 <\- update(lm0, corr = corCAR1(form = ~time | subject/condition))`, I obtain the error:               Error in Initialize.corCAR1(X[[2L]], ...) : covariate must have unique values within groups for "corCAR1" objects      and when trying `lm1 <\- update(lm0, corr = corCAR1(form = ~time | replicate/subject/condition))`, I obtain the error:               Error in lme.formula(response ~ time * condition, data = psdData,  : incompatible formulas for groups in 'random' and 'correlation'      As I didn't manage to use correlation structures, I tried to use variance functions to model heteroscedasticity: `lm1 <\- update(lm0, weights = varPower(form = ~ time))` But when looking at the residuals, it still displayed an increased variance over time. Same thing happened when using "varExp" or "varConstPower". I tried then another model, this time considering trials as units of measurements (and not subject as previously) and including the trial to trial variability within each subject. I created the factor `trialInSub <\- condition:replicate` which list all 24 trials of each subject (12 per condition) and wrote the model:               lm0 <- lme( response ~ time*condition, random = ~time | subject/trialInSub, data = psdData)      giving a random intercept and slope for each subject and for each trial within subject. As this model fits a line for each replicate, I don't have anymore the issue about increasing variance over time. However when looking at the autocorrelation of the models residuals, I would suspect something is wrong with the model: from lag 2 to 7, the autocorrelation decreases from 0.5 to -0.5 and from lag 7 to 10 increases back to 0.04. Additionally I have a clear linear trend when plotting residuals against fitted values of the model for time = 1 (it looks fairly normally distributed for all other times) I tried different correlation structures but it didn't change much my results. I am really not sure that I handle the problem in a proper way and that the models I tried really reflect the design of the experiment and the structure of the data. I hope someone can direct me towards a solution that would be statistically sound.


title: Standard deviation of game results about predictions from a rating system

text: I'm very in to Sports analysis and am keen to look at finessing my analysis models that I have worked up (I don't have a great maths background, I've just done a little bit of reading). Standard deviation of Game Results about prediction from a Ratings system (in this case least squares regression with a few tweaks) equals the stochastic element in a normal distribution function where: x = the margin of victory that you are trying to get a % against (e.g. it could be a point spread e.g. % to cover) the mean = the predicted margin of victory derived in this case from my least squares regression to create power ratings for teams sigma = standard deviation of error on the difference between the predicted margin of victory and actual margin of victory based on past results. Based on the above I would treat in Excel terms a win % for team A as: 1-NORMDIST(0.5,mean,sigma,true) - with the mean and sigma having values obviously. Based on Stern "The Probability of Winning an American Football Game" American Statistician (August 1991) - the sigma value (stochastic element? - sorry this isn't my usual thing) is approximately 13.86. How would you go about refining this for other sports? - e.g. Australian Rules Football. I said this could be done for any sport (and I believe it can within reason unless the results or the sport doesn't lend itself to this type of analysis). Looking at past AFL results from the end of the season point of view the standard deviations are as follows (have only started looking at these in the last two days): 2009: 31.34 (150+ results) 2010: 33.44 (150+ results 2011: 34.91 (150+ results) 2012: 25.62 (only approximately 70 results so far as it is mid season) As you feed more data into the model for 2012 it will build up a more accurate picture but looking at past data I'm putting the sigma value for AFL at around 33.5 (above the 25.62 value). Is there anyway of defining a best fit value? Where I don't have a great background in this kind of thing, if someone could give me some gentle nudges in the right direction in terms of where I should be looking or what I could be applying I would be grateful. I believe this can definitely be applied to other sports as well e.g.: English Premiership: 2008/2009: 1.45 (Season with over 300+ results) 2009/2010: 1.52 2010/2011: 1.54 2011/2012: 1.64 (but when looking over the last 100 results around the 1.50 mark) I've also looked at the percentages that this then gives for home wins, away wins and draws and they are virtually identical to the percentages for bookmakers odds (except they work to 105-106% instead of 100% so they have an edge which accounts for the small differences). NBA: 2008/2009: 11.30 (Season with 1000+ results) 2009/2010: 11.56 2010/2011: 10.88 2011/2012: 11.36 Thanks in advance for any advice,


title: Maximum-entropy joint pdf

text: I have multiple (finitely-many) variables with independent posterior pdfs $f_{x_1}...f_{x_n}$ based on Bayesian updating treating them as independent with independent evidence. If I subsequently discover conditional evidence relating the variables and need a joint prior, is there an analogue to a maximum-entropy prior forming a joint pdf $f_{x_i}(x_1, ...x_n)$ from the individual $f_{x_1}...f_{x_n}$? Sorry if this question is boring, obvious, nonsensical, or underdetermined. If the latter, what extra information is required to make it meaningful?


title: Does event A predict event B?

text: I'm want to show that a customer who buys paper books is more likely to buy an e-reader than a customer who does not buy paper books. Note that the number of e-readers purchased in the past is small (1,000) compared to the total set of events (1,000,000). Here's the historical data that I have:               Customer Name                | Month                       | Purchased Books in Last Month                           | Purchased Books in Last 12 Months                               | Purchased Books in Last 3 Years                                   | Purchased Book this Month                                       | Purchased e-Reader this Month     Joe Smith | 1/12 | 0 | 0 | 0 | 1 | 0     Joe Smith | 2/12 | 1 | 1 | 1 | 0 | 0     Joe Smith | 3/12 | 0 | 1 | 1 | 0 | 1     ...      Using these historical data, I'd like to demonstrate that book-buyers are more likely to buy e-readers, especially if they have recently purchased a book.   1. When I do a linear regression on the binary events for "Purchased Books in Last 12 Months" and "Purchased e-Reader this Month", I get r < 0.1.   2. However, if I compare P(purchased e-reader) to P(purchased e-reader | purchased books in last 12 months), the probability of a previous-book-buyer purchasing an e-reader is 10x higher than a random customer purchasing an e-reader How do I reconcile these two results? Am I using the right tools to answer this question with confidence?


title: Explanation of R diagnostic plot for logistic regression

text: I'm hoping someone can explain this bit of R code for me related to `glm()`. I don't understand the diagnostic plot that has been suggested. It seems a more informative plot would be to plot against the fitted values, but maybe I don't understand something. Here's the code:               result <- glm(survive~age, data=donner, family=binomial)     # Why is this plotted against the respondent index?     plot(residuals(result,type="pearson"), main="pearson residual plot")      Data to reproduce the above example:               > dput(donner)     structure(list(age = c(23L, 40L, 40L, 30L, 28L, 40L, 45L, 62L,      65L, 45L, 25L, 28L, 28L, 23L, 22L, 23L, 28L, 15L, 47L, 57L, 20L,      18L, 25L, 60L, 25L, 20L, 32L, 32L, 24L, 30L, 15L, 50L, 21L, 25L,      46L, 32L, 30L, 25L, 25L, 25L, 30L, 35L, 23L, 24L, 25L), sex = c(1L,      0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L,      0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 0L,      1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L), survive = c(0L,      1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L,      1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L,      0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L)), .Names = c("age",      "sex", "survive"), class = "data.frame", row.names = c(NA, -45L     ))


title: How is it correct to deduce significance level of Pearson's correlation from a table?

text: this is a rather elemental question, but any help would be appreciated anyway. This webpage gives a method of finding the p-value of Pearson's correlation coefficient, using a table. However, the methods listed in Wikipedia's article suggest that _using a table_ is only valid if the variables are special in some sense, or if there is enough data. So, my question is: under no assumptions for the random variables for which Pearson's correlation is calculated, is it correct to use a pre-computed table?


title: How to compute Wilcoxon rank-sum test in Stata?

text: I would like to compute the Wilcoxon rank-sum test in Stata for two countries of my sample with the dummyvalue 3 and 6. Can anyone give the right code for Stata?


title: How to test whether average of ten independent correlations is different from zero?

text: I have calculated the correlation between A and B for 10 different subjects. For averaging purposes, I've converted the r-values to z-values using Fisher's z-transformation.               subj zvalue     s1 -0.04     s2 -0.14     s3 -0.29     s4 -0.20     s5 -0.37     s6 -0.01     s7 -0.09     s8  0.20     s9 -0.15     s10  0.09      I want to test if the correlation between A and B is significantly different from zero. Do I use a paired sample t-test against 0 with the z-values, then report the r-value converted from the mean z-value? Or should I convert each subject's z-value back to an r-value before doing the t-test? Edit: Critically, is a paired t-test against 0 appropriate to test if this correlation is significant across subjects?


title: Generating cumulative incidence curve for competing risk model - assuming changing covariates

text: I am implementing a Fine & Gray competing risk model. The model includes covariates that are static for each person (e.g. location) as well as covariates that change over time (such as local unemployment rates). The issue I have is: how do you generate cumulative incidence curves, assuming covariates that change over time? It is straightforward to use the STATA command "stcurve" and generate a CI curve that assumes fixed values for covariates (e.g. covariate value X for all times), but there is no functionality to generate a CI curve assuming covariate value X in time t, and then covariate value Y in time t + 1. Is there a paper that discusses this, or a package in R/SAS/STATA that lets you generate curves assuming different values of covariates over time? Thanks very much for any help!


title: Baseline measurement using multiple years

text: I would like to know how to move forward in establishing a baseline on multiple observed years of retrospective data. My intent is to first identify the 10% of a given group to include in a study, and then establish a baseline. One concern of the data is that some observations may have an extremely large score in one year due to "bad luck". To attempt to take this into account, what I had done was look back at three years of data, and ranked each observation from lowest to highest (in a year). I then summed the ranks across three years, which would get me a final "ranking". I then took the highest 10%. My question is, would it be wise to use the mean as baseline given that some observations actually had a lower score in the third year OR should I use the third year score as "baseline"? If neither, what would be a better way of establishing a baseline? I've provided some code to replicate what I have done thus far. Thanks in advance. X1-X3: year 1, year 2, year 3   x1.1 - x3.1: year 1 rank, year 2 rank, year 3 rank   total: total sum for year 1, year2 and year 3.   mean: mean of year 1, year2 and year 3. rank: final ranking               x = c('dplyr', 'ggplot2')     lapply(x, require, character.only=T)          props = function(ncol, nrow, var.names=NULL){       if (ncol < 2) stop("ncol must be greater than 1")       p = function(n){         y = 1         z = sapply(seq_len(n), function(i) {           x = sample(seq(0, 10000-y, by = 1), 1)           y = y + x           return(x)         }         )         w = c(z)         return(w)       }       DF = data.frame(t(replicate(nrow, p(n = ncol))))       if (!is.null(var.names)) colnames(DF) = var.names       return(DF)     }          set.seed(9213)     df    = props(ncol = 3, nrow = 100)     df$id = rownames(df)     df    = df[,c(ncol(df),1:(ncol(df)-1))]          df.w = data.frame(df, apply(df[2:4], 2, rank, ties.method = 'min'))          df.rank = df.w %.%       group_by(id, X1, X2, X3, X1.1, X2.1 , X3.1) %.%       summarise(         total = sum(c(X1, X2 , X3)),         mean  = mean(c(X1, X2 , X3)),         rank  = sum(c(X1.1, X2.1 , X3.1))) %.%       arrange(id)          df.rank$rank = rank(df.rank$rank, ties.method = 'min')     df.rank = df.rank[order(df.rank$rank), ]          percent.10 = subset(df.rank, rank > quantile(rank, prob = 1 - 10/100))


title: What can be going wrong when Maximum Likelihood standard errors are high?

text: In maximum likelihood estimation (MLE) a useful result is that the standard errors for some estimated coefficient vector can be computed as the square roots of the diagonal entries of the inverse of the negative expected information matrix. In other words, let the empirical estimate _of the expectation_ of the information matrix be: $$I(\theta) =-\frac{1}{N}\sum_{i=1}^{N}\mathcal{H}(Z_{i}, \theta)$$ for parameter vector $\theta$, and data points $Z_{i}$, where $\mathcal{H}(Z_{i},\theta)$ denotes the Hessian matrix of the log-likelihood function, evaluated at the single datum $Z_{i}$ and the parameter vector estimate $\theta$. Then the matrix given by $I(\theta)^{-1}$ contains the variance/covariance structure of the MLE estimates, and its diagonal entries in particular are the variances of the estimated MLE parameters. In my specific application, the $Z_{i}$ are consumer data about price preferences for different heating options, and the model is a logit probability model. I've done all the work of estimating the MLE and there are standard formulas for the derivative vector and Hessian matrix in this logit scenario. I have about 250 data samples, which is large enough that I expect the MLE estimate to be fairly accurate. However, when I numerically compute the information matrix and the corresponding variances, I am seeing very large numbers. The estimated coefficients range in absolute value from 0.5 to 22, but the smallest standard error is about 7.2, which makes me question what I'm doing. I know that the empirical estimate of the Hessians are correct, as I've used my code on a separate data set to confirm that I'm computing the derivatives and second-derivatives correctly, so I think it's unlikely to be a coding mistake. Does anyone have any advice on what might yield such large standard deviations in practice?


title: Asymmetrical selective sampling for linear classification

text: I've got a online classification problem where I predict a class label `{+1, -1}` for an object and then show it to a user to get a real label. My task is to minimize a number of `-1` objects shown to a user. Obviously, the algorithm will not converge if it learns only on `+1 -> -1` missclassifications, so a certain ammount of `-1` objects needs to be shown to users anyway. The problem is to minimize a number of `-1` objects shown to a user by selecting only the ones most valuable for learning. The above problem is being solved in a field of 'active learning', where a learner gets to decide whether to request a laber for an object from user or not. Searching for a solution that does not imply i.i.d., I've found this work, where the learner decides on whether to require a label for an object by its margin. This is intuitive: the lesser margin means that we're more likely to make an error, which makes an object more valuable for learning. The probability to require an object depends on a margin and some parameter `b`, the optimal value for which is given for a linearly separable case, which applies for my problem. The algorithm is simple and goes as follows:               initialize weight vector w(t) with zeros     for every input vector x(t):         calculate margin p(t) = w(t-1) * x(t)         predict a class with Hebb's rule y(t) = sgn(p(t))         toss a coin with P(heads) = b / (b + |p(t)|)         if Heads:             query the label y' of x             if y' != y:                 update weights w(t+1) = w(t) + y'*x             else:                 w(t+1) = w(t)      Basically, I need to use this rule for querying user only if object prediction is `-1`, and if it is `+1`, I require a label anyway. My question is: would my algorithm converge in that setting? The perceptron algorithm itself and the modified version from the above article are symmetrical in respect to class labels. In my setting there will surely be more errors on `+1` objects, than on `-1`'s. **Does it hurt convergence? If it does, what modification to that algorithm can bring it back to 'symmetry' between positive and negative weight updates?** The modified algorithm querying `+1`s in all cases is as follows:               initialize weight vector w(t) with zeros     for every input vector x(t):         calculate margin p(t) = w(t-1) * x(t)         predict a class with Hebb's rule y(t) = sgn(p(t))         toss a coin with P(heads) = b / (b + |p(t)|)              if Heads or y(t) == +1:                  query the label y' of x             if y' != y:                 update weights w(t+1) = w(t) + y'*x             else:                 w(t+1) = w(t)      Is this algorithm going to converge? If not, what modifications can be made? The first idea was to decrease a learning rate for +1's, but still, I don't know how to prove convergence. Do you have any ideas on how to design an algorithm so that **it will surely find a linear separator in a finite number of iterations**? Thanks. **Update**. If don't know how to answer the question, but have any thoughts on where you would look for a solution or any pointers to publications/books/etc, please, let me know!


title: Record linkage: Weighting matches by estimate of match quality

text: This question is somewhat of a continuation from an earlier posting: Using EM algorithm for record linking I have two data sets of individual some of whom are in both but a prior it is not know which are and there are no identification id's that are common to both. What we know are the individual's first name, last name, birth year, and birth country. Since I still do not understand how to apply EM to link them, I am converting the name string to Soundex or NYSIIS and requiring exact matches. That is, linking individuals only if their NYSIIS-equivalent first names match, last names match, birth years match, and birth countries match, and this match is unique both within and across data sets. Would it make sense to then assess the match quality by say computing the Jaro-Winkler of the names of the match pairs? So say a pair has names Jon Smith and Jonn Smith and NYSIIS(Jon) == NYSIIS(Jonn), does it make sense to compute Jaro-Winkler (Jon Smith, Jonn Smith)? And if so, should I weigh by the corresponding analyses by the Jaro-Winkler estimates?


title: Likelihood-ratio test for different distributions

text: I aim to discriminate between three populations by using a maximum likelihood classifier. The idea is that every population has a unique distribution $\hat{f_i} (x)$. The pdf $\hat{f_i}(x)$ has been estimated on a set of training data. A new set of observations $(x_1,...x_p)$ will be assigned to the population for which $\log(q_i)+\sum_{j=1}^p \log(\hat{f_i}(x_j))$ is maximized. I want to investigate if the difference in the likelihoods are significant. How should I a perform a likelihood ratio test? I could of course perform the "standard" test in which I test $\lambda=\frac{L_0}{L_1} $ and then take $\chi^2= -2 \log \lambda$, but as I understand this test is used when there is a difference in the number of parameters between $L_1$ and $L_0$ . "Let L0 be the maximum value of the likelihood when the parameters are restricted (and reduced in number) based on the assumption. Assume k parameters were lost (i.e., L0 has k less parameters than L1)." Or is this correct? Otherwise how should the test be performed?


title: how to generate samples with the given autocorrelation function

text: I want to generate some spatial data where the points/location in the space form a multivariate gaussian distribution. I want these points to have certain autocorrelation given by the variogram model like spherical with spill 0.025 and range 15. I then want to generate some sample for this How can I do it? I then want to plot the empirical variogram and again fit a model variogram to check if they come out to be similar. How can I do it?


title: Given $n$ ratings from 0 to 1, how to calculate a weighted "average" estimating the "true" rating?

text: I've read How Not to Sort By Average Rating regarding how to average binary positive/negative ratings in a way that takes the number of ratings into account. The author uses the "lower bound of Wilson score confidence interval for a Bernoulli parameter". However, the items I'm dealing with have continuous ratings from 0 to 1. Is there an analogous averaging technique for this case? My ratings collection is long-tailed: the median item has only 2 ratings, but the average one has 80, and the most-rated item has 36,000 ratings. Intuitively ten ratings of 0.8 should "average" higher than one of 0.9, but I'd like a precise formulation of this intuition. (I'm using this to design a recommender system, which has to deal with 50,000 users and 10,000 items. I'm evaluating various known recommenders, like GroupLens and LSI, and have to design one that doesn't perform too much worse than those (and hopefully better). I was reminded of this blog post on averages when using users' average ratings for a baseline RMSE calculation.)


title: How to correctly choose model based on BIC?

text: I have a question about Bayesian Information Criteria. (GARCH models) I have looked for so many hours but still very confused about this BIC especially a negative one. As far as I am concerned it is okay to have a BIC that is negative, but the interpretation of them are different in each book on website. I am not looking for sophisticated answer just a normal explanation as if you were to explain someone who is not math or statisticians. Given same data,length and number of observation which model is better, based on BIC? >   1. -4.98749 > >   2. -4.995782 > >   3. -4.9864 > > Thank you in advance. I am using R software and running 3 models, GARCH-t, GJR model, and simple GARCH (1,1) model. So, I am trying to see which model is better, based only on BIC. I have already concluded what model is better based on other factors but this makes me confused.


title: how to analyze necessary/sufficient relationship?

text: Good day! I'm a second year psych grad student. From my fascination with formal logic, I wondered (just personal; not homework) about how to statistically establish necessary/sufficient relationships between variables. Let's suppose two variables, Intelligence and Creativity, whose relationship is summarized by the conditional, "if creative, then intelligent". So on the high end of Intelligence, we find the full spectrum of Creativity; but on the low end of Intelligence, we find only low levels of Creativity. We can say that Intelligence is a necessary (but not sufficient) condition for Creativity, while Creativity is a sufficient (but not necessary) condition for Intelligence. Two questions: (1) I imagine the scatter plot for Intelligence and Creativity would look like a right triangle whose non-hypotenuse legs are perpendicular to the axes. If Intelligence is on the x-axis, then the right angle would be at the +x,-y corner. Is this correct? intelligence on the x-axis; creativity on y-axis (2) If I imagined it correctly, what formula or procedure do I use to objectively measure how well data follows such a relationship? Is there an already-established way for dealing with this? My psych professor once told me that there is actually such a relationship between Intelligence and Creativity, but when I asked how to measure it, he said something only to the effect of "look at the scatter plot". All I know is that Pearson correlation is inappropriate, as the relationship is not linear. But I don't want to just look at the scatter plot. It's my first time to post here, so please be patient with me. Thank you!


title: can matrix completion work in the presence of many missing values?

text: I have a matrix with about 550k elements (2500 x 220) with 100k values known and the rest are unknown. Would it make sense to use matrix completion in this case, or are there too many values which are missing? What should I be careful of when I use matrix completion with it, if there are enough values there?


title: How to test whether mean and variance is the same in two small samples?

text: I would like to test two relatively small samples against the null hypothesis that both their means and variances are the same. The alternative would be that they in fact differ. I saw a post on this site advocating a ML test but I recall there is also a named test for this case, which I would ideally like to use in R, but whose name I forget. Can anyone help? Normality assumptions might be reasonable but would be difficult to test.


title: Conceptual issues in linear modeling

text: Mathematically an AR model is expressed as: $$ X_t = -\sum_{i=1}^p a_i X_{t-i} + \eta_t$$ Resources say that $\eta_t$ is a white gaussian noise. Is this $\eta_t$ the residual error? When invoking Matlab commands like `ar()` or `arburg()` or `aryule()` do we take into account this term? What is the difference between the $\eta_t$ and the residuals?


title: How do I best estimate population parameters from multiple sets of summary statistics on one sample?

text: Suppose that I have an assortment of summary statistics on a sample, and some beliefs about the underlying distribution, but no access to the sample itself. Each of the sets of summary statistics is sufficient to allow me to estimate the parameters of the distribution that I believe applies. For, instance, I may have (1) the quintiles (quintile break points) (with published sample variance estimates for each), (2) the inter-quintile sample means (with variance estimates), and, (3) for some finer disaggregation into value categories, the number of observations that fall into each category (with variance estimates). I know more or less how to do parameter estimates given any one of these summary statistics – though I am not sure how to make best use of the accompanying variance estimates for each. I could produce 3 distinct parameter estimates with a little work. But suppose I want to produce a single set of parameter estimates that makes the best use of all of the summary statistics that are available?   1. How would I combine the three (or however many I have – this is a simplified example) different types of summary statistics to get one set of parameters?   2. Can I make any further improvement to my parameter estimates by using the variance estimates that accompany each? When I say “best” I am thinking maximum likelihood, but if there is an interestingly different and perhaps better answer for minimum mean square error estimates or MAP estimates or something, I’d love to know that too.


title: Matlab TreeBagger Cost argument not working as it works with similar function fitensemble

text: The cost function of my TreeBagger class and fitensemble (Bag method) are both [0 8;1 0] for binary classification. The confusion matrix on fitensemble shows that the classfication tends to turn in the favor of the costy class (like [100 0; 20 80] favoring false negatives) but the same on TreeBagger does not hold. I have tried it on 3 datasets and TreeBagger seems to ignore the cost as I have omitted the cost argument and the confusion matrix is like when you dont give cost at all. Is there a Problem with my code or is it the TreeBagger? (BTW the problem is not misinterpretation of confusion matrix, it just doesn't work) My Code (TREEBAGGER WHICH DOES NOT WORK)               RF=TreeBagger(150,Xtrain,Ytrain,'oobpred','on','cost',[0 8;1 0])      (FITENSEMBLE WHICH WORKS)               Bag=(Xtrain,Ytrain,'Bag',150,'type','classification','cost',[0 8;1 0]);      Thank you for your help.


title: Can ayone help me in solving the below Problems of probability distribution?

text: Patients arrive randomly and independently at a doctor's surgery from 8·0 A.M. at an average rate of one in Five minutes.The waiting room holds 2 persons. What is the probability that the room will be full when the doctor arrives at 9·0 A.M.? Answer is 53.84% I have tried Poisson distribution with parameter as 12/Hour? and i am getting like 99% as answers.


title: MHT: Pre-Selecting statistical tests without Bias

text: **Summary:** The last formula boxed in red (which is a modified log likelihood from logistic regression) is a special non-differentiable loss function that is adapted to contain a Bonferroni correction in the loss function. The function is non-differentiable with respect to theta which is a problem for gradient descend and thus for using it to train regressors and neural networks. If it could become differentiable the function could optimise a classifier specifically designed to **increase the power of multiple hypothesis testing procedures in the case there is useful additional information that goes with the tests** and would then lead to more significant findings in cancer research. my question thus is how to make the formula boxed in red differentiable with respect to theta (preferably with some theoretical justification! ;P ).      ##  Intro For some research into the **Biology of Cancer** (details excluded for clarity sake) I devised a method to mitigate the issue finding important drivers of the process, which is in essence an issue of Multiple Hypothesis testing.   To do this I combine Machine Learning approaches with Multiple Hypothesis testing correction and aim to end up with valid/significant drivers of the cancer process.      ##  The Challenge Consider the following situation: One is supplied with a bunch of independent statistical tests that each have some additional information that comes with them as displayed in figure 1 below. For these statistical test the goals is to perform Multiple Hypothesis Testing (MHT) and call significance for all the 7 tests, using _Bonferroni correction_ :   ![Figure 1: A bunch of independent statistical tests](http://i.stack.imgur.com/rbNSb.png) _**Figure 1: Seven independent statistical tests, with additional information**_ Confronted with the situation in figure 1 one would perhaps like to make use of the fact that all the tests with additional information 'A' attached to them have very low p-values and in a way _pre-select_ them for Multiple Hypothesis Testing (MHT). In this example case the number of statistical test is small (7), but there are of course situations were one has to correct for many more (>10^6). The use of the additional information could in those cases really increase the overall power of the procedure.   Note that in this case the additional information is comprised of only one bit of information (A/B), but could just as well be a complex feature vector (this is where the **Machine Learning** comes in !)      ## Possible Solution 1 Given that the statistical tests are _independent_ (note: in my case they actually aren't fully independent but that is perhaps something for another post) one could use Machine Learning algorithms to try to predict the -log10(p) value for for instance test 1 by training a linear regression model on tests 2 till 7. One could do this for all the tests, effectively a leave- one-out cross validation and arrive at an unbiased prediction of the -log10(p) and use this information to select the really juicy 'A' type tests over the 'B' type, thereby decreasing the number of test in the MHT, increasing the power.   --> In machine learning terms: The additional information becomes the feature vectors (x), The individual tests are the samples and the -log10(p-values) of the tests become the labels (y). After this the tests/samples can freely be distributed in training and test sets for cross validation because all the samples are stated to be _independent_. There is only one issue with the approach above: It does not optimize for the desired end result being to find as much as possible significant results. In small sets like here this does not constitute a problem, but in large MHT problems with only a few small p-values this leads in effect to regression methods fitting the large 'uninteresting' part of the distribution, leading to poor performance on the small p-values it should actually focus on. A possible solution could be weighting of the samples, but this is problematic since it is not clear what an optimal weight should be.   Therefore another solution is needed that results in a model that actually optimizes the desired end result, being the largest possible number of significant results after MHT correction. This is the target of the next possible solution (2).   ## Possible Solution 2 To have a Machine Learning algorithm try to maximize the number of significant results on a test set, the method should be supplied with an _appropriate cost function_. If we assume for a moment that the dataset essentially consists of 2 classes being the significant tests and non significant test and given that these classes were known one could maximize the following log likelihood (cross entropy) as cost function to arrive at a model: ![Cross Entropy](http://i.stack.imgur.com/gYMhc.png) Which is a nice differentiable cost function that can be used for gradient descend. In this equation n is the total number of samples. yi represents the i'th example in the training set which can have the value 1 for the 'significant class' or 0 for the 'non-significant class' and htheta(xi) represents the output by the model, which has the following structure: ![Hypothesis h_theta\(x\)/Model](http://i.stack.imgur.com/QJhdY.png) With Theta being the model parameters and x being the feature vector. Now,.. If one reckons that the second term of the loss function (1st eq.) is concerned with assigning a cost to wrongly classified samples of the non- significant class removing it makes sense because the classification of non- significant samples is not of interest. One could also take the Bonferroni MHT correction into account in the cost function and one could modify the cost function into: ![Modified loss function](http://i.stack.imgur.com/TfLq4.png) Indeed the strikethrough is not part of the formula anymore. I{A} here is an indicator function that is 1 if A is true and 0 otherwise. In the brackets of the first indicator function above one can see the Bonferroni correction for the i'th p-value, in fact determining its class (sig. or non-sig). This class label also depends on m (the number of statistical test performed) which is determined by the sum of the next indicator function over the sample hypotheses. If htheta(x) reaches over 0.5 it is more likely according to the model that the sample is significant so it should thus be included in the final MHT correction. (note: The value of m/the cost function should eventually be corrected somewhat for the fact that the procedure excluded the one sample in the leave-one-out cross validation).   ## Concluding This last MHT cost function is a function that in fact optimizes the desired goal of finding the largest possible number of significant tests, but because of the indicator functions it is now non-differentiable, which is a big issue of gradient descent optimisation. The **Final Challenge** is to modify the MHT cost function to make it differentiable in such that it can directly be used for (stochastic) gradient descend (SGD). (One of the reason this is desired is because of the fact that that number of test is VERY large, but that SGD implementations on GPUs can crack the problem. Another reason is also the fact that it could be used to train a deep neural network, which has the benefit that it can model the highly non-lin representation that are likely to be found in the data.   I think that this cost function would really help to improve the number of discovered cancer drivers.      Hope you guys can help me out. Really curious though :) Perhaps there is even a better approach to this problem than my suggested path.


title: What's the best way to choose data for Crossvalidation on linear regression settings (PCA, PLS)

text: We are extracting features from EEG, which is a time dependent signal. We have signals of 10,000 datapoints over 64 channels, and we extract 10 features per timestamp per channel, so at the end we have a feature data set of 64x10 features for each time step. We are applying PLS to regress against hand movements (X,Y,Z). We are using cross validation to verify the results, but some of our colleagues are pointing out that using interleaved sampling from the datasets feels like cheating. Their suggestion is to use 2/3rds of contiguous data for training and 1/3rd for testing. I would like to hear the overall opinion on this setting, my intuition tells me that since it is time dependent data, just dividing the datasets like that would yield bad results (which it does)


title: Test to measure the significane of performance enhancement

text: I have 5 different computational tasks. I executed both tasks on two laptops named LA and LB. Out of 5 tasks, 4 tasks were computed in less time on LB, while 1 task took less time on LA. I want to measure, that the achieved performance gain is significant or not. Below is my data example:                          LA        LB     Task1      2.3       4.1 (mean values of 10 executions)     Task2      5.6       2.3     Task3      10.5      3.4     Task4      15.2      4.6     Task5      11.3      3.1      I have the following questions: 1. Which test should I apply to find the significance? 2. Can I apply that test on mean values rather than real 10 execution values? 3. Which parameters (Alpha, P, t) of the test must be presented in table and discussed, as I cannot present/discuss all parameters due to space limitation. Thank you


title: What is an Hypergeometric distribution where the last event is success?

text: I'm trying to find out the name of a distribution that is like negative binomial, only for finite population and without replacement. Or like Hypergeometric distribution where the last event has to be a success. That is: Let's say we have N balls in an urn, where W of them are white balls and B are black balls. I want to know what are the chances of needing to draw exactly n balls so it is the first time I got k of the n drawn balls to be black (e.g: that the last ball drawn was black, and a total of k out of n are black). I can't seem to find the name of such a distribution, so if you could also advise me on where to find names of distributions, that would also be lovely. Thanks.


title: Study- Probabilities problem

text: Can some of you tell me why the following holds true?: $ P(A \& B^C )= P(A)-P(A\&B)$ in which: $ B^C$ is $B$ complement. In which: $A$ and $B$ are independent. thanks..


title: parallelism in data mining softwares

text: I'm working on a data set for order prediction/classification and a close deadline upcoming. Fortunately, my university has a super-computer with restricted access. I was thinking of using a few nodes (each node is a 16-core processor, 2.3 GHz each and 32 GB of memory running on linux) to overcome the time limitation I'm facing. Basic software that I have in mind is WEKA, RapidMiner, MATLAB and maybe R. But as it's my first confrontation to parallel data mining, I need some help on this. I think the software needs to have a parallel nature in order to get most out of multi-core processors and using a lot of cores doesn't make any significant speed-up on serial programs. Best result of my searches was WEKA-parallel that is a modification of WEKA on _distributed_ systems. At last I'm open to any suggestion on getting most out of this computational power specially using WEKA. New software, new method, anything...


title: How do you "control" for a factor/variable?

text: To my understanding, "Control" can have two meanings in statistics.   1. Control group: In an experiment, no treatment is given to the member of the control group. Ex: Placebo vs Drug: You give drugs to one group and not to the other (control), which is also referred as "controlled experiment".    2. Control for a variable: Technique of separating out the effect of a particular independent variable. Some of the other names given to this techniques are, " accounting for", "holding constant", "controlling for", some variable. For example: In a Football viewing study (like or not like), you may want to take out the effect of gender as we think gender causes bias, that is, male may like it more than female.  So, my question is for point (2). Two questions: How do you "control"/ "account for" variables, in general. What techniques are used? (In terms of regression, ANOVA framework). In above example, does choosing male and female _randomly_ constitutes control? That is, is "randomness" one of the techniques for _controlling_ other effects?


title: How to compute the area between two graphs?

text: ![enter image description here](http://i.stack.imgur.com/9JFmv.gif) In the chart above, how can I derive the area between the two graphs, taking into consideration positive and negative areas? ![enter image description here](http://i.stack.imgur.com/lP3EJ.gif)


title: Does the central limit theorem imply the law of large numbers

text: Assuming that the distribution has finite variance (a condition not required for the LLN), then doesn't the LLN follow from the CLT?


title: Can I use a variable which has a non-linear relationship to the dependent variable in logistic regression?

text: Let's say I am building a logistic regression model where the dependent variable is binary and can take the values $0$ or $1$. Let the independent variables be $x_1, x_2, ..., x_m$ - there are $m$ independent variables. Let's say for the $k$th independent variable, the bivariate analysis shows a U-shaped trend - i.e., if I group $x_k$ into $20$ bins each containing roughly equal number of observations and calculate the 'bad rate' for each bin - # observations where y = 0 / total observations in each bin - then I get a U shaped curve. My questions are:   1. Can I directly use $x_k$ as input while estimating the beta parameters? Are any statistical assumptions violated which might cause significant error in estimating the parameters?   2. Is it necessary to 'linearize' this variable through a transformation (log, square, product with itself, etc.)?


title: How can you handle unstable $\beta$ estimates in linear regression with high multi-collinearity without throwing out variables?

text: Beta stability in linear regression with high multi-collinearity? Let's say in a linear regression, the variables $x_1$ and $x_2$ has high multi-collinearity (correlation is around 0.9). We are concerned about the $\beta$ coefficient stability so we have to treat the multi-collinearity. The textbook solution would be to just throw away one of variables. But we don't want to lose useful information by simply throwing away variables. Any suggestions?


title: Collapsing (combining) two variables into one for analysis

text: I came across a topic I am somewhat confused about: the merging of two variables. Assume we have two measurements from the same subjects. The two variables ($x_1$ and $x_2$) are measuring something similar, but not exactly the same thing. The variables (or combined variable, called $x_{12}$) will later be used as an explanatory variable ($X$) of some other variable ($Y$). For example, let's say we want to estimate the IQ of someone, and we only have the IQ of his father and mother (and we don't know the child's gender). What statistical (and non-statistical) issues are relevant for deciding on whether or not to combine the two measurements into one? Some issues to consider:   1. Let us say that we would later fit a linear regression of the type $Y$~$X$ (where $X$ is either $x_1$ and $x_2$ or the combination of the two), is there a time we would rather merge the two variables ($x_1$, $x_2$) into one?   2. How does the association of the two variables ($x_1$ and $x_2$) relevant to the decision on whether or not to merge them?   3. Is there any relation between $x_1$/$x_2$ and $Y$ that might influence the merging decision?   4. What if $x_1$ and $x_2$ are ordinal variables, or forced integer variables, does that make a difference on the value of collapsing them in order to be merged?   5. Is there any other issue to consider on this topic that I didn't mention?


title: What's a good, generic name for chart of things by time of day?

text: We're creating a chart showing traffic by time of day over a given period. So the y-axis is traffic, the x-axis is midnight, 1am, 2am, etc. It could also be days of the week. What's the generic name for this type of chart? I've come up with "cycle chart". Is that the standard? Is there one? ![example](http://i.stack.imgur.com/IDu3k.png) Update: Just to add a bit more clarity, what's being shown in the top chart is _not_ one day, it's an aggregation of many days. E.g. over the last month, 6am has on average been lower than noon. Similarly, in the bottom chart, over the last year, traffic dips on Saturdays.


title: Mean in compositional data

text: I am required to analyse samples a compositional data set. However my samples come from seven different locations. Is there any way for me to check if there is a location/mean effect in the data which I can remove from the Isometric Log Ratio (ILR) transformation of the data before analysing? Can I carry out the test on the ILR data, which is obviously straight-forward, or must the test be on the compositional data?


title: Goodness of Fit Test for Logistic Regression with small n_i

text: I would like to test how well my model fits the data. The response is binary and the Chi-Squared Test cannot be applied for the residual deviance because the $n_i$ are $1$. To use the Chi-Squared GOF Test, the $n_i$ need to be $\geq 5$. What's an alternative method can I use to test for goodness of fit?


title: Standardized Student's-t distribution

text: We know that density for a student-t distribution is given as $$\frac{\Gamma(\frac{\nu + 1}{2})}{\Gamma(\frac{\nu}{2})} \left(\frac{\lambda}{\pi\nu}\right)^{\frac{1}{2}} \left[1+\frac{\lambda(x-\mu)^2}{\nu}\right]^{-\frac{\nu+1}{2}}$$ with $\text{E}(X) = \mu$, $\text{var}(X) = \frac{1}{\lambda}\frac{\nu}{\nu-2}$ where the three parameters are shape factor $\nu$, location $\mu$ and dispersion $\lambda$. Now if I standardize X to Z as $\frac{X-\text{E}(X)}{\text{Std}(X)}$, what is the distribution of $Z$. If I standardize my input data to $Z$ as above, what is the density function for $f_Z(z)$? I came across this question in Stephen Taylor's Asset Dynamics book. I am a bit confused after reading this question. As I understand, $\frac{X-\mu}{\sigma}$ where $\sigma=1/\sqrt(\lambda)$ will give the standardized t distribution with center 0 and disperson 1. But not sure what $Z$ above will look like.


title: feature weights in structured support vector machine

text: I like to find the feature weights in a structured SVM for ranking the features w.r.t. importance. I know that in a binary SVM the weight vector can be written as a linear combination of examples. But how do you compute the same for an SSVM?


title: Confidence Intervals for Holdout R^2?

text: Let's say I'm performing regularized regression and I want to validate the results using holdout. (I'm choosing holdout instead of cross-validation because my dataset is fairly large, so computational power is an issue and the difference between holdout and cross-validation estimates will likely be negligible.) I want to put a confidence interval on the $R^2$ score for the holdout set. Is the following procedure valid?   1. Compute the regression coefficients on the training data.   2. Compute the predictions on the held out data. Treat each prediction as if it were an observed random variable.   3. Compute the Pearson correlation of the predictions with the actual $Y$ values for the held-out data. Compute a confidence interval for this correlation using Fisher's Z-Transformation.   4. If the confidence interval of this correlation is some set of real numbers, then the confidence interval of $R^2$ is just the set of the squares of these numbers. I can't think of any reason why this procedure would be wrong, but I've never heard of it being done before, so I wanted to run it by some more experienced statisticians.


title: For inference of Dirichlet Process Mixture, why the expected value $\int h(x)f(x)$ is desired?

text: Why the expected value $\int h(x)f(x)$ is desired for inference in Dirichlet Process Mixture? What is the intuition for MCMC in Dirichlet Process Mixture? $f(x)$ is the probability density function, and $h(x)$ is the outcome of the Dirichlet Process Mixture (I am not sure). This formula is originally from page 40 of slides.


title: Estimating number of intersections between point and polygons

text: I have a 2D plane (a large rectangle) with a finite size in the x and y direction, which is the field of my problem. The field is covered by $n$ smaller rectangles that are located randomly within its bounds (the smaller rectangles are free to overlap each other). If I "sample" the field at some random location I will intersect some number $i$ of smaller rectangles, which happen to overlap this location. My goal is to estimate the expected number of intersections, with a confidence interval (say 95%). Is there a way to do this analytically? Any help is appreciated. Lots of detail and laymans terms would be appreciated. As you can probably tell, statistics is not my strong suit. Edit: Please see comments below for important additional information.


title: Using Poisson distribution to generate random integers

text: I'm trying to generate random integers which have Poisson distribution. The open source library GSL has one such distribution. > Function: `unsigned int gsl_ran_poisson (const gsl_rng * r, double mu)` This > function returns a random integer from the Poisson distribution with mean > `mu`. The probability distribution for Poisson variates is `p(k) = {mu^k > over k!} exp(-mu)` I want to know what it actually returns. Does it return `k` for a particular value of `p(k)`? For a particular value of `mu`, `k` can take values from `0,1,2,3...` So if I specify a value for `mu`, do I have an upper bound on the values returned by this function?


title: How did Efron imagine the bootstrap?

text: Do you know Bradley Efron? He's a great man.   How did Efron imagine or think about "bootstrap" for the first time?


title: Comparing Hierarchical Linear Model with non-HLM techniques

text: I'm a fairly inexperienced statistician fighting a huge deadline and just need some peace of mind that I'm not making a massive error here. I'd be most grateful for any pointers. I've been working on fitting an `lmer()` model that explores the spatial relationships between crime and socio-economic factors. The ICC demonstrates a 40% attribution of variances to group factors, which is encouraging. The model's evolution has seen a steady improvement of AIC scores. I am now a little confused about how to demonstrate its effectiveness against other non- HLM techniques. Is this purely by standard error residuals, or am I missing something important?! Many thanks


title: Constructing a model from multiple non-independent and unreliable predictors?

text: I have an interesting modelling problem in which I am trying to forecast the occurrence of a type of weather event using an empirical model driven by measurements of a number of different physical quantities in different locations. There are a number of complications. First, the predictor measurements are not independent in the sense that the driving force that generates these weather events (which we don't measure directly) has an impact on all the quantities that we measure, driving them from their normal values and giving us warning of the weather event to follow. The second complication is that the training data is relatively sparse, so it cannot support highly complex models. Using simple linear models is problematic though, because of the stated non-independence of the predictors. We could use an approach like principal component regression to solve this, except for the third complication, which is that for both the training data we have and in operations, data is frequently missing due to sensors being broken or offline at the time. So, what strategies might be useful here? There are a few basics approaches I can see. Say for example we want to forecast the probability of event F occurring, based on the unreliable predictors A and B. We could train 3 completely separate models: $$P(F|AB), P(F|A), P(F|B)$$ using whatever approach works best for each model and then use whichever model is appropriate given the available data. This however does not generalise for many predictors (we have at least 50 different quantities we could use) due to the curse of dimensionality. The next simplest approach I've considered would be to construct separate models for each predictor alone $$ P(F|A), P(F|B) $$ and then learn weights for each model so that we have $$ P(F|AB) = \frac{w_A P(F|A) + w_B P(F|B)}{w_A + w_B} $$ the weights would be learned over the whole training set, setting the relevant weight to zero for each instance where a predictor is missing. This is a much more tractable solution, but I fear it would miss important interactions between the predictors and would sacrifice too much performance for the gain in simplicity. A third approach I can vaguely conceive of but I'm not sure if it is valid, is to try and express $$P(F|AB)$$ in terms of the individual models $$P(F|A), P(F|B) $$ and the conditional probabilites e.g. $$P(A|B)$$ the form of which would need to be learned from the training data. I haven't managed to come up with a useful form for this yet, although maybe it isn't possible (or advisable). The curse of dimensionality is likely to be a factor here as well I guess. I realise this could be phrased as a data imputation problem but I'm looking for a solution that gives a good trade-off between practicality and performance.


title: Forecasting at individual versus grouped level

text: I have monthly usage data (spanning 3 years) for a customer base of around 200K, and I need to generate 1-month ahead forecasts for each of them. There are a couple of exogenous variables that would need to be included too. One way to go about is to build an "Arima with exogenous variables" model for each customer. In which case, I would be looking to manage 200K models. I'm wondering if there is an alternate approach to solving this -- Assuming we are able to segment the customer base into groups that are "homogenous" in terms of usage, then would it make sense to create one model for each group? Would it be a good approximation to have the usage data time series for each group correspond to the mean values across the customers in that group? Any pointers to literature related to this problem is appreciated. Additionally, I'm looking to use R, so pointers to specific packages will be of help as well. Thanks!


title: Test means of populations with lots of zeros, is sampling the way to go?

text: I am trying to test whether means of two populations with lots of zeros are different. Here is the following python code example:               from scipy.stats import mannwhitneyu          import numpy as np          a = np.random.random(100)     b = np.random.random(100) * 2          aa = np.hstack((a, np.zeros(1000)))     bb = np.hstack((b, np.zeros(1000)))          np.random.shuffle(aa)     np.random.shuffle(bb)          mw_stat_1, p_value_1 = mannwhitneyu(a, b)  # 100 obs in a and b     mw_stat_2, p_value_2 = mannwhitneyu(aa, bb) # 1000 zeors added to each a and b          # I take mean of sample of 20 elements, a and b become an array of 55 elements     samp_sum_aa = aa.reshape(-1,20).mean(axis=1)     samp_sum_bb = bb.reshape(-1,20).mean(axis=1)          mw_stat_3, p_value_3 = mannwhitneyu(samp_sum_aa, samp_sum_bb)      Result:               >>> p_value_1       # no zeros     2.5956488654494193e-09     >>> p_value_2     0.42124151395226317     >>> p_value_3       # using sampling     0.0020853586023447269      I find that if I do Mann-Whitney test on raw populations (after zeros are added), my p-value is large; however, if I take random samples, I get a small enough p-value, such that I can reject the null hypothesis for all practical purposes. Is sampling a proper technique here? If so, how do I know what is the right sample size? Are there other methods to address this problem?


title: Is McDonalds (1999) omega a measure of unidimensionality or homogeneity?

text: I recently read a paper in which the authors used Cronbach's $\alpha$ and $\omega$ as measures for unidimensionality and homogeneity. In case of $\alpha$ I know that it is wrong, but I am not so familiar with $\omega$. The literatur I read said it is the proportion of test variance due to a general factor. Does that mean a high value in $\omega$ would imply unidimensionality and homogeneity, because it seems some researchers seem to make a point to differentiate between $\omega$ and unidimensionality and homogeneity, but I do not see why.


title: Proving some properties of expected first order statistics with respect to sample size

text: **Question:** Consider $n$ random variables $x_1, x_2,\cdots x_n\sim \mathcal{N}(0,1)$. The expected value of the $i$th order statistic (the maximum) can be written as $E(\mathcal{O}^n_1)= \displaystyle\int_{-\infty}^{+\infty}nx\Phi(x)^{n-1}\phi(x)\:dx$. I wish to show that for $n_1<n_2<n_3$, $E(\mathcal{O}^{n1}_1)<E(\mathcal{O}^{n2}_1)<E(\mathcal{O}^{n3}_1)$, and $E(\mathcal{O}^{n3}_1)-E(\mathcal{O}^{n2}_1)<E(\mathcal{O}^{n2}_1)-E(\mathcal{O}^{n1}_1)$, where $\mathcal{O}^{n1}_1$ is the first order statistic (maximum) for a sample of $n1$. **Progress so far** I've managed to prove the part for first-order derivative by invoking the concept of first-order stochastic dominance. But still no progress on the second-order... Any help will be greatly appreciated. Thanks!


title: How to use k nearest neighbours for binary classification with unbalanced classes?

text: I have relatively large (100k items) dataset which I need to split in two groups. So far I've tried knn and the results are not good mainly because I have disproportion in my training data: 90% of points belong to the first group. The same proportion is expected to be in test data. Is there a way to improve prediction quality with this kind of data? Performance is not important while quality of prediction is paramount.


title: Definition of a sample: can it include the same object twice?

text: Wikipedia defines a sample as: > a subset of a population. While exploring the reason why we divide by $(n-1)$ instead of $n$ when calculating standard deviation (discussed in this question), I came across this PDF demonstrating why $(n-1)$ is better. When listing all possible samples of $n=2$ from a population of three cards numbered 0, 2, and 4, it includes the samples (0,0), (2,2), and (4,4). I am having trouble reconciling this with the definition of a sample that I thought I knew (and that is given by Wikipedia). A sample of 2 playing cards from a population of 52 would not include the Three of Hearts twice, would it? Likewise, I'd guess a survey of a sample of voters would not include the same voter multiple times. Other sources back the method described in the PDF. What am I misunderstanding here?


title: K successes in Bernoulli trials, or George Lucas movie experiment

text: I'm reading "The Drunkard's Walk" now and cannot understand one story from it. Here it goes: _Imagine that George Lucas makes a new Star Wars film and in one test market decides to perform a crazy experiment. He releases the identical film under two titles: "Star Wars: Episode A" and "Star Wars: Episode B". Each film has its own marketing campaign and distribution schedule, with the corresponding details identical except that the trailers and ads for one film say "Episode A" and those for the other, "Episode B". Now we make a contest out of it. Which film will be more popular? Say we look at the first 20,000 moviegoers and record the film they choose to see (ignoring those die-hard fans who will go to both and then insist there were subtle but meaningful differences between the two). Since the films and their marketing campaigns are identical, we can mathematically model the game this way: Imagine lining up all the viewers in a row and flipping a coin for each viewer in turn. If the coin lands heads up, he or she sees Episode A; if the coin lands tails up, it’s Episode B. Because the coin has an equal chance of coming up either way, you might think that in this experimental box office war each film should be in the lead about half the time. But the mathematics of randomness says otherwise: the most probable number of changes in the lead is 0, and it is 88 times more probable that one of the two films will lead through all 20,000 customers than it is that, say, the lead continuously seesaws"_ I, probably incorrectly, attribute this to a plain Bernoulli trials problem, and must say I fail to see why the leader won't seesaw on average! Can anyone explain?


title: What are the rules / guidelines for downsampling?

text: I have a data set with ~ `7 million` rows, of which ~ `100k` are positives. I'm looking to shrink the data by keeping all the positives and then randomly sampling `several hundred thousand` negative examples to round out the data set. I'm uncertain what the guidelines are though. Are there any rules of thumb for what % needs to be positive / negative? Any suggestions would be appreciated.


title: How to model number of days in the last week smoking cigarettes (0 to 7 - 'U' shaped)?

text: I am currently analysing data where the outcome variable is 'U' shaped. The outcome variable asks 'how many of the last seven days have you smoked'. Most responses to this fall in the first (none) and last (all seven) categories. Because of this I do not think a count data model is appropriate. What would be a good approach to modelling this variable?


title: Singular matrix: eigenvalues perturbation vs Moore-Penrose generalized inverse

text: We often face singular matrices in practice: OLS with singular (X'X), GMM with singular weighting matrix, singular matrix in Wald statistics. I'm wondering how can we overcome this issue. I've seen two solutions used in different contexts: 1) Ridge regression: http://en.wikipedia.org/wiki/Tikhonov_regularization 2) Generalized inverse: http://en.wikipedia.org/wiki/Generalized_inverse What do you think would be preferable to do when the estimator or test statistics depends on the inverse of a matrix which is singular?


title: Fitting an HLM model in lme4

text: I'm a fairly inexperienced statistician fighting a huge deadline and just need some peace of mind that I'm not making a massive error here. I'd be most grateful for pointers. I've been playing around with different model combinations using the `lmer()` function, and have managed to get the AIC scores down from over 10,000 with the most basic combinations to -180 by factoring in various random effects. But I'm concerned I might be creating bias by using variables that already correlate. Is this a valid concern or does the model mitigate for this? Many thanks in advance.


title: Does the vanishing gradient in RNNs present a problem?

text: One of the often cited issues in RNN training is the vanishing gradient problem [1,2,3,4]. However, I came across several papers by Anton Maximilian Schaefer, Steffen Udluft and Hans-Georg Zimmermann (e.g. [5]) in which they claim that the problem doesn't exist even in a simple RNN, if shared weights are used. So, which one is true - does the vanishing gradient problem exist or not? * * *   1. Learning long-term dependencies with gradient descent is difficult  by Y.Bengio et al. (1994)   2. The Vanishing Gradient Problem During Learning Recurrent Neural Nets and Problem Solutions by S.Hochreiter (1997)   3. Gradient Flow in Recurrent Nets: the Difficulty of Learning Long-Term Dependencies by S.Hochreiter et al. (2003)   4. On the difficulty of training Recurrent Neural Networks by R.Pascanu et al. (2012)   5. Learning long-term dependencies with recurrent neural networks by A.M. Schaefer et al. (2008)


title: How many samples do I need to prove that a classification algorithm is better than another?

text: I have two algorithms A and B, used to automatically classify each of N elements into K categories, N and K both being in the millions. Neither A or B is perfect, but it is relatively easy for a human to do the classification, or put another way, to select which one of the algorithms A and B did a better job for a given element. I want to show that algorithm A is better than B, but am not sure how to approach this problem. As A and B most often say the same thing, I want to focus on the cases where they don't. My question is now, how many elements do I need to check to show (with a certain degree of confidence) that A is better than B? Or put another way, given a random sample from where A and B differ where n_a elements were correctly classified by A, and n_b elements were correctly classified by B, how likely is it that A is better than B?


title: How to deal with an SVM with categorical attributes

text: I have a space of 35 dimensions (attributes). My analytic problem is a simple classification one. Out of 35 dimensions, more than 25 are categorical and each attribute takes more than 50+ types of values. In that scenario, introducing a dummy variable also will not work for me. **How can I run an SVM on a space which has a lot of categorical attributes?**


title: Markov blanket vs normal dependency in a Bayesian network

text: While I was reading about Bayesian networks, I run into "Markov blanket" term and got severely confused with its independency in a Bayesian network graph. Markov blanket briefly says that every node is only dependent on its parents, children and children's parents [it is gray area for node A in the picture]. ![Markov blanket](http://upload.wikimedia.org/wikipedia/commons/e/ea/MarkovBlanket.png) What is the joint probability of this BN, $P(M,S,G,I,B,R)$? ![alt text](http://www.aiqus.com/upfiles/IMG_0095_copy.png) If I follow the step parent only independency rule, it is: $$ P(M | S)P(S | G,I)P(I | B)P(R | B)P(G)P(B)$$ However, if I follow the Markov Blanket independency, I end up with this (notice $P(I|\mathbf{G},B)$ is different): $$P(M | S)P(S | G,I)P(I | \mathbf{G},B)P(R | B)P(G)P(B)$$ So which is the correct joint probability of this BN? Update: Crosslink of this question in AIQUS and Respective chapter and diagrams are below: ![alt text](http://img828.imageshack.us/img828/9783/img0103s.png) ![alt text](http://img406.imageshack.us/img406/3788/img0104l.png)


title: Binomial GLM with unknown totals?

text: Essentially, I have some covariate data `X`, and a dependent variable `Y` consisting of proportions of a sample that shown a certain response (i.e. between 0 and 1). I suspect I want to proceed via a GLM approach, but the thing is, I don't know the sizes of each of those samples! My thoughts are to proceed by a quasibinomial methodology, estimating the dispersion parameter. Assuming the sizes of each sample are not too different, I can thus keep the logit link between the linear predictor and the proportion, but disregard the contribution of `n` in the usual binomial variance of `np(1-p)`. Then I can do hypothesis testing the usual way? Does this make any sense? Some R code:               #simulate some data          X = rnorm(500)     Z = rnorm(500)     p = exp(X*0.1 + 2)/(1+exp(X*0.1 + 2))     n = 50     Y = NULL     for (i in 1:length(X)){         Y = c(Y,sum(runif(n) < p[i])/n)     }     Y2 = cbind(Y*n, n-Y*n)          #glm, binomial model is 'true' ?     summary(glm(Y~X+Z, family = "quasibinomial"))     summary(glm(Y2~X+Z, family = "binomial"))     anova(glm(Y~X+Z, family = "quasibinomial"), test = "Chisq")     anova(glm(Y2~X+Z, family = "binomial"), test= "Chisq")      Seems to work, but am I missing something? Surely someone's done something like this before?


title: Explain the distribution of these random variables?

text: Suppose I have 3 eggs, with independent hatching times (in minutes) following an $\text{Exp}(0.05)$ distribution, and A is defined as the earliest hatching time.   If I were to plot A, how would I explain its distribution?   From running a simulation on R, I can see that the histogram of A follows an exponential distribution but I'm not sure why?


title: Is ETA a good measure for computing the efect size between an ordinal and a nominal variable?

text: I have a nominal and an ordinal variable and I would like to look measure the association that is between the two. For now, I have done Chi-square test and non-parametric Kruskall-Wallis Test and I have come to the conclusion that there is an association/relationship between the two variables. How can I further measure the effect of this association? I have read about ETA, which SPSS suggests using when you have a nominal and a quantitative variable, and I have also read that you can use it with ANOVA (which is basically the parametric version of Kruskall-Wallis). Could I use this statistic (ETA) as a measure of size effect and trust the value that Crosstabs option in SPSS provides me? As a secondary question, how can I interpret the values of ETA? I came upon different guidelines on the net. Would you describe a ETA=0.321 as a medium effect? ~Thanks.


title: Functional dependency in Support Vector Machine input

text: I'm currently training a support vector machine with a second order Markov transition matrix. These probabilities are taken from a computer image, and are the probabilities of adjacent pixels having the difference between them being specific values - so each value can be represented as a fraction $a/b$, where a and b are the cardinality of sets. Now, I've come across a curious phenomenon, which I'm trying to explain. I can remove a significant part of my input features with no loss in the accuracy of the support vector machine. My current theory as for why this occurs is because there are functional dependencies between certain statistics. For example, the statistics for one sample might be $(n_1/d_1, n_1/d_2, n_2/d_2, n_2/d_3, n_3/d_3, n_3/d_4, n_4/d_4, \dots, n_n/d_1)$ Each $d_i$ is used seven times in total, whilst each $n_1$ is used either once (all parameters identical) or twice. Now, I'm not as familiar with support vector machines as many, so I'd like to ask - is there a theorem or something which would confidently allow me to state that the above feature set is no better a feature set than: $(n_1/d_1, n_2/d_2, n_3/d_3, n_4/d_4, \dots, n_n, d/n)$ in terms of the constructed hyperplane?


title: What are some examples of real-world processes that are well-described by AR, MA, ARMA, or ARIMA?

text: Subject says it all - AR/MA/ARMA/ARIMA are often described as workhorses of time series analysis. But what are some real-world examples where these methods gave great results, and another more modern method wouldn't have been superior?


title: ICC reduced to 0 in a multi-level model. Should this be suspicious?

text: I calculated a sequence of multi-level models with stata (xtmixed). I have a sample of 800 kids in 46 classes and I want to explain their individual level of prejudice by the proportion of fully prejudiced respondents in the class. This is based on a three categorical index, the level of prejudice is measured by a principal component. In the empty model the ICC is 0.12. In a model where I control for the usual suspects (age, gender, school type, etc), it is reduced to 0.08. In a third model where I include the "proportion of fully prejudiced students in the class" variable it has a significant** positive effect, just as I excepct it to. The ICC is reduced to 0.0000 (23 0-s in a row). My colleagues tell me that this is not usual, they've never seen anything like this before. In order to test for validity I calculated a linear regression with cluster and robust options. To me it seemed reassuring: Model with control variables had R square of 7%, model with the "proportion of fully prejudiced students in the class" variable had 15% of explained variance. This suggested me that the variable "proportion of fully prejudiced students in the class" has an additional explanatory power of 8%. And as the ICC reduced from 8% to 0 in the ML model, this seemed logical to me. Have you ever seen an ICC of 0.0000000000000000? Is it OK? Thanks


title: How to tell if data is filtered well?

text: I have about 3 hours worth of 3-axis accelerometer, gyroscope and compass data taken at subsecond intervals. I'm fully aware that this data needs to be filtered and that there are a shedload of methods for filtering it (FFT, ARMA, et al). How can I tell from my data or a plot thereof, which is the "best" filter to use? For example, if I have 3 hours of a jog on the street and I'd like to track the length of the jogger's stride. I'm thinking I would integrate the gyroscope to get the absolute direction and then double integrate the raw accelerometer data to get distance. Finally, divide the result of this by the number of troughs in the spline encompassing the points (as they would represent mainly downward acceleration). Is this the correct approach? Many thanks! ![](http://analysis.d8u.us/~hdiwan/accelAndGyro.png)


title: Choice of a non-parametric method

text: I want to compare the **delays** of packets belonging to some kind of traffic to the delays of packets belonging to some other (different and larger) traffic, all generated from the same machine. I would like to see if the **distributions** of the two groups of delays **are the same**. Samples for the first group are between **size 20** and **2 000** (bust most of the time between 20 and 75), while the comparison group is always of the order of **a few thousands** (up to 15 000). When the first group has **more than 250** values, I just (randomly) **subsample** it in smaller groups of size 50, since it I noticed that it really takes nothing for large samples to get rejected with Kolmogorov- Smirnov. If the **majority** of subsamples gets accepted by the test, the original sample also does. Now, due to the precision of timestamping, there are lots of duplicate delays (i.e. **ties** ) in my observations. For instance, out of **15000** delays in my second group, I have **only 1000 unique values** , with a mean of 0.7 milliseconds and a standard deviation of around 7 milliseconds. 200 values appear up to 10 times; 90 values between 11 and 100; 30 between 101 and 300; 4 between 301 and 500. What should I take into account in order to **choose a non-parametric test** that best suits my case? Also, is my subsampling correct? I've been using a significance level of 0.01, but I'm getting **20% of false positives** , of which roughly half greater, half smaller than the second group (after the two-sided test, I always used the two one-sided ones just to check). Any ideas?


title: Predict GLM poisson with offset

text: I know this is probably a basic question... But I don't seem to find the answer. I'm fitting a GLM with a Poisson family, and then tried to get a look at the predictions, however the offset does seem to be taken into consideration:               model_glm=glm(cases~rhs(data$year,2003)+lhs(data$year,2003),     offset=(log(population)), data=data, subset=28:36, family=poisson())          predict (model_glm, type="response")      I get cases not rates... I've tried also               model_glm=glm(cases~rhs(data$year,2003)+lhs(data$year,2003)+     offset(log(population)), data=data, subset=28:36, family=poisson())      with the same results. However when I predict from GAM, using mgcv, the predictions consider the offset (I get rates). I'm missing something? I would appreciate any comment, thanks Sandra


title: How to algorithmically find plausible covariates at which to plot predicted responses?

text: I'm trying to plot the predictions of a linear model with multiple numeric predictors as a collection of plots such that $\hat Y$ is on the y-axis, each $X_j$ in turn is on the x-axis, and the remaining $X_{i\ne j}$ are all set to plausible values given $X_j$. Here is some example R code (syntax might look strange, it's that way for conciseness, but it does parse correctly in R):               mdata<-data.frame(x1=x1<-rnorm(100,1.5),                       x2=x2<-x1+2*rnorm(100,0.5),                       x3=x3<-x1*2+rnorm(100,sd=.2),                       y=x1-x2+x3-0.5*x1*x2+0.5*x2*x3-0.5*x1*x3+rnorm(100,sd=.1));          mfit<-lm(y~x1*x2*x3,mdata);      Note that in the above dataset, x2 and x3 are both correlated with X1, and all three are predictors for y. The widespread approach of setting each $X_{i\ne j}=0$ or to the population mean or median is misleading if there is the possibility of multicollinearity among the $X$s because no distinction is made between commonly-occurring combinations of variables and ones that are completely out of the support range of the data. Example using above objects:               predictAtX1<-seq(min(mdata$x1),max(mdata$x1),len=100);          ## Do you trust this? I don't.     plot(predict(mfit,newdata=data.frame(x1=predictAtX1,                                          x2=mean(mdata$x2),                                          x3=mean(mdata$x3)))~predictAtX1,                                          xlab='x1',ylab='y-hat');          ## Ditto.     plot(predict(mfit,newdata=data.frame(x1=predictAtX1,                                          x2=median(mdata$x2),                                          x3=median(mdata$x3)))~predictAtX1,                                          xlab='x1',ylab='y-hat');          ## Ditto, even worse, and yet the most common value at which to fix x_js.     plot(predict(mfit,newdata=data.frame(x1=predictAtX1,                                          x2=0,                                          x3=0))~predictAtX1,                                          xlab='x1',ylab='y-hat');      I thought of fitting each of the $X_j$ to a given $X_i$ and using the predicted values for a given sequence of $X_i$....               plot(predict(mfit,newdata=data.frame(x1=predictAtX1,                                          x2=predict(lm(x2~x1,mdata),                                              newdata=list(x1=predictAtX1)),                                          x3=predict(lm(x3~x1,mdata),                                              newdata=list(x1=predictAtX1))))~predictAtX1,                                          xlab='x1',ylab='y-hat')      ...and that doesn't look as obviously nonsensical as the plots for globally fixed values of $X_2$ and $X_3$ but given that there are additional interactions involving $X_2$ and $X_3$ that are being ignored, these predictor values are not representative of the population either. But. For a given sequence of values for $X_1$ it should be possible to construct a smoothed path through the variable space corresponding to the most likely values of $X_2$ and $X_3$. I can think of brute-force optimization approaches, but given how obvious and universal this problem is for any real- life regression model with multiple predictors, there must be an already- implemented and/or simpler approach to dealing with it. So, my revised question is, does anybody know the name of such a method or the name of the subfield of statistics that deals with it, or the name of an R package that implements it? Thanks.


title: Changepoint Analysis on differenced data

text: I was told to run a changepoint analysis in R on my data. But the data shows a decreasing trend, so the changepoint analysis won't be of help (I don't think?). If I differenced the data, to get rid of the trend, and then perform a changepoint analysis on the differenced data, does that make sense? And how would I interpret the results? Thanks so much


title: MANOVA for mixed model with multiple repeated measurements

text: I struggle performing a MANOVA for a mixed model. I have the following columns: S - subject (1,...,N) B - in between group (B1, B2) W - repeated measurement - within (W1, W2) M - measured dependent variable 1 N - measured dependent variable 2 My data is organized the following way:               S; B; W;   M;   N     1;B1;W1;12.3;14.4     1;B1;W2;35.4;13.5     2;B1;W1;11.1;22.2     2;B1;W2;33.3;11.1     3;B2;W1;15.3;12.3     3;B2;W2;20.3;18.8     ...      If I had only one dependent variable (M), I could obtain the results of an ANOVA in one of these ways:               dat=read.csv(file="data.csv", sep=";")     library(ez)     ezANOVA(dat, dv=.(M), wid=.(S), between=.(B), within=.(W), type=3, detailed=TRUE)          library(lmerTest)     results = lmerTest::lmer(M ~ B*W + (1 | S), data=dat)     anova(results, ddf="Kenward-Roger", type=3, method.grad="Richardson")          library(car)     options(contrasts=c("contr.sum", "contr.poly"))     mod1 <- lme4::lmer(M ~ B*W + (1 | S), data=dat)     car::Anova(mod1, test.statistic="F", type=3)      I wonder, how the code for the MANOVA has to look like? I tried               dat=read.csv(file="vor_actions_new.csv", sep=";")     results = manova(cbind(M, N) ~ B*W, data=dat)     summary.aov(results)     summary(results)      but I doubt, the within-factor is handled correctly. I also tried               mod1 <- lme4::lmer(cbind(M, N) ~ B*W + (1 | S), data=dat)      but, this gave an error, probably because multiple responses are not supported. ezAnova does also not support multiple responses, I think. Is there a way to execute the planned MANOVA?


title: SAS optimal binning

text: In SAS enterprise miner we have the optimal binning feature which allows you to transform continuous variables into an ordered set of bins. The binning, as I read from one of their docs, is done so that the log odds of the predicted categorical variable (good/bad) is monotonically increasing or decreasing. Can we use classification trees to kind of achieve the same result?. Like build a tree with only one variable and how many which splits we can get with sufficiently different log odds.


Is an implementation of a density function for a logit-normal distribution available in R? I have not found one in a package or in the CRAN task view of probability distributions. This is for a MLE fitting of a function to data.


title: What is the distribution of the error around logistic growth data?

text: In ecology, we often use the logistic growth equation: $$ N_t = \frac{ K N_0 e^{rt} }{K + N_0 e^{rt-1}} $$ or $$ N_t = \frac{ K N_0}{N_0 + (K -N_0)e^{-rt}} $$ where $K$ is the carrying capacity (maximum density reached), $N_0$ is the initial density, $r$ is the growth rate, $t$ is time since initial. The value of $N_t$ has a soft upper bound $(K)$ and a lower bound $(N_0)$, with a strong lower bound at $0$. Furthermore, in my specific context, measurements of $N_t$ are done using optical density or fluorescence, both of which have a theoretical maxima, and thus a strong upper bound. The error around $N_t$ is thus probably best described by a bounded distribution. At small values of $N_t$, the distribution probably has a strong positive skew, while at values of $N_t$ approaching K, the distribution probably has a strong negative skew. The distribution thus probably has a shape parameter that can be linked to $N_t$. The variance may also increase with $N_t$. Here is a graphical example ![enter image description here](http://i.stack.imgur.com/EpTp9.jpg) with               K<-0.8     r<-1     N0<-0.01     t<-1:10     max<-1      which can be produced in r with               library(devtools)     source_url("https://raw.github.com/edielivon/Useful-R-functions/master/Growth%20curves/example%20plot.R")        * What would be the theoretical error distribution around $N_t$ (in consideration of both the model and the empirical information provided)?   * How doe the parameters of this distribution relate to the value of $N_t$ or time (if using parameters were the mode can not be directly associated with $N_t$ eg. logis normal)?   * Does this distribution have a density function implemented in $R$? Directions explored so far:   * Assuming normality around $N_t$ (leads to over estimates of $K$)   * Logit normal distribution around $N_t/max$, but difficulty in fitting shape parameters alpha and beta   * Normal distribution around the logic of $N_t/max$


title: Cox regression - Hazard and Survival Estimates?

text: I have a couple questions about Cox survival regression: 1) Is it true that the hazard function h(t) is not available (even WITHOUT time dependent covariates)- and if not, is it because the baseline hazard is not defined (only the cumulative hazard is and you cant go from the cumulative hazard to the hazard)? 2) When you have time dependent covariates is it not possible to estimate the survival curve for a chosen set of covariate values (I'd like to plot prototypical values)? Can anyone explain why not?