title: Taking expectations for Fisher information/Jeffreys' prior

text: I've been working through a few examples of the Jeffreys' prior; however, I'm a little confused by one section. I was hoping that somebody could provide some clarification. If the Jeffreys' prior is defined as: $$ \pi_{j}(\theta) = \sqrt{I(\theta)} $$ where $I(\theta) = {\bf E}\left[\left(\frac{\partial\text{log}p(x_{i}|\theta)}{\partial\theta}\right)^2\right] = -{\bf E}\left[\frac{\partial\text{log}p(x_{i}|\theta)^{2}}{\partial\theta}\right]$ Then, for $p(x_{i}|\theta) \sim Bin(k,\theta)$ with ${\bf E}[X] = k\theta$, Step 1: $$ {\bf E}\left[\frac{\partial\text{log}p(x_{i}|\theta)^{2}}{\partial\theta}\right] = - \frac{k}{\theta} - \frac{k}{(1-\theta)^{2}} +\frac{k\theta}{(1-\theta)^{2}} = \frac{k(\theta-1)}{\theta(1-\theta)^{2}} $$ Now, for the example that I've been looking at, it states that: Step 2: $$ -{\bf E}\left[\frac{\partial\text{log}p(x_{i}|\theta)^{2}}{\partial\theta}\right] = \frac{k}{\theta(1-\theta)} $$ How is this step achieved? Obviously can't we just put a minus sign on the result from Step 1 to yield: Step 1: $$ {\bf E}\left[\frac{\partial\text{log}p(x_{i}|\theta)^{2}}{\partial\theta}\right] = - \frac{k}{\theta} - \frac{k}{(1-\theta)^{2}} +\frac{k\theta}{(1-\theta)^{2}} = \frac{k(\theta-1)}{\theta(1-\theta)^{2}} $$ Step 2: $$ -{\bf E}\left[\frac{\partial\text{log}p(x_{i}|\theta)^{2}}{\partial\theta}\right] = - \left( - \frac{k}{\theta} - \frac{k}{(1-\theta)^{2}} +\frac{k\theta}{(1-\theta)^{2}} \right) = -\frac{k(\theta-1)}{\theta(1-\theta)^{2}} $$


title: Trying to compare a fully fixed and and fully random model in R

text: So I have a data set which I am modeling using the following: `model1 <\- lm(subject ~ .^2, data=hre.train)` This model treats all of the variables as fixed. I want to compare this model to a model which treats all variables as random. Would I use something like: `model2<-lmer(subject~(1|town)+(1|district)+(1|street)+(1|family)+(1|gender)+(1|replicate),data=hre.train)` I can construct a fully nested model by using `model3<-lme(subject~1,random=~1|town/district/street/family/gender/replicate,data=hre.train)` I just keep getting an error when I use `anova(model1,model2,model3)` The goal is to build an optimal mixed model for prediction but I need something to compare to when I start combining fixed and random effects (i.e gender as fixed, others as random etc)


title: Statistical significance between two percentage change

text: I have two groups (test and control) and I wanted to determine the impact of a certain initiative on the test group. The comparison would be 4-month before the initiative and 4-month after. The starting baselines for the groups are slightly different so to calculate the impact use percentages. If the test group improves by 5% after the initiative while the control group changes by 1%. How do I determine whether the 5% change is statistically significance compared to the 1% change in the control group?


title: One-class Classification of multidimensional vectors

text: I have a m x k User-feature matrix (m >> k) obtained by factorizing an original User-websites matrix (m x n) that has #page views as entries. Additionally, there are users (say r) who have been labeled as being part of a particular group. I'm now looking to find a way to see who among the remaining (m-r) users is 'similar' to these pre-classified r users. This seems like a possible scenario for a one-class SVM, but I'd like to know if there are any other methods that I should consider. On a related note, I'm not sure if there is a way I can fit a multivariate distribution to the pre-classified r users, and then for each of the (m-r) users, get a probability of being part of that distribution. One could then use some kind of a cutoff to decide the 'similar' users. I'm using R for my analysis, so any pointers on specific packages to consider would be helpful as well. Thanks in advance!


title: Time series analysis for 2 series where one is dependent on another

text: If I have two time series $A$ and $B$. $A$ is dependent on $B$. I want to forecast future values of $B$. What statistical techniques should I learn and try? As an example consider a time series of sales of Apple ipods every month. Another time series is the number of defects in ipods reported. Clearly if the sales are more the defects will be more too, assuming a constant probability for a defect to occur. How to use both $A$ and $B$ to forecast $B$?


title: Is this statement about the conditional expectation of a sum true?

text: For expectations of random variables (RVs) $X$ and $Y$ it is true that $$E(X+Y)=E(X)+E(Y)$$. My question is whether when conditioning on RV vector $Z_{1...J}$, it is also true that $$E(X+Y|Z_{1...J})=E(X|Z_{1...J})+E(Y|Z_{1...J})$$. Thanks for enlightening answers!


title: Standardizing response variable in shrinkage/regularization

text: I know that I should standardize my predictors before estimating something like `Lasso`. But what about the response variable? Do I standardise this as well? Only concerned with impact on ability for the model to generalise out of sample.


title: Mean residual life with censored observations and covariant eps

text: I have a time variable ( time of claim settlement after occurrence ) that may be censored ( claim still open at the moment of observation) and some additional covariates ( line of business, claim size, ...) I am interested to estimate the mean residual life for claims still open. My question is twofold : 1. If I do not want to use covariates exists any approach to estimate the mean residual time unbiased way ( and an R package (? 2. Is there an approach and corresponding R package to perform such estimation allowing for covariates effect ?


title: Should I exclude predictor variables if used to create a new one?

text: I have a dataset that includes race, gender, income, and family size. In addition, a variable for "sliding fee scale" tier is included, which is determined by income and family size. Should income and family size be excluded when fitting data to a logit model? Or should I first include all variables, examine correlations and then conduct a feature selection procedure to determine whether to include/exclude income and family size? X1: Gender X2: Family Size X3: Income X4: Tier (determined by X2 and X3)


title: Problem constructing model matrix to generate predicted values to plot Bayesian glmmBUGS output

text: My goal is to plot the predicted values generated from a Bayesian model using glmmBUGS run through R. I believe my problem stems from a lack of understanding in 1) how to properly construct a model matrix (R programming problem) because I am not completely sure which values it needs (Bayesian stats problem). I have saved the model output (posterior probability distribution) in the file named larv.r2jags. To get the predicted values, I first generated my matrix of coefficients like this:               larv.mcmc<-as.mcmc(larv.r2jags$BUGSoutput$sims.matrix)      Then made a dataframe to store my predicted values in like this:               newdata<-expand.grid(LARV=levels(larvbays$Mom), SETT=seq(0, 400, 10))      Then made the model matrix:               mm<-model.matrix(~LARV*SETT, data=newdata)      Finally matrix multiplication: aa<-larv.mcmc %*% t(mm) However, at this point, Rstudio throws up the following error:               Error in larv.mcmc %*% t(mm) : non-conformable arguments      Larv.mcmc is an “mcmc” object of dimensions of 9252 X 38 and mm is “matrix” of dimensions of 287 x 14. I believe the object classes are correct (mcmc x matrix), but the dimensions might be incorrect and that is why I am getting the error of “non-conformable arguments.” See http://www.stats.uwo.ca/faculty/braun/SCUT1/notes/complinalgB.pdf for info: "That is, the number of columns of the ﬁrst matrix must match the number of rows of the second matrix." Mine obviously do not match (38 and 287), but I am not sure how to manipulate either newdata or mm to do so. I believe the root of the problem lies in how I have constructed the “newdata” dataframe. The only examples I had for doing this only used factors that were categorical, so I had some trouble adapting my example with my categorical predictor of “LARV” with my continuous predictor of “SETT”. Possibly, I am not constructing my model matrix correctly (because I have constructed the object newdata incorrectly) and that is why I am getting the error message of non- conformable arguments? Does any one have advice on how to deal with the error? I am really looking to understand this process of creating a matrix that represents predicted values so that I can not only make figures of the Bayesian outputs, but also predict values for of hypothetical predictor values. Thanks in advance! If anyone knows of some good resources with dealing with Bayesian interpretation and graphing, that advice would also be much appreciated! cheers, K Edit: I have now changed my model matrix to have the dimensions 35 x 14, which matches the dimensions needed from the mcmc object: 9252 x 35. However, I still get the same error message: "Error in larv.mcmc %*% t(mm) : non- conformable arguments." It doesn't seem to be the dimensions then...


title: How to form a frequency distribution table from different survey results

text: I have a survey I have conducted on a certain market to determine whether or not we should consider reentering that market. After conducting the survey I was hoping to draw some conclusions from this including a frequency distribution table of results as well as measures of central tendency and dispersion. The survey consists of a number of questions with different options for responses. E.g.   Question 1 - the possible responses are Highly likely, Likely, Unlikely.   Question 2 - possible responses are Yes or No. Can someone please help me as to how I am able to analyze this data correctly to draw on results and conclusions. I was hoping to see the frequency of responses per question if that's the right way to do it! Thanks in advance


title: Please help with basic interpretation of an insignificant interaction. (With significant main and conditional effects)

text: I have a very basic understanding of statistics, and have done many correlations and regressions or simple t-tests, but never an interaction and it is confusing me a bit. I only need to understand the basics of it for now and would really appreciate any input from all you statistics experts. I just basically need to understand what an insignificant interaction means, when the main effects and the conditional effects are all significant. I am using modprobe syntax (Hayes, 2010) for SPSS to test an interaction between 3 variables. I have 3 variables they are:               Outcome Variable      InterVio (Inter Group Violence) (scale from 0-13)          Focal Predictor Variable      Ruminati (Rumination) (Scale from 1-4)          Moderator Variable      ColNarc (Narcissism) (Scale from 1-6)      The hypothesis (simple psychology report hypothesis) is: That narcissism and rumination interact at high levels to predict higher aggression and interact to predict lower aggression at lower levels. The overall model shows significant main effects, but insignificant interaction:               Complete Model Regression Summary            R-sq          F        df1        df2          p          n           .3726    28.7066     3.0000   145.0000      .0000   149.0000          ===================================================================                       b         se          t          p     constant     -.0407      .6950     -.0585      .9534     Ruminati      .6450      .2777     2.3229      .0216  <----     ColNarc       .5646      .2258     2.5000      .0135  <----     interact     -.0690      .0838     -.8237      .4114  <----          Interact is defined as:      Ruminati X        ColNarc          =====================================================================      I quite simply want to report this interaction in a psychology research report (basic one). I just do not understand how to interpret the insignificant interaction. I think I could say that "the IVs predict the DV independently but the interaction was not significant (therefore it is by chance that at different levels of one IV the second IV and the DV increase?) But then the output shows 'the conditional effect of focal iv at values of the moderator iv'.               Conditional Effect of Focal Predictor at Values of the Moderator Variable         ColNarc          b         se          t          p    LLCI(b)    ULCI(b)          2.2472      .4900      .1169     4.1920      .0000      .2590      .7210          3.2077      .4237      .0903     4.6905      .0000      .2452      .6023          4.1682      .3574      .1249     2.8614      .0048      .1105      .6043          Alpha level used for confidence intervals:       .05          Moderator values are the sample mean and plus/minus one SD from mean      My problem is that these are significant and therefore makes me wonder what they mean. Because if the entire interaction is insignificant (occurred by chance, as I understand) how can it still interact significantly at different values? (So would I ignore the interaction significance and just report the conditional effects and the slope plot?) Also does this table mean that the slope line plot graph would be statisticall significant? (But how can it be if the interaction is insignificant? Its so confusing for me, how anything can be significant if the interaction isnt...) Basically it would be great if someone could just tell me what the significance of the last table actually is meant to be telling me. Also I was wondering if marginal effects could mean the same thing as conditional effects? Edited to Add: I have supposedly enough literature to support this interaction existing, so maybe I should run a mediation instead after I find the insignificant interaction? Sorry for rambling everyone, but Im just getting more and more confused.... Thank you very much for anyone who bothers to even read this. I feel useless at stats sometimes.


title: How do statisticians determine which distribution is appropriate for different statistical tests?

text: For example, the test statistic calculated for the ANOVA test is compared to an F-distribution, while a t-test comparing means compares the test statistic to a t-distribution.


title: Numerical properties of the logistic growth model for non-linear regression

text: I am using the nls procedure in R to fit a logistic growth model. In their SSlogis function, José Pinheiro and Douglas Bates chose the formulation                Asym / (1 + exp((xmid-input) / scal))      for their model. As I am fairly inexperienced with the numerical properties of such models, I wonder:   * Can somebody explain why the authors chose this formulation instead of possible alternatives? In particular, ecologists seem to prefer a model with initial population, carrying capacity and growth rate. Does the formulation above have favourable numerical properties?   * It seems that when the model is misspecified and the data are actually fairly linear with time, this formulation often fails to converge. Could such a problem be avoided?   * Is parameter orthogonality a key concern here or are other aspects of the model more important?   * Is it trivial to extend this model to allow a flexible intercept? Would the following model provide sensible numerical properties?               Intercept + (Asym - Intercept ) / (1 + exp((xmid-input) / scal))      I am, of course, open for alternatives as long as it allows for some flexibility in intercept, location where 50% of the growth has been achieved and asymptote.


title: How to summarise the ordering pattern of a set of four stimuli by participants?

text: I did an experiment where I asked 13 participants to order 4 stimuli (let's call them stimulus A, B, C and D), so there are 24 possible ways to order 4 items. Now I would like to know how to analyze this data in the correct way. Indeed I would like to say that on average a certain order configuration has been preferred. For example I would like to say that considering all the configurations expressed by a participant, the following represents better the order preference of the participants   1. stimulus D   2. stimulus B   3. stimulus C   4. stimulus A


title: Customizing a Split-Plot Model in SPSS

text: I am having a little trouble with SPSS in a rather complex design. I have to analyze an experiment with 1 within factor and 3 between factors. I don't want to test all possible effects (just the 4 main effects and the 2-way interactions from the between factors), so no interaction effect with the within factor is included. I am using the `GLM Repeated Measures` procedure in SPSS. In SPSS you can build your model any way you want when you're just dealing with between factors, but if a within factor is involved in the model, SPSS automatically builds all possible interactions with the within factor (2-way, 3-way and so on). A full factorial model with all possible effects would be much too complex and could not be interpreted very well, so I just want to use a custom model, despite the fact that the MS Error will increase. (It's worth it...) Maybe someone of the experts here have faced the same problem and knows a solution (e.g. via Syntax)?


title: BRT analysis using count data

text: I have some problems with my BRT analysis. **Introduction to the data:** The dependent variable is count data of a specific palm species in SA, and the predictors consists of nine various kinds of soil related variables(Mg, Na, Ca, %sand etc.) The data consist of two data sets with 86 and 77 plots containing palm count and soil variables. In the data set with 86 plots, only 30 of them had recordings of the palm species, while the species was absent in the remaining 56 plots. **Issues:** I use R with the `gbm` packag , and the expanding codes handed by Elith et al. (2008) "A working guide to boosted regression trees". My first issue deals with deciding what "family" I should use. I don't know if using a Poisson distribution is the right decision because of the amount of 0s in my data set. What are your thoughts? The second problem addresses the terms needed to find the right parameter values in the BRT analysis. There are a lot of parameters you can alter. Normally if you are handling binomial data, you can run a lot of different parameter combinations and select the best model based on the highest-ranked area under the curve of the receiver operating characteristics (ROC) analysis. I am using count data, and since I have so few data points I feel obliged to some sort of cross-validation method. I am pretty convinced that using `gbm.step` (Elith et al 2008) is a good solution to find the best parameter values. I was thinking about running 100 model replications with different combination of parameter values based on what other authors have used. Parameter values found in the literature `tree.depth`/`interaction.depth` from 1–the amount of variables, bagging fraction 0,4–0,75 and `shrinkage.rate`/`learning.rate` from 0,1–0,01. I have problems finding a method to evaluate my models and therefore finding the optimal number of trees for further analysis. The `gbm.step` function supplies you with a CV score. Here is an example from `gbm.step`: > mean total deviance = 292.371   >  mean residual deviance = 206.759 > > estimated CV deviance = 259.41; SE = 69.068 > > training data correlation = 0.57   >  CV correlation = 0.442; SE = 0.063 Could I in some way sum and compare CV deviance values with other combinations of parameter values and thereby obtain the optimal model? Or is there another way to use CV to evaluate the performances of the different replications of the parameter combinations?


title: Are the following GLM models nested?

text: I am new to GLM modeling but committed to learning as much as possible... I have the following situation. Data can fall into one of 4 buckets. It is a GLM model with a poisson distribution and a log link. The three models can be described as:               1) A  A+B  A+B+B   A+B+B+B      The design matrix for this would be with lets say one observation for each bucket:               1 0     1 1     1 2     1 3      or               2) A  A+B  A+B+C   A+B+C+C          1 0 0     1 1 0     1 1 1     1 1 2      or               3) A  A+B  A+B+C   A+B+C+D   (not realistic. for ex)          1 0 0 0     1 1 0 0     1 1 1 0     1 1 1 1      It may not be obvious, but we are saying we need a new parameter to describe the last two buckets in model 2, and we need two new parameters for 3. Would these be considered nested models? How could I compare them to determine which one is best? Thank you for your help.


title: What is a reasonable sample size for correlation analysis for both overall and sub-group analyses?

text: I am running Pearson's correlation on an overall sample of 400 respondents. When I isolate male and female responses, my sample becomes 220 male responses and 180 female responses. If I further isolate male and female responses by (say) age groups, some sample sizes become as low as 35 responses (for example, for females over the age of 65). **My question: How good are these sample sizes for correlation analysis?** (I am looking at the relationship between income levels and overseas travel.) (I think this has something to do with margin of error but how does this apply to inferential analysis which is based on probability. I can understand its role in descriptive statistics such as results of a political poll).


title: Comparing reproducible research strategies: brew or Sweave vs. R2HTML

text: This is a bit more "think about it" question - but I see it as an important one to ask. I have been struggling for the past few days with having a more reproducible- research-like workflow. I am confused with the two different strategies for writing a report. The two strategies are:   1. Sweave or brew. Where there is a report.Rnw or report.brew file that has a mixture of some markup language (either HTML or LaTeX) and R code between special braces (say <<>>= @). This file needs to be run through Sweave or brew in order to create the report file (report.html or report.tex).   2. R2HTML (for HTML) and Hmisc (for LaTeX). Where the .r file uses R functions to construct report.html or report.tex; running the R commands generates the report directly. What is clear to me is that most people online seem to be using option 1. But I do not understand why it is so common, when option 2 seems to me (without too much experimenting) to be less work. When is each of the two strategies better? Any ideas/feedback/thoughts would be welcome, Thanks.


title: How is the proper number of lags for ACF or PACF displaying?

text: How many lags should be used for ACF or PACF displaying if we have $S$ seasonality? For example, for 500 observations I have 25 lags   for 200 observations I have 22 lags It is independent from frequency of seasonality (for $S = 7, 14, 50, 60, ...$ number of lags on the picture is the same). Is it enough for model estimation based on ACF/PACF for high frequencies? We do not need to show lags for few periods; for example, for $S=50$ at least 170 lags (it contains 3 periods)? Sorry, but available materials show different approaches.


title: Flipping identifiable coins in batches

text: I have collected data to estimate a parameter and am now puzzled about how to generate confidence intervals: Setup: 1) We have a bag with $N$ coins.   2) Each coin $i \in N$ has a **known** probability of $q_i$ for being picked from the bag.   3) Each coin $i \in N$ has an **unknown** probability $p_i$ of showing heads when it is flipped.   4) We pick a coin from the bag, flip it $K$ times and note $h_i = $number of times $i$ came up heads. Then we put the coin back into the bag.   5) We repeat step 4) $M$ times. It is possible that a coin is drawn and flipped again, but I would recognise it. I want to estimate the population mean $$ p= \sum_{i \in N} q_i \cdot p_i $$ and give a confidence interval. I know that I could estimate $p_i$ with $\hat{p}_i = \frac{h_i}{K}$, since $h_i$ is binomially distributed with parameters $(K,p_i)$. Of course this is only possible if $i$ was picked from the bag at some point, and if $i$ was picked more than once, I would have to add the $h$'s and divide by $2K$, $3K$, etc. I realize that $$ \hat{p} = \frac{1}{M} \sum_{m=1}^{M} \hat{p}_m $$ is an unbiased estimator for $p$, where $\hat{p}_m$ is the estimated parameter for the coin drawn in the $m$th round. But I don't know how to construct the confidence interval around it. In my dataset,   - $N = 10^{10}$ (really large),   - the $q_i$ are computable, but I would rather not have to compute them if I can avoid it,   - $K=1000$,   - $p_i = 0$ more than half the time,   - $M=10000$.


title: How to test changes in call volume over time?

text: I have carried out a study on change of call volume over a year (oct 07-mar 09). What statistical test on SPSS software can I use to show the difference in call volume on a day by day basis and on a monthly basis?


title: Which are the most effective methods for selecting independent variables?

text: Some clustering algorithms require independence of variables but (especially working with real data) variables are often highly correlated. I have been suggested to apply a Principal Component Analysis to solve the problem (using the resulting linear combination as input for the clustering). But such a solution may affect the resulting clusters since correlated varibles "repeat" the same information in the linear combination. I.T.Jolliffe (1972; 2002) suggest a method of Variables Reductions (called 'method B4') that select one variable for each Component (the one with highest loading). Even if other method have been considered more effective (specially if the goal is to mantaine as much variance of the original variables as possible), Joliffe's method B4 is the only one -- that I have found -- that ensures the independent of the selected variables. **Question:** Are there other more effective methods for selecting independent variables?


title: Difference among bias, systematic bias, and systematic error?

text: Is there any difference among the following terms or they are same?   1. Bias    2. Systematic bias   3. Systematic errors If there exist some differences then, please explain them. Can these errors be reduced when one increase the sample size? **UPDATE:** My field of interest is statistical inference. I mean to say that how we differentiate these term as a statistician.


title: Usable estimators for parameters in Gumbel distribution

text: The Gumbel distribution has the general form: $$F(y)=\exp\left({-\exp{\left(-\frac{y-\mu}{\sigma}\right)}}\right), \quad y\in \mathbb{R}$$ where $\mu \in \mathbb{R}$ and $\sigma >0$. Let $W_1,...,W_n$ be random variables with that distribution and iid. I want to find good estimators for $\sigma$ and $\mu$. What would be the best way to do that?


title: Comparing values over thresholds at several stations -choosing the right method

text: I have a data set containing daily temperatures **T** at several stations and am interested only in whether **T** crosses a certain threshold or not. In the end, I want to answer e.g. the question: "If **T** crosses the threshold at stations 1 and 2, the probability of **T** crossing the threshold at station 3 is XX." I'm not yet interested in the spatial characteristics (I'd assume that you would then need to weigh the stations according to distance somehow), and the threshold is also not high enough that I could think of using extreme value theory. Is a multinomial logistic regression the way to go to model the relationship between the stations? Or would it be preferable to do a simple binary logistic regression for each pair of stations? If I take the actual number of days with **T** over the threshold per year and compare this among the stations, would a poisson regression be the correct choice?


title: Fitting a Poisson distribution with lme4 and nlme

text: I am now looking for a GLMM, which could fitted a Poisson distribution with a log-link. From what I see until now, **lme4** allow to specify the family and the link function for `lmer()` model, but the `lme()` function in **nlme** package doesn't. Is there another way to specify it in **nlme**? May we also do it when fitting a non linear mixed models with both pacakges? Thanks in advance, your help, would be appreciated.


title: Are there any associative multi label classification implementations available?

text: I have seen that it is possible to perform multi label classification using a binary combination of classifiers or reducing a multi label classification to a multi class classification problem by transforming the tags, and in fact this question has been posed around many times. However, I was unable to find anything related to associative multi label classification implementations. Are there any implementations of this? All packages I have seen do not mention about this approach. Thanks!


title: Finding Percentage Point of a Normal Distribution

text: Good evening all, I am doing a self-study exercise, but have been quizzed by a part of the question on finding percentage points of a normal distribution. I fully understand the first part of the question and was able to find the answer, which corresponds to the answer provider. ![enter image description here](http://i.stack.imgur.com/QKuD3.jpg) For the part on "If 80% of the trucks will travel more than k-thoudsand miles, find k", I was unable to understand the answer. My working is as follows: X ~ N(50, 144) P(X > k) = 0.8 P[Z > (k-50)/(12)] = 0.8 This is something that veer far off from what the answer provided, which I wasn't able to concur or understand. The self-study guide uses the following table to get Alpha. I was wondering if the solution provided was because of the limit, thus the need to find 0.2? Even if so, I still do not. understand the need to invert the signs from P[Z > (k-50)/(12)] = 0.8 to P[Z < -(k-50)/(12)] = 0.8 ![enter image description here](http://i.stack.imgur.com/74iZJ.jpg) Appreciate your kind help and advice please.


title: What is the best software for neural network analysis?

text: I need to fit a neural network to my data. I used SPSS 21, but I think it's model summary isn't sufficient. I also used R (nnet & neuralnet), it has a poor report and I can't find a good model summary for it. The last one is NeuroSolutions, I used it. very bad! I can't run a model. What are those graphical icons in it's window after bulding model? You say! What's your advice for Neural Network Analysis?


title: Links between "Interaction Information" and interactions in a linear regression model.

text: I'm trying to get a handle on interactions between three variables in a regression-sense and how they differ from interactions in an information theory sense. Imagine we wish to determine whether two factors X and Z are interacting in the context of Y. For the mutual information based-measure, we use interaction information, expressed as \begin{aligned} I(X; Z; Y) &= I(X; Z|Y)+I(X; Y) \\ &= I(X, Z; Y) - [ I(X; Y)+I(Z; Y)] \end{aligned} where I is mutual information. Here, the interaction information can be interpreted the information expressed by X and Z together in the context of Y minus the information by X and Z separately in the context of Y. The more common approach is to use a regression based approach. Here, we compare a linear model with terms for X and Z and their interaction (denoted here as interaction model) $$ Y=\beta_1X+\beta_2Z+\beta_3XZ $$ to a model with only the X and Z factors as individual terms (denote as additive model) $$ Y=\beta_1X+\beta_2 Z $$ To test for interaction, we may use a likelihood ratio test to determine whether the likelihood of the interaction model is better than that of the additive model. That is, whether Y can be modelled better if we examine the interaction of X and Z as opposed to just modelling the additive effects of X and Y. Is there any link between these two approaches to quantifying interaction? Are they defining interaction in a similar way? To my eyes they appear to be based on two very different underlying intuitions about what constitutes an interaction. Given that the wikipedia page on interactions is pretty clear that interactions are typically defined as deviation from additivity, I'm curious to know whether interaction information can be fit into that framework.


title: Reduced chi square goodness of fit

text: I often use the reduced chi squared as a quick goodness of fit test when fitting histograms. Is there an analogous method that works more generally for interpreting the absolute value of the best fit likelihood when data are not normal? I'm thinking of binomial and poisson with noise here but the more general the better. I guess what I might be fishing for is an answer to this question. Relation between reduced chi square and Akaike criterion


title: Expectation maximization and Gaussian mixtures - bad results

text: I am supposed to find parameters of individual gaussians in a 1D mixture with a known number of components. I use my own implementation of EM algorithm; however, I am not able to find the right solution. ### Data Number of points: 5057   Number of components: 5 ![Approximate components](http://i.stack.imgur.com/i5b4X.png) ### Problem Everytime I launch the algorithm, all resulting component's means are values from 161 to 163 (That is between the first two and largest peaks). Here are my results from individual steps of one run of EM: output   Here is my implementation of EM in C++: EM ## Update Pat, you were right, thank you for your help.   1. I forgot to multiply each item of product by posterior probability   2. I had mistake in log likehood function


title: SPSS from two - to one tailed test

text: If any of you would like to confirm if my understanding is correct or wrong,   I will be most grateful. I have done an independent-samples t-test (two-tailed),   for difference between women/men how they see their own skill sett. In short:   group women,n 76: mean 3,20 - std .980   group men,n 21: mean 3,71 - std .956   df= 95 sig.level set to 0.05 t= - 2.150 sig. (2-tailed) 0.034 (critical value= -1.98 to 1.98) statistics significant difference between men and women ...   I reject the null hypothesis. I am going to find out the the results for a one-tailed   test which SPSS do not check for. H0: There is no difference in drivingskills between women and men.   H1:Women have lower driving skills then men.   I have understood that the principle is the same, but I have to   first be clear on the direction, which is negative (group 2 have a higher mean).   Then I need to divide the significant in 2 which was -2.150/2=-1.075. The area at 0.05 with 95 df one-tailed test in the   "students distribution" shows me about -1.661. I conclude therefore that I must accept the null hypothesis   when doing a one-tailed test (since I have the number - 1,075) Am I very lost here or have I understood it and done it correctly?   (She said with hope). Thank you.


title: Is standardisation before Lasso really necessary?

text: I have read three main reasons for standardising variables before something such as `Lasso` regression: 1) Interpretability of coefficients. 2) Ability to rank the coefficient importance by the relative magnitude of post-shrinkage coefficient estimates. 3) No need for intercept. But I am wondering about the most important point. Do we have reason to think that standardisation would improve the out of sample generalisation of the model? Also I don't care if I don't need an intercept in my model; adding one doesn't hurt me.


title: How do i predict with standard errors using betareg package in R?

text: I'm using 'betareg' package in R to perform beta regression. predict() function with se.fit=T is supposed to return standard errors along with the prediction but it doesn't. Is there any other way I can get the standard error outputs? I'm open to using other packages that can perform beta regression too.


title: Reliability of a network

text: > Consider the communication network shown in the figure below and suppose > that each link can fail with probability **_p_**. Assume that failures of > different links are independent: ![enter image description here](http://i.stack.imgur.com/nC10Q.jpg) > All links are independent. Find the probability that there exists a path > from A to B along which no link has failed. I used this graph for my solution: ![enter image description here](http://i.stack.imgur.com/gghIA.jpg) > Given that exactly one link in the network has failed, find the probability > that there exists a path from A to B along which no link has failed. And for this question, which is my headache, based on my last graph I built this new graph with the solution there indicated. ![enter image description here](http://i.stack.imgur.com/TFkI9.jpg) * * * I've tried to find a book about reliability of networks for more theory and examples but I've had no luck, so I'm developing my solution intuitively but I need to know if I'm in the right path. Anyone can advise me?


title: Comparing the mean square error of a Elastic net and neural network

text: I'm trying to determine the best model for a regression problem between Elastic net regression and a neural network model based on the mean squared error.When I'm fitting the data to the elastic net I do not perform any scaling on the targets, but for the neural network model the target values are scaled to be in the range [-1,1],when I calculate the mean squared error for the elastic net I get errors in the range(order) of say 10^4 and for NNs in the order of 10^-2.   * Is there any measure that I can calculate to bring the two error measure into comparable range (like a percentage or relative error)    * Is is appropriate if I did one of the following: 1.Find the root mean squared error and apply the inverse operation of the scaling operation that was used on the target values so as to get a value comparable to elastic net error 2.Scale the target values of the data to be in the range [-1,1] so the error of the elastic net is comparable to that of NN's.


title: Confusion on the relationship between laplace distribution and exponential distribution

text: If $X \sim Laplace(\mu, b)$ and $\mu = 0$ such that $X \sim Laplace(0,b)$, it becomes a form of the exponential family. Unfortunately, one source tells me that it becomes $|X| \sim exponential(b)$, another source says that $X \sim exponential(\frac{1}{b})$, and my notes say that $X \sim exponential(b)$. Can any of y'all help me out on this?


title: Probability estimation of minesweeper undecidable choices

text: ![enter image description here](http://i.stack.imgur.com/PILnh.png) Someone asked what was his probability to win in this situation? This is an interesting problem and I'll be curious to know an efficient solution. So far, I wanted also to generalize, with a simplification, we consider only:   * pairs of blocks such as the one on top right, and model them with a binary variable which value is the presence of the mine   * twin pairs, such as the one on bottom left, which is a couple of pairs linked together, the knowledge of one value of a pair gives the value of the other pair and vice versa So we omit groups of more than 2 pairs linked together, for simplicity. Now in this simple example, we have 2 single pairs and 2 twin pairs, let's denote twins with binary variables: $A_1, A_2$, and the other one with $B_1, B_2$ and the other single pairs with $C$ and $D$ So we could think the probability is $(\frac 1 2)^4$ but I'm not sure at all The goal is to generalize with m twins pairs and n single pairs


title: How can I test if a marketing effort was successful?

text: I looked at a test group (West Coast pop) versus a control group (non West Coast) and took independent samples pre-message and post message for each group. I measured awareness (top-3 box of a 5-pt scale) for each of the 4 independent samples. I looked at lift (in percentage points) for test group and lift for the control group. What is the statistical test for significance between test and control groups to determine whether or not the marketing messaging got through?


title: Specifying a multilevel model in MCMCglmm (R), that is heteroskedastic at level one

text: I am considering MCMCglmm as an alternative to MLwiN. The former package works perfectly fine, but I cannot figure out how to model heteroskedasticity at level one. For instance, if I have the following model, individuals within countries               Fixed: y = 1 + x1 + x2      Random Level 2: 1 + x1 (country level)     Random Level 1: 1 + x1 (individual level)      For level 2 I can specify the option `random=~country+x1`, how do I specify the level one portion? My guess would be that this has something to do with the `rcov` option, but specifying the formula as `rcov=~indiv.id+x1` doesn't work. Or am I entirely wrong and these type of models are just not fitted in MCMCglmm? Thanks for any tips. EDIT: to provide a clearer idea of what the model is, here is how the model looks in MLwiN. At the top there are fixed effects, at the bottom you have two variance matrices: Ou is variance on the country level, and Oe is variance on the individual level. Ou is easily modeled in MCMCglmm, using the `random` option. I cannot get MCMCglmm to estimate the Oe.


title: How to report the results of Fisher's exact test for larger than 2 by 2 tables

text: What is an acceptable phrase for reporting a Fisher's exact test where the number of rows and the number of columns are both $>2$. An example in R:               > x <- as.matrix(data.frame("A"=c(10,4,6),"B"=c(45,2,6),"C"=c(23,4,20),              row.names=c("df1","df2","df3")))     > x         A  B  C     df1 10 45 23     df2  4  2  4     df3  6  6 20     > fisher.test(x,, simulate.p.value=TRUE, B=1e5)          Fisher's Exact Test for Count Data with simulated p-value      (based on 1e+05 replicates)          data:  x      p-value = 0.00026     alternative hypothesis: two.sided       I'm looking for a phrase / sentence that would encapsulate the above.


title: Association between a categorical and binary variable

text: I have a set of patients that I split into two parts using some obscure algorithm (doesn't matter how): let the $X$ take on values $A$ and $B$. Each patient belongs to exactly one of 5 disease subtypes: let $S$ take on the values $a, b, c, d, e$. Here is the contingency table:                      A    B     a    116  109     b     64   87     c    161   86     d    123  140     e     77   38      Using Fisher's exact test, I get a p-value of $5.274e-07$. So I conclude that $X$ and $S$ are somehow associated (or not independent). What I'd like to ask next is which subtypes are significantly associated with $X$. For example, subtype $a$ doesn't seem to deviate much from the expected number of $A$s and $B$s if $a$ is independent of $X$, while there seems to be an overrepresentation of $A$s within the $e$ subtype. What kind of test would be appropriate here? Could I use Fisher's exact test again for each of the subtypes separately. Or should I use something like logistic regression?


title: Degrees of freedom in my mixed model

text: I am having a serious issue trying to figure out why my degrees of freedom are screwy when running my proc mixed model. I have nested data for a MLM analyses, but SAS isn't providing the correct df and I'm not sure why. I have 50 individuals with 3 or 4 time sensitive datapoints each. Measuring their physical activity per day for 4 days, I have 177 observed instances of physical activity, however SAS keeps outputting 77 df when running ddfm=bw. I have cleaned the data and reimported, but I'm just at a loss for what to do now. I'd really appreciate any help. Thanks!


title: How can I explain the intuition behind ANOVA?

text: I need to explain the intuition behind what ANOVA is doing to a non-technical person. Is there a visual that explains the idea? A visual that illustrates the key idea in the context of a one-way ANOVA with perhaps 3 factor levels might be helpful? Let us suppose that the person has taken some statistical courses as a student in the distant past but has forgotten the details of even performing a z-test. However, he/she remembers that hypothesis testing is used to check if the observed effects are due to random chance or due to a real change in the parameter of interest.


title: Why do coxph() and survreg(, dist="exponential") NOT return the same coefficient (except for expected opposite sign) in R?

text: If I understand correctly, the coefficient of a covariate `X` under a Weibull accelerated failure time (AFT) model is related to the `log(hazard ratio)` of the Cox proportional hazard (PH) model in the following way:               d1 = -c1/shape.gamma      where `d1` is the AFT coefficient, `c1` is the Cox PH coefficient `log(hazard ratio)`, and `shape.gamma` is the shape parameter of the Weibull distribution. See for example slide 36 in this handout. So for `shape.gamma=1` (the Weibull distribution is then equal to the exponential distribution) `d1` and `c1` should be equal (except for sign). My question is: Why, using the survival package in R, do `coxph()` and `survreg(, dist="exponential")` NOT return the same coefficient estiamtes (in absolute value)? Below is the code I used to check this. I first simulate uncensored data using the exponential AFT parametrization and, as expected, `survreg()` returns 'on average' `d1=0.5` (the value I used to simulate the data). So far so good, However, the coefficient returned by coxph is quite different and seems to depend very much on the `scale.lambda` parameter of the exponential distribution. How does this mismatch relate to the theoretical equivalence (in absolute value) of the coefficients of the exponential AFT model and the Cox PH model? What am I missing?               library(survival)       # sample size     n            <- 1000     scale.lambda <- 2     d1           <- 0.5     ph.ests      <- c()     aft.ests     <- c()     X1           <- c(rep(0,n/2),rep(1,n/2))     for(i in 1:100){      e        <- rexp(n,scale.lambda)        # AFT model      logT     <- -d1*X1 + e      T        <- exp(logT)        # no censoring      event    <- rep(1,n)         # cox ph analysis      m.ph     <- coxph(Surv(T,event)~X1)        # exponential AFT analysis      m.aft    <- survreg(Surv(T,event)~X1,dist="exponential")        # collect X1 coefficients (on log scale)      ph.ests  <- c(ph.ests,coef(m.ph))      aft.ests <- c(aft.ests,coef(m.aft)[2])     }       # summary of 100 replications     summary(aft.ests)     summary(ph.ests)


title: Why do ANOVA p and F values change?

text: I am currently conducting a two-way ANOVA on a dataset and have noticed how significance and F values change dependent on the type of ANOVA I use in R. Why is this? For instance, if I run a two-way ANOVA the factor I am analysing becomes outputs F (1,28) = 7.2, P <.05. When I run it by itself however in a one-way ANOVA it becomes F (1, 30) = 21.5, p < .001. If someone could explain to me why the difference, that would be much appreciated.


title: Identifying significantly different groups in a logistic regression model

text: In my experiment, I compare the resistance of 12 different seed types to seed beetles. I have 30 replications (i.e. 30 times a single beetle larva developing in a single seed) for each variety and assess whether an adult beetle emerges from the seed or not. The resulting two-way table I analyzed using logistic regression with the factor "variety" and the binary response "development completed / not completed". The analysis revealed that the groups are significantly different. Now I'd like to do a _post-hoc_ test to identify which varieties are significantly different from one another (similar to a Tukey's HSD for ANOVA). Is there a way to do this? Thank you very much for your support!


title: Test if two measurements come from normal distributions with the same mean?

text: I have two measurements $x_1\pm s_{x_1}$ and $x_2 \pm s_{x_2}$. The error intervals ($\pm$ one standard deviation) do not overlap. Does this necessarily mean that these are two samples from normal distributions with different means? Or should I run a test to exclude (to some degree) the fact that they might still be from the same distribution? Both means have an equal amount of measurements (about 2000). It is given that both are pulled from normal distributions.


title: Which notation and why: $\text{P}()$, $\Pr()$, $\text{Prob}()$, or $\mathbb{P}()$

text: Are these merely stylistic conventions, or are there substantive differences in the meanings of these notations? Are there other notations meaning " _the probability of_ " that should be considered in this question?


title: Disaggregation of spatial autocorrelation parameter

text: I have data aggregated at state level. When I estimate a spatial autoregressive model such as $$y = \rho W y + X\beta + \epsilon$$ on this data, I see that the autoregressive parameter $\rho$ is moderately significant ($\rho = 0.2$, $p$-val$= 0.06$). But I am preparing to apply this estimated model to data aggregated at county level; to be precise, I have ${X, W}$ at both state and county levels, but $y$ only at the state. As I understand spatial autoregression, there should be a significantly higher correlation as I move down in geography. _A priori_ I would expect somewhere around 0.8, based on my prior experience. This project is for a practical model rather than research, so a precise mathematical solution is not necessary (and I imagine doesn't exist). I am interested if others have run into this problem or if what I am doing is recklessly irresponsible...


title: Are HAC estimators used for estimation of regression coefficients?

text: The references I can find on HAC procedures (like Newey-West) in regression focus on the standard error of the estimated regression coefficients and hypothesis testing involving the same. I cannot find any references where the HAC estimate of error covariance is used in a generalized least squares procedure to get 'better' estimates of the regression coefficient. Is this not done? Or are the hypothesis tests just more interesting?


title: Odd error with caret function rfe

text: library(caret)     set.seed(1)      x <- data.frame(runif(10),runif(10))     y <- rnorm(10)     rfeModel <- rfe(x,y,rfeControl = rfeControl(functions = lmFuncs))      returns: > Error in { : task 1 failed - "undefined columns selected" While this code:               data(BloodBrain)          x <- scale(bbbDescr[,-nearZeroVar(bbbDescr)])     x <- x[, -findCorrelation(cor(x), .8)]     x <- as.data.frame(x)          lmProfile <- rfe(x, logBBB,rfeControl = rfeControl(functions = lmFuncs))      Works fine. What am I doing wrong in the 1st example? (I get the same error in R 2.14.0 and 2.13.2)


title: Most efficient way to check answers to exercises when learning from a statistics textbook?

text: There's a community wiki at: Free statistical textbooks that has a nice list of freely available textbooks on statistics. Many of the textbooks that people suggest have exercises at the end of each chapter which is great, but very few of these textbooks have provided solutions to any of the exercises. Sure there may be a variety of reasons for this, but if you are just going through these texts on your own, not in formal course where exercises sets usually have a grader, then that presents an issue as how are you supposed to check your answers. Posting a new question on here to check to see if you got the right answer is not really practical I think. Like that these resources are already free, but, and I'm sure others would agree, it's frustrating for a beginner. e.g. How do i check my answers for @Rob Hyndman book: Forecasting: principles and practice (Hyndman & Athanasopoulos, 2012) List of solution manuals to books at Free statistical textbooks:   1. Elements of Statistical Learning by Hastie, Tibshirani, and Friedman A Solution Manual and Notes for: The Elements of Statistical Learning by JL Weatherwax


title: How to compare the variability of two measures from meta analysis?

text: I have two measures from each of 20 studies. Random-effects meta analysis has been performed for each of the two measures, and for each measure I have two values to indicate data variability, cross-study (random effect) variance $\tau^2$ and its chi-square test statistic $Q$. Now I'd like to compare the cross-study variability between the two measures. $\tau^2$ does not seem to be a good choice because it is dimensional (having different units) and carries the same physical meaning as its corresponding measure. Since the chi-square value $Q$ is dimensionless, I was wondering whether I could construct the ratio of the two $Q$ values as a pseudo $F$-statistic (with the numerator and denominator degrees of freedom being the same) to test for equal variability across subjects between the two measures. Or any better approach?


title: How to remove skewness from the population?

text: Just want to know the ways to remove skewness from the skewed population. I am working on logistic regression (in the insurance domain) for finding out the repurchase propensity of existing users. Income data is right skewed. Will it cause any problem for runing logistic regression? If yes then how to remove it (deletion or any other transformation)?


title: Finding parts of categorical variables that "go together"

text: Good morning all, I'm working on a project related to the occupations that university graduates are likely to go into. I have a large data set (N in the hundreds of thousands) where the unit of analysis is an individual person. For each person, I have two categorical variables - a code representing the field of their university degree and a code representing their occupation. I'm looking for a statistically valid way to find out what fields of study and occupations "go together." In other words, what courses of study prepare people for which jobs? So far, I've considered doing this with simple descriptive statistics... pull, say, the top 10 occupations for every subject area while ruling out occupations like cashiers, fast food workers, etc. But it would be great if there was some sort of more rigorous test that could be used for this. Perhaps one regression per occupation with tons of dummy variables representing each subject area? If a subject area dummy variable is statistically significant, then that subject area is connected with that occupation? I'm open to any and all suggestions.


title: Scores and normal distributions

text: I'm in the process of designing a system. So basically, every day users on my site rack up points based on their activity on my website. So whatever they do on my website, they rack up points. So on 8/3, a user might do multiple things on my website so that they rack up 10 points. Every day, I figure out the number of points for the user for that day and add a row in the database with the number of points for that day. I want to give each person a score. Basically a score from 0 to 10 just like in school, but I want that score to depend on the mean score of all the users. I was thinking of finding the average number of points per day for each user (so user 1 might have an average of 10 points per day) and then finding the average ($X$) of these averages. And using $X$ to create a normal distribution to give each person a score from 0 to 10. But if the person is 2 standard deviations from the mean, how does that translate to a score? I'm not a statistician. Any advice would be greatly appreciated. What's the best way to calculate a score like this? Also I know with normal distributions that the distribution has to be normal. What if the average scores are right skewed? What should I do? Determine whether there is an outlier? Anyway help would be appreciated.


title: Decision trees for advertising data

text: Assuming a dataset with the following attributes:   `Date` (truncated), `f1` ... `fn`, `#impressions`, `#goals`. The problem: I want to grow **$n$ trees** that would find the optimal selection of features and their ranges in each, and that maximize the goal rate (goals / impressions) for future campaigns with **high confidence**. What would be the suggested approach? Thanks.


title: How to interpret the results of ADF test using SAS ARIMA?

text: The SAS System      14:11 Thursday, October 6, 2013   1                                                The ARIMA Procedure                                          Name of Variable = ln_G_S_Index                               Period(s) of Differencing                           1                          Mean of Working Series                       0.094293                          Standard Deviation                           0.316757                          Number of Observations                             15                          Observation(s) eliminated by differencing           1                                                       Autocorrelations           Lag    Covariance    Correlation    -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1      Std Error             0      0.100335        1.00000    |                    |********************|             0        1     0.0026693        0.02660    |          .         |*        .          |      0.258199        2     -0.018517        -.18456    |          .     ****|         .          |      0.258382        3      0.029440        0.29342    |         .          |******    .         |      0.267025                                           "." marks two standard errors                                                   Inverse Autocorrelations                         Lag    Correlation    -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1                           1       -0.14763    |          .      ***|         .          |                      2        0.19526    |          .         |****     .          |                      3       -0.27516    |          .   ******|         .          |                                                   Partial Autocorrelations                         Lag    Correlation    -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1                           1        0.02660    |          .         |*        .          |                      2       -0.18539    |          .     ****|         .          |                      3        0.31522    |          .         |******   .          |                                              Phillips-Perron Unit Root Tests                        Type           Lags         Rho    Pr < Rho        Tau    Pr < Tau                        Zero Mean         0    -11.6883      0.0066      -3.23      0.0033                                     1    -11.4504      0.0074      -3.23      0.0034                   Single Mean       0    -13.7527      0.0129      -3.71      0.0171                                     1    -12.6667      0.0218      -3.76      0.0157                   Trend             0    -14.5288      0.0601      -3.25      0.1144                                     1    -13.1531      0.1022      -3.20      0.1239                                             The SAS System      14:11 Thursday, October 6, 2013   2                                                The ARIMA Procedure                                          Name of Variable = ln_G_S_Index                               Period(s) of Differencing                           1                          Mean of Working Series                       0.094293                          Standard Deviation                           0.316757                          Number of Observations                             15                          Observation(s) eliminated by differencing           1                                                       Autocorrelations           Lag    Covariance    Correlation    -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1      Std Error             0      0.100335        1.00000    |                    |********************|             0        1     0.0026693        0.02660    |          .         |*        .          |      0.258199        2     -0.018517        -.18456    |          .     ****|         .          |      0.258382        3      0.029440        0.29342    |         .          |******    .         |      0.267025                                           "." marks two standard errors                                                   Inverse Autocorrelations                         Lag    Correlation    -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1                           1       -0.14763    |          .      ***|         .          |                      2        0.19526    |          .         |****     .          |                      3       -0.27516    |          .   ******|         .          |                                                   Partial Autocorrelations                         Lag    Correlation    -1 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 1                           1        0.02660    |          .         |*        .          |                      2       -0.18539    |          .     ****|         .          |                      3        0.31522    |          .         |******   .          |                                           Augmented Dickey-Fuller Unit Root Tests              Type           Lags         Rho    Pr < Rho        Tau    Pr < Tau          F    Pr > F              Zero Mean         0    -11.6883      0.0066      -3.23      0.0033                           1    -12.4302      0.0041      -2.42      0.0197         Single Mean       0    -13.7527      0.0129      -3.71      0.0171       6.91    0.0157                           1    -25.2133      <.0001      -3.63      0.0214       6.59    0.0206         Trend             0    -14.5288      0.0601      -3.25      0.1144       6.44    0.0799                           1    -45.0252      <.0001      -3.20      0.1265       6.92    0.0622                                             The SAS System      14:11 Thursday, October 6, 2013   3


title: Analyzing Likert scales

text: The Likert scale I am analyzing has both positive and negative questions, but all the items are answered with 4 - strongly agree, 3 - agree, 2 - disagree and 4 - strongly disagree. For analysis, should I make the weight of the answers the same? For example, if the statement is positive, then 4 must be given to strongly agree, and then for negative statements, a 4 should also be given to strongly disagree. What else can I do aside from Mean and Standard Deviation?


title: p-Value and results of Cramer's V

text: I need a bit help for interpreting this results... I did a correlation analysis with R with the assocstats function. The result is:               $summary     Call: xtabs(formula = ~MH[, i] + MH[, j], data = MH)     Number of cases in table: 2306      Number of factors: 2      Test for independence of all factors:         Chisq = 2806.6, df = 3318, p-value = 1         Chi-squared approximation may be incorrect          $object                         X^2   df P(> X^2)     Likelihood Ratio 1036.1 3318        1     Pearson          2806.6 3318        1          Phi-Coefficient   : 1.103      Contingency Coeff.: 0.741      Cramer's V        : 0.637           $summary     Call: xtabs(formula = ~MH[, i] + MH[, j], data = MH)     Number of cases in table: 2343      Number of factors: 2      Test for independence of all factors:         Chisq = 118.46, df = 73, p-value = 0.000611         Chi-squared approximation may be incorrect          $object                         X^2 df   P(> X^2)     Likelihood Ratio 130.83 73 3.8115e-05     Pearson          118.46 73 6.1100e-04          Phi-Coefficient   : 0.225      Contingency Coeff.: 0.219      Cramer's V        : 0.225       But what does the p-value of 1 in this case mean if Cramer's V is 0.637 and therefor show a strong corelation or with 0.000611 with a Cramer's V of 0.225? And does anyone of you know a good rule of thumb for interpreting Cramer's V? Thanks for your help.


title: Can I subsample a large dataset at every MCMC iteration?

text: I have a large dataset from which I want to perform a bayesian probit regression using Gibbs sampling 1. Since the dataset has one milion rows, and variables from a truncated normal must be sampled [2], the whole process is too slow. I would consider variational or parallel approaches, but before going that far I would like to know whether I could randomly sample (with replacement) from my dataset at every Gibbs iteration, so that I have less instances to learn from at every step. My intuition is that even if I change the samples, I would not be changing the probability density and therefore the Gibbs sample should not notice the trick. Am I right? Are there some references of people having done this? 1 Albert and Chib 1993. Bayesian Analysis of Binary and Polychotomous Response Data. [2] I use python and the function scipy.stats.truncnorm


title: How to statistically and graphically compare distributions for nine groups where group sample sizes are unequal?

text: I conducted a study in which I collected the height, width, and weight of individuals from 9 different races (i.e I collected 10 to about 30 samples per race).   * How can I can display the distribution of the data for each group for a given variable (e.g., height)? I have thought of using histogram but then how do I handle the missing value ranges.   * What other statistical approaches and visualization would be suitable for examining similarity in distribution?   * How could these approaches be implemented in R?


title: Why are two correlations significantly different to zero but not significantly different from one another?

text: Is it 'normal' to have two Pearson correlation coefficients which are statistically significant on their own but **not** when compared with each other. **Example**   * r=.9083 (n=45) for females who buy ice-cream and yoghurt    * r=.8382 (n=60) for males who buy ice-crean and yoghurt (The above correlation is between buying ice-cream and buying yoghurt for the two groups) Based on the above, my interpretation is that both males and females who buy more ice-cream also tend to buy more yoghurt but there is no difference whether one group buys more or less than the other group.


title: Comparison of multiple regression models

text: I have one dependent variable and many predictors, and I need to use a multiple regression model (linear). Now, I performed a stepwise regression to determine which independent variables to include in the final model. However, I see that when results are published, people usually show more than one model. For instance, they include just the first variable and then the first and the second, then the first, the second and the third, and so on. Showing that the adjusted coefficient of determination, $R^2$, improves. Do I always have to do that, in order to show that each variable has a contribution (of course if coefficients are significant)? Is there any reference to a formal procedure that better explains this?


title: How to construct an optimally balanced reference corpus for text categorization?

text: I have a rather large collection of text documents categorized into about 150 categories. While some categories are represented by several thousands of documents, others have only a few hundreds assigned to them. Now I would like to construct a balanced corpus from this data, where each category is represented by the same number of documents and which at the same time maximizes the categorization accuracy. I tried to randomly select documents from the stronger categories but would like to know if there is a more systematic way.


title: What is the best way of weighing cardinal scores and Likert scale scores to create a composite score?

text: I can't thank the experts enough for their clarifications. One final question following my earlier posts on forming coomposite variables here, here, and here. If I have measured some variables as cardinal numbers (e.g. 1, 2, 3, 4, 5 times) and others on a Likert scale, can I then just add the scores together? For example:   1. How many times do you visit the local hospital per month? Answer: 5 times    2. How would you rate the medical service provided to you: Answer 3 on a scale of five going from poor to excellent. Can I add 5 and 3 and say the composite score is 8? What if the likert scale goes from negative to positive (e.g. strongly disagree = -2, disagree = -1, neither disagree or agree = 0, agree = 1 and strong agree = 2)? Would the same rule as above apply if the score is say 0 or -2? My objective is to keep the aggregation process simple and not get 'confused' in complex statistical formulas.


title: Is there an R function that will compute the cosine dissimilarity matrix?

text: I would like to make a heatmap with row clustering based on cosine distances. I'm using R and `heatmap.2()` for making the figure. I can see that there's a `dist` parameter in `heatmap.2` but I cannot find a function to generate the cosine dissimilarity matrix. The builtin `dist` function doesn't support cosine distances, I also found a package called `arules` with a `dissimilarity()` function but it only works on binary data.


title: Best way to select useful features using R software

text: I have a huge matrix (individuals X features with row.names as individuals numbers) and the corresponding segment in another vector of 1D (row.names are the same as in my huge matrix and the vector represent the segments associated). I.E. :               row.names VAR1 VAR2 VAR3 VAR4 … VAR3000         12     4    12    5   18      8         58     6    13    19   3     10      for the huge matrix. and:               row.names  x         12     4         58     2      for the segment representation (where x represent the individual' segment). I have no a priori model and I want to select a subset of variables (variable/feature selection) in order to predict the segment using a minimal subset of variables. I didn't use biclustering technique to detect my classes but a simple-way one. Which technique would you recommend to :   1. select the most discriminative variables (lasso, elastic net) and why?   2. predict the segment from these variables.   3. predict multiple values in another similar matrix (same individuals, few predictors that have been selected). Is it possible in this case to use correlation matrix (or cov) to infer directly the values of predictors that are not known in the other matrix (not using the following method: predict the class, then fill the missing values with medoid values or cluster-mean values)? Thanks in advance.


title: Rolling price returns in a linear regression

text: I want to conduct a linear regression (in matlab) using rolling monthly returns; the aim is to give me a prediction for the next monthly rolling period return. _return calculation_ : $$\mbox{return}(t) = \dfrac{\mbox{Price}(t) - \mbox{Price}(t-30)}{\mbox{Price}(t-30)}.$$ _regression_ : $$\mbox{return}(t+1) = a + b_1f_1 + b_2f_2 +b_3f_3+ e.$$ **My question is** what is the best way to conduct a linear regression using a rolling return with a time horizon greater than $1$ day? Thanks!


title: How can I improve my analysis of the effects of reputation on voting?

text: Recently I had done some analysis of the effects of reputation on upvotes (see the blog-post), and subsequently I had a few questions about possibly more enlightening (or more appropriate) analysis and graphics. So a few questions (and feel free to respond to anyone in particular and ignore the others):   1. In its current in incarnation, I did not mean center the post number. I think what this does is give the false appearance of a negative correlation in the scatterplot, as there are more posts towards the lower end of the post count (you see this doesn't happen in the Jon Skeet panel, only in the mortal users panel). Is it innapropriate to not mean-center the post number (since I mean centered the score per user average score)?    2. It should be obvious from the graphs that score is highly right skewed (and mean centering did not change that any). When fitting a regression line, I fit both linear models and a model using the Huber-White sandwhich errors (via `rlm` in the MASS R package) and it did not make any difference in the slope estimates. Should I have considered a transformation to the data instead of robust regression? Note that any transformation would have to take into account the possibility of 0 and negative scores. Or should I have used some other type of model for count data instead of OLS?   3. I believe the last two graphics, in general, could be improved (and is related to improved modelling strategies as well). In my (jaded) opinion, I would suspect if reputation effects are real they would be realized quite early on in a posters history (I suppose if true, these may be reconsidered "you gave some excellent answers so now I will upvote all of your posts" instead of "reputation by total score" effects). How can I create a graphic to demonstrate whether this is true, while taking into account for the over-plotting? I thought maybe a good way to demonstrate this would be to fit a model of the form; $$Y = \beta_0 + \beta_1(X_1) + \alpha_1(Z_1) + \alpha_2(Z_2) \cdots \alpha_k(Z_k) + \gamma_1(Z_1*X_1) \cdots \gamma_k(Z_k*X_1) + \epsilon $$ where $Y$ is the `score - (mean score per user)` (the same as is in the current scatterplots), $X_1$ is the `post number`, and the $Z_1 \cdots Z_k$ are dummy variables representing some arbitrary range of post numbers (for example $Z_1$ equals `1` if the post number is `1 through 25`, $Z_2$ equals `1` if the post number is `26 through 50` etc.). $\beta_0$ and $\epsilon$ are the grand intercept and error term respectively. Then I would just examine the estimated $\gamma$ slopes to determine if reputation effects appeared early on in a posters history (or graphically display them). Is this a reasonable (and appropriate) approach? It seems popular to fit some type of non-parametric smoothing line to scatterplots like these (such as loess or splines), but my experimentation with splines did not reveal anything enlightening (any evidence of postive effects early on in poster history was slight and tempermental to the number of splines I included). Since I have a hypothesis that the effects happen early on, is my modelling approach above more reasonable than splines? Also note although I've pretty much dredged all of this data, there are still plenty of other communities out there to examine (and some like superuser and serverfault have similarly large samples to draw from), so it is plenty reasonable to suggest in future analysis that I use a hold-out sample to examine any relationship.


title: Finding overall p-value for GLS model

text: Does anyone know how I can find/calculate an overall p-value for a GLS multiple regression model (made with nlme)?


title: Software for drawing bayesian networks (graphical models)

text: I am searching for [free] software that can produce nice looking graphical models, e.g. ![enter image description here](http://i.stack.imgur.com/QXAmU.png) Any suggestions would be appreciated.


title: Questions about how random effects are specified in lmer

text: I recently measured how the meaning of a new word is acquired over repeated exposures (practice: day 1 to day 10) by measuring ERPs (EEGs) when the word was viewed in different contexts. I also controlled properties of the context, for instance, its usefulness for the discovery of new word meaning (high vs. low). I am particularly interested in the effect of practice (days). Because individual ERP recordings are noisy, ERP component values are obtained by averaging over the trials of a particular condition. With the `lmer` function, I applied the following formula:               lmer(ERPindex ~ practice*context + (1|participants), data=base)       and               lmer(ERPindex ~ practice*context + (1+practice|participants), data=base)       I've also seen the equivalent of the following random effects in the literature:               lmer(ERPindex ~ practice*context + (practice|participants) +                      (practice|participants:context), data=base)       What is accomplished by using a random factor of the form `participants:context`? Is there a good source that would allow someone with just cursory knowledge of matrix algebra understand precisely what random factors do in linear mixed models and how they should be selected?


title: Removing outliers from asymmetric data

text: I have a data set that includes the number of visits to a website. Here are some descriptive statistics for my data Median: 4   Mean: 14.1352   SD: 121.8119 Clearly, there are some huge values (individuals who have visited the site thousands of times.) To remove these outliers I considered simply removing any data that falls outside more than 3.5 standard deviations from the mean. The result is that I discovered there is still a significant fat tail with my data. After removal of data that falls outside more than 3.5 standard deviations my descriptive statistics adjust to Median: 4   Mean: 10.2201   SD: 19.7492 I also explored using a winsorized mean but again since my data is asymmetric I feel like my descriptive statistics are biased. Is there a method that I can use to reevaluate my data to provide descriptive statistics that would represent a ‘majority’ of the population? As I understand the concept of bootstrapping, I could sample my population and then resample and generation thousands of populations that may represent my population differently based on the resample of my original sample population. Would this method be appropriate? Any other ideas or direction? Thanks in advance for any help. Any references or examples with R would be very much appreciated as well.


title: How to view GBM package trees?

text: I am trying to view the output from the GBM package for boosted trees in R. Below I am fitting a single tree without any sampling in order to compare the tree to the complete dataset. First, create the data set:                 set.seed(1973)          ############## CREATE DATA#############################################     N <- 1000     X1 <- runif(N)     X2 <- 2*runif(N)     X3 <- ordered(sample(letters[1:4],N,replace=TRUE),levels=letters[4:1])     X4 <- factor(sample(letters[1:6],N,replace=TRUE))     X5 <- factor(sample(letters[1:3],N,replace=TRUE))     X6 <- 3*runif(N)     mu <- c(-1,0,1,2)[as.numeric(X3)]          SNR <- 10 # signal-to-noise ratio     Y <- X1**1.5 + 2 * (X2**.5) + mu     sigma <- sqrt(var(Y)/SNR)     Y <- Y + rnorm(N,0,sigma)          # introduce some missing values     X1[sample(1:N,size=500)] <- NA     X4[sample(1:N,size=300)] <- NA          data <- data.frame(Y=Y,X1=X1,X2=X2,X3=X3,X4=X4,X5=X5,X6=X6)     ########################################################################          #Fit model##############################################################     gbm1 <- gbm(Y~X1+X2+X3+X4+X5+X6, # formula         data=data,                   # dataset         var.monotone=c(0,0,0,0,0,0), # -1: monotone decrease,                                      # +1: monotone increase,                                      #  0: no monotone restrictions         distribution="gaussian",     # bernoulli, adaboost, gaussian,                                      # poisson, coxph, and quantile available         n.trees=1,                   # number of trees         shrinkage=1,                 # shrinkage or learning rate,                                      # 0.001 to 0.1 usually work         interaction.depth=1,         # 1: additive model, 2: two-way interactions, etc.         bag.fraction = 1,            # subsampling fraction, 0.5 is probably best         train.fraction = 1,          # fraction of data for training,                                      # first train.fraction*N used for training         n.minobsinnode = 10,         # minimum total weight needed in each node         keep.data=TRUE,              # keep a copy of the dataset with the object         verbose=TRUE)                # print out progress          ###########################################################################      Next, look at the tree. This suggests, I think, a split on X2 at the value 1.5. However, this suggests 522 records one direction and 478 the other. Looking at the data, this record split does not correspond to the counts. Any insight? Is this a bug?               pretty.gbm.tree(gbm1,i.tree = 1)     length(d<-subset(data, data$X2>1.50,3)[,1])


title: Variable selection with principal components analysis

text: I ran principal components analysis in R on my data. All my regressors are continuous, non categorical variables, except gender which I excluded. I will add it and compare model 1 = PCA to model 2 = PCA + gender and see if its significant. I determined that PC1/2/3 determine 95% of the variance so I only consider the first 3 columns. All the variables have non-zero values and come up in the order I inputted it into the matrix of regressors. I used a spree plot to determine I need 3 regressors to explain most of the variance only (which is rather amazing given I have 30+ variables!). I am confused how I determine which variables I should choose. Do I order the variables per PC column from highest to lowest and see the variables contributing the most? I saw in the Freeway R guide that they actually plot the graphs and compare - similar graphs they take one variable and different looking graphs they subtract the variables from each other. This was rather confusing to me. **Questions:**   1. Is the ANOVA approach correct in this case for this variable?    2. Spree plot - do I look for the first 'kink' or the point where the x axis (no of variables) first goes to 0. Mine jumps dramatically down after 3 variables (to about 0.5) and then curves down to 0 on variable 15. Different websites say different things hence I'm confused.   3. Is the approach above correct to write a linear model? Do I write it as $y={\rm variable}_1 + {\rm variable}_2 + {\rm variable}_3$ or do I multiply them or both?


title: Modeling a birth-death process that is not memoryless

text: How does one approach the problem of modeling a "birth-death process" where the arrivals are dependent on the current state in the following way: if the population is above a certain point, the probability of an arrival decreases. Basically, I'm interested in complicating (slightly) an existing model of "births" that just has Poisson distributed arrivals and thinking of adding in the idea that there's a "saturation point" above which arrivals are less likely (waiting until the population drops back below the point). Should I be reading about nth-order Markov processes? Or should I be looking at queueing theory?


title: Methods of grouping sets of data

text: I have data which has several properties (metadata, as key value pairs, where the keyspace is shared over the whole dataset) per object. I took a sample of objects and divided them in n groups according to an unknown algorithm. What statistical methods or algorithm exists to find the relevant properties and their weight for the division to receive a similar grouping of the data than the "unknown algorithm" did?


title: Why do we make a big fuss about using Fisher scoring when we fit a GLM?

text: I'm curious about why we treat fitting GLMS as though they were some special optimization problem. Are they? It seems to me that they're just maximum likelihood, and that we write down the likelihood and then ... we maximize it! So why do we use Fisher scoring instead of any of the myriad of optimization schemes that has been developed in the applied math literature?


title: Formal test for heteroscedasticity

text: Are there any formal tests for heteroscedasticity for non-normal data? I want to run the test on time series logged returns, so would it be okay to assume a linear relationship? To me it makes intuitive sense that the greater the return (rise or fall in price) the greater the variance.


title: Three-way interaction in LMM: Can I write interpretation in this way?

text: I have a couple of questions regarding to the interpretation of LMM: Question 1: If I have three independent factors and two of them have three categories, I will have to run at least four models to see every angle of it, am I right? For example, I have color (red, blue, green), group (group1, group2, group3), and sex (male, female) as my independent variables. So in order to see every angle, I will have to 1) use red and group1 as baseline, 2) use red and group2 as baseline, 3) use blue and group1 as baseline, and 4) use blue and group2 as baseline. Question 2: If I have red and group1 as baselines and two parameters are significant, i.e. blue:group2 and green:group2, can I write interpretation in this way: The differences between blue and red (β=12, SE=8, p=0.04), and green and red (β=13, SE=6, p=0.03) are higher in group2 compared to group1. As normally we will put this parentheses at the end of the sentence, I am not sure about the case where I combine two parameters together in one sentence. Question 3: I would like to ask about the interpertation of three-way interaction. Suppose that I have red and group1 as baselines, and this parameter: (blue:group2:male) and (blue:group3:male) turn significant. Can I say that: The effect of blue and group2 (β=5, SE=36, p=0.03) and blue and group3 (β=0.02, SE=11, p=0.04) is higher in male than female. Your answers will be helpful for my work.


title: Aggregate ranking from multiple ranked lists

text: I have multiples lists which are ordered. How do I generate a final list with all elements in all lists ranked by their cumulative rank? I.e. the `Final` column here: ![Global ranking of multiple ranked lists](http://i.stack.imgur.com/3eMdD.png) I can't obtain the data with the fruits as a column and other columns with the fruit's numeric ranks. I.e., I can't do: $\begin{array}{c|cc} Fruits&Alice&Bob&...\\ \hline Apple&1&1&...\\ Orange&2&0&...\\ ...&...&...&... \end{array}$


title: Interpreting three forms of a "mixed model"

text: There's a distinction that's tripping me up with mixed models, and I'm wondering if I could get some clarity on it. Let's assume you've got a mixed model of count data. There's a variable you know you want as a fixed effect (A) and another variable for time (T), grouped by say a "Site" variable. As I understand it: `glmer(counts ~ A + T, data=data, family="Poisson")` is a fixed effects model. `glmer(counts ~ (A + T | Site), data=data, family="Poisson")` is a random effect model. My question is when you have something like: `glmer(counts ~ A + T + (T | Site), data=data, family="Poisson")` what is T? Is it a random effect? A fixed effect? What's actually being accomplished by putting T in both places? When should something _only_ appear in the random effects section of the model formula?


title: Standard errors of estimates vary in PROC REG and PROC GENMOD!

text: I am trying to match the outputs of PROC REG with PROC GENMOD. I ran a sample test on the 'iris' dataset of R. The data set is as follows (150 rows in total):                 Sepal.Length Sepal.Width Petal.Length Petal.Width Species     1          5.1         3.5          1.4         0.2  setosa     2          4.9         3.0          1.4         0.2  setosa     3          4.7         3.2          1.3         0.2  setosa     4          4.6         3.1          1.5         0.2  setosa     5          5.0         3.6          1.4         0.2  setosa     6          5.4         3.9          1.7         0.4  setosa      My PROC REG code is:               proc reg data=iris;         model Sepal_Length = Sepal_Width Petal_Length Petal_Width;     run;      My PROC GENMOD code is:               proc genmod data=iris;         model Sepal_Length = Sepal_Width Petal_Length Petal_Width / dist=normal;     run;      The output for PROC REG is:                                       Parameter       Standard     Variable        DF       Estimate          Error    t Value    Pr > |t|          Intercept        1        1.85600        0.25078       7.40      <.0001     Sepal_Width      1        0.65084        0.06665       9.77      <.0001     Petal_Length     1        0.70913        0.05672      12.50      <.0001     Petal_Width      1       -0.55648        0.12755      -4.36      <.0001      The output for PROC GENMOD is:                                              Standard     Wald 95% Confidence      Chi-     Parameter      DF   Estimate      Error           Limits           Square   Pr > ChiSq          Intercept       1     1.8560     0.2474     1.3711       2.3409     56.28       <.0001     Sepal_Width     1     0.6508     0.0658     0.5220       0.7797     97.98       <.0001     Petal_Length    1     0.7091     0.0560     0.5995       0.8188    160.59       <.0001     Petal_Width     1    -0.5565     0.1258    -0.8031      -0.3098     19.56       <.0001     Scale           1     0.3103     0.0179     0.2771       0.3475      According to my understanding, the standard errors of both codes should match as the generalized linear model is run on a normal distribution. Also, I ran both regression on R using `lm()` and `glm(..., family=gaussian)` and the standard error came out equal. Moreover, they are the same as the standard error of PROC REG. Can anyone elaborate on why they are not matching?


What is the difference between discrete data and continuous data?


title: Should sampling happen with or without replacement

text: When doing a simple random sample to estimate population mean for some statistic, how do I know whether sampling happens with or without replacement? It feels wrong to use replacement, because 1) my AP stats teacher never does that and 2) I might use someone's data _twice_ in the average. But on the other hand, the proof that the statistic is an unbiased estimator of the mean is $$E(X)=E(X_1)+\cdots +E(X_n)=\mu+\cdots+\mu$$ which implies $$E\left(\frac{X}{n}\right)=\frac{n\mu}{n}=\mu$$ But doesn't this assume that the $n$ statistics $X_i$ are independent of each other? And isn't that only true if we replace after each sample?


title: Kernel smoothing for Edgeworth expansion

text: Suppose I have an estimator which includes an indicator function in the objective function, then the objective function is not smooth. But if I want to approximate the behavior of this estimator in small samples I would need to have a certain number of higher order derivatives such that I can apply the usual approximation methods like Edgeworth, van Mises or Taylor expansions. I'm particularly interested in the Edgeworth expansion. Would it be possible to smooth out somehow the objective function with a kernel such that I can apply the Edgeworth expansion? I thought of replacing the indicator with a symmetric bound kernel but I'm not sure of whether this would work. Any help would be appreciated!


title: Binary Class Distribution Effects on Probability Scores - (gbm) Boosted Tree Regression Models

text: Any help would be greatly appreciated. **Problem:** I need help to better understand the probability scores that come from the result of a decision tree model. Specifically, I'm using the gbm package from R to create Generalized Boosted Regression Models, but the results I see are common across various ensemble classification models. I am building models with training data that has various levels of unbalanced class levels. In all cases, I'm dealing with two classes, (e.g. Yes/No) with my data. I typically holdout 10% or more of the training data to score and test the model validity. **Findings:** The closer the training data is to a 50/50 distribution (Example: 800 Yes Records and 800 No records), the probability scores assigned to my test data are higher. In fact, the range of scores typically is between 0.01 to .99. This is what I would expect and what I would like for every model. However, as the distribution level becomes more unbalanced (lower amount of "Yes" records) to 40%/60% or even to 10%/90%, the probability scores have a lower range and the maximum score is sometimes below 0.5. As a general rule, if a scored record has a score of 0.5 or greater, then they are predicted to be in the "Yes" class. However, this rule becomes irrelevant if all scored records are below the 0.5 score. The point of building and scoring an unknown universe is to ultimately, RANK the records from most likely to least likely to belong to the "Yes" class. However, even in the cases where all the scores are less than 0.5, the **model does a good job of ranking the records**. Is the probability score that is generated irrelevant and I should re-scale? **I have two main issues:** **Issue #1** - Are these probability scores meaningful as true probabilities? That is, if a model that is built results in zero test records with scores greater than 0.5, is the model poor? Is it common practice to just re-scale the scores from every model to be from 0-1? **Issue #2** - If I re-distribute my classes to be 50/50 in my training file, I sometimes get more "Yes" records than should be expected. For example, a mailing campaign expects 1-5% responses. When building a model using equal classes (50/50), my scored universe returns 20-25%, sometimes even more expected respondents. When is it appropriate to balance the training class and when is it appropriate to leave unchanged and unbalanced? **My Test Results from 3 example Models built from same training data set. The only difference is the class distribution, noted in parenthesis.** Model 1: (50/50)   * No Records = 1,200   * Yes Records = 1,200   * Max Scored record = 0.95   * Median Scored record = 0.48 Model 2: (80/20)   * No Records = 1,200   * Yes Records = 300   * Max Scored record = 0.89   * Median Scored record = 0.16 Model 3: (5/95)   * No Records = 1,200   * Yes Records = 64   * Max Scored record = 0.43   * Median Scored record = 0.039


title: Understanding the hierarchical Dirichlet Process mixture function in DPpackage in R

text: I am going through the `DPpackage` in `R` which provides functions for `non- parametric Bayesian` modeling. I am particularly looking at the function `HDPMdensity()` which relates to `hierarchical Dirichlet Process mixture`. But I am having difficulty in understanding the function and parameter definition provided in the package documentation. The function signature is shown as:               HDPMdensity(y,study,ngrid=100,prior,mcmc,state,status,data=sys.frame(sys.parent()),na.action=na.fail,work.dir=NULL)       Now my major issue in understanding is:   1. What and how this prior is calculated/decided? In one of the example shown in the package documentation it is shown as:              data(calgb)     attach(calgb)     y <- cbind(Z1,Z2,Z3,T1,T2,B0,B1)          prior <- list(pe1=0.1,     pe0=0.1,     ae=1,     be=1,     a0=rep(1,3),     b0=rep(1,3),     nu=9,     tinv=0.25*var(y),     m0=apply(y,2,mean),     S0=var(y),     nub=9,     tbinv=var(y))      The documentation does not specify how the above prior was calculated and thus I am having difficulty in understanding it.   1. What is the `study` parameter? The package documentation says that it is `a (1 by nrec) vector of study indicators` but no where in the examples have they specified the value of this parameter and yet it is passed in the examples.   2. for `y` it is written as `a matrix of responses`. So is this the output responses that you specify for each input while forming the model?   3. If `y` is the response than where is input data specified? I can see that there is one parameter called `data` in function signature but in the examples in documentation this parameter is never used and hence I am unsure if this is what the purpose of the `data` is. **NOTE 1:** The function signature and example shown above are taken from the package documentation for `DPpackage` and are not created by me. **NOTE 2:** I am new to statistics and specially to non-parametric Bayesian modeling and hence request to neglect any lack of understanding that I may have. I am still learning.


title: Cross-validation of a Cox model

text: I have a medical background, with absolutely no knowledge about programming. I am working on survival analysis. I have fitted a multivariate Cox model which I need to cross-validate. My questions are:   1. My dataset is small (79 samples & 37 events). Which is the better method of cross-validation - bootstrapping or leave-one-out method?   2. I am using the `rms` package. I need someone to explain to me in plain English how to enter the code. I have tried reading the help on validation, calibration, but I have not understood much. I tried validation by various methods.  I have indicated below the error message I get. Can someone please explain to me what is going wrong and also give me the correct code to be used for: bootstrapping and cross-validation by leave-one-out method. I have a cox model with four variables: `v1`, `v2`, `v3`, `v4`. `v1` is a factor with three levels; `v2`, `v3` and `v4` are factors with two levels.               data=Model5     fitmodel5 <- cph(Surv(time,event)~v1+v2+v3+v4,na.action=na.omit,data=Model5,x=TRUE,y=TRUE,surv=TRUE)           val1 <- validate(fitmodel5, method="crossvalidation", B = 10, dxy=TRUE)**     Error in measure(xb[test], y[test, , drop = FALSE], strata = stra[test],  :        fit failure in discrim,coxphFit     In addition: Warning message:     In fitter(..., strata = strata, rownames = rownames, offset = offset,  :     Ran out of iterations and did not converge           val1 <- validate(fitmodel5, method=".6", B = 100, dxy=TRUE)     Weights for .632 method (ordinary bootstrap weights 0.01)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      0.00768 0.00925 0.01010 0.01000 0.01070 0.01290      Divergence or singularity in 3 samples     Warning messages:     1: In fitter(..., strata = strata, rownames = rownames, offset = offset,  :       Ran out of iterations and did not converge     2: In fitter(..., strata = strata, rownames = rownames, offset = offset,  :       Ran out of iterations and did not converge           val1 <- validate(fitmodel5, method="boot", B = 100, dxy=TRUE)     Divergence or singularity in 1 samples     Warning message:     In fitter(..., strata = strata, rownames = rownames, offset = offset,  :     Ran out of iterations and did not converge      Thanks