title: How to evaluate a clustering/unsupervised learning problem with massive amounts of data, with labels only for a small fraction of points

text: I'm wondering if anybody can point me to work on the evaluation of unsupervised learning where there are a very large (say hundreds of millions) number of points and manual labelling can only ever be obtained for a small fraction of the points. To elaborate, most clustering evaluations assume gold standard labels for the data points. While I can feasibly obtain this for a small number of points (say in the thousands, or perhaps tens of thousands) this is a small fraction of the total. Furthermore it is obviously not possible to cluster only the labelled points, since the presence of all the unlabelled data would fundamentally change the shape of the solution. The labels are also not well defined (they are a function of the items contained in the cluster) so assigning manual labels to points individually is unlikely to be consistent. One possibility that occurred to me is to label the pairs of points as to whether they should be linked or not: however, chosen at random two points have a vanishingly small chance of being linked, so the true positive rate on the links would have very high variance. I would appreciate if anyone knows of work in a similar area that they could point me to, or can describe an evaluation strategy known to work in such circumstances.


title: How do you handle the situation where the residual variance is very high compared to the other variance parameter estimates?

text: **Context** An experiment in agronomy whose aim is to investigate the possible effect of a treatment, with 13 possible levels, on the height of trees. **Model** $ Y_{ijk} = \mu_{\cdot \cdot \cdot} + \alpha_{i} + \beta_{j} + \gamma_{k(j)} + (\alpha \beta)_{ij} + \epsilon_{ijk} $   * $Y_{ijk}$ is the response for the tree lying in the $k$th row of the $j$th bloc when it has received the $i$th treatment,   * $\mu_{\cdot \cdot \cdot}$ is an overall constant,   * $\alpha_{i}$ are the fixed treatment effects,   * $\beta_{j}$ are the random bloc effects,   * $\gamma_{k(j)}$ are the random row (nested within bloc) effects,   * $(\alpha \beta)_{ij}$ are the random treatment-bloc interaction effects,   * $\epsilon_{ijk}$ are random error terms. **Two important features**   * There is a lot of heterogeneity in response within each treatment.   * The interaction $(\alpha \gamma)_{ik(j)}$ cannot be estimated because there is no replicate. **Partial results** The residual variance is much much higher than the variances of the different random effects. As a consequence, a much simpler model without random effect is selected based on the AIC. EDIT relative to Nick Sabbe's comment: The simpler model I am talking about is $Y_{ijk} = \mu_{\cdot \cdot \cdot} + \alpha_{i} + \epsilon_{ijk} $ **Question** My interpretation is that the residual variance actually contains two parts: the residual variance itself, and the interaction that cannot be estimated. Now, intuitively, I think that that interaction cannot be simply ignored. Hence, I would not compare my model with a simpler model without random effect. **Do you agree with that?**


title: KS, AD and loglike results

text: I'm using R to test some distribution families to my data. I've done KS, AD tests and determined the loglike. For one of the data the indications given by KS and AD do not agree with the ones given by the loglike:               Table: p-values     Test    Normal  Log-normal  Gamma   Logistic    Weibull Gumbel     KS      0,16    0,00        0,00    0,26        0,00    0,49     AD      0,17    0,00        0,00    0,27        0,00    -          Measure Normal  Log-normal  Gamma   Logistic    Weibull Gumbel     loglike 282,86  279,54      308,96  284,41      304,00  291,55      I've read that KS gives more emphasis to the middle part of the distributions and AD to the tails. On the loglike one maximizes the probability of a model fitting the data. The graphical analysis says to me that Log-normal, Gamma and Weibull only fit the data well at the left tail whereas the other distributions fit data quite well all over the domain.... So why does these three dist. have a larger loglike than the others that seem to fit better the data? Thanks


title: Simulating ordinal variables using fitted probit models

text: I have fitted a probit model for an ordinal response and a number of predictors, using polr function in R. Now I want to use this fitted model in order to sample from the conditional distribution of the ordinal variable, given a set of specific values for the predictors. I was thinking of using multinomial distribution, with the probabilities computed based on the coefficient estimates from the fitted probit model. But then I am thinking that this approach ignores any lack of fit for the model. I am not sure how this simulation could be performed. Perhaps I am missing or confusing some key concepts. Any help would be appreciated.


title: Log-Ratio \ Compositional analysis

text: I am not a trained statistician but I am trying to improve on my own MBA thesis which was essentially regression analysis of the factors affecting opening **cinema box-office** in the UK. I am now trying to look at the mix of marketing channels used by distributors, specifically **TV, Press, Outdoor, Radio and Web**. My thesis looked at the absolute levels of marketing spend in each channel, but I am now more **interested in the combinational proportions of spend across the channels** ; (e.g., 35% TV, 25% Press, 25% Outdoor, 10% Radio and 5% Web) which of course sum to 100%. The idea is to try to determine an optimal mix under a range of scenarios of budgets opening dates and screen-counts. I have read that **Log-Ratio compositional regression analysis** can get around the obvious cross-correlations of inputs that sum to a constant, but I do not know where to start. **_Can anyone please point me to a novice's guide to log-ratio analysis_** , especially one that doesn't baffle me with too many strange symbols without definitions, (although I have printed the symbols list from Wikipedia!). I am very competent with **MS.Excel** , (although I would prefer not to have to write macros unless I have to), so basing the analysis in it would please me, (and I don't really have any money for expensive 3rd party software). Would anyone be able to kick-start me please?


title: How to account for different degrees of familial relatedness in linear mixed model? Is heteroscedasticity a problem?

text: (Let me preface this by stating that I am somewhat of an R newb and also a LMM newb.) I have a very, very simple research question. I have a continuous outcome variable ("HAPPY"). I have two factors that divide my sample into 4 cells (Factor1, Factor2), and I want to assess whether there are any main effects or an interactions. A simple 2 x 2 ANOVA should do the trick. There are two issues:   1. The individuals in my sample are family members. Every single individual has a fraternal twin or an identical twin in the dataset. I know that this throws off all the independence assumptions, so I was instructed to use linear mixed models and include the family identification variable (which identifies an individual as being part of a twin pair) ("Family") in my dataset as a random effect. Right now, this is my R syntax:              model <- lme(HAPPY ~ Factor1*Factor2, random= ~1|Family, data=mydata,                  method="REML", na.action=na.omit)     summary(model)     anova(model,type="marginal")      So right now, there is only one random effect. But -- theoretically, we'd expect that identical twins and fraternal twins would differ in their degree of relatedness (the variable is called "twintype"). **Should this be reflected in my model, and if so...how?** Conceptually, I'm struggling with how it would be specified. And practically, I'm also struggling -- putting the above syntax together was actually quite the accomplishment for me!   2. Heteroscedasticity. **What is the best way to assess if heteroscedasticity is affecting my model? Is it something I need to be concerned about and check for?** Right now I'm just eyeballing--`plot(fitted(model), resid(model))`--and looking for cone shapes. **And if heteroscedasticity is an issue, what can I do to fix it?** From what I've gathered from the google machine, I can enter a weighting command that will account for the inconstant variance. Could I modify my syntax above to read --              model <- lme(HAPPY ~ Factor1*Factor2, random= ~1|Family, data=mydata,                   method="REML", na.action=na.omit, weights=varPower())


title: Open-sourced pairwise learning models

text: I am solving classification problem using pairwise-learning training set. We have 2 classes: _bad_ and _good_. We also have pairs of objects $(a_i,b_i)_{i=1}^n$, meaning that object $a_i$ _is better_ than $b_i$. Each object is described through its real coordinates: $a_i = x_i^1, x_i^2, \dots x_i^k$. My aim is to build such an algorithm, that says for test object which is bad or good. The question is what open-sourced algorithms/programs/packages can you suggest for solving that kind of _pairwise-learning_ problem. I have found only svm- rank algorithm from svm-light project that handles such pairwise training data. I wonder, didn't find anything similar in matlab/R/rapidminer.


title: Unit root test for ARIMA models

text: I have a slight confusion regarding **seasonal** models and which polynomial to use for conducting unit root tests. Given a model: $\phi(B)\Phi(B^s)\Delta^d\Delta^D_S X_t = \theta(B)\Theta(B^s)\epsilon $ When I compute the roots for unit root tests, **which polynomials do I calculate the roots of**? Only the basic $AR(p)$ and $MA(q)$ compononents? That is, do I use just $\phi(B)$ and $\theta(B)$? Or instead, do I also include the seasonal polynomials as well? ie, compute the roots by solving for $\phi(B)\Phi(B^s)$ and $\theta(B)\Theta(B^s)$ ?


title: Adjusting for zero mean (standardizing) in a multiple regression model

text: A friend of mine was telling me today about the need to mean adjust input variables to zero in order to "get rid of implicit intercept" or scale terms in the slope coefficients and to make their interpretation easier. He was stressing its use and even standardization by the weighted mean in a weighted regression. Both of these I had never quite heard of other than when the inputs are of a very different scale (like 6 digit dollar values and single digit times or something) for stability of the estimation. Is it sensible what I am hearing? And that this "mean shifting" increases correlation between the slope coefficients? If so where can I read more about this? The argument I was given goes something like this: If i have the model $$y = ax + b + \epsilon$$ and then I shift $x$ by $D$, I will have $y = a' + b(x+D) + \epsilon$ where the estimate for $b$ does not change (and neither do the residuals) but the estimate for the intercept should be shifted accordingly. Thus we can equate $y - \epsilon$ to see that $a' = a - bD$ suggesting that shifting $x$ makes estimates of the intercept more correlated with $b$ than before. It feels to me like no new information has been added at all and there is nothing better about the coefficient estimates with shifts than with no shifts.


title: questions that can be asked for linear regression and logistic regression

text: I am preparing a potential interview, which can involve detailed questions on linear regression and logistic regression. Besides going over the textbook, could you share some possible interview questions on these topics based on your industrial and education experiences? Thanks.


title: What's the simplest way to create a beanplot in matlab?

text: What is the simplest script for generating a beanplot (see article published in the _Journal of Statistical Software_ ) in Matlab?


title: Is the Friedman Test Suitable For an Analysis of Pre/Post-Intervention Choices From Survey Data?

text: We sent out a survey to find if educational material about a type of surgery would potentially influence the decision of subjects. Prior to giving them this educational material we asked them demographic questions (e.g. sex, job, etc.), several questions on a 5 point scale to assess their priorities (e.g. "How important is cosmesis?"). After receiving this educational material we then asked them which type of surgery they potentially would be interested in. For this analysis, I would like to see if   1. the pre and post-educational material preferences changed   2. if demographic information (e.g. age, education, sex) has an influence on choice after receiving the educational material   3. if the pre-educational material priorities are associated with change in choice after receiving This survey was given to a random group of medical professionals, and I am assuming that the distribution of the data collected from this population is non-normal. This makes me think that I need to perform non-parametric tests, but I am having difficulty in figuring out which test would be best. So far I think I may have to perform a Friedman Test since I want to do a pre/post comparison between the same group. However I would appreciate if anyone could clarify the matter for me. Thanks!


title: Reinforcement Learning Value Iteration Explained

text: I am having a hard time to understand the value iteration derived from the Bellman equation: $$V_{k=0}(s) = 0$$ $$\forall s : V_{k+1}(s) = \max_a \Bigr[ R(s,a) + \gamma \sum_{s'} P(s'|s, \pi(s)) V_k(s') \Bigr]$$ I am trying to understand this on a concrete example I have thought of: A maze and a rat in it:                                           +------------------+                                 | s2             s3|                                 +------+ s1 +------+                                        |    |                                        |  R |      There are only three states. $s_1$ is a T-junction with the possible actions, left and right. $s_2$ and $s_3$ is food and something else poision or so. Thus granting rewards $(2,0)$. Everything is determinstic and I would like to use a discounting of $\gamma = 0.5$ How would a concrete example look like: Here are my thoughts:   * $R(s,a)$ is the reward which would be granted immediately (isn't it?), so I guess for $s_1$ this would be $0$?   * the sum $\sum_{s'}$ is the sum over the possible transition states, that would be $s_2$ and $s_3$?   * the transition probability $P(s'|s, \pi(s))$ is not given, but maybe it's just $1/2$ When I try to write it down it does not seem to make a lot sense: \begin{align} V_{k=1}(s_1) = 0 + 0.5 + (0.5 \cdot 0 + 0.5 \cdot 0) = 0 \end{align} This would just become zero again, is that correct? \begin{align} V_{k=1}(s_2) = 3 + 0.5 + (0.5 \cdot 0 + 0.5 \cdot 0) = 3 \end{align} How would the $a$ be evaluated? Its the state of the rat and it should choose the one with the highest result, but what are possible $a$?


title: Longitudinal data where individual id is not available

text: I am associated with a teaching hospital, and medical students are frequently administered tests and measures to monitor their progress. I am not talking about exam results here, but of psychological measures that do not figure on their transcripts. For "ethical" reasons, researchers typically need to annonymize the results. To be specific, I currently have "burnout" measures from a class of first year students. I also have "burnout" measures from a class of third year students ... the same individuals, two years later. I want to compare the results, but I don't have student ID's. A matched pairs analysis would be nice, but I don't have the matching variable. The two samples (year 1 and year 3) are dependent, since they involve the same individuals. A two-sample t-test seems dodgy in this case. I imagine that a permutation test might work, but I'm not sure which one to use. Any ideas?


title: How to best analyze a dataset with only categorical variables?

text: I have a dataset with all multinomial (categorial) variables. The column names look like: religionID, occupationID, sexID, ..., changedID The independent variables are: "religionID" (a factor with 7 levels), "occupationID" (a factor with 3 levels), .... There are 15 total, all factors with multiple levels. The dependent variable is: "changedID" (a factor with 4 levels). In all variables, the ordering of the levels is arbitrary. (E.g., factor level 1 is not "better than" or "greater than" factor level 2) I have about 500 rows of such data. My goals:   1. Find out which independent variables effect the dependent variables (e.g., does religion matter at all to changedID?)   2. Find out which independent variables have more effect than other independent variables (e.g., is religion more important than sex?)   3. Find out which factor levels in each independent variable has the most effect (e.g., which religionID factor leads the most often to changedIDs of factor level "2"?) Of course, I would also like to know if the results above are statistically significant. I've thought of building contingency tables (using R's `table()` function), and then running `chisq.test()` on the result. I guess this could be used to satisfy goal #1, if I build one such table for each independent variable. But it doesn't help with goal #2 and goal #3. Multinomial regression might be the answer, but I am confused as to whether it applies to my situation, and if it does, how to use it and interpret the results. Bonus points if I can have an solution in R.


title: Extremely large class set for support vector machine (SVM) classification

text: For a problem where I have tens of thousands of different classes, wouldn't it be very inefficient to use the typical multiclass methods? If I were to do one vs all, wouldn't that mean that I have would have to run through every individual classes before determining if a new set of data belongs to a certain class? That would be computationally very slow. Is there a way to get around this? Should I even use SVM's when it gets to such a size?


title: Spatially auto-regressive two-stage model

text: I'm working on a project in which I use a 'Generalized Spatial Two-Stage Least Squares' model, mostly known as $y= X \beta + \lambda W y + u$ and $u = \rho M u + \epsilon_n$ where $y$ and $u$ are vectors. The beta's are the general regressionparameters (for instance, from linear regression). $W$ is a weight matrix, which could represent the inverse distance between each two instances, the same holds for $M$ (but not necessarily the same distance-meaning). $\lambda$ and $\rho$ are two parameters which has to be computed by iteration and have both an absolute value strictly smaller than 1. More information about this technique, can for example be found in this paper. By making use of the program R I've already estimated the following variables $\beta_0= -69900$ (intercept), $\beta_1=40630$, $\beta_2=-50$, $\lambda=0.8285$ and $\rho=-0.613$ My problem is that I don't know how to interpret this model as 'one model', ie as one equation. Another problem is that I have no idea how to use this model to make predictions by use of new input data.. Any help will be much appreciated!


title: How can I (Should I?) use logistic regression/the logit function to predict outcome of a tennis match in a simple simulator?

text: I am trying to create a tennis simulator. Specifically I am trying to make a 'random' simulator so that I can see how many times streaks of wins or losses occur, and then compare this to historical data to see if 'The Hot Hand Effect' exists or not. i.e. Do streaks of wins or losses occur more or less often than would be expected by chance. I create a pool of players and assign them a strength score between 0 and 10 from a normal distribution of mean 5 standard deviation 2. When I come to matching up players and deciding the outcome, my supervisor suggested using a logit model/logistic regression. I am not quite sure how this would work practically/mathematically though. I thought to somehow use the difference in strength scores of the two players with the logit function to give me a probability to put into a cointoss, with the probability of a head, i.e. a win for Player A, being equal to the logit model output. 'Somehow' is obviously very vague as I do not know how this would work, or even if it mathematically makes any sense. Does what I intend to do sound valid conceptually? How do I use the strength scores with the logit function to get the probability?


title: Questions regarding Regularity Conditions in mathematical statistics

text: Often times, several theorems in mathematical statistics depend on regularity conditions. For example, there is a theorem that states: "Let $\hat{\theta}_n$ be the maximum likelihood estimator of $\theta$ based on the sample $X_n = (X_1,...,X_n)$, then under **regularity conditions** , $\hat{\theta}_n$ is a consistent estimator of $\theta$. Three of the regularity conditions I typically see are: 1) The set $\mathbb{X} = \{ x_1 \in \mathbb{R}: L_{x_1}(\theta)>0 \}$ doesn't depend on $\theta$. 2) If $L_{x_1}(\theta_1) = L_{x_1}(\theta_2)$ for almost all $x_1 \in \mathbb{X}$, then $\theta_1 = \theta_2$. 3) The parameter space $\phi$ of the unknown parameter $\theta$ is an open subset (though not necessarily a proper subset) of the real line. I do not understand 2) and 3). Would anyone be able to give me an intuitive explanation of number 3) and perhaps why we need 2) to hold? Thank you thank you!


title: Most effective SVM implementation for large classes and datasets?

text: What's the best SVM implementation of SVMs in MATLAB? Right now I'm looking at PEGASOS and SVM-Perf, which both claims to have good performance with large datasets and classes. But there are so many implementation out there, and I definitely missed some. Can I get some other recommendations? My application is in bags of word object recognition that uses the histogram of the words in classification.


title: Logtransforming pre and post data

text: I have pre and post data from a working memory test (following an intervention) that I am going to use in a mixed ANOVA. I have log transformed my post data because it was skewed, adding on 10 because of minus numbers, must I also log transform my pre data, even though it isn't skewed, and add on 10 in the process? Or should I just add on 10 to my pre data or not do anything to it and use it non-logtransformed? Thanks,any help much appreciated.


title: Kmeans plotting on discriminant components

text: When you plot a kmeans model (in R) with the plotcluster() function, it plots the clusters against the axis of the 1st and 2nd discriminant components (dc). In reading about these axis- some state the dc is not that much practical difference from the 1st and 2nd principal components. Other places seem to give a distinct difference. Any clarifications for me? Just want to understand what I am plotting against. Thanks


title: Simulating groups different with respect to the orthogonal complement in R

text: In a paper of Journal of Chemometrics (Naes & Mevik 2001 : Understanding the collinearity problem in regression and discriminant analysis), the authors propose to make simulations by creating two groups which are different with respect to the smallest eigenvector direction. > Here the groups are different with respect to the orthogonal complement to > the ﬁve ‘NIR loadings’. This is achieved in the following way. The constant > 0 ⋅ 18 is multiplied by a sixth loading vector (orthogonal to the other ﬁve) > and added to group 2. Both groups had initially the same means as group 1 How can I compute such a simulation in R? The goal is to obtain group differences which are tied to the "small eigenvectors" space. Then to check whether the LDA results are best when using the whole dataset or better by using only the last components as variables. * * * **Edit**   I am trying to get the same results as these authors by simulation, but unfortunately... I have tried to make two groups in a set of correlated data by adding differences between groups in the last eigen-vector (as it was suggested to me): Here is a try:               require(MASS)     R=matrix(0.9,10,10);     diag(R)=1;     random.normal=mvrnorm(n=20, rep(1,10), Sigma=R) ;     V <- var(random.normal) ## vcv matrix     U <- eigen(V)$vectors ## eigen vectors     Y <- random.normal %*% U     group=rep(1:2,each=10)     Y10 <- Y[,10] ## temporary data     Y10[group==1] <- Y10[group==1]+1 ## add 1 to group 1     Y10 <- Y10*sqrt(var(Y[,10])/var(Y10)) ## scaling to the variance of Y[,10]     Y[,10] <- Y10-mean(Y10)+mean(Y[,10]) ## center          X2 <- Y %*% t(U) #data back transformed          ## Then I try LDA with MASS          test=lda(X2,group,CV=T)     tab=table(test$class,group)     sum(diag(tab))/sum(tab)      But unfortunately I never found the bad classification results that Naes & Mevik 2001 found. when I use the complete dataset instead of using the last PC in the LDA as these authors found. (it's the same with more variables) I found that results are even better with the firsts components (e.g. 1:5 for ten PC)               test=lda(prcomp(X2)$x[,1:5],group,CV=T)         tab=table(test$class,group)     sum(diag(tab))/sum(tab)      Any suggestions to found the effects these authors explain? Thank a lot


title: Extract important features

text: Here is my situation: - A huge amount of data - 600 features - Only one class is provided Now, my question is how can I reduce the number of features to important ones? In another word, all of these features (with data) are intending to predict only one class. but some of features have large impact on the prediction (means their variation come to higher probability).


title: Difference between two mixed effects models

text: I have a question about how to tell two different mixed effects models apart. In the simple case both involve fitting a model with a random group effect and a covariate. I fit the model with `lme4` in `R`. Here is a visualization of the two different scenarios.   ![enter image description here](http://i.stack.imgur.com/GgxYw.png)               library(ggplot2)     library(lme4)     gen_dat2 <- function(group.m,group.v,int, sl,n){           x <- vector()           y <- vector()           g <- vector()              for(i in 1:length(group.m)){              x.t <- rnorm(n,group.m[i],group.v[i])              y.t <- rnorm(n,group.m[i],group.v[i])+int + sl*x.t               x <- c(x,x.t)              y <- c(y,y.t)              g <- c(g,rep(i,n))             }          return(cbind(x,y,g))     }          group.m <- runif(5,1,20)     group.v <- runif(5,1,1.5)          dat2 <- data.frame(gen_dat2(group.m,group.v,1,4,14))     ggplot(dat2,aes(x=x,y=y,colour=as.factor(g),group=g))+geom_point()+stat_smooth(method="lm",se=F)     m2 <- lmer(y~x + (x|g),data=dat2)      Then I can generate and fit the other scenario with similar code: ![enter image description here](http://i.stack.imgur.com/pgGm7.png)                gen_dat <- function(group.m,group.v,int, sl,n){           x <- vector()           y <- vector()           g <- vector()              for(i in 1:length(group.m)){              x.t <- rnorm(n,0,1)              y.t <- rnorm(n,group.m[i],group.v[i])+int + sl*x.t               x <- c(x,x.t)              y <- c(y,y.t)              g <- c(g,rep(i,n))             }          return(cbind(x,y,g))     }          group.m <- runif(5,1,20)     group.v <- runif(5,1,1.5)          dat1 <- data.frame(gen_dat(group.m,group.v,1,4,14))     ggplot(dat1,aes(x=x,y=y,colour=as.factor(g),group=g))+geom_point()+stat_smooth(method="lm",se=F)     m1 <- lmer(y~x + (x|g),data=dat1)      My central question is how do I tell these two models apart? Am I incorrectly fitting the first one, and I need an extra term in there to model the relationships between groups and the x variable as well as y? Both detect substantial between group variation in the intercept and not much in the slope as I would predict. But I need a way to tell these two apart. Any thoughts would be helpful. * * * Edits: This has been helpful in me restating the question. So I want to re-ask the question with an example which I hope will make it clear why I want to be able to tell these two models apart. Let's imagine that Y is the average student test score at a school, and X is spending per student in that school. Our grouping variables are 5 different school districts. Data in the top figure shows that an increase in spending within a district means that test scores increase. It also shows that between districts there are differences is scores, but that's clearly because some districts spend more student than others. Data in the second figure show similarly that within a district student scores increase as spending increases. It also shows that between districts there are differences in test scores. However we don't know what is driving those differences, unlike in the first set of data. This is a pretty common situation I've encountered in building models. The former is not. So what I'm asking is what is the appropriate model that captures the following features from the first dataset:   1. Test scores increase as spending per student does   2. There is also variance between districts in student test scores   3. Part of that difference between districts is because of the underlying relationship between spending and test scores, which also varies with district. More generally stated, how do you handle a scenario where you're building a hierarchical model where the grouping variable is correlated with one of your continuous independent variables (e.g. the first scenario). I feel like the model I've presented get's at points 1. and 2., but not point 3. So I'm really seeking a way to tease these two scenarios apart. Normally I might add an extra level of hierarchy if there was another group level explanatory variable. Continuing our example, maybe in the 2nd dataset there are differences between district because in some districts parents have more time to spend on homework with students. So we would add that as a group level predictor in a hierarchical model. But that wouldn't work in the first scenario.


title: Finding random vector from its covariance matrix

text: I'm having difficulty solving the following problem: > Given the covariance matrix Rxx = [1 0.2 ;0.2 1] for a zero-mean random > vector X, find two different transformations of the form Y = AX such that > the components of Y are uncorrelated. I suppose X cannot be derived back from Rxx by matrix decomposition or any other means since we are dealing with expectations. If this is the case, is there any way to find the required 'Y'?


title: Statistically group chemical batches

text: I have 80 chemical batches, each of which has 8 associated measurements like pH, viscosity, etc. Is there a way to use Principal Components Analysis or Factor Analysis to group the similar batches together? I know I can use cluster analysis, but I want to have a different visual representation of batch groups other than the dendogram. I'd like to have a bi-plot. I tried to perform a PCA, but got an error message that there cannot be more variables than observations. My matrix has the 80 batches as the column names, and the 8 measurements as the row names. When this matrix is transposed, PCA just groups the like measurements, not the batches.


title: Probability of at least one unique outcome

text: We have n actors. Each actor chooses from n*2 actions. How can I calculate the probability that at least one actor will choose a unique outcome? For example, say we have 5 pickup artists in a town with 10 bars. Each PUA chooses a bar at random. What are the odds that at least one PUA will have a bar to himself/herself? This is relevant to some work I'm doing on scheduling in a distributed system. I know that the relevant equation should probably start with $\binom{n + 2n -1}{n}$, but then what?


title: How do I check for variation of a distribution over spatial scales?

text: I have about 30000 particles distributed (not randomly) in space; I have their vector positions and velocities. I'm trying to characterize if and how the (underlying) velocity (magnitude) distribution these particles are assumed to sample changes on different scales. The underlying velocity distribution is of unknown form. Note: My real interest is in the distribution on a scale that I have near zero particles; so I'd like to be able to get a sense of how the scale effects the distributions generally. I can of course look at distributions made from subsets of nearby points; to my eye these look pretty similar near the edges and very discrepant in the centers. I'd like a way to quantify and visualize this change across scales. I'm fairly new to statistics so any/all references are welcome. Feel free to suggest/add tags.


title: Large sample size and partial F-test for multiple regression makes adding a variable always significant

text: I am developing a two variable multiple regression model. ie. $$ Y = b0 - b1 * X1 + b2 * X2 $$ I am using the following formula for partial F-test from http://luna.cas.usf.edu/~mbrannic/files/regression/Reg2IV.html under the section Testing Incremental R2. The F-statistics calculated is supposed to tell me if adding the second variable is significant (more details in that link). $$ F= {\frac{(R_L^2 - R_S^2)/(k_L-k_s)}{(1-R_L^2)/(N-k_L-1)}}$$ My first variable has a strong correlation: regression_coeff_string: b1 = 0.664, b0 = 0.035 R2_val: 0.564 My second variable has a weak correlation: regression_coeff_string: b1 = -25.026, b0 = 0.469, R2_val: 0.027 Adding my seond variable only marginally improves the R2 value regression_coeff_string: b0 = 0.0559, b1 = 0.6633, b2 = -5.2222, R2_val: 0.565 However, because I have a sample size 2949, that With $$ R_L^2 = 0.565, R_S^2 = 0.564$$ $k_L$ the number of predictors in the full set being 2, $k_S$ the number of predictors in the subset being 1 $$ F= {\frac{(0.565 - 0.564)/(2-1)}{(1-0.565)/(2949-2-1)}} = 6.77$$ With F(1,2946) at 0.05 confidence having a F_stat of 4.182, the result is significant. But it seems that it is only because the sample size is large. If I sort the second variable X2 in ascending order in Excel and leave the order of the Y and X1 variables unchanged, I would still get a significant F score. **Question:** How can I do a fair incremental R2 test for the addition of a new variable in multiple regression when the sample size becomes large? Simply looking at the R2 of each X variable individually does not take into account that that they may be cross-correlated, that is why I turned to the incremental R2 test to see how the overall R2 improves relative to adding a new variable. **EDIT1:** The context of my example is predicting solar radiation. The first variable is a solar radiation variable from NWP (numerical weather prediction) software (ie. high correlation). The other variables are other NWP output variables and we are trying to improve our prediction.


title: How to interpret Hidden Markov Model parameters (transition matrix, emission matrix, and pi values)?

text: I am working on channel modeling for cognitive radio using HMM. I've written a MATLAB program for forward, backward and Baum-Welch algorithm for multiple sequences. After given some random input and running the program for 1000 to 4000 iterations I'm getting some results. But I'm not getting how to interpret my results. Will be glad if anyone can talk about inputting matrix/log- likelihood/transition matrix/emission matrix.


title: When to terminate the Bayesian A/B test?

text: I'm trying to do A/B testing the Bayesian way, as in Probabilistic Programming for Hackers and Bayesian A/B tests. Both articles assume that the decision maker decides which of the variants is better based solely on the probability of some criterion, e.g. $P(p_A > p_B) = 0.97 $, therefore, $A$ is better. This probability doesn't provide any information on whether there was sufficient amount of data to draw any conclusions from it. So, it is unclear to me, when to stop the test. Suppose there are two binary RVs, $A$ and $B$, and I want to estimate how likely it is that $ p_A > p_B $, and $ \frac{p_A - p_B}{p_A} > 5\% $ based on the observations of $A$ and $B$. Additionally, suppose that $p_A$ and $p_B$ posteriors are beta-distributed. Since I can find the $\alpha, \beta$ parameters for $p_A\,|\,\text{data} $ and $p_B\,|\,\text{data} $, I can sample the posteriors, and estimate $P(p_A > p_B\ |\ \text{data})$. Example in python:               import numpy as np          samples = {'A': np.random.beta(alpha1, beta1, 1000),                'B': np.random.beta(alpha2, beta2, 1000)}     p = np.mean(samples['A'] > samples['B'])      I could get, for example, $P(p_A > p_B) = 0.95$. Now I would want to have something like $P(p_A > p_B\ |\ \text{data}) = 0.95 \pm 0.03$. I have researched about credible intervals and Bayes factors, but can't understand how to calculate them for this case if they are applicable at all. How can I calculate these additional stats so that I'd have a good termination criterion?


title: Asymptotic Theory in Economics

text: I am interested in deepening my Asymptotic Theory understanding. My current knowledge is that of a typical PhD student (from a decent University), say at the level of Green's textbook. Are there any good book(s) that you would recommend? Much appreciated.


title: Correlational study or ordinal data using 5-point Likert scale

text: How do I calculate the correlation using ordinal data based on a 5-point Likert scale associating perioperative education to patient satisfaction scores? A numerical value (`1`: strongly agree – `5`: strongly disagree) will represent the participant's perception of satisfaction as a patient as related to perioperative educational material. The numerical value is the patient's response to 5 specific statements. The Likert scale is being used during a telephone survey to gather patient responses to 5 statements regarding educational material provided during his or her surgical experience. An example of the statement may appear as, "The day surgery nurse provided clear and easy to understand verbal instruction regarding personal care once at home." The goal is to correlate patient satisfaction with the surgical experience to patient education received by the surgical team of nurses. Question: What is the best way to analyze this data? Would Spearman's rank correlation be appropriate?


title: Random model selection and validity of significance tests

text: Suppose we have some data, $\{y_i, x_{1i}, \dots, x_{ki}\}_{i=1}^n$ and we want to build a linear model of the form $y_i = \beta_0 + \beta_{1'i}x_{1'i} + \dots + \beta_{k'i}x_{k'i} + \epsilon_i$, where the primed indexes are a subset of $\{1,\dots, k\}$. There are $2^k$ possible models to choose from. Now if we randomly (uniformly) draw a model from all $2^k$ , and fit it to the data, once, can we still interpret the standard errors and CIs as usual?


title: Explaining Gaussian Processes

text: I am finding it hard to understand Gaussian Processes. Can someone please explain it here in an accessible way? I do understand what Gaussian distribution is but couldn't understand Gaussian Processes. Or if you could provide some great resources to quickly learn about it.


title: Correlation in Time Series

text: I have two time series $P_i$ and $T_i$ where $i = 0...t$ Now the correlation between the time series ${\rm corr}(P,T)\sim.95$, but the correlation between the derivative of the time series ${\rm corr}(dP, dT)\sim.5$ where $dP = P_i - P_{i-1}$ for all $i$. I'm trying to interpret this but can't come up with any explanations. How can the time series be correlated but the running differences not be?


title: SVD application for a Boolean sparse Matrix

text: Basically, I am trying to have a recommender system based on SVD for a Boolean utility matrix. ie If at all some entries are present in the utility matrix, they will be 1 (I made it pseudo-implicit from explicit ratings available from the Movielens dataset). I used the 'irlba' package in R to compute the partial SVD quickly. My confusion is, whether the SVD is a good technique for boolean user-item preference matrix. I read online that SVD doesn't offer that many 'interpretable' latent factors/features for a Boolean sparse utility matrix. So, will it lead to some kind of an over fitting if the rank-k (the dimension of the latent factor space) is chosen as some value beyond say 20? I used the 10 Million ratings Movielens Data and modified them as 5 Million sparse 1 entries (for ratings= 4 or 5). The problem is, the top-N list from this model for many users do not match at all (at the maximum of 1 or 2 intersections) with the collaborative-user-neighborhood similarity based recommender system. Could someone please tell whether SVD can be applied to my situation ( I need to make the recommendations on boolean sparse utility matrix) and if yes, how should I modify my present approach?


title: How to pick from a set of random variables the one with the highest mean using a fixed number measurements?

text: Suppose we have $N$ approximately normally distributed continuous random variables $X_1, X_2, X_3, \ldots,X_N$, each with an unknown mean and variance. I'd like to find the random variable with the highest mean. Given a fixed budget of measurements $K$, how do I assign these measurements to maximize my chances of correctly picking the $X$ with the highest true mean? Naively one could measure each $X, K/N$ times, then choose the one with the highest sample mean. Intuitively we ought to be able to do better by assigning less measurements to $X$'s with lower sample variance and/or lower sample mean (as they're less likely to be the one with the highest true mean). **Is there an algorithm to produce an optimal sequence of measurement choices?**


title: Fisher classification method for normally distributed classes

text: As a follow-up to my previous question, I want to use Fisher's classification method (I mean, projection method) to project the data given the fact that the two classes are normally distributed. How to approximate the distributions? How to project the data when it's already in a one dimensional space? I remember Fisher's idea (two classes) was to project the data from D-dimensional input space to one dimensional. I don't know why I can't work with real data examples when I know how the actual method works. Please someone guide me through this. The classes are: $C_1=\{3, 3.5, 4, 4.5, 5, 5.5, 7\} \cup \{15, 16, 17\}$ $C_2=\{0,0.5,1,2\}\cup\{8,9,10,11,12,13\}\cup\{20,25,30\}$ And the data will be on the interval $[0,40]$ like in the previous question.


title: Checking homogeneity in an insolation time series

text: I have insolation data for several years on which I want to do some calculations. First, I'd like to determine the level of reliability I can achieve in my posterior calculations. For several references, my first option was use Thom non-parametric test (1966), suggested for data without known distribution, or which is not normally distributed (World Meteorological Organization - http://www.wmo.int/pages/index_en.html - advise). However, I can't get the complete description of the test, and, therefore, I can't implement it correctly in R. Next, I'd like to run tests from the lawstat and tseries packages (the same test I think), but I believe the results only indicate that the data are not aleatory: ------- Console cut ------- Palme.prueba.rachas.lawstat<-lawstat::runs.test(Palme.rad.lintem)               Runs Test - Two sided      data: Palme.rad.lintem Standardized Runs Statistic = -12.1657, p-value < 2.2e-16 ------ End of console cut ------ My question is, how can I check the veracity of the time series?


title: Transform data to a higher dimensional space

text: The classes are as follows: $C_1=\{3, 3.5, 4, 4.5, 5, 5.5, 7\} \cup \{15, 16, 17\}$ $C_2=\{0,0.5,1,2\}\cup\{8,9,10,11,12,13\}\cup\{20,25,30\}$ And we wish to classify numbers from the interval $[0,40]$. I think since they are not linearly separable (the two classes in a one dimensional input space) so they must be transformed to a higher dimensional space where they are linearly separable. Will that idea go anywhere? If yes, how should I do it? I want to use SVM as a classification method to classify a given number on that interval but I don't seem to understand what has to be done in this case even though I am well aware of the theoretical idea behind SVM.


title: Could standardizing an independent variable cause the t-statistic of the OLS estimate to change?

text: Lets say we are looking to estimate the following standard OLS regression: $y_{i} = \beta_{0} + \beta_{1}*X_{i} + \beta_{2}*Z_{i} + \epsilon_{i}$ and that we choose to standardize $X$ as: $X^{*}=((X_{i} - \mu)/ \sigma)$. Suppose we now instead estimate the following regression: $y_{i} = \beta_{0} + \beta^{*}_{1}*X^{*}_{i} + \beta_{2}*Z_{i} + \epsilon_{i}$ While I know that the point estimates for $\beta_{1}$ and $\beta^{*}_{1}$ will be different, I am less certain of how the significance of the point estimate will be affected. My guess is that a simple transformation of an independent variable (or standardization) would not affect the significance of the point estimate (i.e. it would not affect the t-statistic for its $\beta$ estimate)--but I am not sure of this... **My Question:** Is my guess correct? An explanation of why this guess might be right/wrong would be much appreciated. Thanks! PS: Running some simple simulations verified my guess above, but a simple simulation is obviously no proof, so if anyone could provide one that would be great.


title: Why use the parametric bootstrap?

text: I am currently trying to get my head around some things concerning parametric bootstrap. Most things are probably trivial but i still think i may have missed something. Suppose i want to get confidence intervals for data using a parametric bootstrap procedure. So i have this sample and i assume its normally distributed. I would then estimate variance vhat and mean mhat and get my distribution estimate Phat, which is obviously just N(mhat,vhat). Instead of sampling from that distribution i could just calculate the quantiles analytically and be done. a) i conclude: in this trivial case, the parametric bootstrap would be the same as calculating things in a normal-distribution-assumption? So theoretically this would be the case for all parametric bootstrap models, as long as i can handle the calculations. b) i conclude: using the assumption of a certain distribution will bring me extra accuracy in the parametric bootstrap over the nonparametric one (if it is correct of course). but other than that, i just do it because i cant handle the analytic calculations und try to simulate my way out of it? c) i would also use it if the calculations are "usually" done using some approximation because this would perhaps give me more accuracy...? to me the benefit of the (nonparametric) bootstrap seemed to lie in the fact that i dont need to assume any distribution. for the parametric bootstrap that advantage is gone - or are there things ive missed and where the parametric bootstrap provides a benefit over the things mentioned above? thank you very much.


title: Compare averaged GLM with boosted regression trees using cross validation : d2 and RMSE calculation

text: I want to compare BRT and averaged glm models on test sets by calculating the explained deviance and RMSE. How can I calculate d2 and RMSE from predictions?   I use the following functions: > gbm1 <- gbm.step(data=ahe2, gbm.x = c(2,4:5,7:9,11),#15 gbm.y = 10, family = > "poisson", tree.complexity = k, nfolds=10, bag.fraction = 0.50, n.trees=nb, > learning.rate=z, tolerance.method = "fixed", tolerance = 0.01) # GLM poisson > glm1=glm(larve~chla+station+ vit+direct+lune+ > temp,data=ahe2,family="poisson") > > dd=dredge(glm1,rank="QAICc",chat=c_hat(glm1)) > > models.list <- get.models(dd,subset = cumsum(weight) <= .95) > > muminavg.quasi <- model.avg(models.list) # I split data: > kfold1=createFolds(ahe2$larve, k = 10, list = T) After that, I will use predict() on the prediction set. I thought to calculate RMSE with the square root of the mean squared residuals. For deviance, and particularly for gbm1, I don’t know how to proceed (t(residuals(glm1))%*%residuals(glm1) ?). I also wanted to implement a averaged gam but don’t know to plot it.


title: split file in weka

text: I want to run a logistic regression were the file is split on a certain value, namely key (major and minor) and I want to see if the other values are able to predict the ranking of a song. But I have no clue how to split the file and run the analysis. So please help me. Kind regards, Dionne


title: What is the most effective way to combine different rating scales for the same data?

text: I have a collection of data which is inherently the same (would expect it to follow the same distribution) where some of it is rated with one scale of ordinal ratings and the rest of the data is rated with another scale of ordinal ratings. One scale has a larger number of possible ratings than the other. So, if I were to map from one set of ratings to the other by some sort of function (linear transformation for instance) I would get in between ratings and have to round on the output scale to get a legal discrete value. Is there a rigorously justifiable method or at least a heavily used method which allows me to analyze all of the data on a single scale whether it be the larger scale, the smaller scale, or a third new scale which is neither of the first two scales (I was thinking maybe LCM of the size of the two scales if it is not too big)?


title: Compute the user and item features in SVD++

text: I have a sparse matrix. There is lots of missing data. Hence, I can't use SVD naively. I read Koren's SVD++ paper. I'm confused as to how the $q_i$ and $p_u$ vectors are determined. $q_i^Tp_u$ is supposed to capture the interactions between user $u$ and item $i$ (plus some biases). I just don't see how to calculate what $q_i$ and $p_u$ are supposed to be. The most natural things would be to use SVD but you can't since it has missing data.


title: clustering versus projection - what are the best example/scenario to explain their differences

text: I am dealing with non statistician who know are crunching data but don't have a deep understanding of statistics. I am trying to introduce them to non-matrix factorization methods but it has been hard. Since, people avoid leaving their comfort zone ! In my case, my co- workers have been working with hierarchical clustering methods or k-means, A part of the challenge is also, I don't know the answer 100% percent by myself. The way I see it, dimensionality reduction (Non-negative matrix factorization) is distinct from clustering but it can help clustering high- dimensional data by reducing to a "useful" low dimensional representation. I also know, k-means is a relax version of NMF, where W are centers of the cluster and H matrix is a membership matrix. I also know, NMF (non-negative matrix factorization) is a projection methods, where each component is a linear combination of different features; and unlike clustering methods where tries to find the clusters by finding features that maximize the variance; NMF find clusters based on these components. Still, I would like to have more argument to show, NMF is different from clustering methods; and also possibly because of X and Y has more advantages. I appreciate any comment that can help me to understand and explain better the differences between NMF and clustering approaches.


title: using stepAIC of MASS package to select variables with a significance level of 5% in R project

text: First of all, sorry i am new about this and any helps are really welcome. I am reading a reaserch paper where the authors report: _Stepwise forward regression (Zar 1996) was used to select the most informative variables, which were included in a multiple (linear) regression model. A 5% significance level was chosen as a threshold for the inclusion of the model variables._ with a private email the first author told me that the variable selection was performed using stepAIC of MASS library using direction "forward" and they considered only for the final model the variables with a significance level of < 5%. using junk data i tried to rewrite the analysis in order to understand the procedure               state.x77     st = as.data.frame(state.x77) str(st) colnames(st)[4] = "Life.Exp"     colnames(st)[6] = "HS.Grad" st[,9] = st$Population * 1000 / st$Area colnames(st)[9] = "Density"     str(st) model1 = lm(Life.Exp ~ Population + Income + Illiteracy + Murder + + HS.Grad + Frost + Area + Density, data=st)     model1.stepAIC <- stepAIC(model1, direction=c("both"))     summary(model1.stepAIC)               Call:     lm(formula = Life.Exp ~ Population + Murder + HS.Grad + Frost,          data = st)          Residuals:          Min       1Q   Median       3Q      Max      -1.47095 -0.53464 -0.03701  0.57621  1.50683           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  7.103e+01  9.529e-01  74.542  < 2e-16 ***     Population   5.014e-05  2.512e-05   1.996  0.05201 .       Murder      -3.001e-01  3.661e-02  -8.199 1.77e-10 ***     HS.Grad      4.658e-02  1.483e-02   3.142  0.00297 **      Frost       -5.943e-03  2.421e-03  -2.455  0.01802 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.7197 on 45 degrees of freedom     Multiple R-squared:  0.736,     Adjusted R-squared:  0.7126      F-statistic: 31.37 on 4 and 45 DF,  p-value: 1.696e-12      followint the protocol of the paper the final model is               Life.Exp ~ Murder + HS.Grad + Frost (final model)       because Population is > 0.05. I wish to know if this final model approach is correct, and then:               fmodel = lm(Life.Exp ~ Murder + HS.Grad + Frost, data=st)     summary(fmodel)          Call:     lm(formula = Life.Exp ~ Murder + HS.Grad + Frost, data = st)          Residuals:         Min      1Q  Median      3Q     Max      -1.5015 -0.5391  0.1014  0.5921  1.2268           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) 71.036379   0.983262  72.246  < 2e-16 ***     Murder      -0.283065   0.036731  -7.706 8.04e-10 ***     HS.Grad      0.049949   0.015201   3.286  0.00195 **      Frost       -0.006912   0.002447  -2.824  0.00699 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.7427 on 46 degrees of freedom     Multiple R-squared:  0.7127,    Adjusted R-squared:  0.6939      F-statistic: 38.03 on 3 and 46 DF,  p-value: 1.634e-12


title: Best way to test runs

text: I have the following sample:               x <- c(rep(1, 8), 0, 1, 0, rep(1, 5), 0, 0, 1, rep(0, 11), 1, 0, 0, 1, 0)     > x     [1] 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0      The `1` and `0` are actually factors of a nominal variable and the sample is sorted according to a hidden variable. It looks like the sequence is not random, that the 1s tend to rank lower (i.e. appear earlier) in this sequence. I would like, however, to conduct a statistical test to test this hypothesis. So far, my solution was to use the Wald-Wolfowitz test, like so:               library(adehabitat)     wawotest(x)      It yields a p-value of 0.0211, which I consider low enough to reject the hypothesis that the ones and zeroes occur randomly in the given run. However, I've read that the Kolmogorov-Smirnov test is a more powerful solution and I'd like to try it in R, but I'm unsure I'm giving it the right command. This is what I'm running:               ks.test(x, 'punif')      Which gives a p-value of 1,821e-08. Is my command correct? Should I consider any other tests? **[EDIT]** As @Glen_b correctly put it, the KS is not suited for this case, so I'll be sticking with the WW. Here are the two best ways I've found for running this test:               > library(tseries);    runs.test(as.factor(x))     > library(adehabitat); wawotest(x)      I've studied the algorithm behind `runs.test()` and it seems to be conducting the Wald-Wolfowitz, even though the p-values are slightly different. Anyway, WW doesn't look like the most powerful tool out there for this kind of test, so feel free to provide other solutions.


title: Unable to provide random parameter with mlogit

text: I'm having trouble with rpar argument in the mlogit function (package mlogit). My dataset looks like this:               > head(scan.s)        year id.scan day weather sealvl wave lact repro     nb.gr       HS nb.pv pup act     1 2011       1   4       2   0.30    3    1     0 0.6666667 7.600000     3   0   R     2 2011       1   4       2   0.30    3    1     0 0.6666667 7.600000     3   0   R     3 2011       1   4       2   0.30    3    1     0 0.6666667 7.600000     3   1   R     4 2011       2   4       2   0.35    3    1     0 0.6666667 8.100000     2   0   R     5 2011       2   4       2   0.35    3    1     0 0.6666667 8.100000     2   1   R     6 2011       3   4       2   0.40    3    1     0 0.6666667 8.633333     2   0   R          > str(scan.s)     'data.frame':   10140 obs. of  13 variables:      $ year   : int  2011 2011 2011 2011 2011 2011 2011 2011 2011 2011 ...      $ id.scan: Factor w/ 280 levels "1","2","3","4",..: 1 1 1 2 2 3 3 4 4 5 ...      $ day    : int  4 4 4 4 4 4 4 4 4 4 ...      $ weather: Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 2 2 2 ...      $ sealvl : num  0.3 0.3 0.3 0.35 0.35 0.4 0.4 0.5 0.5 0.6 ...      $ wave   : Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 3 3 3 3 3 ...      $ lact   : int  1 1 1 1 1 1 1 1 1 1 ...      $ repro  : int  0 0 0 0 0 0 0 0 0 0 ...      $ nb.gr  : num  0.667 0.667 0.667 0.667 0.667 ...      $ HS     : num  7.6 7.6 7.6 8.1 8.1 ...      $ nb.pv  : int  3 3 3 2 2 2 2 2 2 2 ...      $ pup    : int  0 0 1 0 1 0 1 0 1 0 ...      $ act    : Factor w/ 5 levels "A","C","D","G",..: 5 5 5 5 5 5 5 5 5 5 ...      Then I used mlogit.data to transform my dataset in long shape:               > scan.l<- mlogit.data(scan, varying = NULL, choice = "act", shape = "wide")      There is no variable varying across choices.                   year id.scan day weather sealvl wave lact repro     nb.gr  HS nb.pv pup   act chid alt     1.A 2011       1   4       2    0.3    3    1     0 0.6666667 7.6     3   0 FALSE    1   A     1.C 2011       1   4       2    0.3    3    1     0 0.6666667 7.6     3   0 FALSE    1   C     1.D 2011       1   4       2    0.3    3    1     0 0.6666667 7.6     3   0 FALSE    1   D     1.G 2011       1   4       2    0.3    3    1     0 0.6666667 7.6     3   0 FALSE    1   G     1.R 2011       1   4       2    0.3    3    1     0 0.6666667 7.6     3   0  TRUE    1   R     2.A 2011       1   4       2    0.3    3    1     0 0.6666667 7.6     3   0 FALSE    2   A          > str(scan.l)     Classes ‘mlogit.data’ and 'data.frame': 50700 obs. of  15 variables:      $ year   : int  2011 2011 2011 2011 2011 2011 2011 2011 2011 2011 ...      $ id.scan: Factor w/ 280 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...      $ day    : int  4 4 4 4 4 4 4 4 4 4 ...      $ weather: Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 2 2 2 ...      $ sealvl : num  0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...      $ wave   : Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 3 3 3 3 3 ...      $ lact   : int  1 1 1 1 1 1 1 1 1 1 ...      $ repro  : int  0 0 0 0 0 0 0 0 0 0 ...      $ nb.gr  : num  0.667 0.667 0.667 0.667 0.667 ...      $ HS     : num  7.6 7.6 7.6 7.6 7.6 7.6 7.6 7.6 7.6 7.6 ...      $ nb.pv  : int  3 3 3 3 3 3 3 3 3 3 ...      $ pup    : int  0 0 0 0 0 0 0 0 0 0 ...      $ act    : logi  FALSE FALSE FALSE FALSE TRUE FALSE ...      $ chid   : num  1 1 1 1 1 2 2 2 2 2 ...      $ alt    : chr  "A" "C" "D" "G" ...      - attr(*, "index")='data.frame': 50700 obs. of  2 variables:       ..$ chid: Factor w/ 10140 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...       ..$ alt : Factor w/ 5 levels "A","C","D","G",..: 1 2 3 4 5 1 2 3 4 5 ...      - attr(*, "choice")= chr "act"      Then I ran the model:               mod1 <- mlogit(act ~ 1| nb.gr+nb.pv+sealvl+lact+repro+HS+day+id.scan,data = na.omit(scan.l), rpar=id.scan, format="long", reflevel="R", R=100, halton=NA, print.level=0)      The random parameter here is a factor and I am supposed to specify a distribution for rpar but is it relevant for a factor? (I tried to provide a distribution without any change). And then I get this:               Error in coef(eval(callst, parent.frame())) :        error in evaluating the argument 'object' in selecting a method for function'coef' : Error in solve.default(H, g[!fixed]) : Lapack routine dgesv: system is exactly singular      There is a way to use "HS" and "day" instead, both numerical. But then I get another error:               Error in names (sup.coef) <- names.sup.coef: Attribute 'names' [1] must be the same length as the vector [0]       traceback() did not provide any insight about what happened. I searched for explanations with those errors and found that there could be a problem between one outcome and the random effect so I tried to subset my dataset with every combination of 3 outcomes with the same result. I found nothing relevant about the second error. Maybe it has something to do with the transformation with mlogit.data. I checked the dataset provided with the mlogit package and could not figure out what I did different. I would be grateful if someone could explain what is happening here. Thanks!


title: Distribution of training data for binary classification

text: I'm doing a binary classification using K-nearest neighbor method. However when I look at the training data, 80% of the data is from category 2, and only 20% is from category 1. Is it bad when you try to learn from that data? I feel like most of the prediction go to category 2.


title: Is maximum likelihood a method for estimating parameters or estimating a distribution?

text: I am narrowing down on some more solid idea of what MLE is, and what we use it for. I keep running across different perspectives from different sources. Some places, I read that maximum likelihood is a method for estimating parameters of a distribution, given we have some data. Other places, I read that the maximum likelihood is itself a distribution. What do we get from this process? A set of parameters for a distribution we already are trying to fit? Or do we get a 'new' distribution at the end of the process? Or both? Thanks and I apologize for such a basic question.


title: Measurement error in nonlinear models, when can one ignore it? and what can one do about it?

text: Suppose you know there is measurement error in your system but you don't have any idea what the variance is. In my specific case I have a nonlinear model \begin{align} y_i&=\epsilon_i f(\beta,x_i) \end{align} where $\epsilon_i$ is log-normally distributed with mean 1, and \begin{equation} f(x,\beta)=\frac{\beta_1x}{1+\beta_2x} \end{equation} But we know $x$ and $y$ are measured with some error, which we know to be of the form \begin{align} x_i&=\eta_i x_i^*\\ y_i&=\eta_i y_i^*\\ y^*_i&=\epsilon_i f(\beta,x_i^*) \end{align} where $x^*_i$ and $y^*_i$ are the true variables and $\eta_i$ and $\zeta_i$ are iid lognormal with mean 1 and unknown variance (they have the same variance but not the same as $\epsilon_i$). Currently I have been completely ignoring the measurement error and using nls on the log transformed data. I am only concerned about the parameter estimates and not confidence intervals. $x$ and $y$ should be at least 4 or so orders of magnitude larger than the variance in either $\epsilon$ or $\eta_i$ and $\zeta_i$. I recognize that there is two forms of bias here. First the $\log$ transformation means that the error term in \begin{align} \log(y_i)&=\log(f(\beta,x_i)) + e_i \end{align} does not have mean zero, which creates bias (but possibly negligible because of the 4 or so orders of magnitude condition above), plus there is probably bias from the measurement error in general. My question is (1) how much can this bias affect the parameter estimates and (2) what can I do if this bias is large? My idea was to generate data under this model via simulation under different variances in the measurement error and then see how much it affected the parameter estimates. Is there something better I can do?


title: Random classification forests for extremely sparse response variables

text: I have a response variable that can be $A,B,C$. It is very sparse, meaning 99% of the sample is $B$ and the rest is approximately evenly divided between $A$ and $C$. How do I predict this variable in a random classification forest? I am looking for guidelines:   * Can I use the standard classification splitting criterion with such a sparse response variable?   * Given the asymmetric damage an out of sample misclassification would do (i.e. classifying A or C correctly is most important and B correctly is a lower priority), how do I apply some kind of asymmetric loss function here?    * Are there other special things I need to take into consideration when modelling such a sparse response variable? Related but not duplicated: Random forest implementation for sparse data


title: Achieving a perfect fit using libsvm classification

text: I need to use SVM classification to build a model that predicts a binary code My given data is in the form: `(feature vector, associated bit )`, and contains 5000, 200-dimensional feature vectors. I have tried to use a `libsvm` solution with its associated data file heart_scale, but I would like to try to do this training task perfectly as accuracy is critical in my project. **EDIT** My data ={feature vector , y) where the feature vector is real-valued 265-dimensional vector represent low level visual feature for an image and y ={1,-1}


title: Is calculating "actual coverage probability" the same thing as calculating a "credible interval"?

text: I was reading an entry level statistics textbook. In the chapter on maximum likelihood estimation of the success proportion in data with binomial distribution, it gave a formula for calculating a confidence interval and then nonchalantly mentioned > Consider its actual coverage probability, that is, the probability that the > method produces an interval that captures the true parameter value. This may > be quite a bit less than the nominal value. And goes on with a suggestion to construct an alternative "confidence interval", which presumably contains the actual coverage probability. I was confronted with the idea of nominal and actual coverage probability for the first time. Making my way through old questions here, I think I got an understanding for it: there are two different concepts we call probability, the first being how probable it is that a not-yet-happened event will produce a given result, and the second is how probable is that an observing agent's guess for the result of an already-happened-event is true. It also seemed that confidence intervals only measure the first type of probability, and that something called "credible intervals" measure the second type of probability. I summarily assumed that confidence intervals are the ones which calculate "nominal coverage probability" and credible intervals are the ones which cover "actual coverage probability". But maybe I have misinterpreted the book (it is not entirely clear on whether the different calculation methods it offers are for a confidence interval and a credible interval, or for two different types of confidence interval), or the other sources I used to come to my current understanding. Especially a comment which I got on another question, > Confidence intervals for frequentist, credible for Bayesian made me doubt my conclusions, as the book did not describe a Bayesian method in that chapter. So please clarify if my understanding is correct, or if I have made a logical error on the way.


title: How to model features with NULL (not necessarily missing) values

text: I’d like some advice on my options for modeling features (to be used in R) in the following scenario: I want to make a binary prediction using the positions of the planets as the predictors (the X’s). My issue is how best to model the planet positions to allow them to be used by functions such as randomForest and SVM. The problem is with missing, or more appropriately NULL, values for a feature. As X’s, **for each day** , I want to capture:   1. Whether the moon is new or full (possibly 1 or 2 features – 1 if value = N|F, 2 if each is a Y|N value))   2. Whether a planet is changing signs (e.g., going into Aries – the same as crossing a 30 degree multiple) (one feature per planet, where value = Y|N if planet is changing signs today))   3. Whether there is an important aspect (angular division) between any pair of planets. An aspect is important if it’s 0, 60, 120, 180 degrees. (N(N-1)/2 features where N = 9 … 36 features). The feature will be a factor with likely values: 0, 60, 120, 180 **My big issue is with #3** and how to represent a day when there is no interesting (0,60,120,180) angular difference between a given pair of planets. For example, consider Mercury and Mars. On a day when the angular difference between the Mercury and Mars is 47.52 degrees I really don’t want to record anything, i.e., the value of feature Mercury.Mars would be empty/NULL as 47.52 is not considered important and I don’t want to use it as input to the prediction function (randomForest, SVM, …, etc). The value on such a day is not really missing in the sense that it’s not known or it’s been left out, there just isn’t a value for that feature for that day. Most days will not have an interesting aspect value. So, how best to model this situation? My concerns:   1. If I don’t record a value (i.e., it’s really missing) then randomForest won’t work as it can’t handle missing values. Not sure about SVM and other functions.   2. If I put some generic value to indicate “nothing here today”, for example *, then I’m concerned that the prediction function (rf, SVM) may associate some uncontrolled-by-me meaning to the “nothing here today” value Your thoughts? Thank You.


title: R procedure for comparing multiple categorical variables (similar to anova() followed by t.test() for continuous)?

text: # Big Picture: How can I implement partitioned Chi Square in R? I understand how to perform the overall Chi square, and then how to get individual parameters (observed counts, expected counts, residuals, etc.). However, I don't understand how to get p values for each individual comparison. # Details:               set.seed(200)     alpha = 0.05      I have three groups of patients and I have both categorical and continuous data collected on those patients.               Group <- c(rep('A', 'B', 'C'), 10))     Mass <- c(rnorm(10, mean = 60, sd = 1),               rnorm(10, mean = 70, sd = 1),               rnorm(10, mean = 80, sd = 1))     Sex <- rep(c('Male', 'Female'), 15)     data <- data.frame(Group, Mass, Sex)      I'd like to construct your basic "Table 1" showing whether these groups differ from one another with respect to these variables. For the continuous variables, I would do an anova test. If p were less than alpha, I would do individual t.tests and interpret my p values using a Bonferroni correction.               print('Anova p-value:')     print(anova(lm(Mass ~ Group, data))$'Pr(>F)'[1])          print(t.test(Mass ~ Group, data[data$Group != 'A',])$p.value)     print(t.test(Mass ~ Group, data[data$Group != 'B',])$p.value)     print(t.test(Mass ~ Group, data[data$Group != 'C',])$p.value)                          I'd like an analogous output for my discrete variables. I know that I can do a Chi square test on the data as a whole...               print(chisq.test(x = data$Group, y = data$Sex)[['p.value']])      ...but I'm confused about the correct way to get and interpret individual p values. I found this question, which links to a powerpoint about "partitioning" the chi square test, but I'm having trouble following it. Multiple Chi-Squared Tests


title: Alternatives to Dickey-Fuller test

text: I'm studying unit root tests and therefore the Dickey Fuller test and I can't seem to figure the following out. Is it correctly understood that the Dickey Fuller test (with drift and constant) is designed such that under the alternative the data generating process is $$y_t=c+\lambda t+\alpha y_{t-1}+\varepsilon_t?$$ If yes then i get we want the augmented Dickey Fuller to have a richer class of data generating processes as alternative so we consider the auxilliary regression $$y_t=c+\lambda t+\alpha y_{t-1}+\sum_{i=1}^k \beta_i \Delta y_{t-1} +\varepsilon_t$$ where $k$ is chosen such that the $\varepsilon$'s are white noise. Which now allows for alternatives $AR(k)$ models (is that correct?). Okay, now the way the test was presented to me the alternative is that the process is stationary, but how does the test perform if i.e. the actual data generating process is say $MA(q)$ or $ARMA(p,q)$? Do we somehow beforehand narrow $y$ down to being either $AR$ or a unit root process? Thanks in advance. To the editor: Actually my question was not so much intended to alternatives to Dickey-Fuller as much as to understand the whole "the alternative in the hypothesis is not a general stationary process yet it seems to me that we have that as hypothesis". I feel though my question was answered and of course it is always nice to know generalizations exist. But then again, maybe you just changed the title so that people will find the question if seeking those references.


title: Is the correlation coefficient better than we think?

text: In order to measure precisely the influence of one variable $Y$ over $X$ I wanted to use the mutual information because so far I believed that correlation coefficient (Pearson) was only limited to linear relation. But I realized that even in non-linear cases it gives good results. **Example in Matlab:** $X$ is a uniform random variable and $I$ is the range of $X$ $Y=X^{2},I=[0,100]: \rho=0.96$ $Y=\sqrt{X} ,I=[0,100]: \rho=0.94$ $Y=exp(X),I=[0,10]: \rho=0.70$ $Y=log(X),I=[0,100]: \rho=0.84$ Why the correlation is quite high in all those non-linear cases? I am beginning to think that correlation is only weak when the monotony of $X$ and $Y$ is not the same. Otherwise it is not so bad. Do some people have a point of view on this?


title: Can we estimate the size of a subset X of a set A, by randomly sampling subsets of A?

text: Let $A$ be a finite set and suppose we want to compute the size of some subset $X$. _Motivation_ : If we can generate elements $x$ of $A$ uniformly at random, then we can estimate the size of $A$ by random sampling. That is, we take $n$ random samples from $A$, if $m$ of them belong to $X$, then $|X|/|A| \approx m/n$. Unfortunately, for what I do, usually $|A|$ is massive and $|X|$ (while massive) is quite small with respect to $|A|$. So if I attempt to perform the above estimation, I'm likely to get $m=0$, which, while not useless, is not really all that satisfying. So, I have an idea that I'm hoping will speed up this process. Instead of throwing darts at a massive dart board, why don't I throw balls? That is, instead of sampling elements $x \in A$, we sample subsets of $A$. Surely I should be able to infer something about the density of $X$ in $A$ from this experiment. Suppose $A$ is equipped with a metric $d(x,y)$ (I have in mind the Hamming distance). For any $y \in A$ let $Y(y)=\{x \in A:d(x,y) \leq k\}$ be the closed ball of radius $k$ in $A$ centered at $t$. Since we can sample elements $x \in A$ uniformly at random, we can sample $k$-balls $Y_k(t)$ uniformly at random. Suppose (a) every $x \in A$ belong to exactly the same number of $k$-balls and (b) every $k$-ball has the same size $r$. Now suppose I generate $k$-balls $Y_1,Y_2,\ldots,Y_n$ uniformly at random and suppose $m=\sum_{i=1}^n |Y_i \cap X|$. It seems we can estimate $|A|$ in a similar way, that is $|X|/|A| \approx \frac{m}{rn}$. So my questions are: > Am I correct in that we can approximate $|X|$ this way? If so, I doubt I'm > the first to think of this, so is there a name for this method? I did actually test this on some sets, and it seems to match what I claim. > Are there any drawbacks to this approach? (e.g. is it less accurate? do I > need more samples?)


title: Wavelet auto correlation

text: I have a time serie that I want to analyse through a wavelet decomposition.   I am using the R package **WaveThres**. I am interested in the wavelet autocorrelation, but I struggle to understand what does it mean precisely. I have from the book Wavelet Methods in Statistics with R the following formula $\Psi_j(\tau)=\sum_{k}\phi_{j,k}(0)\phi_{j,k}(\tau)$ $\tau\in\mathbb{Z}$ being the lag of the autocorrelation   and   $\left \{\phi_{j,k}(t)=\phi_{j,k-t} \right \}_{j,k}$ a set of non decimated wavelets I would really appreciate to understand the meaning of this formula, and (why/if) it is different from performing a multi resolution analysis (MRA) and computing the Pearson autocorrelation coefficient on a detail. fRed


title: In a linear regression whose components can also be broken down, is it better to do multi-layered regression, or flatten to final components?

text: Consider a series like CPI (inflation), which I know is composed of a series of component prices (e.g. meat prices, grain prices, non-food prices, etc.), which in turn are also composed of a series of component prices (e.g. average meat prices are a combination of pork, beef, and chicken prices). If I wanted to use a regression to find the components of CPI and their weightings, then is it better to use the final components (regress CPI against pork, beef, and chicken prices), or is it better to create fitted version of the middle components, then regress against CPI (regress average meat prices against pork, beef, and chicken prices, then fit a meat price series, then regress CPI against the fitted meat price series)? Also, I should note that some series which were significant when fitting the interim components in the layered method - the latter method - lose their significance when the regression is flattened. So is is possible that the latter, layered method will return a better result because it includes the effect of more components? I have tried this with actual data and the latter method gives better results.


title: Hinge loss with one-vs-all classifier

text: I'm currently looking at the unconstrained primal form of the one-vs-all classifier $$\sum\limits_{i=1}^{N_I} \sum\limits_{k=1,\atop k \neq y_i}^{N_K} L(1+ \mathbf{w_k}\cdot\mathbf{x_i}-\mathbf{w_{y_i}}\cdot\mathbf{x_i})$$ where $N_I$ is the number of instances,   $N_K$ is the number of classes,   $N_F$ is the number of features,   $X$ is a $N_K \times N_F$ data matrix,   $y$ is a vector of class labels,   $W$ is an $N_K \times N_I$ matrix where each corresponds to the weights for the hyperplane splitting one class from the rest,   $L$ is some arbitrary loss function. My understanding is that the functional above tries to find a hyperplane for each class that maximizes the distance between the samples within the associated class to all other samples. If the hyperplanes are correctly positioned then $\mathbf{w_k}\cdot\mathbf{x_i}$ should always be negative, $\mathbf{w_{y_i}}\cdot\mathbf{x_i}$ should always be positive and our loss function should come back fairly low. I'm trying to implement this using the hinge loss which I believe in this case will end up being $\max(0,1+\mathbf{w_k}\cdot\mathbf{x_i}-\mathbf{w_{y_i}}\cdot\mathbf{x_i}$). However, in the above couldn't we end up with a situation where the hyperplanes classify all samples as belonging to every class. For example, if we are looking at the hyperplane seperating class 1 from all other classes, provided that $1+\mathbf{w_k}\cdot\mathbf{x_i}<\mathbf{w_{y_i}}\cdot\mathbf{x_i}$ then the incurred loss will be 0 despite $\mathbf{x_i}$ being classified as the wrong class. Where have I gone wrong? Or does it not matter whether $\mathbf{w_k}\cdot\mathbf{x_i}$ is negative or positive provided that $\mathbf{w_{y_i}}\cdot\mathbf{x_i}$ ends up with a higher score? I have a feeling that my use of the hinge function as I've described here is incorrect but my use of Google today has only led to more confusion. On a related note, why is there a 1 in the functional above? I would think that it would have little impact.


title: Calculate ROC curve for data

text: So, I have 16 trials in which I am trying to authenticate a person from a biometric trait using Hamming Distance. My threshold is set to 3.5. My data is below and only trial 1 is a True Positive:               Trial   Hamming Distance     1   0.34     2   0.37     3   0.34     4   0.29     5   0.55     6   0.47     7   0.47     8   0.32     9   0.39     10  0.45     11  0.42     12  0.37     13  0.66     14  0.39     15  0.44     16  0.39      My point of confusion is that I am really unsure about how to make an ROC curve (FPR vs. TPR OR FAR vs. FRR) from this data. It doesn't really matter which one, but I'm just really confused about how to go about calculating it. Any help would be appreciated.


title: Theoretical results for cross-validation estimation of classification accuracy?

text: For classification, what theoretical results are between cross-validation estimate of accuracy and generalisation accuracy? I particularly asking about results in a PAC-like framework where no assumptions are made that your function class contains the "true" function. I would love to know if there are theorems of the form: If your leave-one-out cross validation error-rate is $\theta$ on $N$ examples, then your generalisation error rate is lower than $\theta+\varepsilon$ with probability $f(\theta, \varepsilon, N)$. If so, what are the general proof techniques to obtain them? What is the theoretical framework? If a fully general theorem is impossible, what extra conditions, if any, allow you to arrive at this type of conclusion?


title: Model building, p value of GLM for whole model?

text: I have three GLM models (poisson distribution) that each have different factors but the same response variable. I plan to compare the fit of the models using AIC, but first I want to see if any of the models are significant according to a p value. Is is possible to get a p value for a whole model w/glm? Each model has 1 significant factor and one or two non-significant factors, but I was told to leave these non-significant factors in the models instead of removing them with sinful backwards stepwise regression.


title: How to compare results obtained by non-linear regression

text: Disclaimer: Statistics is not exactly something I am particularly good at. So if this is a stupid question, I do apologize I have been using non-linear regression to analyze binding between two biomolecules. Each individual experiment was done in triplicate. The average of these triplicate values (with associated error) were analyzed using regression ([L]=known). > Specific binding = Bmax*[L]/(Kd+[L]) The Bmax and Kd are calculated. Both values have 95% confidence intervals and standard errors reported. I would like to compare the Kd for different pairs of Biomolecules. I.e. I would like to compare the Kd (higher or lower) for A binding to B, A binding to C and A binding to D. In addition, I would like to compare binding between AB and XY. Which test(s) would be appropriate for this analysis?   Any suggestions would be greatly appreciated!!


title: Entropy of noncentral multivariate t distribution

text: I am looking at calculating the entropy of a multivariate noncentral t-distribution. I found in the central case (mean is zero) the entropy is the same as multivariate normal with a correction for degrees of freedom (ie only depends on covariance and $\nu$). But in the noncentral case the extra term is a more complicated infinite sum which depends both on $\mu^T \mathbf{R}^{-1} \mu$ and $\nu$. I guess I was just trying to gain some intuition as to why the entropy would depend on the mean $\mu$. If the distribution is shifted around in space the entropy (which is the integral over the whole thing) should stay the same (as with the normal)? The expression I have for the entropy is (from _Multivariate t Distributions and Their Applications_ , Kotz and Nadarajah): $$ H = \frac{1}{2}\log|R| + \log \left[ \frac{(\nu \pi)^{p/2}}{\Gamma(p/2)} B \left(\frac{p}{2},\frac{\nu}{2}\right) \right] + \frac{\nu + p}{2}M(\nu,p,\Delta) $$ where $\Delta=\mu^T R^{-1}\mu$ and $$ M(\nu,p,\Delta) = \exp\left(-\frac{\Delta}{2}\right)\sum_{j=0}^{\infty}\frac{1}{j!}\left\{ \psi\left( \frac{\nu+p+2j}{2}\right) - \psi\left(\frac{\nu}{2}\right) \right\} $$ A second part of the question would be for any practical tips for calculating this quantity (are there any approximations or tricks for the infinite sum)?


title: Using interaction term for describing non-linear relationship in regression model

text: I have created a Rscript for illustrating my question:               library(epicalc)          temp.df <- data.frame(x=1:100,y=log(1:100))     set.seed(1)     temp.df[,"y"] <- temp.df[,"y"]*sample((10000:10100)/10000,size=100,replace=TRUE)     temp.df[temp.df[,"x"]<=15,"x_gp"] <- "<=15"     temp.df[temp.df[,"x"]>15,"x_gp"] <- ">15"     temp.df[,"x_gp"] <- factor(temp.df[,"x_gp"],levels=c("<=15",">15"))          glm.null <- glm(y~NULL,data=temp.df)     glm.01 <- glm(y~x,data=temp.df)     glm.02 <- glm(y~log(x),data=temp.df)     glm.03 <- glm(y~x*x_gp,data=temp.df)          lrtest(glm.null,glm.01)     lrtest(glm.null,glm.02)     lrtest(glm.null,glm.03)     lrtest(glm.01,glm.02)     lrtest(glm.02,glm.03)     lrtest(glm.01,glm.03)          (logLik(glm.02)*-2)-(logLik(glm.03)*-2)          plot(temp.df[,c("x","y")],col="black",pch=15)     points(temp.df[,"x"],predict(glm.01),col="blue",pch=20)     points(temp.df[,"x"],predict(glm.02),col="red")     points(temp.df[,"x"],predict(glm.03),col="green",pch=16)      The dataframe is describing a log-linear relationship. But in reality, maybe it is not log-linear, even if it is log-linear, it may not be preferred to use `log`, as it may be more difficult to understand. I thought of "bending" the linear regression into two "segments" (`glm.03`) using the interaction term. I am aware that it is not as good as `log` in my example but already a lot better than using a straight line (`glm.02`) alone. I wonder if my method will be too `unorthodox` in the field of public health, so can anyone suggest other ways for describing a log-linear-like relation, without using `log`? Thanks.


title: How to find this conditional probability in sampling without replacement?

text: Here is the complete question: A club consists of $4n$ people of half men and half women. A woman is elected as the chair of a committee from this group. As all members are equally talented, she prefers to select additional $n$ people at random from the remaining $4n-1$ people to work with her for a particular project. If $n$ people chosen happen to be of the same sex, what is the probability that all are (a) male (b) female? Hi Everyone, I am taking 1st year statistics and this is an exam practice question that I am a little stuck on - would greatly appreciate some help! Here is what I have so far: For males, I have thought about it this way: Since a woman is elected as Chair, she would then have $2n$ men and $2n-1$ women to choose from. The probability of selecting $n$ people that are all men for the project from the 1st three picks up to the $nth$ pick, I have: $$\frac{2n}{4n-1} \cdot \frac{2n-1}{4n-2} \cdot \frac{2n-2}{4n-3} \cdot \space ... \cdot \frac{2n-(n-1)}{4n-n}. $$ Similarly, for females I have the probability of selecting all women to be: $$\frac{2n-1}{4n-1} \cdot \frac{2n-2}{4n-2} \cdot \frac{2n-3}{4n-3} \cdot \space ... \cdot \frac{2n-n}{4n-n}. $$ Firstly, am I on the right track, and secondly how can I find a closed formula for this probability? Many thanks in Advance!


title: Problem with ROC curve in R

text: I am trying to plot the ROC curve for a random forest model (`ROCR` package), and I am getting weird results. I have double-checked the code several times, but I can't figure out why the functions give out an AUC of 1 and a perfect ROC curve. If I do the same with the data included in the package (`ROCR.simple`), everything goes fine. As far as I know, an AUC of 1 would mean that the model is predicting with `p = 1` all the `1's` in the data, while the `0's` are predicted as `p = 0`. Right? As the data below shows, this is certainly not the case. Am I missing something? I feel like an idiot with this thing, but I haven't found anything. Thanks! **Code and results**               > pp      [1] 0.679 0.056 0.045 0.194 0.100 0.666 0.028 0.057 0.079 0.042 0.065 0.677      [13] 0.057 0.035 0.679 0.184 0.032 0.127 0.017 0.682 0.030 0.781 0.067 0.281      [25] 0.879 0.009 0.642 0.787 0.215 0.862 0.851 0.095 0.089 0.049 0.158 0.075      [37] 0.669 0.030 0.134 0.641 0.017 0.650 0.027 0.038 0.091 0.011 0.056 0.061      [49] 0.018 0.011 0.040 0.064 0.700 0.208 0.048 0.017 0.127 0.902 0.166 0.110      [61] 0.057 0.027 0.001 0.016 0.121 0.021 0.331 0.044 0.029 0.024 0.021 0.017      [73] 0.005 0.655 0.019 0.148 0.010 0.101 0.233 0.011 0.121 0.768 0.091 0.010      [85] 0.096 0.308 0.047 0.067 0.654 0.014 0.018 0.018 0.083 0.218 0.100 0.765      [97] 0.069 0.779 0.084 0.063 0.134 0.155 0.133 0.088 0.819 0.044 0.251 0.704      [109] 0.141 0.664 0.066 0.029 0.870 0.020 0.050 0.189 0.041 0.041 0.029 0.055      [121] 0.075 0.022 0.006 0.009 0.012 0.646 0.057 0.012 0.001 0.055 0.645 0.073      [133] 0.915 0.130 0.166 0.142 0.700 0.005 0.014 0.063 0.008 0.670 0.003 0.042      [145] 0.065 0.006 0.027 0.741 0.907 0.951 0.759 0.762 0.266 0.002 0.021 0.074      [157] 0.029 0.896 0.126 0.155 0.000 0.000 0.690 0.140 0.019 0.006 0.028 0.037      [169] 0.036 0.033 0.695 0.032 0.038 0.067 0.778 0.804 0.076 0.097 0.055 0.031      [181] 0.058 0.120 0.134 0.035 0.112 0.963 0.126 0.036 0.016 0.063 0.035 0.009      [193] 0.082 0.805 0.800 0.175 0.024 0.040 0.032 0.014 0.088 0.203 0.082 0.023      [205] 0.131 0.727 0.925 0.044 0.046     > obs      [1] 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 1 0 1 1 0 1 1 0 0 0 0 0 1      [38] 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1      [75] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0      [112] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1      [149] 1 1 1 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0      [186] 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0     > require(ROCR)     > pred = prediction(predictions = pp, labels = obs)     > perf = performance(pred, "tpr", "fpr")     > plot(perf)     > performance(pred, "auc")@y.values[[1]]     [1] 1      ![My data](http://i.stack.imgur.com/uJov3.png) And with the `ROCR.simple` data:               > data(ROCR.simple)     > performance(prediction(ROCR.simple$predictions, ROCR.simple$labels), "tpr", "fpr")     > plot(perf)     > performance(prediction(ROCR.simple$predictions, ROCR.simple$labels), "auc")@y.values[[1]]     [1] 0.8341875      ![ROCR.simple data](http://i.stack.imgur.com/15ayM.png)


title: Sampling from dataset without sacrificing power of multivariate logistic regression

text: I'm developing a multivariate logistic regression model for a large healthcare dataset (100 million+ rows of demographic and diagnosis columns) where the event in question, an inpatient admission day, occurs infrequently - say for every 10,000 patient days there may be only 100 IP admit days. The sheer size of the dataset is going to tax my desktop machine's resources, so I'm looking for ways to improve my computer performance. Given that the power of detecting a rare binary event is determined by the number of events rather than the proportion itself, can I safely reduce my dataset size by randomly sampling from the non-event records without sacrificing power? For example, instead of 1 million records where Y=1 and 100 million records where Y=0, take a random 5% sample from the 100 million non- event records.


title: Neural Networks: weight change momentum and weight decay

text: Momentum $\alpha$ is used to diminish the fluctuations in weight changes over consecutive iterations: $$\Delta\omega_i(t+1) = \omega_i - \eta\frac{\partial E}{\partial w_i} + \alpha \Delta \omega_i(t),$$ where $E({\bf w})$ is the error function, ${\bf w}$ - the vector of weights, $\eta$ - learning rate. Weight decay $\lambda$ penalizes the weight changes: $$\Delta\omega_i(t+1) = \omega_i - \eta\frac{\partial E}{\partial w_i} - \lambda\eta\omega_i$$ **The question is if it makes sense to combine both tricks during the back- propagation and what effect it would have?** $$\Delta\omega_i(t+1) = \omega_i - \eta\frac{\partial E}{\partial w_i} + \alpha \Delta \omega_i(t) - \lambda\eta\omega_i$$


title: Multinomial Logistic Regression: IIA violated - What to do?

text: I have a data set with 561 observations. I want to compare two groups on an outcome measure with 3 categories. I have a couple of covariates. I used the following syntax (Stata SE 12.1):               mlogit outcome_measure covariate_1 covariate_2 i.covariate_3 i.covariate_4 group      To check for Independence from Irrelevant Attributes (IIA) I used               mlogtest , haus      This gave me the following output: ![output of "mlogtest , haus"](http://i.stack.imgur.com/6FbpS.png) I understand that this means that the data do not meet the IIA assumption. Stata then suggests to use the suest (Seemingly Unrelated Estimation) command. I could try this, but I am not sure if this is a valid alternative. What do you guys think? **EDIT (2013-09-18):** @GUNG: I asked one of my statistics professors. He said that the additional tests won't be performed if the assumption is not violated for one of the outcomes (in the case of three outcome categories). He therefore thinks that I can proceed with the mlogit model. @Dimitriy V Masterov I ran the Small-Hsiao IIA test. These are the outcomes: ![Outcomes Stata Small-Hsiao IIA test](http://i.stack.imgur.com/kpPvn.png) I guess this means IIA is not violated. Do you agree? **EDIT (2013-09-22):** As Dimitriy V. Masterov suggested, I ran the new syntax for Hausman. Because I have 3 categories in my outcome variable, I compared the Full model to all models with one category of my outcome variable omitted (Full vs. 1+2, Full vs. 1+3, Full vs. 2+3). These are the outcomes: ![Hausman - Ignoring Category 1](http://i.stack.imgur.com/RVuL2.png) ![Hausman - Ignoring Category 2](http://i.stack.imgur.com/ZUdRT.png) ![Hausman - Ignoring Category 3](http://i.stack.imgur.com/Zls6w.png) It looks like leaving out the third category changes the model significantly. Does this mean that I cannot use Multinomial Logistic Regression and that I should move to (the suggested) suest (Seemingly Unrelated Estimation)? suest is giving the following outputs: ![suest output](http://i.stack.imgur.com/5Fwto.png) The output looks fine to me and it supports my hypothesis, but I am not sure if suest is valid, what it assumptions are and how I can test these. **Any comments are welcome!**


title: Fisher's exact test in R - 2x4 table

text: I have a question on categorical data - I have a 2x4 matrix and I want to test the difference in the number of males or females in each group, and so I entered the data into R. So I input this:               A=c(31,7)     B=c(8,1)     C=c(39,16)     D=c(2,6)     tab=as.table(cbind(A,B,C,D))     row.names(tab)=c('males','females')     fisher.test(tab)      This is the output that I get -               Fisher's Exact Test for Count Data          data:  tab      p-value = 0.01077     alternative hypothesis: two.sided      The only statistical information I get is the p-value - how do I know _how_ the groups are different? Thanks!


title: Trajectory Clustering: Which Clustering Method?

text: As a newbie in Machine Learning, I have a set of trajectories that may be of different lengths. I wish to cluster them, because some of them are **actually the same path** and they just **SEEM different** due to the noise. In addition, **not** all of them are of the **same lengths**. So maybe although Trajectory A is not the same as Trajectory B, yet it is **part of** Trajectory B. I wish to **present this** property after the clustering as well. I have only a bit knowledge of `K-means Clustering` and `Fuzzy N-means Clustering`. **How may I choose between them two? Or should I adopt other methods?** **Any method that takes the "belongness" into consideration?** (e.g. After the clustering, I have 3 clusters `A, B and C`. One particular `trajectory X` belongs to `cluster A`. And a shorter `trajectory Y`, although is not clustered in `A`, is identified as part of `trajectory B`.) =================== UPDATE ====================== The aforementioned trajectories are the pedestrians' trajectories. They can be either presented as a series of (x, y) points or a series of step vectors (length, direction). The presentation form is under my control.


title: How to implement bagging with under sampled majority class in decision trees (e-miner)?

text: I have a dataset where the target variable = Y to N ratio is highly skewed around 1% (events/successes) and 99% (non events/failures). I need to detect patterns using the independent variables which distinguish the events from non-events, in my target variable. As of now from the original input data which has this skewed proportion, I have been able to extract a sample which is balanced in nature (i.e 50% proportion for both events and non events). In this sample I have selected all the records for events and randomly selected same number of records from the non event class. I further ran a decision tree and noted the results. I have few questions regarding how to implement the same in E-Miner 7.1: 1) I now want to create multiple trees with a process flow which looks like: Input Data --> Sample --> Data Partioning --> Decision Tree. How can I create multiple samples from the larger set of non-events and compare them with the same "event" records for say 10 to 15 iterations? This is because of the huge number of non event records, which necessarily cannot be ignored. Can someone let me know how can bagging be implemented using Start Group node in the process flow? Where should it be placed (right after Sample?)? 2) How do I combine the results/patterns from multiple iterations post running various decision tree nodes? I understand that "Ensemble" node can be handy in such a scenario, but again it would do a voting/max of predicted probabilities for every records and combine results, however, I am interested in the patterns/business rules? Not exactly sure on how would patterns/rules get combined in E-Miner? Any help on these would be highly appreciated. Thanks.