title: Which tests can I use to analyze dependent likert-type data?

text: In a survey, where students were presented various stimuli, after each stimulus they were asked to rate the following:   * **Quality** , on a five-point scale from "Excellent" to "Bad", let's say $Q$   * **Confidence** in their quality ratings, on a four-point scale from "Very Confident" to "Not confident", i.e. $C$ Here are my assumptions:   * Obviously these samples are _dependent_ , since they belong to exactly one stimulus and student pair (am I right assuming so?)   * Also, the questions are very similar to the Likert scale, therefore I'm unsure of whether to classify them as _ordinal_ or _interval-based_ (and others seem too). I'd say they're ordinal, but many researchers in my field seem to ignore that and treat them as interval-based. Basically, what I'd like to find out is whether these ratings are dependent on each other. So, my questions are:   * Would that data be considered interval-based or ordinal?   * Which tests can I apply here? $\chi^2$? Wilcoxon signed rank test? I've already come up with three dimensional plots that show the counts for each pair, e.g. to say "In ten cases, users chose _Bad_ and were _Not Confident_ about it". Or, $count(Q_1,C_4) = 10$ … But there's nothing I can statistically prove from that alone.


title: Linearity Assumption in OLS with Dummy Variables

text: Dumb question but let's say that I have a continuous response variable and have constructed a regression model with multiple predictors. Most of my predictors are continuous but I have one which is a dummy variable. How do I determine whether the ols assumption of linearity is met or violated? And if the linearity assumption is violated, how would I go about transforming the dummy variable?


title: Odds ratios question

text: If you have a categorical variable with $k$ levels (so $k-1$ indicator variables), is there an easy way to calculate the odds ratio? Suppose the first level is the reference level. The model would be $$ \frac{p(x_{2},\dots, x_{k-1})}{1-p(x_{2}, \dots, x_{k-1})} = \text{exp}(b_0+b_{2}x_2 + \dots + b_{k-1}x_{k-1})$$ So if we want to compare the odds between two groups, just plug into the RHS and divide? For example, comparing the odds between $x_2$ and the reference: $\exp(b_0+b_2)$ and $\exp(b_0)$ are the odds? So then just divide them to get the odds ratio?


title: Choosing the proper statistical approach for glm

text: I have data on `alcohol use by drivers` as a binary response variable along with many other variables such as:   1. Driver's age, gender and race    2. Passenger's (if present) age, gender   3. Location variables   4. Road characteristics Here is sample of the data:               data1 <- structure(list(drunk = c(1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0,      1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0,      1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1,      0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0,      0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,      0, 0, 1, 1), DriverAge = structure(c(1L, 3L, 2L, 1L, 2L, 2L,      1L, 3L, 3L, 1L, 1L, 2L, 1L, 3L, 1L, 3L, 1L, 1L, 2L, 3L, 2L, 2L,      3L, 2L, 3L, 2L, 2L, 1L, 3L, 2L, 1L, 3L, 2L, 3L, 2L, 3L, 1L, 2L,      2L, 1L, 1L, 2L, 2L, 2L, 1L, 3L, 1L, 3L, 1L, 1L, 2L, 3L, 2L, 1L,      2L, 3L, 3L, 2L, 2L, 2L, 1L, 3L, 2L, 1L, 3L, 1L, 1L, 2L, 1L, 3L,      3L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 2L, 3L, 1L, 3L, 2L, 2L, 3L, 3L,      2L, 2L, 3L, 3L, 3L, 1L, 3L, 2L, 3L, 1L, 2L, 2L, 2L, 2L), .Label = c("MidAge",      "Old", "Young"), class = "factor"), DGender = structure(c(1L,      2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L,      1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L,      1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,      1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L,      2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L,      1L, 1L, 2L), .Label = c("F", "M"), class = "factor"), location = structure(c(1L,      1L, 2L, 3L, 2L, 2L, 3L, 2L, 3L, 2L, 3L, 3L, 3L, 2L, 2L, 2L, 3L,      3L, 2L, 1L, 2L, 3L, 3L, 2L, 2L, 1L, 3L, 2L, 2L, 3L, 3L, 3L, 2L,      3L, 2L, 3L, 3L, 1L, 2L, 3L, 1L, 3L, 2L, 2L, 3L, 3L, 2L, 2L, 2L,      2L, 3L, 3L, 3L, 2L, 1L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 2L,      3L, 2L, 1L, 1L, 2L, 3L, 3L, 1L, 3L, 3L, 2L, 2L, 3L, 2L, 3L, 3L,      3L, 2L, 2L, 2L, 3L, 1L, 3L, 1L, 3L, 2L, 2L, 3L, 2L, 3L, 2L, 2L,      2L, 3L, 2L), .Label = c("other", "rural", "urban"), class = "factor"),          PassAge = structure(c(4L, 3L, 1L, 4L, 1L, 1L, 3L, 3L, 4L,          2L, 2L, 2L, 3L, 2L, 3L, 3L, 4L, 1L, 2L, 3L, 3L, 3L, 1L, 3L,          1L, 3L, 2L, 2L, 4L, 1L, 3L, 3L, 3L, 2L, 4L, 2L, 1L, 3L, 2L,          2L, 4L, 3L, 1L, 4L, 4L, 4L, 2L, 3L, 1L, 1L, 3L, 4L, 4L, 2L,          4L, 1L, 1L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 4L, 1L, 4L, 4L, 2L,          2L, 1L, 1L, 1L, 1L, 3L, 3L, 4L, 4L, 3L, 3L, 2L, 2L, 4L, 2L,          3L, 2L, 3L, 4L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 4L, 3L, 3L, 2L,          2L), .Label = c("MidAge", "None", "Old", "Young"), class = "factor"),          PGender = c("M", "M", "M", "M", "M", "M", "M", "M", "M",          "None", "None", "None", "M", "None", "M", "M", "M", "M",          "None", "M", "M", "M", "M", "M", "M", "M", "None", "None",          "M", "M", "M", "M", "M", "None", "M", "None", "M", "M", "None",          "None", "M", "M", "M", "M", "M", "M", "None", "M", "M", "M",          "M", "M", "M", "None", "M", "M", "M", "None", "M", "M", "M",          "M", "M", "M", "M", "M", "M", "M", "None", "None", "M", "M",          "M", "M", "M", "M", "M", "M", "M", "M", "None", "None", "M",          "None", "M", "None", "M", "M", "M", "M", "M", "M", "M", "M",          "M", "M", "M", "M", "None", "None")), .Names = c("drunk",      "DriverAge", "DGender", "location", "PassAge", "PGender"), row.names = c(NA,      -100L), class = "data.frame")        1. Now, how do I go about analyzing the data. Will my method depend on what I want to model / analyze?   2. My question is more general than this particular data-set. Should we always take all the available data together. For example, suppose I want to see how passenger's age and gender affect driver's being drunk. One thing I can do is do an ANOVA with respect to passenger's gender. One can also do a binary logistic regression with all data and study the marginal effects/odds-ratios for passenger's gender. How about spiting the data into two sets based on passenger's gender and then getting two separate models and comparing the co-efficient? To me it seems wrong since i have not seen anything like that. Can anyone explain to me why is it so?    3. You may also give a hint on how you would go about analyzing this data from scratch.


title: Is the theory of minimum variance unbiased estimation overemphasized in graduate school?

text: Recently I was very embarrassed when I gave an off the cuff answer about minimum variance unbiased estimates for parameters of a uniform distribution that was completely wrong. Fortunately I was immediately corrected by cardinal and Henry with Henry providing the correct answers for the OP. This got me thinking though. I learned the theory of best unbiased estimators in my graduate math stat class at Stanford some 37 years ago. I have recollections of the Rao-Blackwell theorem, the Cramer - Rao lower bound and the Lehmann-Scheffe Theorem. But as an applied statistician I don't think very much about UMVUEs in my daily life whereas maximum likelihood estimation comes up a lot. Why is that? Do we overemphasize the UMVUE theory too much in graduate school? I think so. First of all unbiasedness is not a crucial property. Many perfectly good MLEs are biased. Stein shrinkage estimators are biased but dominate the unbiased MLE in terms of mean square error loss. It is a very beautiful theory (UMVUE estimation), but very incomplete and I think not very useful. What do others think?


title: Confidence Interval for ratio (response time 2 / response time 1)

text: We have a population of n=1000, nine binary responses (0 symptom absent, 1 present), and 2 timepoints. Note that some of the responses are heavily skewed, only 1% of the subjects are in the "symptom present" category at baseline for the most skewed response. We calculate the change over time as ratio: participants with symptoms present at time2 / participants with symptoms present at time1. This leads to increases between 100% and 600%. Now we want to examine the confidence intervals around these changes. How would you go about this?


title: Is it plausible to get a Cronbach's alpha of .85 with only two Likert-type items?

text: I was reading an article that reported a Cronbach's alpha of .85 for a two- item scale. The scale was measuring self-perceived competence. Each item was on a 7-point Likert-type scale. Intuitively, that seemed like quite a high Cronbach's alpha for a scale with only two items.   * Is it plausible to get an alpha that high on a two-item psychological scale based on Likert-type items?   * How through equations or simulations can I better understand what degree of consistency would be required between the items to achieve a reliability coefficient of .85 with only two 7-point items?


title: Posterior simulations of the variances with the mcmcsamp function

text: I would like to get the posterior simulations of the variance components of a lmer() model with the mcmcsamp() function. How to do ? For instance, below is the result of a lmer() fitting :               > fit     Linear mixed model fit by REML     Formula: y ~ 1 + (1 | Part) + (1 | Operator) + (1 | Part:Operator)        Data: dat        AIC   BIC logLik deviance REMLdev      97.55 103.6 -43.78    89.18   87.55     Random effects:      Groups        Name        Variance Std.Dev.      Part:Operator (Intercept) 2.25724  1.50241      Part          (Intercept) 3.30398  1.81769      Operator      (Intercept) 0.00000  0.00000      Residual                  0.42305  0.65043     Number of obs: 25, groups: Part:Operator, 15; Part, 5; Operator, 3      Now I run mcmcsamp() :               > mm <- mcmcsamp(fit, n=15000)       I expected that the simulations of the residual variance are stored in the "sigma" node but this does not seem to fit the results of lmer() :               > sigmasims <- mm@sigma[1,-(1:5000)]  # discard first 5000 simulations (burn-in)     > summary(sigmasims)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      0.8647  1.4960  1.7040  1.7460  1.9480  3.7920       Similarly I expected that the simulations of the other variance components are stored in the "ST" node but I get a similar observation.


title: Ordered logit with (too many?) categorical independent variables

text: I am doing some inferential analysis on some ordered dependent variables (about 70). Their scales range from 4-10 possible (ordered) responses. To give some context, it's all social data--happiness, feelings towards minorities, etc.--and my inferential question is on the social returns to university degrees. For my binary and continuous dependent variables (not a part of this problem), I have a nice set of controls. Quite a few of these are categorical (did you go to a private high-school? What is your marriage status?), and others are continuous or dummy. My problem is that when I use the same controls for the ordered variables, I get plenty of outcome level/categorical independent variable combinations which are empty. For example, there are no widowers who give themselves a happiness score of 2. As there are many such cases like this, the R functions polr and lrm don't work (if this is not the case, please let me know). The choice, then, seems (to me) to be either: a) reduce the set of controls to eliminate the empty cells, or b) do it in ols. I'm aware that for my categorical data, the assumptions of ols are not met. But given my task, I'm wondering how much would be lost from either solution?


title: ROC graph interpretation

text: I'm reading Fawcett's 2004 paper on ROC graphs for machine learning algorithms, which can be found here. On page 7-8 he shows a simple ROC example and makes some interpretations that I don't understand. Below is the ROC graph: ![example ROC graph](http://i.stack.imgur.com/vE5u4.png) And here is what he wrote: > Although the test set is very small, we can make some tentative observations > about the classifier. It appears to perform better in the more conservative > region of the graph; the ROC point at (0.1,0.5) produces its highest > accuracy (70%). This is equivalent to saying that the classifier is better > at identifying likely positives than at identifying likely negatives. Note > also that the classifier’s best accuracy occurs at a threshold of ≥ .54, > rather than at ≥ .5 as we might expect with a balanced distribution. I don't understand how he derived his interpretations.   1. `the ROC point at (0.1,0.5) produces its highest accuracy (70%)` How is the highest accuracy of 70% for point (0.1, 0.5) found from that graph, and how do we know it's the highest accuracy?   2. `This is equivalent to saying that the classifier is better at identifying likely positives than at identifying likely negatives.` I don't see how that interpretation is determined.   3. `Note also that the classifier’s best accuracy occurs at a threshold of ≥ .54` How was this found?   4. `rather than at ≥ .5 as we might expect with a balanced distribution` Why would we expect that? Thank you for any help.


title: Famous statistician quotes

text: What is your favorite statistician quote? This is community wiki, so please one quote per answer.


title: What are the good algorithms for feature extraction for large dataset?

text: I have KDD dataset for detecting fraud actions on networks but it has millions of lines and >20 feature columns. Thus it is not viable to process all these on my personal computer. I am thinking about the random sampling on data to reduce the rows to some legible numbers than I want also to reduce the number of features. I know PCA can be used to reduce the feature number but first I also want to get rid of the some dirty features that are not helpful enough to classification purpose. What are the algorithms can be used for that feature elimination and extraction purpose?


title: How can I test whether my clustering of binary data is significant

text: I'm doing shopping cart analyses my dataset is set of transaction vectors, with the items the products being bought. When applying k-means on the transactions, I will always get _some_ result. A random matrix would probably also show some clusters. Is there a way to test whether the clustering I find is a significant one, or that is can be very well be a coincidence. If yes, how can I do it.


title: Expectation of a minimum

text: Assume that $c_{ij}(W)$ is a simple function of a random variable W. (This comes is Dynamic Programming where we are currently at node i and we want to go to node j in the shortest path problem and W is the random variable denoting the time that we wish to minimize) What does the following expression mean? $$\mathbb{E} \{\min_{j} \quad c_{ij}(W)\}$$ Intuitively, I understand that we are minimizing the expectation of the minimum possible of all possible functions of $W$ but is there more?


title: I'm stuck: How to 'bin' multiple binary outcome variables and then analyze within subject?

text: I have what feels like it should be a very easy analysis to work out, but I've been trying for hours and I can't figure it out. I'm using SPSS by the way. I have 10 (5 of each group) different binary response variables (yes/no) that I need to 'bin' into two different groups for comparison. Each subject answered each of these ten questions so I need to analyze this within subjects. However, since there was no independent variable, I want to figure out a way to bin these variables and compare the two bins against each other. I tried creating multiple response sets in SPSS but I can't seem to figure out how to use them after the fact. I also thought maybe I would need to restructure the dataset but again I'm not sure how to optimally do that. Finally, I thought maybe I could create new variables based on the counts for each response within each group and subject (so how many said 'yes' in group 1 variables, how many said yes for group two etc) so ending up with 4 new variables. Or create 2 new variables with proportions of yes to no in each.but I'm not sure how to create proportion variables for each subject in SPSS I feel like there must be an easier way to do this analysis that I'm missing, and that takes into account all of the response variables, while analyzing group differences. Any thoughts? I would really appreciate it!


title: How to determine which method is the most reasonable clustering results?

text: * * * **Method 1: Cluster by K-means with initial centroid {27, 67.5}** * * * **Method 2: Cluster by K-means with initial centroid {22.5, 60}** * * * **Method3: Agglomerative Clustering** * * * How can I know which method gives a more reasonable clustering results?


title: Non significant Pearson correlations included in hierarchical regression?

text: I would like to perform hierarchical regression in which all variables are based on previous research/theory. But when I perform Pearson correlations, I found that some variables did not correlate to the DV. Should I still include those non-significant variables in the hierarchical regression as suggested by previous research even though they were not significant in Pearson correlations? Thank you very much for your help.


title: Interpreting the change of significance when adding an independent dummy variable

text: I have a regression with my dependent variable being the number of visitors and my independent variable being the presence or absence of a certain kind of ad. Now, when I run the regression my independent variable is significant. When instead, I add another independent variable (the presence of a big exhibition at the museum) the same variable becomes non significant. Can someone enlighten me on what to make of this result? Thank you very much!


title: Understanding Dempster et al. on expectation maximization

text: I'm reading about expectation maximization from Dempster, Laird and Rubin's original paper which can be found from the following link: http://web.mit.edu/6.435/www/Dempster77.pdf My questions are from pages 1 and 2. I have added those particular parts into the image below and tried to be as clear as possible with my questions (4 questions in the image, highlighted with red numbers 1, 2, 3, 4). Hope someone could make things clearer for me :) To summarize my questions in the picture are: 1) Have I understood the relationship between the variables correctly? 2) What does the text in the green box mean exactly? 3) How to interpret these formulas 4) WHAT IS GOING ON HERE?!...this doesn't seem clear at all! x) ![enter image description here](http://i.stack.imgur.com/aqX8Q.png)


title: Regression to predict probability - what transformation to use?

text: I am working on a project where the training labels are given to me as a probability value in the range [0,1]. My first approach was to fit a simple linear ridge regression to predict the probability. This isn't ideal as: 1) Predictions from this model end up with values outside the range of 0 and 1 2) I don't think linear regression works that well if you want to make it predict something in a fixed interval. I try the logit transformation, but since I have probabilities that are exactly 0 and 1, I perturb them slightly so I get 0.00001 and .999999 instead to avoid +/- infinity. I train my model on the transformed labels, and then make a prediction on the test set and undo the transformation with the inverse of the logit function (logistic function). Frustratingly though, this gives me even worse results than the naive linear regression! Any suggestions on other transformations I can try or what I am doing wrong? Thank you.


title: Confusion related to derivation of the log likelihood

text: I was reading the Hastie, Friedman, Tibshirani paper "Sparse Inverse Covariance Estimation with the Graphical Lasso" and it had the following ![enter image description here](http://i.stack.imgur.com/915iR.png) I couldn't get how the following expression was derived $\hat\Sigma^{-1} = \displaystyle\arg \max_{X \succ 0} \log \mbox{det}X - \mbox{trace}(SX) - \lambda||X||_1$ I know that they are calculating the likelihood of the data and maximizing with respect to the mean $\mu$ as given in the paper and it results in the expression above. How is the expression obtained?


title: Multivariate distance function in Excel

text: I want to do a simple nearest neighbour calculation in Excel over a multivariate space to get an idea of how my data clusters. I have a set of data points $\{X_1, X_2, \ldots X_n\}$ and a set of possible cluster centers $\{M_1, M_2, \ldots, M_m\}$. Initially, I propose to simply use a weight vector $w$, and compute the scalar product of $w$ and $[X_i - M_j]$ to get the distance of point $i$ to cluster centre $j$. This is all very simple... but completely unwieldy if $n$ and/or $m$ and/or the dimensionality of my data points is large. So is there any (simple) way of hiding all the matrix calculations and producing a (matrix/vector) formula that produces the (scalar) distance without using temporary arrays?


title: Matching covariates between more than 2 groups

text: I have three groups that I would like to match on age and gender. They don't need to be matched one-to-one, I just need group differences p < .05 for both variables. I can find packages that will match TWO groups, but not more than that. Does such a package or code snippet exist? Secondarily, because the available matching packages only work for two groups, the treatment parameter is required as a logical vector, whereas right now mine are factors. What is the most efficient way to convert a factor to a logical when there are only 2 levels? Muchos gracias!


title: Basic Bayesian MCMC to estimate two parameters from binomial distributions given unknown number of trials

text: This is a very basic question about Bayesian inference. I'm not grasping one or more fundamental concepts. Let's say I have two observed outcomes, _X_ and _Y_. I want to infer the probabilities ( _px_ and _py_ , respectively) of each occurring given _X_ and _Y_. I do not know _N_ , the total number of trials. I'm assuming that _X_ and _Y_ are binomially distributed. How do I calculate the likelihood without _N_? **What I ultimately want is to show the bivariate posterior distribution of _px_ and _py_ via MCMC. I do not care about estimating _N_ --I just want to show the chain in the plane of _px_ and _py_. No convergence necessary.** * * * _Clarification:_ _X_ and _Y_ are drawn from the same _N_ : _X_ ~Binom( _N_ , _px_ ) and _Y_ ~Binom( _N_ , _py_ ). We have no other information about _px_ or _py_ , although I'll use a beta prior to start. Also assume that _X_ and _Y_ are independent.


title: Cluster evaluation based on different similarity algorithms

text: At the moment I'm researching clustering of single words. The input of this research is a list of words (unigrams). During the research I want to compare several similarity algorithms to see how they impact the cluster result. At the moment I’m using a 70word data set, where the words belong to 7 clusters (i.e., animals, music instruments, vehicles etc). How do I evaluate the different algorithms? I can use measures like: Rand measure, F-measure or Jaccard index. But I’m struggling with two questions. The first one is the cluster count and the second how do you measure accuracy? In order to compare you need the same measure, otherwise what you are comparing to each other doesn’t say anything. Thus, I need to keep the cluster count equal, or can I compare while using different cluster counts? The silhouette function in Matlab shows me that the optimal cluster count differ between the similarity algorithms, when I look at the data I see a lot of errors and this error count will go down when I choose more clusters (this is no surprise). Second is the accuracy measure. How do I determine if the word is in the right cluster? Because it is more a data exploration tool. For example with the data set from above and I choose to cluster the words into seven clusters I got the following result: six vehicles and one animal. Is my precision and recall then .86? Or do I need to calculate this over all clusters? If so, how do I cope with clusters that have six vehicles and six animals? What is the right cluster in order to determine if the right word is inside? Or doesn’t it matter because the math will stay the same?


title: Computing Image Similarity based on Color Distribution

text: ## Image Similarity based on Color Palette Distribution I am trying to compute similarity between two images based on their color palette distribution, let's say I have two sets of key value pairs as follows, Img1: `{'Brown': 14, 'White': 13, 'Black': 40, 'Gray': 31}` Img2: `{'Pink': 82, 'Brown': 8, 'White': 7}` Where the numbers denote the % of that color present in the image. What would be the best way to compute similarity on a scale of 0-100 between the two images?


title: IRR on two raters’ integer scores

text: I am writing with questions about a study I am helping with: Two “raters” will score athletes subjectively from 0 to 3 (3>2>1>0) on each of 10 skill performances. Scores will be summed to a maximum of 30, and the total bracketed into approximately three categories for comparison with regard to knee injury incidence. Reliability/agreement between the two raters’ scores must be determined. How should the inter-rater reliability testing be set up? Kappa and/or Kendall? Same data for individual skills and total score? How much data is needed for this irr testing? The lead investigator anticipates highly consistent scores between raters and plans to divide their testing of subjects accordingly, for the main study.


title: Change-point in Cox survival model

text: I am looking to fit a Cox proportional hazard survival model. Looking at the K-M curve (below) for one variable (with 2 categories) it appears there is a change in hazard ratios at around day 110. I was thinking of modeling this with a change-point model. ![KM-Curve](http://i.stack.imgur.com/xZBjJ.png) I'm having trouble implementing it. I have defined days_ind as 1 if days>=110 and 0 otherwise. Then I run the model:               coxph(Surv(time=days,event=event2)~x*days_ind)      I get several warning messages about convergence and the results don't seem to make any sense. Am I approaching this in the correct way? I thought of bringing in the interaction of x*days instead but this too does not converge and also leads to strange estimates.


title: Estimating Pearson III density using R's lmomco package

text: I am trying to estimate parameters for a Pearson III distribution using package "lmomco" (session information given at the end). I have 40 precipitation values in vector "val.x". I estimate P3 parameters using the sequence of pwm.ub(), pwm2lmom() and parpe3() function calls (results for each are shown below). Although estimated parameters are checked as part of parpe3(), I use is.pe3() and are.parpe3.valid() to check for any errors. Once parameters are estimated, I try to compute Pearson III densities for the val.x vector. I then plot a histogram of values and try to overlay the line of densities to see whether the P3 provides a good fit. PROBLEM: For THIS particular set of x-values, the computed densities are way too high (by maybe two orders of magnitude). As a result, in my diagnostic plot the densities fall much above the histogram. BUT, if I generate a slightly different data set by generating a sample with replacement (see below), then the densities are OK... Any clues of what may be happening? Suggestions/comments are most welcome! Thanks in advance!!! Reproducible example (results shown in commented lines):               # Vector of precipitation values          val.x <- c(71.6, 80.1, 87.7, 91.1, 117.3, 125.5, 128.4, 129.3, 135.9, 137.3, 143.5,       150.3, 162.2, 178.8, 186.3, 189.5, 192.1, 192.2, 210.4, 217.4, 217.6, 218.5,       219.6, 227.0, 235.6, 241.2, 247.3, 251.6, 260.7, 267.1, 267.3, 286.3, 286.8, 288.7,       289.2, 302.2, 312.2, 328.1, 335.4, 338.5)          pppwm <- pwm.ub(val.x)          # $betas     # [1] 208.69500 126.26442  91.32573  71.75488  59.19168     # $source     # [1] "pwm.ub"          lmom <- pwm2lmom(pppwm)          # $lambdas     # [1] 208.6950000  43.8338462  -0.9371457   1.8036514   0.4569509     # $ratios     # [1]  NA  0.21003784 -0.02137950  0.04114746  0.01042461     # $source     # [1] "pwm2lmom"          p3.pars <- parpe3(lmom, checklmom=TRUE)          # $type     # [1] "pe3"     # $para     #         mu       sigma       gamma      # 208.6950000  77.7352993  -0.1312404      # $source     # [1] "parpe3"          # Check the Pearson III parameters          is.pe3(p3.pars)     are.parpe3.valid(p3.pars, nowarn=FALSE)          #[1] TRUE     #[1] TRUE          # Compute Pearson III density values for val.x          den.val.x <- pdfpe3(val.x, p3.pars)          #> den.val.x     # 0.08423841 0.10154640 0.11880345 0.12705149 0.19986800 0.22500346     # 0.23400587 0.23680527 0.25732763 0.26166189 0.28065471 0.30086561 0.33360020     # 0.37050541 0.38272547 0.38696339 0.38995436 0.39006114     # 0.39884633 0.39644872 0.39633320 0.39578137 0.39503603 0.38803349 0.37574868     # 0.36554600 0.35265859 0.34257802 0.31895037 0.30083653 0.30025457 0.24237631     # 0.24081995 0.23490770 0.23335295 0.19351901 0.16440994 0.12262096 0.10568267     # 0.09895299                   # Plot histogram of val.x and overlay fitted Pearson III density              hist(val.x, freq = FALSE)         lines(val.x, den.val.x)              # Code used to generate an alternative data set for which the estimation works...               set.seed(1234)         val.x <- sort(sample(val.x, length(val.x), replace = TRUE))      R session info: > sessionInfo(package = NULL) R version 3.1.1 (2014-07-10) Platform: > x86_64-w64-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252   [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C   [5] LC_TIME=English_United States.1252 attached base packages: [1] parallel grid splines stats graphics grDevices utils datasets methods base other attached packages: [1] PearsonDS_0.97 yaml_2.1.13 xts_0.9-7 zoo_1.7-11 weathermetrics_1.0 [6] stringr_0.6.2 SPEI_1.6 spdep_0.5-74 Matrix_1.1-4 reshape2_1.4   [11] RColorBrewer_1.0-5 plyr_1.8.1 nortest_1.0-2 mice_2.22 Rcpp_0.11.2   [16] mgcv_1.8-1 nlme_3.1-117 maps_2.3-7 MASS_7.3-33 lubridate_1.3.3   [21] locfit_1.5-9.1 lmomco_2.1.1 insol_1.1.1 Hmisc_3.14-4 Formula_1.1-2   [26] lattice_0.20-29 geosphere_1.3-8 futile.logger_1.3.7 fitdistrplus_1.0-2 survival_2.37-7   [31] dplyr_0.2 Cairo_1.5-6 boot_1.3-11 ADGofTest_0.3 sp_1.0-15   [36] RPostgreSQL_0.4 DBI_0.2-7 loaded via a namespace (and not attached): [1] assertthat_0.1 class_7.3-10 cluster_1.15.2 coda_0.16-1   [5] deldir_0.1-5 digest_0.6.4 e1071_1.6-3 futile.options_1.0.0 [9] lambda.r_1.1.6 latticeExtra_0.6-26 LearnBayes_2.15 memoise_0.2.1   [13] nnet_7.3-8 randomForest_4.6-7 rpart_4.1-8 tools_3.1.1 >


title: Conditional Inference Trees in R

text: Consider the following data set **test** with binary outcome variable **z** and predictor variables **a,b,c**.                   z   a   b   c         0   1   40  185         0   1   25  128         0   0   32  100         0   0   29  100         1   1   30  107         0   0   30  133         1   1   38  132         1   1   37  127         1   0   30  184         1   0   40  199         1   1   26  185         0   1   21  185         0   0   21  134         0   0   20  137         1   1   22  135         0   0   23  189         1   0   32  109         1   0   31  152         1   0   38  130         1   1   37  191         0   1   39  168         1   0   28  183         0   1   26  171         1   1   23  164         0   1   32  111         0   0   34  131         1   0   30  121         1   0   27  195         1   1   29  117         1   0   26  187         1   0   34  183         0   0   28  189         0   1   34  150         0   1   34  176         0   1   24  140         1   0   37  181         0   1   36  109         1   0   39  198         0   0   32  164      Suppose I use to caret package to train a conditional inference tree:               library(caret)     set.seed(825)     fitcontrol4 <- trainControl(method="repeatedcv", number=10, repeats=10, classProbs=TRUE, savePred=T)     RFfit5 <- train(z~a+b+c, data=test, method="ctree", trControl=fitcontrol4)      I then get the accuracies for each repetition using ` RFfit5$resample `. Is this method of cross validation good to use for conditional inference trees (e.g. is oob a better method)? Are there any suggestions about how else I could display a conditional inference tree?


title: One random variable bigger/smaller than another independent

text: Suppose $X$ and $Y$ are independent. To be concrete, suppose $X$ is $N(0,a)$ and $Y$ is $N(0,b)$. For any scalar c,d with c < d, is there a way to bound the following probability from above in terms of a,b,c,d? $P(cX<Y<dX)$


title: Can i use rotated PCA factors to make models and then subsitute these back to my original variables?

text: I have been using PCA results and variable loadings on the factors generated to select original uncorrelated input variables for model making I was putting what i thought were **rotated** factors from PCA into models and seeing what factors were significant. Then i would go back and see if i could make sense of the results in regards to the original input variables. For example if rotated PC #4 was highly significant in my model and had a single loaded input variable (say variable 8 loaded onto PC #4 to 0.9) i would substitute PC#4 for variable #8. This did not work for one data set and i found that the factors i was using were generated raw **unrotated** PCA factors. These are obviously much less distinguishable in terms of the original x variables. Im not very familiar with the pure mathematical side of PCA but i had assumed the factors i generated related to the rotated factors not the unrotated factors. The initial data set i worked on appeared to be successful as the unrotated and rotated factors were fairly similar thus i didnt come across the issue Can i use rotated factors after PCA to make a model and then substitute these factors back to original X variables?


title: Structural equation model: Bootstrap, Bayesian statistic or rescaling binary data to deal with small sample size?

text: I’m trying to solve a SEM. Mi model includes 4 latent variables and 11 observed variables. My data is binary. **Question 1:** I just have 47 observations and in the literature a sample size over 150-200 is recommended. I understand that bootstrapping or Bayesian statistic may be used to deal with small sample sizes. Am I right? I’m using the R package `lavaan` that has implemented bootstrap, but as my data is binary I’m not sure if I can use bootstrap. For binary data it is recommended to use diagonally weighted least squares (DWLS). **Question 2:** If I can’t use bootstrap because my data is binary, what if I rescale my data? There are methods to rescale categorical and ordinal data based on correlational aspects of data. Variance covariance and correlation matrices may be calculated with these methods. Is it a reasonable solution? **Question 3** : Regarding a Bayesian approach, is it a good approach?


title: Regression with repeated measures in Matlab

text: Is there a way to perform multiple logistic regression on repeated measures data using Matlab? I have a data set containing a daily measurement recorded from 20 participants for 60 days. I am calculating 18 features for each participant for each day. I wish to determine if there is a relationship between each of these features and a 8 binary response variables, also recorded for each participant each day. What is the best way to do this in Matlab? My question is related to a previous question


title: How to compare multiple groups of unblanced repeated measures non normal data?

text: I'm trying to compare three groups data. But the data set is about a new drug trial. The data set has these characteristics:   1. Follow-up set. That is, after administration of the drug, a series of parameters were collected in the following date, day0, day1,day2,...,day28.   2. Unbalanced set. Because some patients died and some recovered, not all patients were followed to 28 days; for instance, some were followed for only 20 days.    3. **Non normally distributed parameter. Especially, day25-day28, sometimes sd >> mean.** So, I try to use the linear mixed-effect model with `group`, `time` and `group*time` for the comparison. ANOVA is not suit for the unbalanced set. However, I wonder if it's correct model for **non normal data set**? Or is there another method suited for this kind of set? I use SPSS 19. Some of my friends use STATA or JMP.


title: Benjamini-Hochberg dependency assumptions justified?

text: I have a data set where I test for significant differences between three populations with respect to some 50 different variables. I do this using Kruskal-Wallis tests, on one hand, and by likelihood ratio tests of nested GLM model fits (with and without population as an independent variable), on the other. As a result, I have a list of Kruskal-Wallis $p$-values on one hand, and what I think are Chi square $p$-values from the LRT comparisons, on the other. I need to do some form of multiple testing correction since there are >50 different tests, and Benjamini-Hochberg FDR seems like it is the most sensible choice. However, the variables are probably not independent, with several "clans" of them being correlated. The question is then: how can I tell if the set of underlying statistics for my $p$-values meet the requirements of positive dependence that are needed for the Benjamini-Hochberg procedure to still be bound to the FDR? The Benjamini-Hochberg-Yekutieli paper from 2001 states that the PRDS condition holds for multivariate normal and studentized distribution. What about my likelihood ratio test Chi square values for the model comparison? What about the $p$-values I have for the Kruskal-Wallis tests? I can use the Benjamini-Hochberg-Yekutieli worst-case FDR correction that assumes nothing on the dependency, but I think it may be too conservative in this case and miss some relevant signals.


title: Selecting the best indicator of disease progression

text: Quite a while ago, I asked a question for which Peter Ellis provided a very interesting answer. Now I'd like to follow on that and have your comments and ideas on how to actually put it to use. I try to restate the problem, hopefully this time in a more clear way using a synthetic example. **Problem** : Two groups of subjects participate in a longitudinal study. One is a control group and the another is a group of subjects with a degenerative disease. That is to say, from all other sources we know and expect to see a decline in the conditions of the subjects as time goes by. Over the course of the study, we measure all subject using a series of performance measures $M_{1}...M_{n}$. Our research question is: which of these measures can pickup the decline in patient's conditions the best. **Suggested solution** : Peter's insight was to build a _theoretical_ response variable and use model selection methods to see which of $M_{1}...M_{n}$ could predict it the best. Here is a simple simulation with 2 measures to demonstrate the problem the best:               library(ggplot2)     library(nlme)          set.seed(1234)          # Followups in the longitudinal study     t <- 1:5          # Number of subjects in each group and number of followups     N <- 20     M <- length(t)          # Subject IDs     id <- c(paste('C', 1:N, sep=''), paste('P', 1:N, sep=''))          # The first measure. Each subject has a unique value at the baseline.     # There is also a difference between the groups at the baseline.     m1.c.start <- 150 + rnorm(N, 0, 5)      m1.p.start <- 140 + rnorm(N, 0, 5)          # The control group does not change over the time.     # The patient group, however, decline over the time.     m1.c <- rep(m1.c.start, each=M) + rep(0 * (t-1)/M, N) + rnorm(N*M)     m1.p <- rep(m1.p.start, each=M) + rep(-5 * (t-1)/M, N) + rnorm(N*M)          # Generate the data for the second measure.     # The separation of the two groups and rate of progress are weaker.     m2.c.start <- 70 + rnorm(N, 0, 3)      m2.p.start <- 71 + rnorm(N, 0, 3)     m2.c <- rep(m2.c.start, each=M) + rep(0 * (t-1)/M, N) + rnorm(N*M)     m2.p <- rep(m2.p.start, each=M) + rep(1 * (t-1)/M, N) + rnorm(N*M)          # Now put everything in a dataframe     data <- data.frame(       m1 = c(m1.c, m1.p),        m2 = c(m2.c, m2.p),       group=c(rep('C', M*N), rep('P', M*N)),        id = rep(id, each=M),       t=rep(t, N*2))      Here are the plots of measures over the time: ![enter image description here](http://i.stack.imgur.com/bFsxv.png) ![enter image description here](http://i.stack.imgur.com/0abJA.png) Now to build a theoretical response: In the case, we'd like to see a constant response over the time for the control group and a change of performance for the patients. Also, at the baseline groups are typically already a bit different.               # The desired, theoretical response     data$theory <- 2 + (data$t-1)/M     data$theory[which(data$group=='C')] <- 1      Here is the plot of the theoretical response: ![enter image description here](http://i.stack.imgur.com/2nFTS.png) Notice that the exact value of the slope and the differences between the two groups in this theoretical curve should not be important as long as model we fit could consider them. Finally, we fit a linear mixed model for each measure as well as a null model:               fit.0 <- lme(theory ~ 1, random=~1|id, data=data, method='ML')     fit.1 <- lme(theory ~ m1, random=~1|id, data=data, method='ML')     fit.2 <- lme(theory ~ m2, random=~1|id, data=data, method='ML')      I fit the models with _ML_ option as I'm interested in comparing the fixed effects. The final step is to compare the models to pickup the best measure of the disease progress. I use $AIC_{c}$ for this purpose. The _AICcmodavg_ has a convenient function to build an $AIC$ table for us:               library(AICcmodavg)          f <- list(fit.0, fit.1, fit.2)     aictab(f, modnames=c('Null', 'm1', 'm2'))          Model selection based on AICc :               K   AICc Delta_AICc AICcWt Cum.Wt     LL     m1   4 -45.15       0.00      1      1  26.68     m2   4 119.51     164.66      0      1 -55.65     Null 3 130.22     175.37      0      1 -62.05      The results suggest that $M_{1}$ is the best predictor of the theoretical response here. Now what I'd like to know is:   * Do you think the whole idea of this approach to pickup the best measure related to disease progression make sense? Have you seen something similar? Is this the right way to do it?   * Do you think the theoretical model I built is sound? Are there cases that this model would mislead us?   * Can we expand the same kind of approach for picking up the best measures in other types of designs (like e.g. before/after interventions, comparing two different interventions)?   * Any other comments or suggestions? **EDIT** : I changed _bio-marker_ to _measure_ in the text as I suspected that saying what we deal with are bio-markers draws attention to issues that are not pertinent to this narrow definition of the problem. **EDIT** : With this approach looking at combination of measures is easy. For example:               fit.3 <- lme(theory ~ m1 + m2, random=~1|id, data=data, method='ML')      In this example, the $AIC_{C}$ of fit.3 will be slightly worse than the model with only $M_{1}$ (i.e. fit.1).


title: How to summarize categorical data?

text: I've been struggling with the following problem with hopefully is an easy one for statisticians (I'm a programmer with some exposure to statistics). I need to summarize the responses to a survey (for management). The survey has 100+ questions, grouped in different areas (with about 5 to 10 questions per area). All answers are categorical (on an ordinal scale, they are like "not at all", "rarely" ... "daily or more frequently"). Management would like to get a summary for each area and this is my problem: **how to aggregate categorical answers within the related question?**. The questions are too many to make a graph or even a lattice plot for each area. I favor a visual approach if possible, compared to, say, tables with numbers (alas, they won't read them). The only thing I can come up with is to count the number of answers in each area, then plot the histogram. Is there any thing else available for categorical data? I use R, but not sure if it's relevant, I feel this is more of a general statistics question. Thanks for any suggestions.


title: SAS PROC LOGISTIC: Hosmer and Lemeshow test is good but Gini is bad?

text: I am using `PROC LOGISTIC` along with `Class` statements to do binary logit model(default=1,non-default=0) on a bank loan dataset where I have 7 numeric independent variables and 1 categorical independent variable. Hosmer and Lemeshow test yields a `p = 0.6757`, and `Percent Concordant = 73.3` so the model is fine I think. To my dismay Gini is coming up as `0.2305`, and ks is only `35.1`. I there any way I can do anything about the results or am I missing anything?               %let Keepvar= Age_in_years Level_of_education Years_with_Current_Employer         Years_at_Current_Address Household_income_in_thousands Debt_to_income_ratio         Credit_card_debt__in_thousands Other_debt_in__in_thousands;      PROC LOGISTIC DATA = work.bankloan descending outest=estimates;         class Level_of_education (ref='1') / param=ref;         MODEL Defaulted = &keepvar/selection=backword rsquare lackfit;     RUN;          %macro rankorder(data=, est=, dependent= );         proc score data = &data score = &est out = mlscore type = parms;         var &keepvar;     run;          data mlscore;         set mlscore;         prob = exp(&dependent.2) /(1+exp(&dependent.2));     run;          proc rank data = mlscore out = mlrankscore groups =10;         var prob;         ranks probrank;     run;          proc sql;     select probrank, avg(prob) as probability, count(*) as total,avg(&dependent) as        Act_Avg_Response, sum(&dependent) as sum_of_response from mlrankscore group by        probrank order by probrank descending;     quit;     %mend;     %rankorder(data=bankloan, est=estimates, dependent=Defaulted);


title: Stats 101 basic hypothesis testing questions

text: > The annual profit for your organisation last year was £1203, per employee. > The average for the sector was £1228 per employee. You need to know if there > is evidence that your organisation is not as profitable, on average, as your > competitors. The standard deviation for the sector is £104 per employee. **Questions** :   1. State $H_0$ and $H_1$ for the test to check  Assuming that annual per employee profit follows a Normal distribution, the resulting P-value for the test is, P= 0.405 (3 d.p.):   2. Based entirely on this P-value, what would you conclude from the test?    3. What is the probability of getting a profit per employee of £1203, or less, when there is no difference between that for your company and the industry average?    4. If your annual per employee profit was such that the test had resulted in a P-value given by P=0.032 and you were testing at the α=0.01 significance level, what would you conclude and why?    5. If you chose to reject $H_0$ but $H_0$ was actually true, would you be making a Type I or Type II error?    6. What assumptions have you made about the annual per employee profit in this test?


title: Fitting a dataset into another dataset

text: I had heard that there is a body of literature devoted to the following problem: You have a dataset and you produced a good predictive model for it. Now you have a different dataset, derived from different instruments, or different data sources, but similar enough that you can hope to scale one to the other so that you can use the predictive model on it. Alternately, you do not have enough data to produce a model and hope that combining different datasets by scaling one to the other will help you reach the quantity of data you need. First off what is this problem called so that I can search it more effectively? And also does anyone know of a good recent survey of techniques for this, either in book or paper form? I currently have access to most academic journals so those links work for me as well. Say I have dataset A that I have a model for, dataset B occupies the same database schema but is from a difference source with different factors that are not included in my feature set. Initially my intuition was to construct a QQplot and fit lines (or curves) to features that I thought should be similar. If the difference in the way feature 1 from A increases is similar to the way feature 1 from B increases but with a constant factor then fitting a line can reveal this factor. If the difference was exponential or logarithmic then I could scale using a fitted function. In this way I could constrain the way one variable increased to fit how another variable increased. However this is just my intuition. I can certainly test it for overfitting but when I had heard that there was a lot of literature devoted to this subject then it seemed as though I should learn a few ways in which I could question my assumptions. It would probably be good for me to review the literature. Does anyone know what the tag for this literature might be?


title: Effect on a rate

text: I've got one independent variable and two dependent variables. As a concrete example, consider the independent variable to be _browser_ (categorical, what Web browser a visitor to a Web site uses) and the dependent variables to be _visitors_ (number of visitors) and _enrollers_ (number of visitors who sign up for a newsletter). Note that _enrollers_ depends on _visitors_ (because every enroller is also a visitor). I care about the rate _enrollers/visitors_ only. How can I determine the effect size of _browser_ on the rate _enrollers/visitors_? My specific concern is that I don't know to study _enrollers/visitors_ where power and presumably effect size depend (also) on _visitors_ itself.


title: Simulation of Logistic Regression Power Analysis - Designed Experiments

text: This question is in response to an answer given by @Greg Snow in regards to a question I asked concerning power analysis with logistic regression and SAS Proc GLMPOWER. If I am designing an experiment and will analze the results in a factorial logistic regression, how can I use simulation ( and here ) to conduct a power analysis? Here is a simple example where there are two variables, the first takes on three possible values {0.03, 0.06, 0.09} and the second is a dummy indicator {0,1}. For each we estimate the response rate for each combination (# of responders / number of people marketed to). Further, we wish to have 3 times as many of the first combination of factors as the others (which can be considered equal) because this first combination is our tried and true version. This is a setup like given in the SAS course mentioned in the linked question. ![enter image description here](http://i.stack.imgur.com/4LIvh.jpg) The model that will be used to analyze the results will be a logistic regression, with main effects and interaction (response is 0 or 1).               mod<-glm(response ~ Var1 + Var2 + I(Var1*Var2))      How can I simulate a data set to use with this model to conduct a power analysis? When I run this through SAS Proc GLMPOWER (using STDDEV =0.05486016 which corresponds to `sqrt(p(1-p))` where p is the weighted average of the shown response rates):               data exemplar;     input Var1 $ Var2 $ response weight;     datalines;     3 0 0.0025 3     3 1 0.00395 1     6 0 0.003 1     6 1 0.0042 1     9 0 0.0035 1     9 1 0.002 1     ;run;          proc glmpower data=exemplar;     weight weight;     class Var1 Var2;     model response = Var1 | Var2;     power         power=0.8         ntotal=.         stddev=0.05486016;     run;      Note: GLMPOWER only will use class (nominal) variables so 3, 6, 9 above are treated as characters and could have been low, mid and high or any other three strings. When the real analysis is conducted, Var1 will be used a numeric (and we will include a polynomial term Var1*Var1) to account for any curvature. The output from SAS is ![enter image description here](http://i.stack.imgur.com/T44tM.png) So we see that we need 762,112 as our sample size (Var2 main effect is the hardest to estimate) with power equal to 0.80 and alpha equal to 0.05. We would allocate these so that 3 times as many were the baseline combination (i.e. 0.375 * 762112) and the remainder just fall equally into the other 5 combinations.


title: Effect of log transform on Pearson correlation

text: Is there a generic answer as to how the log transform affects the Pearson correlation between two variables, i.e., is $\rho(\log(x),\log(y)) < \rho(x,y)$ always or the other way around and why?


title: How should I interpret the Kappa Fleiss test?

text: I noticed that when I have tables which the values are only 0 and 1, I get a kappa of 1 when the table is completely full of one, and when I have a table of zeros I get NaN as result using the irr package and the kappa fleiss function. I would expect the kappa to be also equals to 1 when the table is full of 0, as the 0 would represent agreement on 'No' and the full of 1s would represent agreement of 'Yes'. What is the right way to interpret the results?


title: Simulate censored data cox model

text: I would like to simulate interval-censored data in a Cox model. In the R package `intcox` I found the following code:               sim.weibull.intcox.rfc <-function(N=200, beta.0=0.1, beta.cov=c(0.5,-0.5,0.5,0.5),                                       alpha=0.75, p.cov=c(0.5,0.75), grid=10)     {       x.design <- cbind(rbinom(N,1,p.cov[1]), rbinom(N,1,p.cov[2]), runif(N,-1,1),                          rnorm(N,0,1))       colnames(x.design) <- paste("x.",1:4,sep="")       lambda  <- exp(beta.0+x.design%*%matrix(beta.cov,ncol=1))       scale   <- (1/lambda)^(1/alpha)       t.true  <- rweibull(N,alpha,scale)       T.max   <- max(qweibull(0.9,alpha,median(scale)))       t.left  <- NULL       t.right <- NULL       for (i in 1:N) {         tt <- unique(c(0,sort(runif(grid,0,T.max)),T.max))         if(t.true[i]>=T.max) {           x.left  <- T.max           x.right <- NA         } else {           x.left  <- max(tt[t.true[i]>tt])           x.right <- min(tt[t.true[i]<tt])         }         t.left  <- c(t.left,x.left)         t.right <- c(t.right,x.right)       }     return(data.frame(ID=1:N, left=t.left, right=t.right, x.design))     }      They applied Weibull distributed random variables with `shape=0.75`, while `scale` is derived from $(1/λ)^{(1/{\rm shape})}$ with $λ = \exp(β_0 + β'X)$ where $β = 0.5, −0.5, 0.5, 0.5$ and design matrix $X$ which is formed by the four covariates. My question is:   1. Why did they use $λ = \exp(β_0 + β'X)$ for calculating the scale parameter beta? Also, is their baseline hazard function decreasing given that `shape = 0.75`?


title: Expected agreement if random

text: I am looking to measure agreement between participants in choosing which member of a group is most like certain attributes. I want to calculate the expected agreement if they were to choose by chance. e.g. assume we have 3 participants (P1, P2, P3) and they can choose from a list of 3 people (A,B,C) to associate with a characteristic (Cx). The results could be:               characteristic   P1   P2   P3     C1               A    B    A     C2               C    B    B     C3               A    A    A     ...      If, in this case, the agreement is an average of P1vP2, P1vP3, and P2vP3 and agreement was {0,1} then we would have:               C1  1/3     C2  1/3     C3  3/3      How could I calculate the expected agreement if the choices were made at random?


title: What are appropriate tests for goodness of fit on glm with a small sample size?

text: I've thought quite a lot on large sample size inference where the strong law of large numbers is easily validated. In my case however, I'm trying to infer the sign and magnitude of an outcome where the noise-to-signal is quite large in comparison to my available sample size. **How do I properly test for goodness of fit? I possess a very small sample size ( _n_ = 16).** I cannot possibly get more data as it simply doesn't exist. I'm fitting a generalized linear model with Gaussian errors for which I obtained the following p-values with `summary(my.glm.fit)` using R's `glm()` function for fitting.               intercept   0.34172       slope       0.00734        * Using Cook's _D_ gives one data point that is problematic. If anyone knows about its validity with small samples, please speak out!   * I don't think a Durbin–Watson test statistic is appropriate for _n_ = 16. Even if my data has almost surely no auto-correlation, I have very sound reasons to think it doesn't. What could I use instead?   * Normality tests for small sample size...is Lilliefors OK? or should I go with a Bayesian test? Anything else you can point out? Anything else missing in my maybe too-broad question? I could provide more details if that's of any interest to you guys. Partial answers are also welcome. Thanks.


title: document similarity with documents using synonyms

text: I have a bunch of documents where some of the documents are a copy of other documents with their text jumbled up and some of the words replaced by their synonyms. Mentioned below is one such example of a sentence: > **_Article 1 (original)_** : I caught up with John Snow in town making > purchases at Kingslanding Hardware store to repair a broken tractor. Snow > has farmed soybeans his entire life, as did his father and their fathers. I > asked him about his life on the farm. > > **_Article 2 (duplicate)_** : I obtained John Snow which in city in > purchases make rise of the hardware at Kingslanding to repair a broken motor > tractor. Snow have soya broad beans complete life have been treated, such as > its father and their fathers. I asked him concerning its life on the > agriculture company. > > **_Article 3 (duplicate)_** : I took for above with John Snow in the city > that made purchases in the warehouse of the hardware of Kingslanding to > repair an broken tractor. Snow has cultivated the soybeans its whole life, > like its father and his parents. I asked to him about its life in the farm. > > **_Article 4 (duplicate)_** : I caught up with myself compared to John Snow > downtown making of the purchases to the kingslanding store of material to > repair a broken tractor. Snow cultivated soya its life whole, just as > his/her father and their fathers. I questioned it about his life with the > farm. I want to do a document similarity which ends up tagging all these documents in the same group. Any suggestions along with examples or tutorials will be greatly appreciated.


title: What is the probability that two independent random vectors with a given euclidean distance $r$ fall in the same orthant?

text: Consider two independent and identically distributed random vectors of dimensionality $N$, $\mathbf{x}$ and $\mathbf{y}$, where their elements are iid generated from a Gaussian with zero mean and variance equal to $\sigma^2$; i.e., $\mathbf{x} , \mathbf{y} \sim \mathcal{N}(0, \sigma^2)$. What is the probability that $\mathrm{sign}(\mathbf{x})=(\mathrm{sign(x_1)},\mathrm{sign(x_2)}, \ldots, \mathrm{sign(x_N)})$ be equal to $\mathrm{sign}(\mathbf{y})=(\mathrm{sign(y_1)},\mathrm{sign(y_2)}, \ldots, \mathrm{sign(y_N)})$ if their Euclidean distance is smaller than $r$, $r \geq 0$: $\Pr\left[ \mathrm{sign}(\mathbf{x})=\mathrm{sign}(\mathbf{y})\ \vert\ \Vert \mathbf{x} - \mathbf{y} \Vert \leq r\right]$. Note: $\mathbf{x}$ and $\mathbf{y}$ are independent.


title: Distribution of reciprocal of regression coefficient

text: Suppose that we have a linear model $y_i = \beta_0 + \beta_1 x_i + \epsilon_i$ that meets all the standard regression (Gauss-Markov) assumptions. We are interested in $\theta = 1/\beta_1$. **Question 1:** What assumptions are necessary for the distribution of $\hat{\theta}$ to be well defined? $\beta_1 \neq 0$ would be important---any others? **Question 2:** Add the assumption that the errors follow a normal distribution. We know that, if $\hat{\beta}_1$ is the MLE and $g(\cdot)$ is a monotonic function, then $g\left(\hat{\beta}_1\right)$ is the MLE for $g(\beta_1)$. Is monotonicity only necessary in the neighborhood of $\beta_1$? In other words, is $\hat{\theta} = 1/\hat{\beta}$ the MLE? The continuous mapping theorem at least tells us that this parameter is consistent. **Question 3:** Are both the Delta Method and the bootstrap both appropriate means for finding the distribution of $\hat{\theta}$? **Question 4:** How do these answer changes for the parameter $\gamma = \beta_0 / \beta_1$? **Aside:** We might consider rearranging the problem to give $$\begin{align*} x_i &= \frac{\beta_0}{\beta_1} + \frac{1}{\beta_1} y_i + \frac{1}{\beta_1} \epsilon_i \\ &= \gamma + \theta y_i + \frac{1}{\beta_1} \epsilon_i \end{align*}$$ to estimate the parameters directly. This doesn't seem to work to me as the Gauss-Markov assumptions no longer make sense here; we can't talk about $\text{E}[\epsilon \mid y]$, for example. Is this interpretation correct?


title: Heuristic for bagging in random forests

text: Is there a good heuristic for what fraction of the data to use when training a random forest with bagging? I imagine this fraction should depend on the number of trees in the forest, like 2/3 1/ _N_ trees, where _N_ trees is the number of decision trees in the forest. Please point me to a reference! Thanks!


title: Why does a time series have to be stationary?

text: I understand that a stationary time series is one whose mean and variance is constant over time. Can someone please explain why we have to make sure our data set is stationary before we can run different ARIMA or ARM models on it? Does this also apply to normal regression models where autocorrelation and/or time is not a factor?


title: Verifying that measured time does fit O(n^3)

text: I'm dealing with a homework. We were given a software, with inputs, that we should measure and then verify that the time complexity is O(n^3). I measured the software, picked best times from each batch (since that's the least affected one), and ended up with this:               304   70000     574   480000     775   1190000     850   1580000     1070  3130000     1557  9740000     1965  19540000      Now I'm really lost. I know how to calculate a model using least squares, but I have no idea how to do that for O(n^3). Plus I don't know how I would verify that the model is correct.


title: What statistical test do I use for differences in proportions?

text: I have two groups of samples with a proportion of mutations in each sample being a particular type (the proportions are represented as a fraction of 1). How to I check if that type of mutation is significantly likelier to happen in one group over the other?


title: Constructing an interval estimate for a multivariate output

text: I am building a multivariate Gaussian Process model that predicts output $X_i$ and $Y_i$ jointly from input location $z_i$. At the end of it all I want to be able to plot my best estimate of $X_i$ against $Y_i$ where I know that $$\begin{pmatrix} X_i\\ Y_i \end{pmatrix} \sim N_2 \begin{pmatrix} \begin{pmatrix} \mu_{X_i}\\ \mu_{Y_i} \end{pmatrix}, \begin{pmatrix} \sigma_{11i} & \sigma_{12i}\\ \sigma_{21i} & \sigma_{22i}\\ \end{pmatrix} \end{pmatrix}$$ So as you can see, each $X_i,Y_i$ pair has their own mean and covariance matrix. Now, in order to plot my best estimate of $X_i$ against $Y_i$ I can simply plot my my predicted values of $X_i$ against $Y_i$ (the red line shown below). Now what my question is whether or not it is appropriate (or if there is a better way) to construct a 95% interval estimate in the following fashion. So for each point $X_i,Y_i$ I can plots its 95% confidence ellipse as seen below (the blue ellipse) ![enter image description here](http://i.stack.imgur.com/rcSNL.png) And for another point I can do similarly. ![enter image description here](http://i.stack.imgur.com/pVcO3.png) And for a third point I can do the same again where we see that there is a very apparent difference in covariance matrix here for this set of $X_i$ and $Y_i$ ![enter image description here](http://i.stack.imgur.com/gGmEE.png) ![enter image description here](http://i.stack.imgur.com/ch4g9.png) ![enter image description here](http://i.stack.imgur.com/I6wno.png) And so finally we see that we could plot a 95% interval estimate for every single point $X_i$ and $Y_i$ (the last two plots, where the second to last plot has the same covariance matrix for all the points and the last plot has varying covariance matrix for all the points). Now, can I consider this a 95% joint interval for all points $X_i$ and $Y_i$ or is there a better way to go about this?


title: Bernoulli random variable parameter estimation

text: Suppose $\theta$ is the probability that a Bernoulli random variable is one (therefore $1-\theta$ is the probability that it's zero). I have a sequence of $n$ of these i.i.d. Bernoulli random variables, $m$ out of which are ones. It is a well-known fact that $\hat{\theta}=m/n$ is a maximum likelihood estimator (MLE) for $\theta$ (and it is also minimum-variance unbiased estimator). Now, in my case, I know that $\theta=\frac{1}{2}+\frac{1}{2}e^{-x}$ for some unknown $x$, where $x$ satisfies $0<x<\infty$. This implies that $\theta$ satisfies $\frac{1}{2}<\theta<1$. I am trying to find an estimator for $x$ given $m$ and $n$ (and, perhaps the full vector of observed random variables). So, I try to find an MLE by first writing down the log-likelihood function: $$\log\mathcal{L}(x;n,m)=-n\log 2+n\log(1-e^{-x})+m(\log(1+e^{-x})-\log(1-e^{-x}))$$ Then I take the derivative with respect to $x$ and set it to zero: $$\frac{\partial\log\mathcal{L}(x;n,m)}{\partial x}=\frac{n}{e^x-1}-\frac{2m}{e^{-x}-e^{x}}=0$$ Solving for $x$ yields the estimator $\hat{x}=-\log(2m/n -1)$. However, if less than half the random variables are one (i.e. $m/n<1/2$ -- and this certainly can happen in the scenario that I described), then the estimator yields an imaginary number (due to a negative in the logarithm). What did I do wrong? How does one estimate $x$ in this case? **EDIT** John A. Ramey pointed out that how for certain values of $\theta$, MLE leads to a degenerate answer. Perhaps there are methods other than MLE?


title: If k-means clustering is a form of Gaussian mixture modeling, can it be used when the data are not normal?

text: I'm reading Bishop on EM algorithm for GMM and the relationship between GMM and k-means. In this book it says that k-means is a hard assign version of GMM. I'm wondering does that imply that if the data I'm trying to cluster are not Gaussian, I can't use k-means (or at least it's not suitable to use)? For example, what if the data are images of handwritten digits, consisting of 8*8 pixels each with value 0 or 1 (and assume they are independent thus it should be mixture of Bernoulli)? I'm a little bit confused on this and will appreciate any thoughts.


title: Qvalue.R results

text: I have some p-values that I want to convert to q-values. The issue I am running into is that the qvalues being generated by storey's qvalue package are all in a very tight range (i.e. .25-.35). I am not sure what this means. I am working with a ridiculously small dataset, given the number of hypotheses, could that be the problem? I could understand if I had all high q-values, but that all of them having a ~.25 q-value... I'm just lost. I would really appreciate it if someone could explain what I am not understanding here about q-values. Thanks. Edit: A histogram of the p-values ![enter image description here](http://i.stack.imgur.com/tvUQK.png) Descriptive Statistics   * count 249.000000   * mean 0.395376   * std 0.263631   * min 0.006288   * 25% 0.173421   * 50% 0.351731   * 75% 0.588190   * max 1.000000 Dispersion   * p-val < .2 0.026   * p-val > .2 0.092


title: Mean absolute error vs sum absolute error

text: I'm using the Matlab Neural Network toolbox for a regression problem (twelve inputs, one target). My target has a high dynamic range and is strongly non- gaussian. So far I've solved this by a log-transformation, but as an alternative I'm exploring if a different error function than the standard mean square error would be of use. The Matlab Neural Network toolbox comes with four built-in "performance" functions:               >> help nnperformance       Neural Network Toolbox Performance Functions.              mae - Mean absolute error performance function.         mse - Mean squared error performance function.         sae - Sum absolute error performance function.         sse - Sum squared error performance function.      In my understanding, such a performance function is a cost function. But as a cost function is simply something to be minimised, what is the difference between _Mean absolute error_ and _Sum absolute error_ in practice? Similarly, what is the difference between _Mean squared error_ and _Sum squared error_? Shouldn't those be the same? The documentation for the different functions is of no use. In fact, the documentation for _mse_ and _sse_ is identical apart except for one word.


title: How to add weight to data in descriptive statistics?

text: I have a class with a set of descriptive statistic functions (mean, median, kurtosis, etc...). Now I need to include weight (array) into my equations. My first thought was to just create a weighted version of the functions where the additional weight array is passed. However, I was wondering if there is a way to alter the data (array) so that the same functions can be used. Note: This is my first Statistics project so I am learning as I go. Sorry if this is a stupid question. **Question:** Is there a way to weight the data before calculating the descriptive statistics - so that the existing functions can be used. If this is possible then I could add an additional function to weight the data before it is passed into the descriptive functions.


title: What is the difference between Exploratory Factor Analysis and Principal Components Analysis (PCA)?

text: I know what you're thinking, this is a duplicate of "What are the differences between Factor Analysis and Principal Component Analysis", but it isn't really. That other question deals with Confirmatory Factor Analysis. Either way, I would like to know what the difference is :) Thanks!


title: Using fast99 sensitivity method in R with normal distributions

text: If I use the fast99 method (from the sensitivity package) with qnorm, it generates parameter sets with -Inf and Inf in, which are not much use for then running models:               > f <-fast99(model=NULL,factors=c("a","b"),n=100,M=4,q=c("qunif","qnorm"),q.arg=list(a=list(min=0,max=1),b=list(mean=0,sd=5)))     > summary(f$X)            a              b             Min.   :0.00   Min.   :-Inf        1st Qu.:0.26   1st Qu.:  -3        Median :0.50   Median :   0        Mean   :0.50   Mean   : NaN        3rd Qu.:0.74   3rd Qu.:   3        Max.   :1.00   Max.   : Inf        I guess I could use qtnorm instead, but I'm a bit suspicious about making up upper and lower bounds, and it leads to some weird histograms:               > f <-fast99(model=NULL,factors=c("a","b"),n=100,M=4,q=c("qunif","qtnorm"),q.arg=list(a=list(min=0,max=1),b=list(mean=0,sd=0.1,upper=100,lower=-100))); summary(f$X)                a              b                       Min.   :0.00   Min.   :-100.00000            1st Qu.:0.26   1st Qu.:  -0.06433            Median :0.50   Median :   0.00000            Mean   :0.50   Mean   :   0.00000            3rd Qu.:0.74   3rd Qu.:   0.06433            Max.   :1.00   Max.   : 100.00000           > hist(f$X$b)      ![Histogram of qtnorm sample parameters](http://i.stack.imgur.com/99RQ3.png) So, is there any way to use fast99 with infinitely wide distributions, or is this a fundamental issue? And if I use a truncated normal, how should I set the upper and lower bounds so as to not bias the result unduly? I've generally preferred using fast99 to the Sobol monte-carlo based methods, but if it's unworkable, I can switch over.


title: Test for multinomial distribution

text: I have a sample of animals and a categorical variable of $L$ behaviors $i=1,...L$. I computed the proportions of this sample corresponding to each behavior category $p_{i}$, $\sum p_{i}=1$.   1. Can I use a common method of confidence interval for a proportion (e.g. Clopper-Pearson) to compute the CI for each of the $p_{i}$?   2. Is it possible to test that proportion $p_{i}$ is significantly different than $p_{j}$? for instance if their confidence intervals do not overlap?   3. Is it possible to test that $p_{i}$, $i=1,...n$ are equal? for instance conduct a chi-squared test of goodness of fit: `chisq.test(c(p1*n,p2*n,p3*n), c(1/3,1/3,1/3))`, if L=3 and my sample size is n. Thank you


title: Mutual information with outliers

text: I am using Mutual Information on a data set with some outliers. The univariate variables don't have a known distribution, so I cannot remove them automatically with any neat method, like Grubbs or Studentized tests. Hence, As I must work with a dirty data set, the methods that I apply must be robust against outliers. I want to apply Mutual Information, but the outliers makes a long tail in the distribution that hides the density information in the interesting range of values. Thus, I am thinking about unevenly dividing the bins of the histogram with any method, like Freedman–Diaconis rule. However, I couldn't find any similar approach on-line, so I would appreciate a confirmation that I am not incurring into any problem by doing so. Therefore, as a sum-up, using any non-equal width method to estimate the bin sizes of the 2D histogram will hide any dependency or increase the bias in the Normalized Mutual Information? Thank you,


title: Density estimation using knn

text: I am dealing with sentiment results of web articles. Sentiment is represented with two int values, +ve and -ve. For given article "xyz" I am getting different sentiments while testing at different time period. Each result has +Ve and -ve value and value are close to each other. I represent them using 2-D Graph on x-y axis. values are like this: x=10, y=100 - 10 For my given point (x,y) `Euclidean distance` for each point can be calculated but how it is useful? Because mypurpose is not to get nearest neightbour. But getting most relevant and eligible point from all available. My purpose is to get sentiment result with highest density. Is KNN also is suitable for this kind of scenario? Or any othe approach if someone can suggest .. I tried to understand this doc: http://www.statsoft.com/textbook/k-nearest- neighbors but not clear whether it can do what I want or not.


title: Applying LASSO with cross-sectional/time-series data

text: Even after reading the original paper on LASSO, I am still confused as to how to apply LASSO for variable selection when the response variable _y_ is a matrix (and not a vector) of _n_ dependent variables with _t_ observations each. As an example, the response variable _y_ is monthly household consumption for 50 households over 12 months, while the predictors are price of basket of goods, income... Should LASSO be done for each household over the 12 months, and the models for all households be averaged after that? Or should LASSO be applied across households for each month, then the models be averaged for all months?


title: multiple imputation

text: I am running multiple imputation for a set of variables including clinical data. I am wondering if I can use (or should use) outcome variable (follow-up is 99%) to predict missing clinical data. There is about 12 % of non complete cases due to mostly one variable. I intend to use Amelia package (R)? I would be most grateful for answer.


title: How to compare features and classifiers which achieve perfect accuracy?

text: So I'm looking to compare different combinations of features and classifiers. But I'm getting a lot of combinations that achieve 100% cross validation accuracy. I'm trying to figure out how I would compare the usefulness of each combination. For example I can both train an SVM using Features 1, 10, 15 to get 100% accuracy. But at the same time I can train a logistic regression classifier only using Feature 7 to get 100% accuracy. Also this is a binary classification problem.


title: Random sample before or after data transformation

text: I have a very basic model development question. The objective is to develop a predictive model. Let's say I have a dataset with 50,000 observations. Take a random split of this dataset into train (N=25000) and validation (N=25000) datasets (we know some folks prefer a test dataset but assume we are not going to do that for now).   * Should I transform (binning, any scaling, creating dummy variables, declaring factors etc.) the original raw dataset and then take the random sample for train and validation datasets?  **OR**   * Should I apply data transformation to the train set first and then apply the data transformation schemes learned from train set to the validation set? Binning can differ if applied to the raw dataset as compared to train set if we are trying to use supervised binning algorithms like chi-merge or other optimization algorithms.


title: Developing a prediction model for bus stops

text: I need to predict if a bus will stop at a particular scheduled route stop. I have GPS coordinates, I can check door signals, the stops have radii, and I also have historical data about if the bus skipped a stop or not. I want to use this information to train a model predicting whether a bus will stop at a particular stop or not. I was looking at Bayesian rule, and I was thinking it might be possible to use a likelihood ratio to predict the bus reaching a stop. But I'm not sure if this is the right way to look at this problem. I'm not a statistician, so all the jargon to choose the right way to model this problem makes me nervous. I would like to get your advice on this. I'm willing to put in time to learn, but I don't want to spend a few days to find out this might not be the way to look at it in the first place. I think my problem is limited understanding of various probability based modelling techniques. If my statement of my problem needs any clarification, or any section needs more detail, I'm happy to add it. Any help would be much appreciated. Cheers! UPDATE: Sorry guys—these are buses that take employees to work and the other way around. So these buses have the liberty to avoid a stop if there are no passengers; the buses can even avoid the last two stops, for example, when taking the passengers from work to home, they usually go back to a different place to pick people from. I have GPS data from which I can check if a bus reached a stop, and there's also data for route configuration. In each route there is data about if a stop is mandatory as well. With this data in hand I'd like to predict if a bus will reach a certain stop. Please let me know if you require further information. Thanks!


title: Linear regression with strongly non-normal response variable

text: I have carried out a linear regression. The plot below shows the distribution of the response variable: ![enter image description here](http://i.stack.imgur.com/4RtjY.jpg) I believe the response variable is beta distributed, therefore virtually the exact opposite of normally distributed. However, when including all my predictor variables in the linear regression, the residuals turn out to be quite normally distributed, as shown in this plot: ![enter image description here](http://i.stack.imgur.com/kFJSJ.jpg) Has my model satisfied the assumptions of linear regression? Might there be a better model to use?


title: Non linear maximization problem

text: ![enter image description here](http://i.stack.imgur.com/QcPfC.jpg) I am not sure how I can solve this optimization problem. I have a (MATLAB) code which can solve this but only for very low dimensions, that is, low values of $n$ (around $n = 4$ or $5$), once $n$ exceeds $10$, the code suffers from severe curse of dimensionality and takes days to solve. What's more is that I need to solve the above maximization problem for many different sets of $P_i$'s (around 2000 sets), that is, I have 2000 different sets of $\{P_1, P_2, \cdots, P_n\}$. Numerical approximations to this would suffice. This problem is derived from a much larger problem (which I have not shown here) and the lambdas represent Lagrange multipliers of the original problem. If anyone needs more context for this problem (or even papers that hints how to solve this problem), please let me know and I will edit/email you any further details.


title: Memory problems using pROC library in R to estimate confidence interval for area under the curve

text: I am trying to use the R pROC package to estimate the confidence interval for the area under the curve (AUC). I am running into a memory limitation. My R script looks like the following:               library(pROC)     data <- read.csv("data.csv",header=TRUE)     rocobj <- roc(data$observed, data$predicted)     ci(rocobj)      I have modified my 64-bit R v3.0.2 shortcut with the following > "C:\Program Files\R\R-3.0.2\bin\x64\Rgui.exe" --max-mem-size=8192M still, when I run the ci(rocobj) command, I get the following message:               Error: cannot allocate vector of size 256 Kb     In addition: There were 11 warnings (use warnings() to see them)     warnings()     Warning messages:     1: In m * n : NAs produced by integer overflow     2: In lapply(X = X, FUN = FUN, ...) :       Reached total allocation of 8192Mb: see help(memory.size)     3: In unique.default(unlist(lapply(x, length))) :       Reached total allocation of 8192Mb: see help(memory.size)     4: In unique.default(unlist(lapply(x, length))) :       Reached total allocation of 8192Mb: see help(memory.size)     5: In lapply(X = X, FUN = FUN, ...) :       Reached total allocation of 8192Mb: see help(memory.size)     6: In lapply(X = X, FUN = FUN, ...) :       Reached total allocation of 8192Mb: see help(memory.size)     7: In unlist(x, recursive = FALSE) :       Reached total allocation of 8192Mb: see help(memory.size)     8: In unique.default(unlist(lapply(x, length))) :       Reached total allocation of 8192Mb: see help(memory.size)     9: In unique.default(unlist(lapply(x, length))) :       Reached total allocation of 8192Mb: see help(memory.size)     10: In unique.default(unlist(lapply(x, length))) :       Reached total allocation of 8192Mb: see help(memory.size)     11: In unique.default(unlist(lapply(x, length))) :       Reached total allocation of 8192Mb: see help(memory.size)     Warning messages:     1: Reached total allocation of 8192Mb: see help(memory.size)      2: Reached total allocation of 8192Mb: see help(memory.size)       my laptop is windows 7 64-bit with 16 GB of RAM. But I am running other stuff, and so I can only really use 8 GB (which is already pushing it). I tried to limit the data size and it works.               n <- 10000     o <- head(data$observed,n=n)         p <- head(data$predicted,n=n)     rocobj <- roc(o, p)     ci(rocobj)      But obviously, I want more than 10,000 data points. Any help on alternative packages or how to resolve this memory issue is appreciated. By the way, I have 165,000 predicted points.               > head(data)       observed   predicted     1        0 0.005101049     2        0 0.042003527     3        0 0.045466933     4        0 0.017828838     5        0 0.092385962     6        0 0.053715844     > summary(data)         observed        predicted              Min.   :0.0000   Min.   :        0        1st Qu.:0.0000   1st Qu.:        0        Median :0.0000   Median :        0        Mean   :0.1714   Mean   :    11071        3rd Qu.:0.0000   3rd Qu.:        3        Max.   :1.0000   Max.   :277494082       >       using the caTools package, the script to compute the confidence interval of AUC is as follows.               r <- 1000     s <- 120000     bAuroc <- replicate(r, {      m <- sample(s, replace=TRUE)      colAUC(data$predicted[m], data$observed[m])     });     quantile(bAuroc, c(0.05, 0.95))      note that changing the value of r (e.g. r = 1000, 2000, 3000, ...) does not really produce any problems. however, changing the value of s (e.g. s = 100000, 110000, 120000, 130000) does produce problems. it seems that s = 120000 is the upper limit, and going anything above that produces the integer overflow problem. the solution for my problem is as follows. note that there is no more limit on the number of samples, s < nrow(data) but s = nrow(data), and also, i had to specify alg="ROC".               r <- 1000     s <- nrow(data)     bAuroc <- replicate(r, {      m <- sample(s, replace=TRUE)      colAUC(data$predicted[m], data$observed[m], alg="ROC")     });     quantile(bAuroc, c(0.05, 0.95))


title: Test statistic and null rejection of independence in normalized redundancy, as implemented and automated in Python

text: I currently have an output function that roughly follows a $\chi^2$ distribution for random data with one degree of freedom. As a cursory analysis that has since proven mathematically unrigorous, this appeared to be bounded by $[0,1]$ on the x-axis, with a mode of 0, and a steep falloff such that my PDF offers highly significant results after a value of roughly 0.1 on said x-axis. My input data set is derived from the normalized redudancy of categorical bins of data. My null hypothesis is that the data coming in for two arbitrarily selected datasets are independent (that is, $p(x,y) = p(x)\ p(y)$ for all $x \in X$ and $y \in Y$ ). My test statistic is the measure from which this test for independence deviates from the observed data. At present, the closed mathematical model for my test statistic remains undevised. Unsurprisingly because of this, after slogging through the documentation for various tools to provide a statistical measure of significance in the context of this distribution, including SciPy (the framework I am using for other tests, as all of this work is currently in Python), I am unable to come to a good, working definition for my test statistic and significance thereof. Even if I do, I am unsure how to implement it effectively, as tools like SciPy are sparsely and inadequately documented. I understand, for example, how to integrate an arbitrary distribution, but not how to construct one of my devising in line with these requirements. To this end, I've played with the input parameters for `scipy.stats.chi2` and can get a non-robust, 0th level approximation that roughly aligns with my observed output. However, I am unable to see how my model aligns with the constructs provided by SciPy, as I fail to grasp SciPy's object model in the context of the statistical and mathematical definitions of, for example, a $\chi^2$ distribution. I am also not cognizant on the correct terminology to further specify this question, and as evidenced, I have altered its title and content severely because of this. So, as a further specified question, this becomes two pieces:   * Given these requirements, how might I construct a good, statistically rigorous test statistic and rejection criteria for null? This will be used in an automated manner, so showing the closed form for the approach is greatly appreciated, albeit not required, to solve this problem.   * How might I implement this in Python? Because the documentation is so tersely insufficient, links to further sourcing and crisp definitions of all constants, variables, and non-obvious operators (to a high school graduate with a fundamental course in probabilistic statistics) would be greatly appreciated.


title: Overcoming memory constraints in rpart?

text: I have 4.4 million observations, 160 binary features, and a binary response. Using rpart on Windows (64-bit, with the 64-bit R v2.13.0 build), I run out of memory on a machine with 64GB RAM. My memory limit seems reasonable:               > memory.limit()     [1] 61440      But on running rpart, I still receive this error after watching memory usage spike:               > mytree = rpart(formula=fmla,data=mydata,method="class",xval=0)     Error in rpart(formula = fmla, data = mydata, method = "class",  :        cannot allocate memory block of size 5.0 Gb      This r-help post suggest that memory fragmentation may be an issue. Rebooting the machine and starting with a fresh workspace has no effect, nor does --max- mem-size=60G as recommended here. I'm hoping to avoid sampling from my population. I've tried using logical and factor representations of my data. I've also tried `tree` and looked at `gbm`, `randomforest` and `party` packages, none of which appear to offer more efficient multi-class CART implementations. Can anyone suggest an alternate package or different run parameters? Barring that, I'm looking for a package to streamline running it through rpart in pieces. I could run rpart on subsets of features and choose the best split from all my runs at each node, but reassembling the results into an rpart object for plotting and pruning sounds like a chore. Can anyone suggest a cleaner solution or a pre-written package for this? Thanks.


title: Converting time to readable form in R

text: 10/16/2011  9:33:46 AM     10/16/2011  9:54:48 AM      The above are a few rows from my G$time column. I wanted to convert them to readable form when I graphed Y over Time, so I tried               Time<-strptime(G$Time,format="%D  %H:$M:%S %r")      and then               Time<-strptime(G$Time,format="%m/%d/%y  %H:$M:%S %r")      but it's not working. Can anyone help?


title: Beta binomial Bayesian updating over many iterations

text: I'm using a beta binomial updating model for a piece of code that I am writing. The software is real time updating - meaning that data is continually being gathered and after N data points are gathered, the bayesian model is updated using the N data points. Under this logic, I am using the posterior output as my prior for the next iteration. My problem is that over billions/trillions/maybe more of iterations, the bayesian beta parameters (alpha and beta) will grow very large. I am worried that eventually the parameters will become so large that they will cause an integer overflow in memory. So my question is twofold -   1. Is it reasonable to be worried about this integer overflow. I understand that $2^{32}$ is an extremely large number, but I'm building this software for an internet service that will be running 24/7, 365 days a year and I don't want it to crash. For example if I was updating it with 1,000,000 data points a day then the model would only last ~4000 days before an integer overflow.   2. Is it possible to transform a Beta(x,y) r.v., where x and y are extremely large, to a Beta(x*,y*) r.v. where x* and y* are relatively smaller? The transformed Beta doesn't have to be exact, just similar.


title: How to form a meaningful statistical indicator to reflect user interaction with a website

text: A given user can interact in multiple ways with a website. Let's simplify a bit and say say a user can:   * Post a message   * Comment a message   * "like" something on the website via Facebook _(after that we could add, following the site on twitter, buying something on the site & so on, but for readability's sake let's stick to these 3 cases)_ I'm trying to find a formula that could give me a number between 0 and 100 that reflects accurately the user interaction with the given website. **It has to take the following into account:**   * A user with 300 posts and a one with 400 should have almost the same score, very close to the maximum   * A user should see his number increase faster at the beginning. For instance a user with 1 post would have 5/100, a user with 2 would have 9/100, one with 3 would have 12/100 and so on.   * Each of these interactions have a different weight because they do not imply the same level of involvement. It would go this way: `Post > Comment > Like`   * In the end, the repartition of data should be a bit like the following, meaning a lot of user around 0-50, and then users really interacting with the website. ![enter image description here](http://i.stack.imgur.com/65JVm.jpg) * * * This is quite specific and data-dependent, but I'm not looking for the perfect formula but more for how to approach this problem.


title: Can statistics based on crude units of measure be improved?

text: I'm collecting online recipes in order to calculate average ingredient proportions for, say crepes or pastry. Most volume based units of measure in english recipes (cups, tablespoons, teaspoons) differ from country to country and sometimes differ within the same country depending on whether traditional, metric or government standardised units are used. I have used my own guesses to convert from these units to milliliters. So, for example, if it is unclear exactly which _cup_ measure is indicated, I use a metric cup (250ml). I wonder if there are tools that can be used to infer the most likely conversion for these crude units given that a percentage of the data is recorded using more accurate units (grams, milliliters). Or, failing that, to account for the uncertainty in some way (larger confidence intervals?). **Update:** I see that some questions have comments asking what the stats are used for. In case it is helpful, I'm using them to answer questions like "How do crepes differ from Swedish pancakes?" and "How do English recipes for crepes differ from French recipes?". For example... ![Ingredient proportions for pancake batter](http://i.stack.imgur.com/iEAMA.png)


title: Next step after principal component analysis

text: I am looking for improving a certain process. Several variables are measured about the performance of this process. My endeavour is to identify some critical variables so that I can run DOE to optimize this process. I ran a PCA. I find that too many variables are entering each PC. What can be done?


title: Confidence interval for proportions

text: I have some data like this:               id    pop   var     1     593   51     2     592   31     3     346   20     4    1214   70     5    1063   66     6    1370   71      each `id` represents a school. There are around 1847 of them, of which 1833 are unique . `pop` is the number of pupils in the school. `var` is the number of pupils with a particular trait. Interest in centered on the ratio of `var/pop` . I have used a binomial regression:               fit.bn <- glm(cbind(var,pop-var)~1, data=dt, family=binomial(link=logit))      which produced the results:                            Estimate Std. Error z value Pr(>|z|)         (Intercept) -2.748516   0.003295  -834.2   <2e-16 ***      And to get the estimate for the proportion and it's confidence interval I used:               invlogit(coef(fit.bn))     invlogit(confint(fit.bn,level = 0.95))      where               invlogit <- function (x) { return(exp(x)/(1+exp(x))) }      which produced:               (Intercept)       0.06017054                2.5 %     97.5 %      0.05980603 0.06053643       Is this a valid confidence interval for these data ? I was surprised that it was so narrow, but I guess this is because the total number observations is high (altogether, `sum(var)= 98008` and `sum(pop)=1628981`). Nevertheless, I am still unsure about this, due to:               sd(var/pop)     0.01373888      In other words the standard deviation of these proportions is around 20 times larger than their confidence interval. Am I doing something stupid ? It just doesn't seem quite right. I wondered about using a variance components model to try to account for variability between schools and within schools. I did this with               require(lme4)     fit.lme4 <- glmer(cbind(cbind(var,pop-var)~1+(1|id), data=dt, family=binomial(link=logit))      This produced:               Random effects:      Groups   Name        Variance Std.Dev.      id      (Intercept) 0.024125 0.15532      Number of obs: 1847, groups: id, 1814          Fixed effects:                  Estimate Std. Error z value Pr(>|z|)         (Intercept) -1.776882   0.005423  -327.7   <2e-16 ***      So I computed another confidence interval by:               invlogit(-1.776882  + (c(-1.96,1.96)* 0.005423))     0.1433781 0.1460089      which is wider than the first one, but the estimate is very different. Is this a valid procedure ? How do I interpret this estimate and confidence interval as compared to the first one ? And how do I interpret the variance of the random effects ? Am I barking up the wrong tree ? **Update:** Here is a histogram of `var/pop` ![enter image description here](http://i.stack.imgur.com/Mxj6O.png)


title: Independence of mutually exclusive dummies in same regression?

text: I have been wondering about this for a while. When I have mutually exclusive dummies in the same regression, are their coefficients independent? For example, suppose I have the following: $y_{i}=\mathbf{\mathbf{\alpha'}}\mathbf{x_{i}}+\beta_{0-5}I\left\{ z_{i}\in\left[0,5\right]\right\} +\beta_{6-10}I\left\{ z_{i}\in\left[6,10\right]\right\} +\beta_{11-15}I\left\{ z_{i}\in\left[11,15\right]\right\} +\epsilon_{i}$ where $\mathbf{x_{i}}$ is a vector of controls and $z_{i}\in\left\{ 0,1,2,3...15\right\}$ . In the end, I would like to test a non-linear combination of the coefficients $\beta_{0-5}$ , $\beta_{11-15}$ , and $\beta_{6-10}$ . Here is my guess: $\hat{\beta}_{0-5}=f\left(Y,X|Z\in\left[0,5\right]\right)$ and $\hat{\beta}_{6-10}=f\left(Y,X|Z\in\left[6,10\right]\right)$ and so the random variables $\hat{\beta}_{0-5}$ and $\hat{\beta}_{6-10}$ are functions of different random variables, which leads me to think they are independent. But something doesn't seem right to me. Don't the $x_i$ and $y_i$ of all $z_{i}$ -categories affect the estimates of $\alpha$ ? And don't the estimates of $\alpha$ affect the estimates of the $\beta$ 's (e.g. the more $\alpha$ explains about the variance of $y_{i}$ , the less the $\beta$ 's will)? Does this depend on the type of regression? I am interested in linear regression and ordered logit regression (where $y_{i}\in\left\{ 0,1,2,3,4,5\right\}$) . I don't have a specific application in mind but I can come up with an example if that helps.