title: How to visualize Bayesian goodness of fit for logistic regression

text: For a Bayesian logistic regression problem, I have created a posterior predictive distribution. I sample from the predictive distribution and receive thousands of samples of (0,1) for each observation I have. Visualizing the goodness-of-fit is less than interesting, for example: ![enter image description here](http://i.stack.imgur.com/bDczC.png) This plot shows the 10 000 samples + the observed datum point (way in the left one can make out a red line: yea that's the observation). The problem is is that this plot is hardly informative, and I'll have 23 of them, one for each data point. Is there a better way to visualize the 23 data points plus there posterior samples. * * * Another attempt: ![enter image description here](http://i.stack.imgur.com/pQ42U.png) * * * Another attempt based on the paper here ![enter image description here](http://i.stack.imgur.com/Wy8rE.png)


title: How to determine multivariate confidence in laymans terms

text: I am trying to determine colors to best grab a web-users attention. I have a simple system which randomly chooses one of 4 foreground colors, and one of 6 background colors. I then save the stats in a table like this:               fore_color     back_color    shown    clicked     ----------     ----------    -----    -------     red            green         16       1     blue           green         15       2      At the moment, I require every combination to be shown 250 times before I come to a conclusion about which combination is best... For a test size of ~6000 views. I would like to know how I can reduce the size of the test while maintaining high confidence that we are making the correct choice (>90% admittedly, I don't know if my current method meets that standard). Because I am only interested in 90% confidence for the best result, and not at all interested in stats for the losers, we can obviously disqualify early-bad- performers more quickly than final-contenders. To be clear, my goal is to attain 90% confidence of the best results while showing each combination the least amount of times possible. I will be running this test over 1000 times on many different web pages, and would like to eventually include more variables (like font size). If the math is simple, I could implement this directly in my own programing. If it's more complex, I'd be more interested in a (preferably free) package that I could feed the data through, and return confidence levels that a particular combination is or is not our winner. UPDATE: Here is a help page from google website optimizer that does exactly what I am looking for: http://support.google.com/adwords/bin/answer.py?hl=en&answer=61146#results I would like to be able to reproduce this type of report locally. Using google website optimizer is not a good choice for me because I would need to enter in every combination for >1000 tests by hand.


title: Estimating distribution from censored data

text: $X$ is a positive variable with known support (assume discrete support, if that simplifies solution). $Y$ is another variable with the same support. $X$ and $Y$ are independent. $Z$ is equal to $X$ if $X < Y$, and equal to $0$ otherwise. $Y$ and $Z$ are observed, $X$ is not. How to estimate distribution of $X$? I realize that there may be no single objectively optimal answer. Assume prior estimates, if needed for Bayesian inference.


title: I have a problem in bayesian networks get p(E|A)

text: I'm doing this book "Modeling and reasoning with Bayesian Networks" and I have this problem:               ConsidertheDAGinFigure4.14:     (a) List the Markovian assumptions asserted by the DAG.       (b) Express Pr(a, b, c, d, e, f, g, h) in terms of network parameters.      (c) Compute Pr(A = 0, B = 0) and Pr(E = 1 | A = 1). Justify your answers.      (d) True or false? Why?     - dsep(A, BH, E)      - dsep(G, D, E)      - dsep(AB,F,GH)      I have done almost all the question, but I'm missing c) p(E=1|A=1), I have tried p(E and A)/p(A), but this does not seem correct and for p(A=0, B=0) is this correct p(A) p(B)=(0.8)(0.3), any help?? Thanks. ![Network](http://i.stack.imgur.com/vkMnQ.png) ![Parameters](http://i.stack.imgur.com/cGYpA.png)


title: Guess several distribution's parameters from a list of interarrivals

text: Just a small introduction to the setting: I have traffic that is generated using several layers. Layer 3 consists of the basic packets, while Layer 2 is a more high level grouping of packets, and Layer 1 is even higher. I have now this list of interarrival times between packets/groups for each layer. These interrarrival times are taken from distributions, with each layer having its own distribution and parameters. Now I want to guess the used distribution and parameters for each of these layers. My first approach to just make a histogram of all interarrival times, and extract from this histogram 3 distributions. However, as these different layers can generate values close to each other, this proved to be a useless approach. My second idea is to iterate over each number, and see if it belongs to Layer 3, by first estimating the parameters using the observations we think belong to Layer 3 (e.g. for an exponential distribution we just take the sample mean), and then see if our new observation belongs to this distribution by checking the CDF against the lower and upperbound of the confidence interval. If it seems the value does not belong to Layer 3, we repeat the procedure at Layer 2 and so on. When all numbers are processed, I perform a goodness of fit test to select the most probable combination of distributions. As I am still new to this, I was wondering if this last idea would be a good approach, or do I have some wrong ideas on how to use these tools, like the confidence interval? Or maybe there has been some research done on this subject, that I haven't found? For the geometrical distribution there doesn't seem to be a closed formula for the confidence interval for a small amount of observations. Are there other means to determine if a value would would fit a distribution, given its parameters that are guessed from a limited number of observations?


title: Expectation of correlated variables

text: I'm looking to compare effects $ \delta = \frac{\mu_T - \mu_C}{\sigma}$ for two studies, compared with the same control group. In order to find the covariance of effects for treatment A and treatment B, I need to find $$ E \left[ \left( \frac{\bar{A} - \bar{C}}{s} \right) \left( \frac{\bar{B} - \bar{C}}{s} \right) \right] $$ I can assume $\bar{A}, \bar{B}$ and $\bar{C}$ are sample means and approximately normal, $s$ is pooled standard deviation with the usual relation to chi square. I believe all individual terms should be independent as well. My instinct is to simply expand $$ E \left[ \left( \bar{A} - \bar{C} \right) \left( \bar{B} - \bar{C} \right) (s^{-2}) \right]$$ $$E \left[ \bar{A} \bar{B} - \bar{A} \bar{C} - \bar{B} \bar{C} - \bar{C}^2 \right] E[(s^{-2})] $$ But this loses the difference groupings that I'm interested in. My advisor also said that this covariance should be related to the variance of the common components, but I'm failing to see how this is the case. Any advice on how to proceed would be very much appreciated.


title: Conditional Logistic Regression, Rare Events and R

text: I have a panel data set organized around 35 provincial units, with 110 binary "positive" outcomes in the dependent variable. This is for about 2600 observations in that data set. I have a couple questions. 1) Would it even be reasonable to run a logit at the national level, combining each province count? 2) Is it reasonable to try for a fixed-effects model, using condionoanl logit/clogit() in R? 3) Would a rare events model be a better alternative and is it possible to do it with fixed effects in R? Thanks.


title: Assessing multicollinearity of factors

text: I have a multifactor model (with 7 factors currently) and 754018 observations. In order to check for multicollinearity issues as the model grows I wrote an R script to compute a correlation matrix from the factors (so with 7 factors it's a 7x7 matrix). I then apply the function: $$\frac{r}{(1-r^2) / (N-2)}$$ where $N = 754018$ and $r$ is the sample correlation in order to get a test statistic according to http://faculty.vassar.edu/lowry/ch4apx.html. Then I get a corresponding p-value, and display those factor pairs (and their sample correlation) whose correlation p-value is less than $0.05$. After running this I get 12 pairs displayed! With 7 factors the total number of possible pairs is 21 so this is pretty bad. Out of these 12, however, only 4 of them have correlations above .1 and the rest have sample correlations of around .02 or .01...from a practical standpoint should I worry about any nonzero correlation with very small pvalue (as in all 12 pairs) or only those with small pvalue AND high sample correlation? As in maybe only those 4? If the latter, are there empirical ways of choosing a threshold sample correlation? Thanks


title: Multivariate normality test from R's energy

text: Could someone explain me which is the null hypothesis of the energy test of multivariate normality (R package `energy`)? By using this test I obtained a very different p-value from that I obtained with `mshapiro.test`.


title: What to do when the intersection of some time series doesn't contain enough data point?

text: I have to do some studies on time series, which involve computing the covariance matrix. However, my time series span on different time intervals and their intersection is too short to provide me a good estimator. I guess I have to remove some time series but I don't know how to select the ones to keep. I can come up with my own methodology (it could be something like keeping the data that have the most recent points) but I would like to know first if there is a standard procedure for this kind of things.


title: HMM library, different length sequences training

text: I'm using the Kevin Murphy's HMM library in MATLAB(http://www.cs.ubc.ca/~murphyk/Software/HMM/hmm.html) There is a section called 'How to use the toolbox'. There is this example for GMM ouputs:               Let us generate nex=50 vector-valued sequences of length T=50; each vector has size O=2.     O = 2;     T = 50;     nex = 50;     data = randn(O,T,nex);      But in my case, The length T of the training sequences is variable. I could solve the problem for the discrete case using cells, but I don't know how to do it here because 'data' is a three-dimensional variable.


title: Drawing random vectors from multivariate normal distribution with limited deviation

text: In the one-dimensional case it is easy to obtain random numbers that are not too far away from the mean by checking if they are within one or two standard deviations. Is there some similar way to obtain random numbers that are not to far away from the mean vector of a multivariate normal distribution? Thank you very much


title: Backtesting/cross-validation for time-series and prediction intervals

text: Suppose I carry out the following exercise with my trusty statistics software:   * Fit some time-series model to the data $y_1,\dots,y_t$ and calculate $\hat{y}_{t+1}$, the forecast of the next observation, and the error $e_{t+1}^*=y_{t+1}-\hat{y}_{t+1}$ for that forecast observation.    * Repeat the previous step for $t=m,\dots,n-1$ where $m$ is the minimum number of observations needed for fitting my model.   * I plot the distribution of the errors $e_{m+1}^*,\dots,e_{n}^*$ or calculate its percentiles **What relationship does that distribution have to the analytical prediction interval for $\hat y_{t+1}$?** My intuition is that the distribution of the errors from this iterative CV process does not tell you very much about the variability of the prediction from the forecast made using the _final_ version of the model. As the model is trained on more data, the errors will tend to decline with each step. So the large errors will be from early versions of the model and the small errors will come from the later versions. The final version of the model is more like the small-error late version, so it makes no sense to consider the early large errors as coming from that final model's error distribution. Even if there's no improvement in the model as it is fed more data, many time series model produce analytical prediction intervals. That will tell you if the difference in the actual and the prediction you observe is an outlier or not.


title: How to best account for population sizes in disease outbreaks

text: Say one has the fatality and survivor rates of a pathogen in terms of geographical location (e.g., the latest data on the MERS-CoV).               France                      2 (1)     Italy                       3 (0)     Jordan                      2 (2)     Qatar                       2 (0)     Saudi Arabia               49 (32)     Tunisia                     2 (0)     United Kingdom (UK)         3 (2)     United Arab Emirates (UAE)  1 (1)     ______________________________________     Total                      64 (38)      The total fatality rate is ~59%. However, I am not sure that this is a good measure: the rate for Saudi Arabia is ~65%, while the smaller cases have rates of 0 to 100%. My thoughts are that the size of the outbreak should weight the outbreak's impact on the overall rate. However, I don't think averaging the individual rates would be correct. Is there a more appropriate measure that would allow me to take the size of the outbreak into account when calculating the overall impact?


title: Anomaly Detection with Dummy Features (and other Discrete/Categorical Features)

text: ## Intro This is my first time posting on here, so please, if anything doesn't seem technically correct, either in the formatting, or the use of correct definitions, I'm interested to know what should've been used instead. Onwards. I've recently been taking part of the Machine Learning class by Andrew Ng For anomaly detection we've been taught to determine what the Normal/Gaussian distribution parameters are for a given feature/variable, ${x_i}$ within a data set, and then determine the probability of a chosen set of training example's/observation's value given that particular Gaussian distribution, and then taking the product of the probabilities of the features. ## Method Choose $x_i$ features/variables that we think explain the activity in question: $$\{x_1, x_2,\dots,x_i\}$$ Fit the parameters of the Gaussian for each feature: $$\mu_j = \frac{1}{m}\sum_{i = 1}^m x_j^{(i)}$$ $$\sigma^2 = \frac{1}{m}\sum_{i = 1}^m (x_j^{(i)} - \mu_j)^2$$ For each training example, $x$, compute: $$p(x) = \prod_{j = 1}^n \ p(x_j; \mu_j, \sigma_j^2)$$ We then flag as an anomaly ($y = 1$), given: $$y = \left\{ \begin{array}{l l} 1 & \quad p(x) < \epsilon\\ 0 & \quad p(x) \geq \epsilon \end{array} \right.$$ This gives us the method with which to determine if an example requires further inspection. ## My Question(s) _This seems fine for continuous variables/features,_ but discrete data is not addressed. What about dummy variables, e.g. a gender flag feature, possibly called `[IsMale]` that can be of the value $0, 1$? To take a dummy feature into account would we use the binomial distribution instead to calculate $p(x)$? What about categorical data such as car colour? While we could map colours to numerical values, e.g. $red \to 1, blue \to 2$, the distribution of such a categorical feature could be close to uniform (i.e. equally likely chance to be any of the colours), and further, as any numerical mapping that occurs (i.e. $red$ having the value $1$, etc) is not ordinal, does it make sense to try and transform any non-normal distribution of frequencies for colours to be normally distributed (does it even matter that it is not ordinal??)? For example, to me, it wouldn't make sense to do a $log()$ transform as the data is neither continuous nor ordinal. So perhaps it would be best to find a discrete distribution that fits the feature, as opposed to "torturing" the data to fit the Gaussian? **Questions:**   * Can binary variables be modeled with a binomial probability distribution and become another factor in the $p(x)$ calculation?   * Should categorical variables should be modeled with a discrete probability distribution instead of a Gaussian, and become another factor in the $p(x)$ calculation?   * Is there another method altogether that takes into account what I'm asking here that I can further research/learn about?


title: Finding average number of rolls for recursive rolls on success

text: If there is more appropriate terminology to better phrase my question, please edit. Imagine I have a 3 sided die: two faces marked 0, one face marked 2 Whenever a two is rolled, two more rolls are to occur. How do I determine the average number of rolls that will occur? I tried figuring this out based on an explanation of exploding dice (http://axiscity.hexamon.net/users/isomage/rpgmath/explode/) but I'm not sure I'm even headed in the right direction.               p(number of rolls)          p(1) = 1/3     p(3) = 1/3 * 2/3 * 2/3     p(5) = 1/3 * 1/3 * 2/3 * 2/3 * 2/3     p(7) = 1/3 * 1/3 * 1/3 * 2/3 * 2/3 * 2/3 * 2/3     p(n) = n+1/3^n      Solve for the sum of p(n) ??? Are those probabilities correct or do i have to include the alternate ordering (ex: p(5) = 1/3 * 2/3 * 1/3 * 2/3 * 2/3) what about a 3 sided die: 0,1,2 where 1 = +1 roll and 2 = +2 rolls? what about a 20 sided die: 15 zero faces, 3 two faces, 1 three face, and 1 five face? Curious why in the heck I care? In an effort to learn python i thought it'd be fun to generate a random dungeon. Now i'm curious about the probabilities of the random (or should I say psuedo-random) generation. I tried digesting the question/answers linked below ... snow crash! If you do help me, please, please, feed me via spoon. I am a statistics baby. How to easily determine the results distribution for multiple dice? Help with probabilities on a game I am making Now i'm going to read http://www.diku.dk/hjemmesider/ansatte/torbenm/Troll/RPGdice.pdf while I wait for enlightenment Thanks in advance for said enlightenment.


title: Why is this result not significant?

text: I am reading a research paper about the relationship between service quality and consumer behavior. One hypothesis is: customers experiencing no service problems have the strongest levels of favorable behaviors and the smallest levels of unfavorable behaviors. ANOVA was operated, results are: ![enter image description here](http://i.stack.imgur.com/00Hut.png) The author intepreted: "The pay more intentions of customers experiencing no service problems are not significantly higher than those of customers experiencing service problems that were resolved." But the significance level for this factor was ,036, less than ,05. It should be significant, right? (The author used ,05 as standard for other analysis )


title: Estimate variance of sub-sets from overall variance

text: I am looking for a way to estimate the variance of a summed sub-set based on the variance of those sums.               Si = sum( Ai )     S = { S0...Sn }     V = variance( S )      That is, each sample value `Si` is actually the sum of a set `Ai`. `V` is then the normal variance of the set `S`. Is there any way to estimate the variance of the values in sets `Ai` using the final variance, or some result from the above? _Currently I do`variance / (len(Ai)^2)` but the value is then an averaged value of the variance. I'm wondering if there is a way to improve upon this._ * * * _Background:_ The reason the I don't calculate variance on individual samples from the population is that the measurement is then too small. I am measuring elapsed time, and so individual samples the elapses time is too small for the timing device to measure, so instead I do a series of operations and calculation the time for the whole group. _Note:_ In the application I have no ability to store samples and am using an online algorithm (Knuth).


title: Feature selection for one class SVM

text: I have around 300 features, i need to choose features for one class svm. can some one tell me the ideal algorithm for this use case. I know about that for feature selection regularised random trees/SVM/correlation can be used. But I am not sure which specific algorithm would be useful for one class svm. Please guide. Thanks a lot.


title: Are there any tutorials on Bayesian probability theory or graphical models by example?

text: I've seen references to learning Bayesian probability theory in R, and I was wondering if there is more like this, perhaps specifically in Python? Geared towards learning Bayesian probability theory, inference, maximum likelihood estimation, graphical models and the sort?


title: How to interpret results if a reference category of a categorical variable in multivariable logistic regression is not significant?

text: I am trying to do a multivariable logistic regression and using a normal binomial logistic regression, using binomial variable `X` (coded `0`/`1`) and `Y` (4 sub-categories). Unfortunately, the reference category was shown to be not significant. Should I then reject the null hypothesis for this variable?               Cat.,   Sig.,   Exp(B)     Y,         1,   .161     Y(1),   .032,   .226     Y(2),   .106,   .370     Y(3),   .061,   .309


title: Kohonen SOM: Maximum number of neurons for any given dataset

text: I have an unknown data set consisting of 1024 vectors (integers with `1024 rows`, `30 columns`). What is the `maximum` number of neurons (`x_dimension * y_dimension = number_of_neurons`) for a 2-dimensional Kohonen Self Organizing Map for this dataset (without having any knowlege what this dataset is about)? Any explanation _why_ would be appreciated. My guess is that the maximum appropriate map size would be 1024 neurons (32×32 neurons). Much larger than that would create a sparse SOM with many empty regions between the classified objects, whereas a too-small map would force objects to share the neurons to make groups. Any ideas? One paper describes an empirical guideline for the optimum map size (no. of neurons) for any given dataset example: `number_of_neurons = 5 * sqrt(1024*30)`. Also, is there a mathematical relationship between the increase in number of neurons and the decrease in quantization error? Thanks for any help on this.


title: Why naive bayes classifier cannot perform correctly on 1 conditional probability?

text: I use e1071 package for naive bayes classification. My dataset is having 2 condition/feature (FF1 and FF2) to predict a class of signal (ELO or FP). When I used the following line of code:               A_classifier <- naiveBayes(A_train[,c(1:2)], A_train$Label)      it works fine and the conditional probability display               Conditional probabilities:                  FF1     A_train$Label     [,1]     [,2]               ELO 506.8857 95.28520               FP  618.1716 64.53944                       FF2     A_train$Label     [,1]     [,2]               ELO 1533.629 324.8894               FP  1361.693 144.5700      However I want to test and train the classifier with only FF1 as single feature. But the result is quite weird.               Conditional probabilities:                  x     A_train$Label     [,1]     [,2]               ELO 506.8857 95.28520               FP  618.1716 64.53944      Isn't the conditional probability should display FF1 to 'x'? Where am I doing wrong.?I am a new beginner in R and NB. Please help me..TQ


title: Factor analysis with binary variables using Stata 13

text: I am trying to do confirmatory factor analysis on data that is coded binary (`0` no, `1` yes). UCLA suggests using a tetrachoric correlation matrix, which, however, assumes that binary variables reflect underlying continuous variables. Is there a way to relax this assumption and use logistic procedures instead? I am using Stata Version 13. Could generalized SEM provide a solution?


title: Regression with exponentially-smoothed errors

text: I'm just starting to look into exponential smoothing models. Is there a way to fit a linear regression with exponentially-smoothed errors, similar to the standard technique of fitting a regression with ARMA errors? Or is this generally considered a bad idea? Of course, I can come up with heuristics, like fit with AR(1) errors, then exp-smooth the residuals (possibly on the differences, then look for a trend) but that feels very ad hoc.


title: Mixture model and Pymix (python package for mixture models)

text: I have a data set that behaves approx. Standard normal. It is an image where each observation is a pixel intensity. I want to cluster this into three different sets by fitting a 3-gaussian mixture model. I want one cluster to be "significantly" smaller than the mean (0) one to be around the mean and one "significantly" larger than the mean. I find the initial moments and mixture coefficients by determining intial sets and calculatins the corresponding sample means and std's. This is:   * Set1:=(x|x smaller than -2std)   * Set2:=(x|x in (-0.5std,0.5std))   * Set3:=(x|x larger than 2std) Taking the sample means and std's of these sets I can estimate initial moments and mixture coefficients to later use the EM algorithm to get 'optimum' parameters for the 3-gaussian mixture. Using the pymix package I do:               from pymix import mixture          # standard deviation of whole data set     standev= np.std(x,ddof=1)          # Initial Sets for Negative-Change, No-Change and Positve Change         CN = x[x<(-2*standev)]     NC = x[ (x > (-1*standev)/2) & (x < (standev/2))]      CP = x[x>(2*standev)]               # Initial Means and Variances        CNimean = np.mean(CN)     CNistd = np.std(CN)          NCimean = np.mean(NC)     NCistd = np.std(NC)          CPimean = np.mean(CP)     CPistd = np.std(CP)          # Load data for mixture model     data = mixture.DataSet()          data.fromArray(np.squeeze(np.asarray(x)))          # Initialize mixture model          n1 = mixture.NormalDistribution(CNimean,CNistd)     n2 = mixture.NormalDistribution(NCimean,NCistd)     n3 = mixture.NormalDistribution(CPimean,CPistd)          m = mixture.MixtureModel(3,[0.01,0.98,0.01], [n1,n2,n3])          # Resolver el Mixture model          c = m.EM(data,50,0.1)      What I'm no sure on how to estimate is the proportion. For pymix's em to solve this mixture model you also need to supply an initial guessing of the proportion of observations that belong to each of the three gaussians (sum = 1). And after trying some judicious guessing, for example 1% 98% 1% as the extremal sets are before -2stds and after 2stds I end up with a mixture having for example 0.03 0.02 and 0.95 in that order, and with parameters that dont make a lot of sense. Does anyone know what I might be doing wrong? or is my approach wrong all over? Any help would be deeply appreciatd! Greetings!


How do I interpret hierarchical regression results?


title: Implementing MCMC

text: I'm writing an MCMC algorithm in R and I'm wondering about the following: say we have two parameters, $\theta_1$ and $\theta_2$. I want to update each one at a time from the corresponding posterior conditional distributions. Say $\theta_1^{(0)}$ and $\theta_2^{(0)}$ are the initial values. Then at iteration 1 I update first $\theta_1$ from $\theta_1^{(1)} \sim f(\theta_1 | \theta_2^{(0)}, \text{Data})$ Now when updating $\theta_2$, I should use $\theta_2^{(1)} \sim f(\theta_2 | \theta_1^{(1)}, \text{Data})$ Is there a theoretical problem if when updating $\theta_2$ I use $\theta_1^{(0)}$ I instead of $\theta_1^{(1)}$? That is $\theta_2^{(1)} \sim f(\theta_2 | \theta_1^{(0)}, \text{Data})$?


title: Comparing the ratios of occurrences

text: I am trying to find a statistical test for comparing the ratios of occurrences and would appreciate any help. There are two groups; treatment and control. Each group has three categories of outcomes; High, medium and low blood pressure. The number of subjects is 100 for each group. Now, let's say, for control group, the number of observations is 40, 40, 20 for high, medium, and low blood pressure, respectively. And, for treatment group, the number of observations is 10, 50, 40 for high, medium and low. What I want to test is whether the ratio of occurrences (High/Low) ignoring the medium is different between these two groups. 2 (=40/20) for control vs. 0.25 (=10/40) for treatment. What statistical test should I use? Many thanks.


title: Estimating State Space Model in R with MARSS package and shared parameters between Q and R

text: I am trying to estimate the following unobserved components model using the MARSS package $y_t = \mu_t + \varepsilon_t $ $\mu_t = \mu_{t-1} + \beta_{t-1}$ $\beta_t = \beta_{t-1} + \zeta_{t}$ with a restriction on the variance of the errors: $\sigma^2_{\zeta} = \frac{1}{\lambda}\sigma^2_{\varepsilon}$ for a given $\lambda$. $y_t$ are the observed series and $\mu_t$ and $\beta_t$ are the states. Those of you familiar with these models will recognize this as the Hodrick-Prescott filter. MARSS handles models in this form $y_t = Z x_t + v_t , v_t \sim MVN(A,R) $ $x_t = B x_{t-1} + w_t, w_t \sim MVN(U,Q) $ $x_0 \sim MVN(x0,V0)$ So I tried this               lambda=1600     B1 = matrix(c(1,0,1,1),2,2)     U1 = matrix(0,2,1)     Q1= matrix(list(),2,2)     Q1[1,1]=0     Q1[1,2]=0     Q1[2,1]=0     Q1[2,2]="q11*1/lambda"     Z1 = matrix(c(1,0),1,2)     A1 = as.matrix(0)     R1 = as.matrix("q11")     pi1 = matrix(0,2,1)     V1 = diag(1,2)     model.list = list(B = B1, U = U1, Q = Q1, Z = Z1, A = A1, R = R1, x0 = pi1, V0 = V1)     fit = MARSS(data, model = model.list, control = list(kf.x0 = "x00"))      But the specification for Q[2,2] is wrong and MARSS reads it as a new parameter to estimate. I know I can restrict parameters to be the same within the Q matrix and R matrix, but I don't know how to impose the variance restriction that affects parameters of both matrices. Does anyone know a way around this? Maybe using substitute() ? I know there are packages to implement the HP filter as well as other state space estimation packages but I want to use MARSS since it uses the EM algorithm to estimate the parameters and I want to extend this model to others where convergence is harder to achieve. Thank you.


title: How to test whether linear models fit separately to two groups are better than a single model applied to both groups?

text: My question is how to tell if two regressions explain the data better than one. Let me be more concrete with an example (which I'm making up as I go, it's not meant to be plausible). Say I'm interested in the relationship between income and happiness. Then I ask a lot of people about their income and happiness. But say some of the people I ask are men and some are women. Now, I perform a linear regression for income and happiness and see that income does predict happiness according to some coefficients, it's significant, the model explains a lot of the variance, etc. My question is, how do I test whether there's a unique relationship between happiness and income, that applies across men and women, or the data are better explained by two models, each obtained by regressing separately on men and women? I'm thinking of calculating residuals for the data with respect to each model, then comparing the sum of the residuals, or something along those lines. Is that a reasonable thing to do? How would I go about comparing the residuals? So, in summary, how do I know if the better hypothesis is that of one unique linear model for all the data, or the one where each subset (subsets defined a priori) of the data has its own linear model?


title: How to interpret a multi-variable regression model?

text: How do I interpret this predictive model? > S&P500 Next Week = Earnings Per Share * 0.25 + Dividends Per Share * -1.519 > + Current S&P500* 1.018 + 5.409


title: Generating likely populations given a subsample and control totals

text: **Context** In transportation planning, agent-based microsimulation is a method to deal with the complexity of the problem. Instead of computing aggregate flows (as in the classical four-step model), agents (persons, households, businesses, ...) are simulated on an individual basis. Results are obtained by aggregating after the simulation. For the US, data sources for the simulation include the PUMS (Public-Use Micro Sample), an anonimyzed 5% sample from the full population census, and aggregated control toals (also spatially). The population census contains mostly demographic data. For other countries, similar data sources exist. This data cannot be fed straight away into a microsimulation model, as it requires   * Ideally, a full sample with demographic variables   * Spatial information for each agent   * Adherence to control totals This problem has received some attention in the recent literature, starting with Beckman et al.'s 1996 paper (reference below). However, my feeling is that this is often seen as an engineering rather than a statistical analysis task. I am aware of this related question, but the answer doesn't seem to fit my needs. **Questions**   1. While inferring from a sample is a common statistical task, inferring a _likely full population_ doesn't seem to be that common. In which scientific fields could this (or a similar) problem be also relevant?   2. Which methodology/techniques does one need to know in order to perform a maximum-likelihood estimation of a population given sample? I have peeked into entropy maximization, bootstrapping and Gibbs sampling/MCMC, but I'm missing a feeling for which technique would work best, and I don't have a deep statistics background. (Actually, near-maximum likelihood should be fine, too. Eventually, fractional weights will have to become integers during the process, and this is computationally expensive if it has to be done exactly.)   3. Which approach is best suitable to validate the synthetic population, other than running the microsimulation model? Should I reserve a validation data set, or would some sort of cross-validation be suitable in this case? Again, I am missing some background here, so I'd appreciate any hints.   4. Do I really need to carry out statistical equivalence testing for the validation (as described by Wellek), or "will it do" to fail to detect a difference? I am aware that the latter is incorrect from a statistical point of view, but I also was unable to apply Wellek's methods to my data. I appreciate any hints, remarks, and literature references. * * * Beckman, R.J., Baggerly, K.A., & McKay, M.D. (1996). Creating synthetic baseline populations. Transportation Research-A, 30(6):415-429. Wellek, S (2010). Testing Statistical Hypotheses of Equivalence. Chapman and Hall/CRC.


title: How to visualise "equal in distribution" in the context of stochastic dominance?

text: I am reading about stochastic dominance in wikipedia. I don't comprehend the meaning of "equal in distribution". $x_B \overset {d}{=} (x_A+y)$ **What does $\overset {d}{=}$ mean?** I understand that equal in distribution means that they have the same distribution functions. $P[X \le x] = P[Y \le x]$ and two random variables having equal moment generating functions have the same distribution. But I have difficulty visualizing this. **Is there a nice way to visualize this?**


title: Error term in repeated measures ANOVA in R

text: I am having problems with defining the error terms for a two way repeated measures ANOVA in R. My data consists of wood density estimates for three radial positions (inner, middle, and outer) along a core extracted from a tree. There are a total of 20 tree species, 6 individuals of each species, and two cores from each tree. To test the effect of radial position on wood density, I use the following two way ANOVA model with an error term that accounts for variability within individuals:               radpos.aov<-aov(WD~Species*Radialposition + Error(Individual), data=Radpos)      However, I am not sure if my specification of the error term is adequate. Am I also supposed to account for the variability within a core? To me this variability is the same as that due to radial position which is the main factor I am interested in. Although I have devoted some time to reading about specifying the error term in Repeated measures ANOVA, I still have problems with practically specifying the error term. I will appreciate some help with this.


title: What model should I use in this situation?

text: I have income based on zipcode from US census data, I think they also provide income based on gender for whole of US, income based on age for whole of US, and income based on zip code. Now I need to find out what the income of a person from a certain zipcode, of a certain age and gender. What model should I use? How should I go about it?


title: Log transformation in CART analysis

text: I'm working as a liaison between a researcher and a stats team at a university. I'm a database admin who is working on using business intelligence tools to offer the option of (as of right now) offer CART analyses in .PDF form over the internet. I'm using Microsoft Business Intelligence Development Studio to do this. Recently I received an email saying that the measure (or independent variable I guess it should be called here) was highly skewed and the original CART analysis didn't make sense. Fair enough, I know that sometimes these things don't work out like you had hoped. Then, another analysis was done with a log transformation. Can someone point me to a good quick and dirty resource so I can read up on this? I'm not a complete stats novice, but ANOVAs and logistic regressions are about as far as my knowledge extends.


title: Why would we want to use average cost as response but not total cost?

text: Assume I have a data set as follows:               Total_Cost  Average_Cost  Quantity         12           6            2         21           7            3         28           7            4         30           6            5         36           6            6      The linear relationship between Total_Cost and Quantity is nearly perfect, while between Average_Cost and Quantity is very poor. However, I saw many researches that basically use average cost as a response but not total cost. My question here is therefore: Are there any particular reasons for us to study the linear relationship between Average_Cost and whatever independent variables? (For instance, the use of Total_Cost as a response is way too straightforward.)


title: Exogeneous regressors in auto.arima and using them in forecast function in R

text: I'm trying to forecast a seasonal time series based on its historical values, and also two more time series (that are seasonal themselves.) I'm trying to use an **auto.arima** , and I'm going to input the other two time series (the exogeneous regressors) as a contatenated list of dummy variables, in auto.arima's **xreg** parameter. I am having difficulty how to use the forecast function after this point. I've written up the following code, but I don't understand what I should put in the **xreg** and **newxreg** parameters of the forecast function.               tempfit<-auto.arima(dnew, xreg=dExt)     plot(forecast(tempfit, xreg=dnew1,newxreg=dExt1))      Also, my data points for these three series were all values per day that had a seven day seasonality. In order to let auto.arima calculate the (p,q,d) for seasonality, I converted them to time series with a frequency of 7. Now, after forecasting is done, the plot shows one unit for every seven days. How can I covert this back to one unit per day? Further, do you happen to know how we can input a set of external regressors to an ETS model? I would greatly appreciate your inputs! Thank you. **EDIT** : I just saw the following page from Dr. Hyndman: Time series modeling with dynamic regressors in SAS vs. in R Is it safe to assume that I don't need to enter a newxreg parameter for my forecast? Also, I really want to know if it's statistically correct to use the two external regressors in xreg, but then also use a number of dummy variables in xreg that will represent the seasonality of these two variables.


title: Reporting percentages in preference type survey questions

text: I came a cross some survey results that show user preferences, asked in the form of the most and least preferred choices. The result reported was something like "x % of people chose option A as their most preferred" and "y % of people chose option B as their least preferred". It seems to me that there is some dependency between what someone would prefer to what someone would least prefer, so instead they should show the breakdown of least preferred options for people who chose option A as their most preferred option (and vice versa). I want to ask whether presenting the percentages as a total only is misleading.


title: Error function of noisy input and target variables

text: Why is the sum of squares error function of noisy input and noisy target variables very similar to the error function for only noisy input? This is the relevant part in Bishop's book: > Another viewpoint on kernel regression comes from a consideration of > regression problems in which the input variables as well as the target > variables are corrupted with additive noise. Suppose each target value > $t_{n}$ is generated as usual by taking a function $y(z_{n})$ evaluated at a > point $z_{n}$, and adding Gaussian noise. The value of $z_{n}$ is not > directly observed, however, but only a noise corrupted version > $x_{n}=z_{n}+\xi_{n}$ where the random variable $\xi$ is governed by some > distribution $g(\xi)$. > > $$E=\frac{1}{2}\sum_{n=1}^{N}\int > \{y(x_{n}-\xi_{n})-t_{n}\}^{2}g(\xi_{n})d\xi_{n}$$ However, the error function when we consider a noisy input variable is extremely similar: $$E=\frac{1}{2}\sum_{n=1}^{N}\int \{y(x_{n}+\xi_{n})-t_{n}\}^{2}\nu(\xi_{n})d\xi_{n}$$ Why?


title: Why are drug trials designed with discrete treatment values?

text: I'm not involved in drug testing, but it is my understanding that most trials are performed using a control and a treatment group (or perhaps a few discrete treatment groups). This data is then used to determine effectiveness of a drug, and hopefully uncover any deleterious side-effects before the drug is released on the population. I'm wondering why it would not be more sensible to use a continuous distribution of drug dosages, rather than a discrete distribution? For example, give 20% of patients the control, and the give the other 80% a randomised (uniform?) dosage between 0% and ~110% of the expected optimal dosage. It would seem to me that this would provide more information than data that is inherently discretised, as well as potentially making it easier to more accurately fit models like logistic functions for both effects and multiple side-effects that occur at different rates? The obvious draw-back would be the increased cost of manufacture for the trial, but that would be fairly minimal in comparison to potential benefits, especially for large trials. Also, this method could just as easily be applied to any "treatment" experiment (e.g. in agriculture, metalurgy, etc.). So why aren't more trials conducted like this?


title: How to adjust machine learning training data set with time

text: I'm using machine learning to do text classification right now, I first use a training data to train my classifier, then use this classifier to classify text document into different classes. With the time goes on, some new words will arise and I need to adjust my training vocabulary to correctly classify these documents, could you please tell me how to do this?


title: Measuring the performance of Logistic Regression

text: Being quite new to the field, it occurs to me that there are multiple and fundamentally different ways of assessing the quality of a logistic regression:   * One can evaluate it by looking at the confusion matrix and count the misclassifications (when using some probability value as the cutoff) or      * One can evaluate it by looking at statistical tests such as the Deviance or individual $Z$-scores.  Although I more or less understand the aforementioned tests, I am curious to know when it is appropriate to use the misclassification count and when one should apply more statistical tests?


title: Likelihood for negative binomial distribution

text: One of the parameterization of the negative binomial distribution is NB($m,r$), $$\Pr(X = k) = \left(\frac{r}{r+m}\right)^r \frac{\Gamma(r+k)}{k! \, \Gamma(r)} \left(\frac{m}{r+m}\right)^k \quad\text{for }k = 0, 1, 2, \dots. $$ I would like to consider the parameterization NB($m$, $\phi$) where $\phi=1+\frac{m}{r}$. I would like to find the maximum likelihood estimate and its standard error of $m $and $\phi$. Can I just use existing R package which finds the MLE for $m, r$ and just set $\hat{\phi}=1+\hat{m}/\hat{r}$. Or do I have to use some profile likelihood to fix one of the parameters.


title: Calculating the probability of 31 of 628 items are sampled (no replacement) more than 10x amongst 150 participants drawing 50 items each

text: I am sorry for having to ask a simple probability question, but I have been thinking about it for weeks and an extensive google search has given no answers. I have a group of 628 questions. 31 Questions belong to category A. 150 participants will each be given 50 randomly sampled questions from the total 628 questions (no replacement). What is the probability that questions belonging to category A will be answered more than 10 times by the group of participants? I looked at a binomial distribution using R's binomial density function, but the best I could come up with was the probability of 1 item from category A being administered to the 150 participants if only 1 question was sampled from the total of 628 questions.                   x <- seq(1,50,by=1)         high.biDen <-dbinom(x,size=150,prob=((31/628)))         round(high.biDen,2)              round(sum(high.biDen[10:50]),2)      The answer obtained is .208 or a 20.8% probability that a question from category A will be administered to the participants more than 10 times if 1 question is sampled from the pool of 628 questions. I would like to know the probability if 50 questions are sampled instead of 1. Is the answer obtained by some sort of bayesian technique? It seems to rely heavily on conditional probability seeing how if one item is sampled only 30 remain to be randomly selected. Thanks so much for sharing your knowledge and expertise! -Xander


title: Dummies only apply to a subgroup of the data

text: I have a dataset covering small and large companies. The three dummy variables are "large manufacturing", "large service", and "large governmental" companies. Roughly 90% of the firms can be classified by these categories (the remaining firms are small and are not classified in manufacturing, service, and government). Now, I would like to include the dummies for manufacturing and service in my regression models. The problem is, the base category is not "large government firms" but rather large government firms and small firms. Would it be ok to just term the base category other and compare against it? Thanks!


title: Multiple Imputation Using Different Data Sets

text: I realize that this is contrary to "standard" multiple imputation using Rubin's rules, but I have a data set were a variable with a lot of missing values is most likely MNAR. It represents a large portion of the data (40%), so I would like to fill it in, but doing so using the mi command in Stata rests on faulty assumptions. My thought is to build an imputation model on the census data of the population at hand, and then use that model to impute values in my dataset. It might work better because the census represents the same population but with a more accurate distribution of values. Thoughts?


title: How to calculate the first principal component from multi-data sets (using WGCNA)

text: I have a bunch of datasets (gene expression profiles).   From the datasets, I am gonna analyze co-expression profiles using WGCNA (http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/).   Actually, briefly speaking, I   1) got specific gene-centric PCC from each dataset,   2) convert it into fisher-Z statistic,   3) calculate the mean effect size of each gene from all datasets, and   4) re-converted the value into correlation coefficient of each gene.   This correlation coefficient serves as input value for topological overlap (such as adjacency value in WGCNA tutorial).   My question is about how to calculate first PC from multi-data sets.   FYI, in WGCNA, eigengenes is calculated like this.               PC = moduleEigengenes(datExpr)     ME_1A = PC1A$eigengenes


title: Sparse variable selection algorithms that account for parameter changes

text: I am variable selecting for a time-series forecasting model that has parameters sampled from a high variance sampling distribution centred near zero and that undergo changes over time. Each predictor has a tiny but non- zero correlation (in population terms) with the response. Lasso et al. are problematic here since the correlation sampling distribution is so wide, and correlations change over time; features that are crappy out of sample get selected. What algorithms are there that I can use (apart from manual hacks)?


title: Variance of sample mean of bootstrap sample

text: Let $X_{1},...,X_{n}$be distinct observations (no ties). Let $X_{1}^{*},...,X_{n}^{*}$denote a bootstrap sample (a sample from the empirical CDF) and let $\bar{X}_{n}^{*}=\frac{1}{n}\sum_{i=1}^{n}X_{i}^{*}$. Find $E(\bar{X}_{n}^{*})$ and $\mathrm{Var}(\bar{X}_{n}^{*})$. What I have so far is that $X_{i}^{*}$ is $X_{1},...,X_{n}$ each with probability $\frac{1}{n}$ so $$ E(X_{i}^{*})=\frac{1}{n}E(X_{1})+...+\frac{1}{n}E(X_{n})=\frac{n\mu}{n}=\mu $$ and $$E(X_{i}^{*2})=\frac{1}{n}E(X_{1}^{2})+...+\frac{1}{n}E(X_{n}^{2})=\frac{n(\mu^{2}+\sigma^{2})}{n}=\mu^{2}+\sigma^{2}\>, $$ which gives $$ \mathrm{Var}(X_{i}^{*})=E(X_{i}^{*2})-(E(X_{i}^{*}))^{2}=\mu^{2}+\sigma^{2}-\mu^{2}=\sigma^{2} \>. $$ Then, $$E(\bar{X}_{n}^{*})=E(\frac{1}{n}\sum_{i=1}^{n}X_{i}^{*})=\frac{1}{n}\sum_{i=1}^{n}E(X_{i}^{*})=\frac{n\mu}{n}=\mu $$ and $$ \mathrm{Var}(\bar{X}_{n}^{*})=\mathrm{Var}(\frac{1}{n}\sum_{i=1}^{n}X_{i}^{*})=\frac{1}{n^{2}}\sum_{i=1}^{n}\mathrm{Var}(X_{i}^{*})$$ since the $X_{i}^{*}$'s are independent. This gives $\mathrm{Var}(\bar{X}_{n}^{*})=\frac{n\sigma^{2}}{n^{2}}=\frac{\sigma^{2}}{n}$ However, I don't get the same answer when I condition on $X_{1},\ldots,X_{n}$ and use the formula for conditional variance: $$ \mathrm{Var}(\bar{X}_{n}^{*})=E(\mathrm{Var}(\bar{X}_{n}^{*}|X_{1},...,X_{n}))+\mathrm{Var}(E(\bar{X}_{n}^{*}|X_{1},\ldots,X_{n})) \>. $$ $E(\bar{X}_{n}^{*}|X_{1},\ldots,X_{n})=\bar{X}_{n}$ and $\mathrm{Var}(\bar{X}_{n}^{*}|X_{1},\ldots,X_{n})=\frac{1}{n^{2}}(\sum X_{i}^{2}-n\bar{X}_{n}^{2})$ so plugging these into the formula above gives (after some algebra) $\mathrm{Var}(\bar{X}_{n}^{*})=\frac{(2n-1)\sigma^{2}}{n^{2}}$. Am I doing something wrong here? My feeling is that I am not using the conditional variance formula correctly but I'm not sure. Any help would be appreciated.


title: Design a Data Model to predict value of sum of Function

text: I am working on a data mining projects and I have to design following model. I have given 4 feature x1, x2, x3 and x4 and four function defined on these feature such that each function depend upon some subset of available feature. e.g. F1(x1, x2) =x1^2+2x2^2 F2(x2, x3) =2x2^2+3x3^3 F3(x3, x4) =3x3^3+4x4^4 This implies F1 is some function which depend on feature x1, x2. F2 is some feature which depend upon x2, x3 and so on Now I have a training data set is available where value of x1,x2,x3,x4 is known and sum(F1+F2+F3) { I know total sum but not individual sum of function) Now using these training data I have to make a model which which can predict total sum of all the function correctly i.e. (F1+F2+F3) I am new to the data mining and Machine learning field .So I apologizes in advance if this question is too trivial or wrong. I have tried to model it many way but I am not getting any clear thought about it. I will appreciate any help regarding this.


title: How to do model selection in dynamic linear model?

text: I am trying to use DLM to model a time series. Candiate model includes local level, local trend and local trend with seasonal part. I do not know how to do model selection. Can AIC be calculated? I found no function in the R package `dlm`.


title: Overfit by removing misclassified objects?

text: Actually this question may be simple for you, but I need to learn the correct answer.   If I remove misclassified instances from data set with Naive Bayes (it gives minimum FP rate) and then train logistic classifier on this filtered data set, will it overfit or not? Thanks in advance.


title: Selecting individuals from a population using a binary classifier

text: I have a dataset consisting of around 200 individuals, whose outcome is either of state $0$ or $1$. I am able to make binary classifiers and predictors on this set and build ROC-curves for them just fine, but now I'm interested in a problem which has me kind of stumped. Suppose I have some even larger population of individuals, of size say $N = 500$ or $1000$. Assuming that they have the same distribution as my original population and given the classification models I developed, I'm interested in selecting the $100$ "best" individuals, i.e. who are deemed as most likely to be of state $1$. Finally, I want to calculate the resulting precision (the proportion of selected individuals who are predicted to indeed turn out to be of state $1$). What I'm uncertain about is how such a selection process would look like computationally. I'm thinking I could maybe bootstrap populations of size $N$ from my data, run my models on those and study the resulting confusion matrices given a threshold value (a false positive is "worse" than a false negative in my case, so I will want a threshold value that minimizes FPR to a reasonable extent). This is unlikely to give me exactly $100$ selected individuals, however. Can I force the model to pick exactly $100$ somehow? Is there a better and more precise method available in this case than approximating through simulation like this? Any help would be appreciated.


title: Data imputation question

text: I have a variable with some missing values               a <- rnorm(100);     a[sample(1:100,10)] <- NA;     a;      How can I fill missing values with previous non missing value? for example if I have sequence: a<- (3, 2, 1, 6, 3, NA, 23, 23, NA); first NA should be replaced by first previous non NA number 3, second NA should be replaced with 23 etc. Thanks


title: How to deal with a specific case of unbalanced within-subjects design?

text: I have a completely within-subjects design with 3 independent variables:   * Trial type (3 levels)   * Task order (2 levels)   * Modality (3 levels) However, for one of my levels of Trial type, the Task order and Modality levels are redundant (because it is essentially a baseline measurement). Ideally, I'd like to run a 3 (trial type: 1, 2, or 3) x 2 (task order: a or b) x 3 (modality: vm, av, hm) repeated measures ANOVA. Do you know of a way around this? If I were to duplicate the Task 3 data across 'modalities' and 'task orders' I would have the correct number of levels for each within-subjects factor (i.e. six identical columns of data to represent imaginary modalities and task orders for Task 3). I'm assuming this violates an assumption but I'd be worried about looking at a similar analysis with a huge number of paired sample t-tests. Is lots of t-tests the way to go? Thanks for your help!


title: In Bayesian statistics, does variance parameter actually mean 1/variance (i.e. precision)?

text: Say, $\psi | Y, {Y}^{*}, {\sigma}^{2}$~$N(\mu, {\Lambda}^{-1})$ this is a bivariate normal Does the ${\Lambda}^{-1}$ actually mean precision or variance. I was told in some Bayesian papers, the second parameter actually means precision. I didn't think so, but maybe someone else has read more paper would know? The paper I am reading www.jds-online.com/file_download/353/JDS-746.pdf


title: Finding confidence intervals for the click-through-rate of a website

text: This is a very basic question, so please bear with me. I've been learning about AB Testing, which is largely used in internet marketing to examine the effectiveness of certain aspects of ads, websites, etc. Here's a couple links to people who want to know more about AB Testing: http://visualwebsiteoptimizer.com/split-testing-blog/what-you-really-need-to- know-about-mathematics-of-ab-split-testing/ http://20bits.com/articles/statistical-analysis-and-ab-testing/ http://elem.com/~btilly/effective-ab-testing/ Let's say that I have a website that registers users for a forum. I want to know if Headline 1 or Headline 2 is more effective at getting visitors on the web site to register for the forum. So I have the following data.               dat = data.frame(Headline=c("Headline 1", "Headline 2"),                      Visitors=c("1000", "1300"),                     Clicks=c("500", "600"),                     Conversions=c("100", "150"))       And here are the click through rates and conversion rates for each of the headlines.               ctr1 = (500/1000)*100  # for headline 1     ctr2 = (600/1300)*100  # for headline 2     ctr1; ctr2          conv1 = (100/1000)*100  # for headline 1     conv2 = (150/1300)*100  # for headline 2     conv1; conv2      According to the sites above, I'm really interested in determining the confidence intervals for the conversion rates for each headline. While 95% confidence would be ideal, I'm really open to anything 80% and up, so I need to calculate confidence intervals where I am 80%/85%/90%/95% confident that the conversion rate for a headline is within a certain range. I'm really not sure how to go about this. Are there specific tests and/or functions in R that with provide me with the appropriate information? (confint, chisquare, gtest, etc.?) Thanks for your patience and help. EDIT: So I tried the following, and I'm not sure if I'm doing it properly or making the right conclusions. Furthermore, there has to be a more efficient way to perform this task in R. For a given conversion rate (p) and number of trials (n):               p1 = 0.1     n1 = 1000     se1 = sqrt( p1  * (1-p1) / n1 )     se1     se1 * 1.96     (p1 + 1.96*se1) * 100     (p1 - 1.96*se1) * 100                      p2 = 0.11     n2 = 1300                 se2 = sqrt( p2 * (1-p2) / n2 )     se2     se2 * 1.96     (p2 + 1.96*se2) * 100     (p2 - 1.96*se2) * 100          (8.1, 11.8) # headline 1     (9.2, 12.7) # headline 2     # these confidence intervals for the two headlines overlap.     # therefore, the variation (headline 2) isn't more effective     # than the control headline      Thanks again.


title: Whether to include $x$ and $x^2$ in regression model examining diminishing returns when only $x^2$ is significant?

text: I have a data set with sales as $Y$ and total retail shelf space ($x$). I want to investigate the diminishing return, i.e. whether additional shelf space would contribute to more sales. I did regression or Arima model where I found that the squared shelf space ($x^2$) was significant. Thus there is a diminishing return. My question is how to proceed from here? Should I transform my data to $$\log Y = \log x\times \text{shelf_space} + 2\log x \times \text{shelf_space}^2$$ or should I simply ignore the $\text{shelf_space}$ (the one which was not significant i.e. the non squared one) and only do my analysis with shelf space squared? How to do regression which takes diminishing return into account?


title: How to draw a scree plot in python?

text: I am using singular vector decomposition on a matrix and obtaining the U, S and Vt matrices. At this point, I am trying to choose a threshold for the number of dimensions to retain. I was suggested to look at a scree plot but am wondering how to go about plotting it in numpy. Currently, I am doing the following using numpy and scipy libraries in python:               U, S, Vt = svd(A)      Any suggestions?


title: Change of measures with Wiener process

text: I am trying to test importance sampling for a simple a Wiener process $W_t$ in R:               set.seed(123)     Z <- matrix(rnorm(12*1000),12,1000)     W <- apply(Z,2,cumsum) #Wiener process simulated 1000 times for 12 periods     B <- W-(1:12) #Brownian motion with drift -1     w <- exp(-W[12,]) #Radon-Nikodym derivative (actually w <- exp(-W[12,]-12/2))     mean(B[12,]) #-11.973     sum(W[12,]*(w/sum(w))) #-7.881 to compare with -12!      I think that the estimate with the derivative is not really "precise"! Any suggestions to improve the accuracy, besides increasing the number of simulations?


title: Updating/ adjusting theoretical growth curves when raw data becomes available

text: my stats question is as follows: The research group I work for have developed a theoretical growth model for a particular species of fish. The idea is that if you provide some initial starting values for the model you then generate an expected growth curve along with 95% confidence bands. To extend the model we would like to be able to update the model and recalculate the curve when real data becomes available. For example, imagine that at age 6 the model predicts an average weight of 34g, but in a random sample of fish we find that the mean age is, say, 24g, we would like to be able to use this new data to 'tweak' our previously estimated curve. At the moment I am not sure quite how to address this problem so any suggestions would be much appreciated. Thanks


title: How do I get coefficients of a random forest model?

text: I am using `randomForest` to generate a model, and at the end I don't know how I can get the final coefficients that the model is fitting. I know that for linear regression, you just type `summary(lm)` where `lm` is your model, and you get the coefficients.


title: k-means + linear regression: How to split the data for validation

text: I want to cluster my data first using k-means and then determine a regression model for each cluster. Then I want to evaluate the performance of this approach using split validation. I can think of two alternative approaches for that: a) clustering the complete data set and then splitting each cluster into a testing and training data set for training and evaluating the linear models b) splitting the complete data set first into a training and testing data set, then clustering and learning the linear models on the training data, while testing the results with the test data So, in other words, when is it more appropriate to split the data: before or after the clustering step? Which approach is more reasonable?


title: SMOTE implementation in R

text: I am working on an imbalanced classification problem (about 90 thousand records total, of which 2% are target = Y and remaining are target = N). My dataset has been prepared in SAS Enterprise Guide and has a mixture of binary flags (coded as 0 or 1) and interval/numeric variables. I want to implement SMOTE function from DMwR package in R to re-balance this data. I have few questions for which I would seek suggestions/replies: 1) How do I preserve the data type of each variable (be either binary or numeric) to come up with a dataframe which can be given as input to SMOTE? There are around 100s of variables in my dataset? 2) Also, please confirm if SMOTE converts all the variables to factors and does the artificial instance generation procedure? Would really appreciate a quick response from the community folks.


title: Grouping predictor factor levels based on response variable

text: I've read that it's bad to do this, but am looking for details as to why. Suppose we're trying to fit the linear model $Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \epsilon_i$ where $Y$ is continuous and $X$s are categorical with many levels. We don't have any opinion on how the levels should be grouped. We then calculate the average response for each level of the predictors, then collapse levels with similar averages. How badly does this screw up parameter estimates and inferences on them, and why?


title: Derive squared exponential covariance function

text: In Gaussian Processes, SVMs, kernels are used (as to my understanding) as similarity measure. However, they have the constraint that any kernel has to be represented as a dot product. i.e. $k(x_1,x_2)=\langle f(x_1), f(x_2)\rangle$. Note that $f(x)$ could map x to a higher dimension e.g. $f(x)=[x^2\,\,\,\, \sqrt{2}x]^T$ My question is how can you derive the square exponential function $k_{SE}=\exp(-(x_1-x_2)^2)$ as a dot product.


title: mlogit query with multiple options: response rate not good for third variable

text: I am stuck with one query on mlogit. I have three options (0,1,2) for prediction of probabilities. The model equation is:               Model1 = mlogit(A ~ 1 | B +  C + D, data = XYZ, reflevel = "0")     Predicted = predict(Model1, newdata = XYZ, probability = TRUE)      Or               Model2 = mlogit(A ~ 1 | B +  C + D, data = XYZ, reflevel = "1")     Predicted = predict(Model2, newdata = XYZ, probability = TRUE)      Or               Model3 = mlogit(A ~ 1 | B +  C + D, data = XYZ, reflevel = "2")     Predicted = predict(Model3, newdata = XYZ, probability = TRUE)      The above models are running fine. But when I see the probabilities for each response variable (0,1,2)…The confusion matrix doesn’t hold well for response variable 2. The frequencies in the data are: $$\begin{array}{ccc}0&1&2\\\hline61\%&28\%&10\%\end{array}$$ This model is biased for 0 and 1, but not predicting with 2 as an option. Response rate for 2 is coming out to be less than 1%, but for 0 and 1, it's good $(>90\%)$. Can you help me with this?


title: How can I get a "whole treatment effect" for a 3 level categorical variable in a GLMM (glmmPQL in R)

text: I work with habitat use data (summarized binomial data: visited / did not visit) which I fitted with a mixed model of the binomial family with a random factor accounting for repeated sampling of the same individual. I simply want basic descriptions of the data - effects of sex and year on the proportion of time spent in each habitat. The code is as follow (for a basic model without interaction):               Colony.glmmSexYear <- glmmPQL(cbind(ColoYes.allnoF24, ColoNo.allnoF24) ~                                    year.coded + sex, random= ~1| BirdID,                                    family=binomial(), data=mydata)     summary(Colony.glmmSexYear)      Here is the fixed effects section of my output:               Fixed effects: cbind(ColoYes.allnoF24, ColoNo.allnoF24) ~ year.coded + sex                               Value Std.Error DF   t-value p-value     (Intercept)      1.4442270 0.1824173 74  7.917161  0.0000     year.codedY2011 -0.1733713 0.1864744 74 -0.929732  0.3555     year.codedY2012 -0.3004284 0.2027411 74 -1.481833  0.1426     sexM             0.4403108 0.1507471 74  2.920857  0.0046      You were supposed to see that 2011-2012 is also not significant (not shown). In traditional linear models I would report the p value for the whole treatment "year" from an ANOVA associated with a F test (i.e, just one p value), not the individual p values of the various year comparisons (3 p-values). But my understanding is that producing an ANOVA from `glmmPQL` is not possible (?). Reporting those three individual p values of my multiple year comparisons really seems clumsy and out of place... What is the proper way to do this? Should I not be using `glmmPQL`? Specification: I cut into the background description of the study for simplicity, but this is not a resource selection study - please do not suggest using resource selection models.


title: Can I pull a confidence interval out of a single sample by dividing it into sub-samples?

text: Let's say I have taken a sample $S$ of a population. I am trying to figure out the population mean. Because I have only made one sample, the best I can do is assume that $\bar S$ is the mean of the population. I can't get a confidence interval for the population mean, because I can't calculate the standard deviation of a single sample mean. But I could break up my initial set $S$ into disjoint subsets $S_i$ such that $S=S_1 \cup S_2 \cup S_3 \ \cup ... \cup \ S_i$. Then, I pretend that each $S_i$ was a separate sample of the population. I estimate the population mean by averaging $\mu_i = \bar S_i$. I calculate the error of the mean by calculating the standard deviation of $\sigma_i = std(\mu_i)$. I now have an estimate of the sample mean, and a confidence interval for it too, even though in reality I only sampled the population once. This isn't right, is it? Why/why not? You can assume that the population is very large, and $S$ is very small, if that helps.


title: Difference between effect size (partial $R^2$) and coefficients

text: I am working with spoken language data and use linear models do determine the relationship between different phonological processes in my data. **Background** Measures of the regularity of syllable durations are influenced by how _syllable_ is defined. I used two methods to determine syllable boundaries. Method 1 and 2 give different results for regularity of syllable durations, but the difference between the two methods is not the same for every speaker. **Linear model** I think the frequency of certain phonological processes (vowel elision, insertion of glottal stops and others) might be responisble for these differences. I want to find out if that's true and if yes how much each of these phonological processes contributes to explaining the differences. **Measuring effect size with partial $R^2$/sum of squares** To measure effect size, I could use partial $R^2$. The data below is the ANOVA of my model, and from the variance explained by each factor (`Sum Sq`) I can derive how much of the total variance each factor explains.               Analysis of Variance Table          Response: varcos                               Df Sum Sq Mean Sq F value    Pr(>F)         rel_gs                     1 395.48  395.48 46.9135 7.022e-07 ***     rel_all2_gs                1   9.47    9.47  1.1236    0.3006         rel_gs:rel_all_gs          1   1.37    1.37  0.1622    0.6910         rel_gs:rel_all2_gs         1  13.63   13.63  1.6168    0.2168         rel_all_gs:relDur.mean.y   1  15.94   15.94  1.8910    0.1829         rel_all2_gs:relDur.mean.y  1  69.19   69.19  8.2079    0.0090 **      Residuals                 22 185.46    8.43                            **Coeffcients** Beta coefficients can also be interpreted as a measure of effect size, as an answer to this question points out. I was expecting that factors that explain a large part of the variance would also have large coeffcients, and the other way around. But that's not the case. `rel_gs:rel_all2_gs` has a huge negative coefficient, but explains only a small part of total variance.               Call:     lm(formula = varcos ~ rel_gs + rel_all2_gs + rel_all_gs:rel_gs +          rel_all2_gs:rel_gs + relDur.mean.y:rel_all_gs + relDur.mean.y:rel_all2_gs,          data = read.wwb3.diff_syll.perc)          Residuals:         Min      1Q  Median      3Q     Max      -5.7929 -1.5933 -0.5208  1.5569  5.4834           Coefficients:                                Estimate Std. Error t value Pr(>|t|)         (Intercept)                  15.539      2.725   5.702  9.8e-06 ***     rel_gs                     -336.020    118.306  -2.840  0.00952 **      rel_all2_gs                 492.997    183.440   2.688  0.01345 *       rel_gs:rel_all_gs          2533.220    810.012   3.127  0.00490 **      rel_gs:rel_all2_gs        -2292.357    716.625  -3.199  0.00414 **      rel_all_gs:relDur.mean.y   -955.239    315.827  -3.025  0.00623 **      rel_all2_gs:relDur.mean.y   691.578    241.394   2.865  0.00900 **       **Question** Should I only trust partial $R^2$ for effect size and ignore differences between the coefficients? P.S.: I asked a related question some time ago, and it hasn't received much attention and left some issues open. I suppose that might be because the original question wasn't clear enough or was based on toy data. In this question I'm using actual data. If my question could be imporoved in some way I'd be happy to consider any suggestions.


title: Zipf's vs Self-similar: are they really the same

text: Recently I ran into a test using both zipf's and self-similar generated datasets. I followed the description from Jim Gray's paper on generating such datasets (Quickly Generating Billon-Record Synthetic Databases). In that paper it mentions: > It is commonly thought that self-similarity and the Zipf distribution are > the same, or at least close. This misconception ... In particular, other > statistics can yield very different results for the self-similar and Zipf > distributions. So my question is, from a statistic point of view, what is the actual difference between these two distributions intuitively? I am not a statistician, but (probably wrong) my thought is that they are all some power- law shaped distribution... Just clarify: when I am talking about "the same", I am meaning that there is a (mathematic or statistic) way to bridge these two distributions... or one can be converted (mapped) into another...


title: Comparing means - 2 repeated measures variables, controlling for covariate with heterogeneous slope

text: I have a 2 level repeated measures DV of accuracy, and a covariate of response bias. As response bias increases, accuracy level 1 increases while level 2 decreases. I want to see if there is a difference in the means of the groups after controlling for response bias. I can't do it with an ANCOVA. Can I just manually calculate expected values based on a regression equation, and run an ANOVA? I'd like to do it in R or SPSS, so specifics for either would be welcome, but not necessary.


title: Choice of neural net hidden activation function

text: I have read elsewhere that one's choice of hidden layer activation function in a NN should be based on one's _need_ , i.e. if you need values in the range -1 to 1 use tanh and use sigmoid for the range 0 to 1. My question is how does one know what one's _need_ is? Is it based on the range of the input layer, e.g. use the function that can encompass the input layer's full range of values, or somehow reflects the input layer's distribution (Gaussian function)? Or is the _need_ problem/domain specific and one's experience/judgement is required to make this choice? Or is it simply "use that which gives the best cross-validated minimum training error?"


title: Success of a logistic regression model

text: Say I have a model $y=f_n(x_1,x_2,x_3)$. Here say $y$ is categorical and binomial response. i.e. $y$ can be only 0 or 1. Data shows 87% 1 and 13% 0 values. I fit a multinomial logit on a test dataset to try to predict $y$. Validation against a validation dataset showed 70% success. (Success meaning Model prediction matches $y$ in the validation set). Normally I'd be happy with the 70% success rate. But I'm a bit confused because: Say, I had a blackbox model that always answered "$y=1$" no matter what $x_1$, $x_2$ and $x_3$ are. Wouldn't this model anyways achieve 87% success? What's the point behind my multinomial logit unless I can beat 87%? What gives? Am I making a blunder?


title: How to describe the differences in skewed data with same median but statistically different distribution?

text: I am comparing length of stay after laparoscopic and open appendectomy in over 160000 patients. LOS is typically a skewed variable so I use the median and interquartile range and ranksum test to describe the result, which is 2, IQR 2-3 and 2, IQR 2-4, p<0.0001. So there is a highly significant difference but it does not show that well in the median and IQR. The difference in length of stay is very small. This finding is a negative finding but still interesting as it goes against what most people think about LOS after laparoscopic appendectomy. That is why I want to report on the difference in LOS which is not evident if I only report median and IQR with p-value. The arithmetic mean (3.1 vs 2.8 and a difference of 0.3 days) is not recommended as the distribution is skewed. I wonder if any other measure could be used like, geometric mean (2.5 vs 2.3, difference 0.2 days) or harmonic mean (2.1 vs 1.9 days, difference 0.2 days)? When can I use geometric mean or harmonic mean? I am inclined to use the harmonic mean which comes closest to the median and gives a difference of 0.2 days. Is this incorrect?


title: How to deal with an unavoidable correlation between two independent variables?

text: In a simple linear regression model with two independent variables, if there is a strong correlation found between the variables, it is suggested that we should include only one of them in the model. I am building a model where the dependent variable is _Lunch Cost_ and the independent variables are _Number of Class 1 Students Buying Lunch_ and _Number of Class 2 Students Buying Lunch_. However, I found there was a strong correlation between the two independent variables; when _Number of Class 1 Students Buying Lunch_ increases, _Number of Class 2 Students Buying Lunch_ increases as well. In this case, removing one of them from the model seems not ideal since my goal is to know how the two independent variables would describe _Lunch Cost_. What should I do here? Thanks!


title: method of p value correction for a set of dependent p values

text: Here is my situation, I just got a list of 100 p-values, some of which are associated. For details, these 100 pvalues could be seperated into 7 groups, each group was calculated from a expression dataset. However, these datasets had 10%~60% overlap among each other, that is the pvalues inside a dataset were independent, but were not across the dataset, which bring the whole 100 p-values not independent. So is there any solutions you know to directly and simply solve this? Some ref. pointed out to use the BH adjustment inside the R package, while choosing a loose criteria, e.g. 0.2, to define the significance level. Any help will be appreaciated. Thanks! tmp$V7 = p.adjust(tmp$V6,method="BY") tmpV6[1:5] [1] 0.040461040 0.001250268 0.037407030 0.009361665 0.006866937 tmpV7[1:5] [1] 0.31423224 0.08659423 0.29570166 0.14456543 0.13611251 tmp$V7 = p.adjust(tmp$V6,method="BH") tmp$V6[1:5] [1] 0.040461040 0.001250268 0.037407030 0.009361665 0.006866937 tmp$V7[1:5] [1] 0.05252837 0.01447545 0.04943072 0.02416616 0.02275313 plz take a look at this, tmp$V6 is the original pvalue, while tmp$V7 is after correction. However, it seems that after BY correction, I get even worse result than BH(or fdr) method


title: Transforming variables

text: Let $W = X_1/X_2$ where $X_1$ and $X_2$ are normally distributed with means $\mu_1,\mu_2$, variances $\sigma_{1}^{2}, \sigma_{2}^{2}$ and correlation coefficient $\rho$. The exact distribution of $W$ and the standard approximation assuming $X_2 > 0$ is wanted. Why is studying the distribution of $W$ the same as studying the distribution of $$Z = \frac{a+Y_1}{b+Y_2}$$ where $Y_1$ and $Y_2$ are $N(0,1)$? Hinkley mentions that $W$ can be expressed as $$\frac{\sigma_1}{\sigma_2} \sqrt{1-\rho^{2}} \left[\frac{\frac{\mu_{1}}{\sigma_{1}}+Y_1}{\frac{\mu_{2}}{\sigma_{2}}+Y_2} + \frac{\rho}{\sqrt{1-\rho^{2}}} \right]$$ Simplifying I get $$W = \frac{\sigma_1}{\sigma_2} \left[\frac{\sqrt{1-\rho^{2}} \left(\frac{\mu_1}{\sigma_1}+Y_1 \right)}{\frac{\mu_2}{\sigma_2}+Y_2}+ \rho \right]$$ $$= \frac{\sigma_1}{\sigma_2} \left[\frac{\sqrt{1-\rho^{2}} \left(\frac{\mu_1}{\sigma_1}+Y_1 \right)+\rho \left(\frac{\mu_2}{\sigma_2}+Y_2 \right)}{\frac{\mu_2}{\sigma_2}+Y_2} \right]$$ How is this equivalent to $W = X_1/X_2$?


title: How should outliers be dealt with in latent growth curve/GMM modeling?

text: When doing analysis of growth trajectories you often find natural outliers. Often these are not because of any "errors" in the dataset, they are simply there because certain cases grew at an explosive rate compared to the majority of the cases in the dataset. A case in point is open source repositories where most seem to exhibit moderate growth, whereas some (Linux, Twitter Bootstrap, etc.) exhibits very high growth rates. This can lead to difficulties when estimating latent growth curve models and growth mixture models. For example, in the latter case it is often difficult to perform statistical tests (such as the LRT significance test for K and K-1 class models) when outliers are present, since they tend to be placed in a class of their own with a very small class membership (low N) therefore preventing reliable significance testing. Some packages, such as Mplus even refuse to calculate the LRT statistic if a class has fewer than 5% of total N in a single class. **What statistical tests or rules of thumb can be used as a basis for excluding outliers in latent growth curve/GMM analysis?** **Are there any special considerations for growth mixture modeling?**


title: Pre, Post, and Change predictors all in the same model

text: We have a variable, $X$, measured at pre-study and post-study and are studying the effects of $X$ across changes of several outcome variables (post-pre). Currently we are using regression models to investigate the change in $X$ (post-pre) effect while including pre-study $X$ to control for the initial level of $X$. We are also interested in if the relationship between change in $X$ and the outcome variable differs by post-$X$ level. Thus I am considering adding post-$X$ level and the interaction between change in $X$ and post-$X$ level to the model. Therefore the model would include pre-$X$, post-$X$, change-$X$, and change-$X$*post-$X$ as predictors. This seems like a problem to me, especially since change-$X$ is a function of pre-$X$ and post-$X$. However we would still like to control for pre-$X$ level. Is there a better model to address the question of if the change-$X$ effect differs by post-$X$ level?


title: Pattern Recognition Problem

text: I have the following questions: Suppose we have the following chart, which represents the relationship between faculties and types of objects; for example CTE represents Computer Technologies, and the intersection between CTE and Wood says that at the CTE faculty there are 56 types of wood. My question is: Can I use this table as a tool for PCA analysis? Or can I represent these data in matrix form and perform a principal component analysis from this data table? I need it for pattern recognition rate problems. I know Matlab well and linear algebra as well; I just need advice. ![enter image description here](http://i.stack.imgur.com/QonjI.png)


title: In a Monte Carlo approximation of a product of expectations, can the same samples be used for both expectations?

text: I am trying to approximate a product of expectations: $\operatorname{E}[f(x)]\operatorname{E}[g(x)]=\sum_x P(x) f(x) \sum_x P(x) g(x)$ with $N$ Monte Carlo samples $(x_1,x_2,...,x_N)$ from $P(x)$: $\operatorname{E}[f(x)]\operatorname{E}[g(x)] \approx \left( \frac{1}{N} \sum_i f(x_i) \right) \left( \frac{1}{N} \sum_i g(x_i) \right)$ The Monte Carlo samples are very expensive to generate. What are the consequences of using only $N$ samples total, reusing the same samples for each expectation (vs. using $2N$ samples so that each expectation uses a different set of $N$ samples)? Will that produce an incorrect/biased estimate, higher variance, etc.? ### More Context (added Aug 2, 2013) Specifically, the application that inspired this question concerns gradient descent optimization of entropy in Boltzmann machines. Assume $M$ visible/observable binary neurons (Bernoulli variables), and no hidden neurons for simplicity. The state of each neuron is $X_i \in \{0,1\}$, and the joint state of all neurons is $X \in \{0,1\}^M$. These neurons are fully connected with a set of ${M \choose 2}$ weights $\theta_{ij}$ between neurons $i$ and $j$ with $i \neq j$ (a symmetric $M \times M$ parameter matrix with $0$s on the diagonal). Assume the following relationships: $$ p(\left. X \right| \theta) = \frac{1}{Z_{\theta}} \exp \left(-\frac{E(X,\theta)}{T} \right)$$ $$ Z_{\theta} = \sum_X \exp \left( -\frac{E(X,\theta)}{T} \right) $$ $$ E(X,\theta) = -\sum_{ij} \theta_{ij} X_i X_j$$ where $Z_{\theta}$ is the partition function, $E(X,\theta)$ is the energy function, and $T$ is the temperature (a global weight scalar which is usually, but not always, assumed to be $1$ for Boltzmann machines; note that Boltzmann's constant $k_B$ is assumed to be $1$ here and is not shown). For this model, say we want to maximize the entropy of $X$ by adjusting the parameters $\theta_{ij}$. We can compute the gradient of the entropy expression and apply a gradient descent procedure to maximize it. The entropy expression is: $$ S_{\theta} = -\sum_X p(\left. X \right| \theta) \log p(\left. X \right| \theta)$$ The partial derivative of the entropy wrt each weight $\theta_{ij}$ is (derivation not shown): $$ \frac{\partial}{\partial \theta_{ij}} S_{\theta} = \frac{1}{T^2} \left[ \sum_X p(\left. X \right| \theta) \, E(X,\theta) X_i X_j - \sum_X p(\left. X \right| \theta) \, E(X,\theta) \sum_X p(\left. X \right| \theta) \, X_i X_j \right] $$ These expectations cannot generally be computed exactly if $M$ is large because they require summing over $2^M$ possible states. Instead, we just approximate them with $N$ Monte Carlo samples from $p(\left. X \right| \theta)$: $$ \frac{\partial}{\partial \theta_{ij}} S_{\theta} \approx \frac{1}{T^2} \left[ \left( \frac{1}{N} \sum_k^N E(X^k,\theta) X_i^k X_j^k \right) - \left( \frac{1}{N} \sum_k^N E(X^k,\theta) \right) \left( \frac{1}{N} \sum_k^N X_i^k X_j^k \right) \right] $$ where $X^k$ and $X_i^k$ refer to the $k$th Monte Carlo sample. Note that optimizing entropy alone would simply make all $\theta_{ij}=0$ because that maximizes $S_{\theta}$; however, in practice we may want to combine this objective with something else, like a maximum likelihood objective, and optimize a weighted combination of the two objectives. In that case, the max likelihood objective would make the model fit the data, and the max entropy objective would regularize the model to avoid overfitting. Regarding the original question posted here, I was curious if this gradient approximation could use only $N$ Monte Carlo samples total (i.e. can the product of the two expectations use the same $N$ samples?), or whether it requires $2N$ samples (one set of $N$ samples for each expectation in the product) to remain an unbiased estimator of the gradient. Note that, in any case, the expectation on the left side of the gradient expression can reuse the $N$ samples from either one of the two expectations on the right side without bias because we can simply rearrange things to put a summation outside the outer brackets. Specifically, we can do either this: $$ \frac{\partial}{\partial \theta_{ij}} S_{\theta} \approx \frac{1}{T^2} \frac{1}{N} \sum_k^N E(X^k,\theta) \left( X_i^k X_j^k - \frac{1}{N} \sum_l^N X_i^l X_j^l \right) $$ or this: $$ \frac{\partial}{\partial \theta_{ij}} S_{\theta} \approx \frac{1}{T^2} \frac{1}{N} \sum_k^N X_i^k X_j^k \left( E(X^k,\theta) - \frac{1}{N} \sum_l^N E(X^l,\theta) \right) $$ ...reducing the problem from 3 to 2 expectations. Again, the question is (essentially) whether it is unwise to try to reuse the same set of $N$ Monte Carlo samples for these two expectations.


title: colinearity between variables

text: I am running a multiple regression of Y~a+b+c+d etc... I want to do a quick check to see whether my different explanatory variables are colinear (they're a mix of categorical and continuous). There seems to be a whelm of complicated statistics behind all this -- would looking at the R-squared value of a simple regression between each variable in turn (a~b, a~c, a~d etc) be a satisfactory coarse estimate of colinearity? Or would I have to use VIF statistics and other more complicated methods? I hope that this is not off-topic... :-) Edit: Also, if a simple look at the R-squared value is sufficient, what is an acceptable amount of correlation? I have two variables that are 54% correlated for example - this seems high to me...


title: How to associate the probability output of a Random Forest classifier with reality

text: I am currently working with the national forest inventory of Mexico. A sample sceheme that collected data concerning the state of forests on points spread out all over the country. One of these variables is a cualitative one concerning ground damage/degradation due to many causes. E.g. pasture activities. This variable also has associated an ordinal degree of "damage" if there is any, 1: 1%-20%, 2: 20%-40% ... (I dont recall the size of the intervals exactly). As a quick excercise I took this pasture variable and ignored the degree of damage. I just made a binary variable that indicates damage is present / not present. Then I trained a Random Forest classifier using many remote sensing, topographic and climatological variables. These variables are available wall to wall in the country so I can later use the model to generate a map. If I train the classifier using a balanced sample between 1's and 0's and do a hard-classification I get a precision with CV of around 75%. I'm more interested in generating a probability of degradation map so I didnt balance the sample and generated a soft-classification. This would be the probability of class 1 (degradation is present): ![enter image description here](http://i.stack.imgur.com/pXzqq.jpg) Does this excercise make sense to anyone? My true concern is the interpretation. How would one associate these probabilities with what is really going on in the country? Or does anyone have an example of when one would want to associate a soft-classification with a real life phenomenon? To me it makes much more sense to think fuzzily. But I am lost when it comes to describing what is going on.


title: Eigenvectors of a shuffled correlation matrix

text: A question about a shuffled vs. unshuffled correlation matrix I took the correlation matrix and shuffled between its values symmetrically. (shuffled only the left lower triangle of the matrix and copied the values to the upper right triangle. saving the diagonal of ones untouched) The eigenvalues spectrum of both the correlation matrix and the shuffled matrix have one large eigenvalue. As I understand from my results, this eigenvalue reflect a uniform collective mode of the fluctuations (because most correlation values are positive and are the same for the shuffled matrix) However, the eigenvector of these first eigenvalues is different: For the correlation matrix, most are positive and very distributed. For the shuffled matrix, all of them are positive and are centered around one positive value. (Of course that if by chance the shuffling was the same as the correlation matrix I would have reached similar results... Can somebody please explain to me this phenomenon? EDIT: Let's say that I have a correlation matrix between the neural response of 10 neurons -> a 10X10 correlation matrix with a diagonal of ones. I extract all elements from the matrix (I have 45 elements which appear twice Cij=Cji). Then I randomly insert these elements into a new matrix while saving the diagonal of ones and the symatry of the matrix: Cij=Cji. Is that clear now?


title: Text mining "power calculations"?

text: I have a list of 6,500 or so medical treatments, of which I have classified 700 or so as involving a physician or not. I am interested in both the specific question of how to calculate whether 700 observations is sufficient, and the general question about how one does power calculations in text mining. A typical treatment might have 3-10 words describing the intervention (e.g. "Coronary artery bypass graft surgery for left main artery disease"). My training dataset adds to that a binary variable indicating whether a physician is involved in the process or not. I want to predict this binary variable on all the remaining observations. I can classify more if necessary, but obviously the time and money cost increases to do so. I am curious how I might know when to stop creating a training dataset. I assume that a "power calculation" would involve:   * Characteristics of the text (number of words per observation, prevalence of sparse words, overlap of the words between each of the binary categories, etc.)   * Desired prediction accuracy   * Algorithm used Is there a formal method for this? I'm using RTextTools if that helps.


title: What does the following message mean in SAS?

text: I received the following message in SAS: > Measures of association between the observed and predicted values were not > calculated because the mean scored predicted probabilities are > indistinguishable when they are classified into intervals of length 0.008.


title: Kullback-Leibler divergence

text: Suppose we seek to approximate an arbitrary distribution $p_1(x)$ by a normal $p_2(x) \sim \mathcal N(\mu, \Sigma)$. How can I show that the values that lead to the smallest Kullback-Leibler divergence are: $$ \mu_1 = \mathbb E_1[X] $$ and $$ \Sigma_1 = \mathbb E_1[(X − \mu)(X − \mu)^T], $$ where the notation $\mathbb E_1(\cdot)$ indicates the expectation is taken over the density $p_1(x)$? For reference, the definition of Kullback-Leibler divergence is $$ D( p_1 \|\; p_2) = \int p_1 \log(p_1/p_2) \text{d}\lambda \> . $$


title: Can I ignore coefficients for non-significant levels of factors in a linear model?

text: After seeking clarification about linear model coefficients over here I have a follow up question concerning non-signficant (high p value) for coefficients of factor levels. Example: If my linear model includes a factor with 10 levels, and only 3 of those levels have significant p values associated with them, when using the model to predict Y can I choose to not include the coefficient term if the subject falls in one of the non-signficant level? More drastically, would it be wrong to lump the 7 non-significant levels into one level and re-analyze?


title: What is a practically good data analysis process?

text: Maybe this is a dumb question as my experience has been more on statistical computing algorithms than on data analysis. But I would like to know, or have references on, analysis process most of statistical data analysts go through for each data analysis project. If I make a "list", to complete data analysis project, an analyst has to   1. first collect requirements for the project,    2. plan/design his data analysis based on those requirements before    3. actually pre-processing data,    4. executing the data analysis and    5. writing a report based on his analysis results. For this question, I am interested in more details of Step 2. But I understand this is not practically clear cut as the analyst might have to change his plan or design according to data analysis output. Is there any reference on this subject?


title: Quick and simple cluster analyses for univariate data

text: Can you suggest any quick and simple clustering analyses, for univariate real- valued data? In other words, I have $n$ real numbers, $x_1,\dots,x_n$ where $x_i \in \mathbb{R}^+$, and I want to cluster them. I don't know a priori the best number of clusters, so that's something the method would need to discover as well. It'd be nice if it were simple to code up in Python. Something quick and dirty -- say, easy to understand, easy to implement, and pretty effective-- beats something complex but optimal. * * * _My motivation:_ As mentioned elsewhere, in the application in front of me now, a reasonable model would be to say that the points were generated from a mixture of Gaussians. I don't know the parameters of the mixture model, but if it helps, I can reasonably assume some lower bound on the probability of each component: for instance, if you like concrete numbers, you could imagine I have $n=40,000$ samples and each component of the mixture model is guaranteed to have proability at least $0.0001$. A twist is that there may be a few outliers thrown in as well, and I want to detect the outliers. @whuber suggested that a good approach to outlier detection would be to start by clustering the points, so I'm looking for quick-and-dirty clustering methods. That's my motivation at the moment -- but I expect the broader question is of general, independent interest, so feel free to ignore this specific motivation.


title: Is there a difference between the "maximum probability" and the "mode" of a parameter?

text: I am reading a manuscript that provides the "maximum posterior probability" in a Bayesian context as a statistical summary of a parameter. Is the term "maximum posterior probability" equivalent to (i.e., "exchangeable" with) "posterior mode"?


title: Fit a logistic regression code in R

text: If I have 10 Variables (Q,W,E,R,T,Y,U,I,P,A) and I want Q to be my response variable and other 9 to be my predictors variable. Do I write it in R like this `EXAMPLE<-glm(Q~W+E+R+T+Y+U+I+P+A,family=binomial)` Furthermore, what if Q is Binary (goes from 1 to 0) and all the other 9 variables are categorically numbered. is it still the same or do i need to write it diff in R


title: Normalization of data before using a Kolmogorov-Smirnov test

text: I am working with two sets of data and would like to determine if they come from the same distribution -- hence the two population KS-test. At this point, I have binned both sets of data, counting the number of entries in each bin, to create relative frequency histograms. These binned arrays look like:               dis1 = [0.0, 4.0, 9.0, 10.0, 23.0, 26.0, 43.0, 44.0, 38.0, 4.0]     dis2 = [116.0, 160.0, 209.0, 247.0, 415.0, 455.0, 382.0, 288.0, 161.0, 44.0].      When I perform the KS-test for 2 datasets, I get:               [h,p]=kstest2(dis1,dis2)     h = 1     p = 1.7012e-04      meaning that the null hypothesis was rejected and the samples do not come from the same distribution. Now, if I normalize the binned data with respect to the maximum in each array, I have:               norm_dis1=[0.0, 0.090909090909090912, 0.20454545454545456, 0.22727272727272727, 0.52272727272727271, 0.59090909090909094, 0.97727272727272729, 1.0, 0.86363636363636365, 0.090909090909090912]     norm_dis2 = [0.25494505494505493, 0.35164835164835168, 0.45934065934065932, 0.54285714285714282, 0.91208791208791207, 1.0, 0.83956043956043958, 0.63296703296703294, 0.35384615384615387, 0.096703296703296707],      such that now the norm_dis arrays vary from [0,1]. The 2 population KS-test for these arrays gives:               [h,p]=kstest2(norm_dis1,norm_dis2)     h = 0     p = 0.3129,      such that the null hypothesis is accepted and the populations come from the same distribution. These two results are seemingly discrepant. What is the statistical difference between normalizing and not normalizing? And which is the correct statistical analysis of these two populations? Thanks for your help.