title: Random forest algorithm

text: I have a question about `Random Forest algorithm`: 1.Let the number of training cases be N, and the number of variables in the classifier be M. 2.We are told the number m of input variables to be used to determine the decision at a node of the tree; m should be much less than M. 3.Choose a training set for this tree by choosing n times with replacement from all N available training cases (i.e., take a bootstrap sample). Use the rest of the cases to estimate the error of the tree, by predicting their classes. 4.For each node of the tree, randomly choose m variables on which to base the decision at that node. Calculate the best split based on these m variables in the training set. 5.Each tree is fully grown and not pruned (as may be done in constructing a normal tree classifier). I think I understand the most steps, except `step 3`. How do we in practice predict the classes? Does anyone have a simple explanation? Because I cant find any explanation online... I understand that : `Finally a response class is predicted in each terminal node of the tree (or each rectangular section in the partition respectively) by means of deriving from all observations in this node either the average response value in regression or the most frequent response class in classification trees. Note that this means that a regression tree creates a piecewise (or rectangle-wise for two dimensions and cuboid-wise in higher dimensions) constant prediction function.` But I am trying to understand `How do we use Out Of Bag data` later for `finding prediction error`? For example I am trying to understand `In an ensemble of trees the predictions of all individual trees need to be combined. This is usually accomplished by means of (weighted or unweighted) averaging in regression or voting in classification. The term “voting” can be taken literally here: Each subject with given values of the predictor variables is “dropped through” every tree in the ensemble, so that each single tree returns a predicted class for the subject. The class that most trees “vote” for is returned as the prediction of the ensemble. This democratic voting process is the reason why ensemble methods are also called “committee” methods. Note, however, that there is no diagnostic for the unanimity of the vote. For regression and for predicting probabilities, i.e. relative class frequencies, the results of the single trees are averaged; some algorithms also employ weighted averages.` Best Regards!


title: In SPSS, how do I analyze the similarity of multiple scores, differentiated by another variable?

text: I could use some help figuring out how best to analyze my data with SPSS. I'm new to statistical analysis and SPSS, and trying to solve these issues is melting my brain. Sorry if this question is confusing or too nonspecific. Some of the details of my research have been changed for privacy reasons. I created and distributed a survey to two schools--let's call them Alpha College and Beta University. I'm studying how important students think it is to follow rules with regards to plagiarism, drugs, and sex, and I've hypothesized that, due to a difference in campus culture, students at Alpha College are going to be more consistent in how they view the rules--that is to say, Alpha students are going to find the rules for all three subjects similarly important. On the contrary, I hypothesize that the campus culture at Beta University is a lot less salient when it comes to rules, and so I expect that Beta students will experience more variation in importance across the three subjects. I measured "rule importance" with three different multiple-item Likert scales according to subject--so I have three variables, `plagiarismImportanceScore`, `drugsImportanceScore`, and `sexImportanceScore`. Each of these Likert scales contained a different number of items, so where a maximum score on `plagiarismImportanceScore` might be 25, the maximum on `drugsImportanceScore` might be 15 instead. I've also measured `school`\-- a value of 1 means "Alpha College," and a value of 2 means "Beta University." How, in SPSS, do I test how similar the scores for `plagiarismImportanceScore`, `drugsImportanceScore`, and `sexImportanceScore` are, differentiated between Alpha students and Beta students? Ultimately, what I'm looking for is data that essentially says, "Alpha students have more consistent attitudes toward rules as a whole than Beta students." Any and all help will be appreciated. Thank you.


title: Interpreting graphs

text: I'm not sure what to ask so forgive me if this is a bit all over the place. I've written a program to analyze the time it takes to query a data structure and generated a graph which shows on the Y axis, the time a query took and on the X axis the number of items in the structure. I didn't know what to expect at all but I'm trying interpret the results and objectively describe what I've recorded. Probably from my description you can guess I'm a complete noob in the stats arena. What I'm getting is something like this ![random graph](http://i.stack.imgur.com/t6Mms.png) Couple of questions come to mind: How do I go about inferring anything from a graph like this? Given the erratic Y values is there a common approach to somehow normalizing these? It also doesn't seem like a situation where you'd just kinda draw a line of best fit and go with that... Is that even the right thing to do? From my research so far people have mentioned standard deviation, log distribution and a host of other things but not being trained in this I don't understand when to use these things and what purpose they serve as it would apply to my case. I could upload the real graph but I'd prefer some pointers as to how to approach this generally so that I can learn from it and have a better understanding to be able to do this next time. Graph shown is taken from http://openclipart.org/detail/170148


title: How to aggregate human estimates?

text: I'm working on a program that involves people estimating numbers. For instance, when estimating $X$, one person says that $X$ is the product of two numbers over which he has probability distributions $\mathcal N(2000,50)$ and $\exp(0.05)$, and another person says that he has probability distribution $\mathcal N(30000, 200)$ over $X$. If the people are equally trustworthy, what should our posterior be? Actually I don't need the entire posterior, I just need its expectation and its variance. Is there a way to do this for arbitrary products/sums of distributions without using a symbolic math engine? I know that we can compute, say, the expectation of the product of independent random variables by just taking the product of their expectations, but these sorts of estimates aren't independent, are they?


title: Recalculate log-likelihood from a simple R lm model

text: I'm simply trying to recalculate with dnorm() the log-likelihood provided by the logLik function from a lm model (in R). It works (almost perfectly) for high number of data (eg n=1000) :               > n <- 1000     > x <- 1:n     > set.seed(1)     > y <- 10 + 2*x + rnorm(n, 0, 2)     > mod <- glm(y ~ x, family = gaussian)     > logLik(mod)     'log Lik.' -2145.562 (df=3)     > sigma <- sqrt(summary(mod)$dispersion)     > sum(log(dnorm(x = y, mean = predict(mod), sd = sigma)))     [1] -2145.563     > sum(log(dnorm(x = resid(mod), mean = 0, sd = sigma)))     [1] -2145.563      but for small datasets there are clear differences :               > n <- 5     > x <- 1:n     > set.seed(1)     > y <- 10 + 2*x + rnorm(n, 0, 2)     >      > mod <- glm(y ~ x, family = gaussian)     > logLik(mod)     'log Lik.' -8.915768 (df=3)     > sigma <- sqrt(summary(mod)$dispersion)     > sum(log(dnorm(x = y, mean = predict(mod), sd = sigma)))     [1] -9.192832     > sum(log(dnorm(x = resid(mod), mean = 0, sd = sigma)))     [1] -9.192832      Because of small dataset effect I thought it could be due to the differences in residual variance estimates between lm and glm but using lm provides the same result as glm :               > modlm <- lm(y ~ x)     > logLik(modlm)     'log Lik.' -8.915768 (df=3)     >      > sigma <- summary(modlm)$sigma     > sum(log(dnorm(x = y, mean = predict(modlm), sd = sigma)))     [1] -9.192832     > sum(log(dnorm(x = resid(modlm), mean = 0, sd = sigma)))     [1] -9.192832      Where am I wrong ?


title: How to interpret the significance of logistic regression variable? Something analogous to a partial F test in logistic regression

text: For ordinary linear regression with Gaussian noise, it is easy to interpret the significance of a variable. This is consistent with a partial F test. The square of the t-test for the second variable equals to the partial F-test statistic, and their p-values are the same. I wrote simple R codes to confirm this. Is there something like this for logistic regression? I thought/hoped that the likelihood ratio test would correspond to this, but no. What should I do if the variable and the likelihood ratio test (of adding that particular variable) do not have the same (in)significant effect? I appreciate your time and help,               rm(list=ls(all=TRUE))      n = 100   ;       x1 = runif(n,-4,4)   ;       x2 = runif(n,6,10)     y = 3*x1 + 8*x2 + rnorm(n,2,4)     l1 = lm(y~x1)  ;  l2 = lm(y~x1+x2)  ;  a = anova(l1,l2)          summary(l1)$coeff                  Estimate Std. Error   t value     Pr(>|t|)     (Intercept) 66.093853  1.0123131 65.289929 1.385202e-82     x1           3.199212  0.4292828  7.452458 3.664499e-11          summary(l2)$coeff                 Estimate Std. Error    t value     Pr(>|t|)     (Intercept) 2.767750  2.7871368  0.9930441 3.231592e-01     x1          2.870897  0.1707022 16.8181610 1.648852e-30     x2          7.871545  0.3428392 22.9598753 5.370614e-41          (summary(l2)$coeff[3,3])^2     527.1559     >     a          Analysis of Variance Table          Model 1: y ~ x1     Model 2: y ~ x1 + x2       Res.Df    RSS Df Sum of Sq      F    Pr(>F)         1     98 9899.1                                       2     97 1538.4  1    8360.6 527.16 < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          >     a$F ; a$Pr        [1]       NA 527.1559     [1]           NA 5.370614e-41     >      >      >      > rm(list=ls(all=TRUE))      > n = 100     > x1 = runif(n,-4,4)     > x2 = runif(n,6,10)     >      > y = rbinom(n,1,1/(1+exp(-3*x1 - 2*x2 + 20)))     >      > l1 = glm(y~x1,family=binomial)     > l2 = glm(y~x1+x2,family=binomial)     >      > a = anova(l1,l2)     >      > summary(l1)$coeff                      Estimate Std. Error   z value     Pr(>|z|)         (Intercept) -2.988069   0.812041 -3.679702 2.335068e-04         x1           2.115333   0.498431  4.243984 2.195858e-05         > summary(l2)$coeff                   Estimate Std. Error   z value     Pr(>|z|)     (Intercept) -17.215960  5.5710699 -3.090243 0.0019999276     x1            3.048657  0.8618367  3.537395 0.0004040949     x2            1.675323  0.5976386  2.803238 0.0050592272     >      > (summary(l2)$coeff[3,3])^2         [1] 7.858145         >          > l1$deviance -  l2$deviance         [1] 13.65371         > pchisq(l1$deviance -  l2$deviance,df=1)     [1] 0.9997802     >      > a     Analysis of Deviance Table          Model 1: y ~ x1     Model 2: y ~ x1 + x2       Resid. Df Resid. Dev Df Deviance     1        98     45.534                 2        97     31.880  1   13.654     > a$F         NULL         > a$Pr         NULL


title: How to summarize GAM model result from multiple imputation data in R

text: I am very new to R and not very experienced in statistics. I have this general question regarding applying Generalized Additive Models (GAM) in multiple imputation dataset. I used R package _mice_ for multiple imputation and _mgcv_ for GAM. But I have no idea how to combine the results from each imputed data set. Does anyone have any good idea about this? Thanks a lot!


title: Permutation test: Exactness

text: I am currently reading about permutation/randomization tests and have some difficulties to understand why they are exact. More precisely, I consider two groups of random variables with means $\mu_1$ and $\mu_2$ and variances $\sigma^2_1$ and $\sigma^2_2$, which are assumed to be equal if $ \mu_1 = \mu_2$ holds. To test the one-sided hypothesis $H_0: \mu_1 \leq \mu_2$ versus $H_1: \mu_1 > \mu_2$, I apply a permutation test with the same test statistic as for a two sample t-test with unequal variances and unequal sample sizes. It was no big deal to prove and to understand that the level of significance of the permutation test is equal to $\alpha$ for $\mu_1 = \mu_2$. However, I don't understand why the level of significance is less than $\alpha$ for $ \mu_1 \leq \mu_2$!? Can someone give me a hint?


title: Evaluating the relationship between two attributes in a gene study

text: I have a data set for thousands of genes. Each gene has two attributes represented as a percentage range(0-100). One attribute (V3) has a normal distribution but the other (V2) is highly negatively skewed (e.g 99% of the values are greater than 99%). I want to know if the out liars of V2 have any relationship with the V3 attribute? I thought about partitioning the data into percentiles and treat it as categorical which would allow me to use chi- squared test but I am not sure if this is the best method or not. I also would like to get your input as to if a numerical method can be used with this data? ![scatter plot of the two attributes](http://i.stack.imgur.com/NARjE.png) I am very new at asking questions here so please let me know if there is anything I can do to improve the question. Thank you


title: Multilevel model with an ordinal outcome

text: I have a data set from a survey of health problems and annoyance by traffic noise in people living in residential buildings, with "5-point scale answer" questions. I want to perform a multilevel model for an ordinal response (health problem or noise annoyance), and I tried with the `gllamm` command in Stata but it's not allowed the option `link(ologit)`. Does anyone know another way to perform it? If not in Stata, in R or SPSS would be useful anyway.


title: What is the proper way to estimate the probability (proportion of time) a rare event occurs?

text: Often, I need to estimate the probability (proportion of time) a rare event occurs. The standard MLE estimate often gives me extreme estimates since the denominator is usually 1, and the numerator is either 0 or 1, giving me either 100% or 0%. For example, I am trying to estimate the proportion of web referrals as a result of my email campaign for each of my users. Since the events are rare, most of my users usually have only 1 web referral, and they have either 0 email referral or 1 email referral. In such cases, the MLE estimate is quite unreliable. Are there standard tricks to correct this over-under estimation? Perhaps something like the laplace smoothing? If yes, how should I go about it?


title: Which algorithm would be most appropriate for this

text: Can someone tell me which algorithm would be best to solve this problem? I am thinking this belongs in the category of 'Market Basket Analysis' aka ('Affinity Analysis'), but I am not sure. Here's the problem: A user belongs to certain types of groups. When our client gives us a set of users that belong to a particular group, we need to tell them which other groups those users might be most interested in as well as most disinterested in. To me this sounds similar to the classic Market Basket Analysis example: 'customers who buy diapers also buy beer'. So I was thinking - if I pre- compute all the affinities within groups then I will be able to answer our client's question. No? Would greatly appreciate any thoughts on this? What kind of algorithm would be most appropriate for this? Thanks for your help.


title: Are ratio, interval, ordinal and nominal variables nested?

text: It seems to me that ratio, interval, ordinal and nominal variables nested in the following sense:   * A ratio variable is also interval, ordinal and nominal.   * An interval variable is also ordinal and nominal.   * An ordinal variable is also nominal Are the above statements accurate? If not, are there any counter-examples where one or more of the above statements will not hold?


title: Why do we examine the rate of convergence when we find a confidence interval?

text: > ![enter image description here](http://i.stack.imgur.com/KEcP0.jpg) I would like to help me with this. I don't understand why do we examine the rate of convergence. Also what do we mean by saying "the error is usually dominated by the variance, not the bias" Thank you very much.


title: What is the problem with my lagged variables?

text: There are variables a,b,c, and d.   Variables a and b is from firm's financial statements and they have time- lag.(financial statements take time to be announced)   Variable d is a dummy variable indicating 2008 financial crisis and it should be significant and positive. I tested the normal one with fixed effect model,   and then using lagged variables of a and b tested again.   R square increased 0.02 but dummy variable has a negative insignificant effect!   So I reject the hypothesis they have a time-lag effect but how can I say to confirm this? Since R square is increased a little I don't have any reason to ignore that lagged model...


title: Prove Confidence of Rule

text: Using Association. I want to prove the confidence of the rule $B'->(A\B')$ can't be more than the confidence of $B->(A\B)$ Given:   * frequent itemset $A$ and subsets $B$ of $A$.    * $B'$ is a subset of $B$. "\" is the Set Difference Operator Sometimes written like this: $(X-Y)$ ![Dataset](http://i.stack.imgur.com/yB2cg.jpg) I found the exact question on this book: Data Mining: _Concepts and Tenchniques by Jiawei Han_ online, under Mining Frequent Patterns, Association, and Correlations chapter 6.3 (c) on this book. Doesn't say anything on how to solve it though.


title: Select a formula-based model by AIC

text: Let's say I want to investigate if there is any correlation between the response (continuous) variable `individual_fish_size` and the three explanatory (continuous) variables `Depth`, `substrate` and `Temperature`. But as my fishes come from several different species and the distribution of these species along the three variables is certainly not random, I'd like to include the categorical variable `Fish.species` to get rid of the variance due to difference of size along `fish.species`. Once I ran such a model, does it make sense to do a selection (a backward selection for example) of the best formula by AIC (This can be realized with the function `step` in R)? Or would we do better to avoid such a thing because we might lose the part of the variance that is explained by `fish.species`.


title: Likelihood of Conditional Grouped Continuous Model

text: ![Likelihood](http://i.stack.imgur.com/J3Lhq.jpg) I would like to find MLE of the likelihood above by using optim function in $R$. However, I couldn't understand the terms. I couldn't write the likelihood in $R$. I have the data given, some of them are cdf denotions but what is the 2nd log, or 4th log term. I would be so glad if anyone help me. Thanks


title: standard score or p value

text: I have two groups of data, each contains five sets. Each of them has a standard score. I would like to compare them. However I don't know if it is good to directly compare standard scores. As you know, the difference in standard score -2 and -1 is different from -1 and 0. So, should I use p value instead? Sorry for this silly question.


title: What is a loss function in decision theory?

text: My notes define a loss function as the 'cost' incurred when the true value of $\theta$ is estimated by $\hat\theta$. What kind of cost is it talking about? monetary cost? or is it something related to errors?


title: Joint distribution of two multivariate normal distributions

text: If we define 2 independent variables $Y_1$ and $Y_2$ as follows: \begin{align} Y_1 &= (Y_{11},Y_{12},Y_{13})^T \sim\mathcal N_3(\mu_1,\Sigma_{11}), \\\ Y_2 &= (Y_{21},Y_{22})^T \sim\mathcal N_2(\mu_2,\Sigma_{22}) \end{align} where, \begin{align} \mu_1 &= (2, 2, 2)^T &\Sigma_{11} &= \left[\begin{array}{ccc} 3 &1 &0 \\\ 1 &2 &0 \\\ 0 &0 &3 \end{array}\right] \\\ \mu_2 &= (3, 4)^T &\Sigma_{22} &= \left[\begin{array}{cc} 4 &2 \\\ 2 &4 \end{array}\right] \end{align} Then how can I find the joint distribution of $Y_{11}-Y_{13}+Y_{22}$ and $Y_{21}-Y_{12}$? I know its a simple question but I could find if it was asked for $Y_1-Y_2$ or something. How am I supposed to solve it when it is like that?


title: On FDA guidance about Bayesian practice

text: US FDA authorizes the use of Bayesian statistics with informative priors (in certain contexts): http://www.outsourcing-pharma.com/Clinical-Development/US-FDA-says-Bayesian- analysis-could-cut-device-trial-costs http://www.stat.rutgers.edu/iob/bioconf09/slides/Campbell.pdf FDA Bayesian guidance In the slide "Importance of Simulation" of Campbell's slides, it is written: _"So simulate to show that Type 1 error (or some analog of it) is well- controlled"._ However Type 1 error is not well-controlled with an informative prior. For instance consider a simple binomial model $x \sim \text{Bin}(n,\theta)$ with a Beta prior $B(a,b)$ on the unknown proportion parameter $\theta$, and consider a credibility interval $I(x)$ for some given credibility level $100(1-\alpha)\%$. Then the frequentist coverage function $\theta \mapsto \Pr(I(x) \ni \theta \mid \theta)$ is close to $100(1-\alpha)\%$ for the Jeffreys prior $B(a=\frac12,b=\frac12)$, but when the sample size is not large and $(a,b)$ is far from $(\frac12,\frac12)$ then the frequentist coverage is far from $100(1-\alpha)\%$ (possibly except for some very particular values of $\theta$, but anyway the coverage is not "controlled"). The Type 1 error of Bayes factor tests is not controlled too. So according to which point of view could one achieve good frequentist properties with Bayesian inference under an informative prior ?


title: Why do we need sparsity for auto-encoders?

text: Auto-encoders are used for dimensionality reduction, and as a tool for unsupervised feature learning. Further, auto-encoders have been used for building and training multi-layer neural networks. When we talk about auto-encoders, we introduce a sparsity and I want to know, what is the purpose of introducing an sparsity term for auto-encoders?


title: Auto.arima vs autobox do they differ?

text: From reading posts on this site I know there is an R function `auto.arima` (in the `forecast` package). I also know that IrishStat, a member of of this site built the commercial package autobox in the early 1980s. As these two packages exist today and automatically select arima models for given data sets what do they do differently? Will they possibly produce different models for the same data set?


title: Performing multiple linear regressions, in Excel, that have a common x-intercept?

text: I was plotting some linear data sets in Excel, including linear trend-lines: ![enter image description here](http://i.stack.imgur.com/KRRJu.png) I was about to perform 5 separate linear regressions, so I could get the `slope` and `y-intercept` of each _"independent"_ data set. But then, in a flash, I realized that data might have a common point: ![enter image description here](http://i.stack.imgur.com/OH5IF.png) In fact, it might even be: the **x-intercept** itself. So what I need now is way to run multiple simultaneous linear regressions, with the assumption that all of the lines intersect at a common point. Does such a linear regression analysis method exist? Does it have a name? Does it exist in Excel? * * * The best I've been able to muster so far is to run five independent linear regressions, getting the slope and intercept of each data set:               Slope (m)  Intercept (b)     =========  =============      1.15287     11484.8        0.86301      7173.5      0.43212      4306.4      0.25894      2853.6      Plotting `slope` against `y-intercept` you see some kind of correlation: ![enter image description here](http://i.stack.imgur.com/pwgAb.png) If I **_assume_** that my data sets share an `x-intercept`, then I can find that x-value through:                y = mx + b      0 = mx + b     -b = mx      x = m / -b      Which gives:               Slope (m)  Intercept (b)  Common x-intercept (assuming their is one)     =========  =============  ============================      1.15287     11484.8        -9961.9      0.86301      7173.5        -8312.2      0.43212      4306.4        -9965.6      0.25894      2853.6       -11020.2      Which, aside from one really wonky point, converges pretty well.


title: Advice needed for time series correlation

text: I have two time series data sets that I wish to compare. One is from an automated process of laboratory generated medical test results. The positive test results occur every day over N years. The second is from manual collection of clinical information about the patients from on which the medical tests were performed. My null hypothesis is that there is no correlation between the aggregated monthly counts of positive tests results over time and the aggregated monthly number of clinical follow-ups on those patients during the same period. So, I am not looking for validation of laboratory test results. I am looking for correlation between two data generating processes over time. Can someone advise if this sounds feasible / appropriate and suggest a strategy to undertake the analysis. Any point me in the direction of a useful resource. Thanks.


title: Dealing with Big Data and Lots of Variables

text: What is a good technique to use on data that has many categorical variables with many possible values? For example, let's say you are trying to determine what kind of people are more likely to purchase again from your online store and you have E-mail, Country, Browser. Each variable may have 10+ possible values (e.g. Email: yahoo, gmail, hotmail; Country: USA, Canada, Australia, etc.). Plus you also have continuous variables such as the shoppers age and how much they have spent so far. I have tried using a logit regression but with so many categorical variables, it gets too big and unwieldy. A multiple regression also has the same problems (too many dummy variables). A decision tree seems to work best, but it requires turning the continuous variables into categorical variables. Just wondering what kind of solutions/techniques people have applied to similar situations. Note: I can use R and other stat software.


title: Build corpus with phrases

text: I have my documents as:                doc1 = beautifull, very good, very bad, you are great      doc2 = very bad, good restaurent, nice place to visit      I want to make my corpus separated with `,` so that my final `DocumentTermMatrix` becomes:                     terms      docs       very good      very bad        you are great   good restaurent   nice place to visit       doc1       tf-idf          tf-idf         tf-idf          0                    0       doc2       0                tf-idf         0                tf-idf             tf-idf      I know, how to calculate `DocumentTermMatrix` of individual words but don't know how to make the corpus `separated for each phrase` in R. Though I would prefer a solution in `R` but solution in `Python` is also welcome. What I have tried is:               > library(tm)     > library(RWeka)     > BigramTokenizer <- function(x) NGramTokenizer(x, Weka_control(min = 1, max = 3))     > options(mc.cores=1)     > texts <- c("very good, very bad, you are great","very bad, good restaurent, nice place to visit")     > corpus <- Corpus(VectorSource(texts))     > a <- TermDocumentMatrix(corpus, control = list(tokenize = BigramTokenizer))     > as.matrix(a)      I am getting:                                        Docs       Terms                   1 2        are                   1 0        are great             1 0        bad                   1 1        bad good              0 1        bad good restaurent   0 1        bad you               1 0        bad you are           1 0        beautiful             1 0        .......      What I want is not combination of words but only the phrases that I showed in my matrix.


title: Comparing proportions from the (somewhat) same sample

text: I’m looking for some help on testing whether there are statistical differences between proportions when the data is dependent. By dependent, I mean that all the data is from the same participants though I’m not talking about a repeated measures experimental design but instead just different questions from the same participants. Usually, I’d go straight for the McNemar test. However, the problem I’m facing is that... although the data is all from the same group (1000 people), there are slight differences in the numbers of people answering the questions that I’m using to calculate the proportions. I asked all 1000 people whether they regularly used Windows (all of them said yes; i.e. 1000 said yes), Mac OSX (500 said yes), a Linux distribution (100 said yes) and a variety of other operating systems (10 said yes). Participants could combine their answers - they may have used both Windows, Mac OSX… or just Windows…. or Windows and Linux, or all three, etc. For each operating system that a person indicated using, I asked them whether they used certain software on those operating systems (for instance, whether they used a Word processor). So if a person used OSX and Windows, they’d be asked whether they used a word processor for WIndows and separately asked the same question for OSX. Another participant might use Windows and Linux, and would have been asked whether he/she used a word processor for WIndows and separately asked the same question for Linux. Another person might have said he/she used only Windows and so was only asked about his/her word processor usage on Windows. That is where the problem arises. If I work out that 85% of 1000 Windows users regularly access a Word processors, how do I compare this to 65% of 500 OSX users regularly accessing a Word processor. I had thought about using the McNemar test but I’m not sure how to deal with the discordant pairs - there is no concordant/discordant pair for somebody who uses Windows but doesn’t use OSX, for instance. I had thought about using a Z-test for a proportions but that is for comparing proportions that are strictly independent - I’m fairly sure my data wouldn’t satisfy that condition (but stand to be corrected). Any ideas? Thanks everybody! I hope somebody can help me! Simon :) **Addition of extra information as requested** I've just mocked up some raw data. So this is a representation of the data that I have for each participant (-1 = participants who did not answer that question because they don't use that operating system). ![Mock Participant Data](http://i.stack.imgur.com/hfSKw.jpg) I'm looking to see whether the proportion of people using a word processors in Windows is similar to the proportion of people using a word processor in OSX (where the data is from the same sample 'somewhat' - 'somewhat' given the issue above about participants only answering the questions relevant to them meaning I'll be comparing slightly different carvings of the same sample). So yes, I am trying to compare the proportion of people who use word processors given their OS. However, I'm trying to break it down to compare between individuals operating systems. So in the summary table below... I'm trying to compare whether 85% (Windows) is significantly different from 65% (OSX). The reason why I've been avoiding a standard Chi-square (Pearson's) is that that data is somewhat dependent - the questions do not involve different participants. ![Mock Data Overview/Summary](http://i.stack.imgur.com/MKsH4.jpg) Does that make it slightly clearer? If not then let me know and I can expand :)


title: What is the minimum training set size required for a given number of features for document classification?

text: For document classification problems, is there a rule of thumb for the number of training instances required for the number of terms in the vocabulary? I am using a logistic regression classifier with TF-IDF weighted features. After stop-word filtering, stemming, and filtering by minimum and maximum document frequencies, I have a vocabulary of ~13,000 terms for a training set with ~20,000 documents. I have attempted using LDA for dimensionality reduction by adding topic probabilities as features, but this did not significantly affect performance. The performance of a classifier trained only on LDA topic probability features was inferior to the performance of the classifiers trained on TF-IDF features and TF-IDF+LDA topic probability features.


title: Statistical test to assess differences between two sequencing methods

text: To asses differences between two DNA sequencing methods (454 and MiSeq) I have made artificial mixtures of different bacteria in the lab (mock communities). An Example might be (numbers are percentages of bacterial cells in the sequenced samples):   * Bacterial strain A: 25%   * Bacterial strain B: 25%   * Bacterial strain C: 25%   * Bacterial strain D: 25% I have also more complex mixtures (up to 20 different bacteria and other ratios than 1:1, for exp. logarithmic) If you sequence this you get for example the following results: ## 454:   * A: 20%   * B: 30%   * C: 24%   * D: 26% ## Illumina:   * A: 22%   * B: 28%   * C: 20%   * D: 30% From every sample (like the one above) I have triplicates. What statistical test would you suggest to test if there is a significant difference between the two methods (I guess analysing the difference in ratios is the right thing?). Is Chi-Square the right thing if I have three different samples (One is the one mentioned above)? Would even some bootstrapping approach be an idea? I am familiar with R, so if there is any package you can suggest, that would be nice as well... * * * An image to explain the workflow. I hove this helps to understand my problem better: ![Workflow](http://i.stack.imgur.com/bShvA.png) I also have mock communities that are more complex (10 or 19 bacterial strains in different ratios). But the principle is the same. For the sequencing process you decide how many DNA sequences you want to sequence! The number of reads (Sequences you get) are not (only) a matter of DNA concentration (as long as you have enough DNA to have a "satturation"). So what you can compare at the end between 454 and MiSeq are rather ratios (after clustering) than complete numbers, as far as I understand...


title: Correct methodology to repeat testing of classifier to get good estimate of performance

text: I'm having trouble with a basic machine learning methodology question. I understand the concept of not using the same data to both train and evaluate a classifier, and furthermore when there are parameters in an algorithm to be optimized, you should use an independent third test set to get the final reportable performance figures (e.g. recall rate). However, using a _single_ test set to measure performance seems to be problematic because the measures of performance will likely differ depending on how the data is partitioned into training (plus validation) and test sets, especially for small datasets. It would be better to average the results of N different partitions. For the training stage, this is exactly why people normally use k-fold cross- validation, but if you were to repeat the final test stage with different test sets (i.e. k-fold cross validation), then you are mixing data used for parameter tuning with the test data. This is what all the machine learning references I've seen say not to do. But how to get a good performance measure that is not dependent on partitioning of the data set? * * * ## EDITED Related to cbeleites's comment: I'm not sure whether or not I understand the concept of nested cross validation. Here is a check of my understanding. Please comment if I am not correct. Assume I have 1200 labeled samples, and I want to compare 3 different models (e.g. a decision tree algorithm with three different parameter settings). For each model, I would have two cross-validation loops--an outer loop that performs testing and an inner loop that generates the model and also tests it. For the "outer" loop, if I were to choose 4-fold I first make a 75/25 split to create a training/validation set (900 samples) and test set (300 samples). With the 75% portion, I perform k-fold cross validation (in my example below I choose 3-fold, in order to keep the figure simple, so there are 300 samples per fold) and calculate the average performance of the k-folds. I would then re-train the same model (i.e. same parameter settings) using the full 75% portion and then test with the remaining 25%. I would then repeat this procedure three times using different folds from the outer loop (see figure). The final part is where I'm a bit unsure what is the right procedure. Do I average the four tests from each model (i.e. the folds from the outer loop). What actually is the purpose then of the inner loop? Or am I not doing the nesting correctly at all? ![enter image description here](http://i.stack.imgur.com/oPxJF.png)


title: Unequal variances t-test or U Mann-Whitney test?

text: I have two independent samples and I need to compare their values, to say if values in the first group are higher than values in the second group. I performed the Levene's test and the variances are very unequal. Is it better to use the SPSS version of the t-test (equal variances not assumed) or the U Mann-Whitney test? Do you have any formal reference (a paper I can cite) for the SPSS t-test version when equal variances are not assumed?


title: Multiple linear regression in MATLAB

text: Which is the easier way to perform multiple linear regression in MATLAB given that my dataset consists of 384 explanatory variables and 1 dependent variable? In fact, I need to compute coefficients, corresponding residuals, as well as the generalization error when testing this linear model using, 30 say, examples not in the training set. I think that this function works fine for me, but I am not sure how to make it work with a big number of explanatory variables. I know that there should be a more efficient way to construct such a model, but this is just a preliminary task of my research. All ideas are welcome. EDIT: I think I found it. Could you please confirm that the following is correct?               X = DATA(1:101,1:99);     [M,N] = size(X);     y = DATA(1:101,100);     X = [ones(M,1) X];     b = regress(y,X);      BUT, what about the residuals and the generalization error?


title: MLE/Likelihood of lognormally distributed interval

text: I have a variable set of responses that are expressed as an interval such as the sample below.               > head(left)     [1]  860  516  430 1118  860  602     > head(right)     [1]  946  602  516 1204  946  688      where left is the lower bound and right is the upper bound of the response. I want to estimate the parameters according to the lognormal distribution. For a while when I was trying to calculate the likelihoods directly I was struggling with the fact that since the two bounds are distributed along different set of paramaters, I was getting some negative values like below:               > Pr_high=plnorm(wta_high,meanlog_high,sdlog_high)     > Pr_low=plnorm(wta_low, meanlog_low,sdlog_low)     > Pr=Pr_high-Pr_low     >      > head(Pr)     [1] -0.0079951419  0.0001207749  0.0008002343 -0.0009705125 -0.0079951419 -0.0022395514      I couldn't really figure out how to resolve it and decided to use the mid- point of the interval instead which is a good compromise until I found mledist function which extracts the loglikelihood of an interval response, this is the summary I get:               > mledist(int, distr="lnorm")     $estimate     meanlog     sdlog      6.9092257 0.3120138           $convergence     [1] 0          $loglik     [1] -152.1236          $hessian              meanlog       sdlog     meanlog 570.760358    7.183723     sdlog     7.183723 1112.098031          $optim.function     [1] "optim"          $fix.arg     NULL          Warning messages:     1: In plnorm(q = c(946L, 602L, 516L, 1204L, 946L, 688L, 1376L, 1376L,  :     NaNs produced     2: In plnorm(q = c(860L, 516L, 430L, 1118L, 860L, 602L, 1290L, 1290L,  :     NaNs produced      The parameter values seem to make sense and the loglikelihood is greater than any other method I have used (mid-point distribution or distribution of either one of the bounds). There is a warning message which I don't understand so could anyone tell me if I am doing the right thing and what this message means? Appreciate the help!


title: Ordinal regression: proportional odds assumption

text: I'm analysing survey results with most responses coming in the form of Likert scales. Many of these scales have either very few or 0 responses in the bottom categories. As you can imagine, this is leading to some complications when I try to run my models. I’m still getting encouraging results, but want to be sure that my information is valid before reporting it. Below is a scenario I’ve been working with over the last couple of days: I’m using SPSS to run an ordinal regression with two predictors. In this case, the predictors themselves are actually responses on a Likert scale (but entered into the model as nominal variables). My DV is, of course, also an ordinal scale. My two predictors each have five categories (levels on the scale). My dependent variable also has 5 levels. Just like with my predictor variables, the dependent variable has very few observations in the bottom categories. In fact when I run the regression, it says that 47.5% of cells have frequencies of 0. Yet all my coefficients are significant, the overall model fit (-2 log likelihood) is significant at .000, and the odds ratios (exponentiated form of my coefficients) all seem reasonable. The model looks like a good one other than these cells with frequencies of 0. My test is failing the proportional odds assumption, which says that the coefficients for each predictor category must be equal across all DV levels. I know this based upon the results of the Test of Parallel lines, which SPSS reports as part of the ordinal regression output. So, on the recommendation of an article I found online, I've done two things to explore further. First, I've run separate **logistic** regressions with new dependent variables, each one representing a cutpoint in my original DV— in other words, they indicate whether Y is less than each of my original DV categories (excluding the bottom one). So my new DV's are level 2 or above vs. not; level 3 or above vs. not; etc.. These did not yield significant relationships for most IV-DV combinations (cells). The idea is to compare the odds ratios across the different cutpoints to see if they’re fairly constant. In my case, since few are significant, they’re not. The second thing I’ve done is estimated separate **ordinal** regressions using my original dependent variable— I did one model for each category in my predictors, coded as dummies. So, in 10 separate models (2 predictor variables with 5 categories each), my single predictor would be: a 1 for level 2 and a 0 for all other levels; or a 1 for level 3 and a 0 for all other levels; etc.. For most of these categories, the parallel lines is failed (i.e. the null that the proportional odds assumption is upheld is proven true – a good thing). However, in a couple of these categories, I have no observations (nobody responded Very Poor or Poor on one of my predictor Likert scales). Thus I cannot get a parallel lines P-value for these categories. **My question has two parts**. One is whether it is the bottom levels of the predictor variables that are causing the parallel lines test to fail… and if the reason is that there are no observations in these categories, whether I can still use the overall odd's ratios from my full model. I think this should be no problem since these categories automatically drop out of the model. The second question is whether instead, it might be the low/0 frequencies in the bottom levels of my DV that is causing the parallel lines test to fail. I don’t think it is based on the fact that the test is passed for all predictor variable categories that have observations in them. I have tried combining the bottom categories of my DV, and this decreases the % of cells with frequencies of 0, but does not totally eliminate the problem. Many thanks for taking the time to consider my question. I would be tremendously grateful for any guidance that you can provide.


title: Predicting chemical property (Boiling Point) from a SMILES string

text: I was trying to develop a model for predicting Boiling Points (BP) given a chemical name. One good and unique (ok, almost) way to encode a name is the SMILES notation string. The details of the notation are a bit complex ( see here)                e.g.       Name     Smiles  BP      Benzene c1ccccc1 80      Propane CCC      -42      Ethanol CCO      78      Phenol  c1ccc(cc1)O  181      ....      I was wondering if I might be able to extract features from the Smiles String and use it to predict BP's. A large training data set (~10,000 compounds) is already available. One tricky point may be that the notation is case sensitive and various parenthesis have a topological meaning (branching) in there and also certain symbols (= for a double bond, # for a triple bond) Traditionally the chemical literature uses group contribution methods but those depend on either a human or some smart code to parse the SMILES strings into its constituent entities (e.g. 3 OH groups, 1 double bond etc.) There are good thumbrules already known e.g. Among all alkanes every extra CH2 group adds so many degrees to the BP etc. OR Alcohols are higher boiling due to hydrogen bonding. But I was wondering if it'd be nicer to let the machine learning algorithm figure these out by itself. Any ideas about what methods to use? Or is this too much to ask of an automated algorithm?


title: Dealing with 'Don't Know' answers for a categorical outcome variable

text: I have a survey data with categorical outcome variable (yes, no, don't know) which reflects the acceptance of some situation by respondents. My concern is how to deal with Don't know answers, I really doubt I should drop these observations, because:   1. it shrinks my dataset from around 14400 to around 13000, which is considerable;   2. I have intuition that DK answer carry some info and thus not random. So my questions are:   1. One suggested that non-randomness influences the estimated probability and I should check for it, but how do we check for randomness in Stata?   2. If keeping DK answers is desired then multiple imputation (for example) is the way to deal with this issue. Is there any source/links that I could use to make myself familiar with what multiple imputation is and how it is done in Stata?   3. Almost all papers I read on my topic use logistic regression, I wonder what is the justification behind it. Is there any links/source that compare different probabilistic approaches for not-binary outcome variable (in my case it will be three-asnwers categorical outcome variable) and how we choose between them?


title: Question regarding Bonferroni correction

text: Prove the following version of the Bonferroni inequality- $$P\left(\bigcap_{i=1}^kA_i\right)\ge1-\sum_{i=1}^kP(A_i^c)$$ When creating simultaneous confidence interval, what are $A_i$ and $A_i^c$? Would appreciate any help.


title: Comparing two logit/probit curves with a single parameter

text: I've conducted a psychological experiment on the same subject, under two different condition. For each condition I've collected the number of correct and wrong answer for each stimulus (number of trials per stimulus = 10; number of stimulus = 15). From the data collected I fitted a psychometric curve (probit or logit curve), and now I want to compare the results. Consider **subject 1** : His answers are collected in the following matrices               **COND_1**         correct     cnt   Y  N       1   0 10       2   0 10       3   0 10       4   0 10       5   0 10       6   0 10       7   1  9       8   4  6       9   5  5       10  7  3       11  8  2       12  7  3       13 10  0       14 10  0       15 10  0          **COND_2**         correct     cnt   Y  N       1   0 10       2   2  8       3   2  8       4   3  7       5   2  8       6   4  6       7   8  2       8  10  0       9   8  2       10 10  0       11  9  1       12 10  0       13 10  0       14 10  0       15 10  0          COND_1.1 <- matrix(c(0, 2, 2, 3, 2, 4, 8, 10, 8, 10, 9, 10, 10, 10, 10, 10, 8, 8, 7, 8, 6, 2, 0, 2, 0, 1, 0, 0, 0, 0), byrow=F, ncol=2)      COND_1.2 <- matrix(c(0, 0, 0, 0, 0, 0, 1, 4, 5, 7, 8, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, 6, 5, 3, 2, 3, 0, 0, 0), byrow=F, ncol=2)          cnt <- seq(from=0, to=1.4, by=0.1)          ddprob.1.1 <- glm(COND_1.1 ~ cnt, family = binomial(link = "probit"))      ddprob.2.1 <- glm(COND_2.1 ~ cnt, family = binomial(link = "probit"))      ![enter image description here](http://i.stack.imgur.com/gDGhi.jpg) Observing the plot, you can surely say that performance under condition 2 is better then the performance under condition 1. In fact if you compute some parameters like AUC, p25, p50, p75 (the value of the stimulus at which the number of correct answers are 25%, 50% and 75% each) [threshold parameters], you'll see that:               COND_1      COND_2     AUC    <    AUC     p25    >    p25     p50    >    p50     p75    >    p75      * * * Now consider **subject 2** :               COND_1         correct     cnt   Y  N       1   0 10       2   0 10       3   0 10       4   0 10       5   3  7       6   3  7       7   4  6       8   5  5       9   9  1       10  9  1       11 10  0       12 10  0       13 10  0       14 10  0       15 10  0          COND_2         correct     cnt   Y  N       1   0 10       2   1  9       3   0 10       4   4  6       5   2  8       6   6  4       7   4  6       8   7  3       9   4  6       10  7  3       11  7  3       12 10  0       13  9  1       14 10  0       15 10  0               COND_1.2 <- matrix(c(0, 0, 0, 0, 3, 3, 4, 5, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 7, 7, 6, 5, 1, 1, 0, 0, 0, 0, 0), byrow=F, ncol=2)     COND_2.2 <- matrix(c(0, 1, 0, 4, 2, 6, 4, 7, 4, 7, 7, 10, 9, 10, 10, 10, 9, 10, 6, 8, 4, 6, 3, 6, 3, 3, 0, 1, 0, 0), byrow=F, ncol=2)     cnt <- seq(from=0, to=1.4, by=0.1)          ddprob.1.1 <- glm(COND_1.2 ~ cnt, family = binomial(link = "probit"))     ddprob.2.1 <- glm(COND_2.2 ~ cnt, family = binomial(link = "probit"))      The plotted curves are shown in the image below. ![enter image description here](http://i.stack.imgur.com/7Qkct.jpg) Can you say what performance is better? Under condition 1 or under condition 2? Extrapolating AUC, p25, p50, p75, you have:               COND_1      COND_2     AUC     >   AUC     p25     >   p25     p50     =   p50     p75     <   p75      So my question is: Is there a method to say that a generic curve (logit, probit or whatelse) is HIGHER then another? Is there a method to compare (a single number of "performance") that describe the differences of the two curves? My example shows that AUC, p25, p50, p75 are not good parameters. I'd like to compute a single numeric parameter for each curves, to make a simple t-paired test, extrapolated from the distribution seen.


title: How to calculate the expected number of distinct items when drawing pairs?

text: Suppose I have a set $\mathcal{S}$ of $N$ distinct items. Now consider the set $\mathcal{P}$ of all possible pairs that I can draw from $S$. Naturally, $|\mathcal{P}| = \binom{N}{2}$. Now when I draw $k$ items (pairs) from $\mathcal{P}$ with a uniform distribution, what is the expected number of distinct items from $S$ in those $k$ pairs?


title: Social network dataset

text: I'am looking for social network datasets for determining the strength of (acquaintance, friendship, family ,professional,...) relationships can any one help me to find it thanks in advance


title: What is a good tool and format for representing and communicating the design content of a survey?

text: Are there any decent tools for writing/designing questionnaires before providing to programmers? Currently Microsoft Word is being used and tracking changes and keeping it standardized has become a headache. **Update** : I think I'm being a little misunderstood here. Here's a scenario: A client speaks to a statistician/expert in survey design so they can design a survey to receive feedback from their customers. Once that survey is designed (not in the fashion sense) they might give it to an artist and programmer to actually create the survey (for the web or whatever). I am wondering what some good ways to communicate the logical design of the survey to the programmer might be.


title: Why semi/nonparametric models?

text: Increasing the flexibility of models makes it prone to overfitting. On the other hand, it looks to me that, if the space function classes $\mathcal{F}$ is too big, it is hard to prove bounds on empirical risks bounds and that sort of stuff. That's why I am questioning the necessity/importance/applicability of non-parametric models. Here by nonparametrics I mostly mean, Dirichlet Processes and Beta Processes (and related family). Any comments?


title: Integral in Stata when upper limit is infinity

text: How can we calculate the integral of the integrand with the lower limit 0 and upper limit infinity in Stata? I am aware of the `integ` command, but I am not sure whether I can use that when the upper limit is infinity.


title: Why do they call it "sampling distribution of the sample mean"?

text: Ok, I understand that there is a true population mean and one that I get from the sample. It is different for every sample and, thus, I can build the distribution of the sample means. I arrive at a _distribution of sample means_. But why is it a _sampling_ distribution? Is the whole point that the moniker must be longer than necessary? What do I lose if I omit the extra qualifier and get away with _distribution of sample mean_ alone?


title: Supervised dimensionality reduction-applications

text: What are the applications or advantages of dimension reduction regression (DRR) or supervised dimensionality reduction (SDR) techniques over traditional regression techniques (without any dimensionality reduction)? These class of techniques find a low-dimensional representation of the feature set for the regression problem. At times there are assumptions that the joint distribution of the features $X$ and the response $Y$ lies on a manifold. It makes sense to learn the manifold from the observed sample in this context for solving a regression problem. But the DRR, SDR techniques do not mention any assumption of this kind which makes me wonder why they are useful. Are they just space and time complexity reducing procedures-for regression? If that is the primary advantage-some resources on complexity reduction for high-dimensional datasets when this techniques are put to use would be helpful.


title: How to get started with neural networks

text: I'm completely new to neural networks but highly interested in understanding them. However it's not easy at all to get started.   Could anyone recommend a good book or any other kind of resource? Is there a must-read?   I'm thankful for any kind of tip.


title: How to get 1, 3 and 5 year KM survival for large database using SPSS?

text: I have large database of hundreds of thousands observations. When I used spss to analyze KM survival, it gave me mean and median survivals with 95 % confidence interval. However, I need 1-, 3- and 5- year survival too. It's difficult to figure out using the curves. I am missing something? I also tried life tables, but I am not sure about it's accuracy? Can anyone please throw light here?


title: Factor scale as a dependent variable in repeated measures ANOVA

text: I have to analyse data from an experiment in which participants performed two tasks. Each task had a different content and also a different "space condition". Content, space condition and order of presentation of content and condition were counter-balanced between participants. When participant finished the tasks, they answered a questionnaire about the difficulty of each one, containing a set of 18 questions (for each task). I ran a factor analysis for the 18 questions referred to task with space condition 1, showing that there was just a single component, and another factor analysis with the questions related to space condition 2, that also resulted in only a single component. My question is: Is it correct for me to run a repeated measures ANOVA comparing the two resulting factor scales (DV) as within subject IV ("space condition"), and with "combination of content and space condition" as between subjects IV? Thank you for your help.


title: How to best analyze these data?

text: I'm trying to figure out a way to analyze these data that I have: These are data that I pulled from various studies, and each study has a different sample size, and each study reported an obliteration rate and seizure-free rate. I.e. study A has 20 people with obliteration rate of 70% and seizure-free rate of 65%. Study B has 15 people with obliteration rate of 65% and seizure-free rate of 50%, and so on. I am trying to demonstrate correlation between obliteration rate and seizure- free rate, what is the best way to go about doing this? I was plotting each study in a scatterplot with obliteration/seizure-free rates as the _x_ and _y_ axes and calculating Pearson's $r$... however each study has a different number of people, so should be weighted differently... Please help me with this...


title: Discretization of the split in random forest estimation

text: If the dataset I'm working with has too many independent variable dimensions, I may need to choose a proper subset of all possible locations in my feature vector to determine whether I am going to use it to split a node. This is to manage the computational complexity of the task. I am wondering what results there are in the literature surrounding this, or what practical guidelines I should follow. Should I use a uniform probability mass function to pick random points, or a more complicated approach? Are there any guidelines about the ratio of points I choose to the total number of possible points, and how this ratio will impact model generalisation?


title: How to show operations on two random variables (each Bernoulli) are dependent but not correlated?

text: I was looking at the following question from "One Thousand Exercises in Probability" by Grimmett, page 25, question 16 (not homework just self-study): > Let $X$ and $Y$ be independent Bernoulli random variables with $p = 1/2$. > Show that $X + Y$ and $|X − Y|$ are dependent though uncorrelated. Now, the solution given on page 176 goes as follows: $\mathrm{cov}(X+Y, |X-Y|) = E[(X+Y)\cdot(|X-Y|)] - E[X+Y]\cdot E[|X-Y|]$ I can work out that $E[X+Y] = (0)(1/4) + (1)(1/2) + (2)(1/4) = 1$ and $E[|X-Y|] = (0)(1/2) + (1)(1/2) = 1/2$ But I can't work out how to do $E[(X+Y)\cdot(|X-Y|)]= ?$ The solution in the book for the covariance is: $1/4 + 1/4 - 1(1/2) = 0$ To show dependence, the book provides the following solution: $P(X+Y=0, |X-Y|=0) = 1/4$ is not the same as $P(X+Y=0)\cdot P(|X-Y|=0) = 1/8$ Thank you kindly in advance, CC P.S. First time I've ever used latex, yay! Couldn't work out how to put a table in for the values of the joint mass functions.


title: How to interpret the output of choicemodelr (rhierMnlRwMixture) in R

text: Can someone help me with this one? ## My Problem I just started using the R library `choicemodelr` and succeded in getting some beta values as a solution. But I wonder how do I assign these values to the specific attribute-levels. As a result I only get values for A1B1, A1B2, A1B3,... etc. How does this generic output generally connect to my Design? Didn't find a hint in the documentation. Neither for the `choicemodelr` library, nor the `bayesm` library (`rhierMnlRwMixture`) to which it is connected to. I hope you can help me with this one. Thanks in advance, Phil to illustrate this, some code and output:   ## my design **Attribute**     **Attribute-Levels**   Color          red, blue, yellow, green   shape         square, round   ...                ... ## my data-input               respondent  choice-set  stimulus   Attr. Color  Attr. Shape  Choice     1           1           1          blue         square       1     1           1           2          green        round        0     1           1           3          red          round        0     1           2           1          blue         square       4     1           2           2          green        round        0     1           2           3          red          round        0     1           3           1          yellow       round        3     1           3           2          blue         square       0     1           3           3          green        round        0     1           4           1          red          round        2     ...        ...         ...        ...          ...          ...      ## my code in R               # loading neccesary librarys     library(bayesm)     library(MASS)     library(lattice)     library(Matrix)     library(ChoiceModelR)     library(XLConnect)          #DATENSATZ:     setwd("C:/DATA/CBC/")   # set workingdirectory     .Workbook <- loadWorkbook("DataCBC-R2.xls")     data <- data.frame(readWorksheet(.Workbook, "DataCBC-R"))     remove(.Workbook)          #set parameter for calculation     R = 50000 #Total Iterations of the Markov Chain Monte Carlo     use = 100 #Iterations for Paramerter-Estimation          #Parameter of datainput     none = TRUE #TRUE, if the questionaire has a none-Option but is not coded in the data     xcoding = c(0,0,0,0,0,0,1,1) #0=nominal scale; 1=metric scale          #Parameter dataoutput     save = TRUE #TRUE saves the calculated parameters     keep = 500 #number of random parameter draws to save (thinnig Parameter)          mcmc = list (R=R, use=use)     options = list(none=none, save=save, keep=keep, restart=restart)          #final calculation of the betas     out = choicemodelr(data, xcoding, mcmc = mcmc, demos = demos, options = options, constraints = constraints)         ## I get the following Output (excerpt):               RespID  A1B1     A1B2     A1B3     A1B4      A2B1     A2B2    ...     NONE     001    -2,56    -6,54    -18,49    27,59    -1,74     1,74    ...    -1,94     002    -3,18    -6,52    -19,79    29,49     0,50    -0,50    ...    -0,58            So, which of the A1 betas refeers to what color? Hope that piece of information helps. Let me know, if you need further information.


title: How to calculate the statistical relevance for the hidden nodes?

text: I am implementing ELM (Extreme Learning Algorithm) **in MATLAB**. I separate the training set into two non-overlapping subsets for learning and validation and calculate the hidden layer output matrix H using the learning subset. Now I want to calculate the statistical relevance of the hidden nodes using some method like Chi-squared (w2) , Information Gain etc. So, I need some Matlab code to identify irrelevant/relevant hidden nodes.


title: Multinomial goodness of fit

text: Suppose I generate a multinomial distribution with probabilities $p_i$, i from 1 to k. Now suppose I test it back using goodness of fit (chi square) with the probabilities we know. I read that this works well when the sample size, N is large enough. Questions   1\. What are the factors that influence how large N should be so that the test works well? Hypotheses: k or the way the probabilities are spread out   2\. What if we change one of probabilities by a small amount (say $10^{-3}$). How big N should be to be confident the test rejects the false distribution? Notes: 1.the question also works if instead of generating, we have real data, the sample size would need to be bigger due to noise   2\. this can be seen as a test that the generation was done without error   3\. a program could be exemplary to test all these


title: Within and between group covariates in mixed model

text: I am modeling solute concentration downstream from a source as a function of distance and flow. This was measured during nine unique sampling sessions within the same stream. Concentration was measured at nine evenly spaced distances downstream from the source. For each sampling session, there was a unique stream flow (measured in l/sec). There was also a unique mass for the solute source for each session. This results in three continuous explanatory variables, one within session (distance) and two between session (biomass and flow). So, this gives nine sampling sessions, with nine distance measurements within each sampling session. To model this, I used a mixed effects model allowing for a random intercept for each sampling session and a random slope by distance for each sampling session, each sampling session being treated as a group. The code for the model looks like... Model <\- glmmadmb(Concentration ~ Distance + Biomass + Flow + Biomass*Flow + Distance*Flow + Distance*Biomass + (Distance | Sampling Session), data = Stream, family = "nbinom", zeroInflation = FALSE) The paper containing this analysis was recently reviewed. One of the reviewers said, “The most problematic part of the manuscript is the statistical analyses, which are confusing, unnecessarily complex, and have pseudo- replication issues… the way models are constructed, there is pseudo- replication in that if the random effect is session (which I believe is the case) and there is only one biomass and only one flow measure per session, distance is the only thing that varies. The nested nature of the repeated (and spatially autocorrelated) distances within session while all other variables were collected at the session level violates the assumptions of the model and could have led to erroneous results.” My understanding is that the random effect on the repeated distance measures by session essentially treats the session as the sampling unit and takes care of the pseudoreplication issue. The between groups covariates of biomass and flow are then evaluated at one distance per session (in this case the mean because I z-score standardized these). Am I wrong about this? If not, are there any recommendations on how I could respond to this reviewers comments?


title: Conditional Distribution of Poisson Variables, given $\sum X_i$

text: An exercise asks me to show that if $X_1,X_2,\ldots,X_n$ is a random sample from a Poisson distribution with parameter $\theta$, the conditional distribution of $X_1,X_2,\ldots, X_{n-1}$, given $Y=\sum_{i=1}^n X_i$, is multinomial. Normally, this would have been an easy one but unfortunately the joint distribution here throws me off a bit. This is where I would appreciate some help. Thank you in advance.


title: Calculate residual after SPSS runs ANOVA

text: I have a task something like compute residuals in dollars for one observed person (with gender code, a class they're in, their age, and they make let's say 50,000 a year-the Dependent Var). After SPSS runs and outputs ANOVA with the regression value, the residuals, and the total (and I'm referring to the sum of the squares). Then there is the degrees of freedom, and the mean square. Is is possible from these values to back into the residual value for a given observation? The (square root of the mean-square residuals) divided by the (square root of the mean-square regression) calculates to the same value as the standard error of the estimate where the R-square for that model is .203. So is it correct to say that for that 50,000 income I mentioned above, the residuals are 20.3% or $10,150. My brain is melted here - any help is greatly appreciated. As opposed to Bob John's post ( Finding value of residual ) I don't have the best-fit line formula in this case.


title: Combining normal distributions

text: Imagine that I take two separate measures and I get two separate normal distributions $\mathcal{N}_1(m_1, s_1^2)$ and $\mathcal{N}_2(m_2, s_2^2)$   1. How can I find a single normal distribution $\mathcal{N}_3$ which is a type of average of these two normal curves?   2. As an extension, how could I find a distribution which represents both normals? Meaning, it would be a single distribution and it would have two peaks and would be centered at $\frac{m_1 + m_2}{2}$?


title: Bootstrapping first, then data mine?

text: I previously posted about bootstrapping my data. I'm now wondering about the whole process of testing/data-mining for a strategy and THEN testing on bootstraped data. Does it make sense to bootstrap your data first and then data-mine for the best result? This came to me because my current strategy had significant decay when I tested it against the bootstrapped data. The initial strategy is one that re- balances between the SPY and the TLT on a monthly basis. Back-tested on 6 years worth of data (or 72 total trades), the strategy has a Sharpe of 2.06. But when I test it against the bootstrapped data, the Sharpe ratio drops significantly to 1.04. It seems like I'm ultimately looking to raise the 1.04 number, and to do that I'll need to #1) start the data-mining process over again to find a new back- tested strategy, and *#2) test again against the bootstrapped data. So can I skip step #1? Obviously this will require some significant computational power. Thanks!


title: Financial Time series prediction/ SV Regression

text: I'm working with R software (Lib e1071) and I'm trying to get predictions using Support Vector Regression. The way I'm doing it is the following: I'm windowizing the raw closing prices using N=3 yields:               s[t-3]    s[t-2]     s[t-1]   ->   s[t]            1.2350    1.2358     1.2354        1.2360            .         .          .             .            .         .          .             .    etc...      The value I want to predict is y=s[t]. The SVM Type is "eps-regression" and the Kernel is "radial". Also, I perform a 10-fold cross validation to get the best parameters, gamma and cost.      But I have a problem:   The prediction is always a very close value of the last value s[k-1], for instance:      Last Vector:                 s[t-3]    s[t-2]     s[t-1]   ->   s[t]            1.2350    1.2358     1.2354        1.2355            The predicted values are going to be very close values of the last s[t-1] value. I have tried increasing the number of learned vectors (10K) and increasing the N yields (up to 7), but the results are the same.      Can anybody tell me why is this happening and how could I get real predictions?


title: Statistical method to map a set of values to another

text: I am completely newbie to stats. And I want to know which statistical method would be applicable for following scenario. I have a bunch of sentences extracted from a text document. Each sentence has a specific integer value say X which ranges from 1-5. Each sentence also have another logarithmic value Y which can be positive or negative. There is not relation (direct or inverse proportion) between X and Y. I would like to decide in a mapping of each Y value to a specific X value using those bunch of statements. e.g. Lets say I have 12 sentences with X and Y values as follows:                                X    |    Y          Sentence 1  :    1    |  0.00     Sentence 2  :    1    |  1.05     Sentence 3  :    2    |  0.00     Sentence 4  :    3    |  12.0     Sentence 5  :    5    |  2.0     Sentence 6  :    4    |  -1.0     Sentence 7  :    5    |  -9.0     Sentence 8  :    2    |  0.0     Sentence 9  :    3    |  15.0     Sentence 10 :    1    |  0.0      Then I want to decide what ranges of Y would fit in for 1,2,3,4,5 values of X. I guess some form of regression would be right option here. But I am not sure the exact method will be applicable here.


title: Dinamic changing of parameter in a formula during optimization

text: I have two functions, the first one calulates the result of a generalized linear model for a series of values (`Dose`) respect to a given series of `outcome`:                   glm.Vdose <- function(dvh.matrix, InputDose, outcome) {           Dose<-Vdose(dvh.matrix=dvh.matrix, Dose=InputDose)           model <- glm(formula=outcome ~ Dose, family=binomial(link="logit"))           P<-summary(model)$coefficients[, "Pr(>|z|)"][2]           return(P*1e12)         }      `Vdose` is another function that extract from some input matrices the value of `Dose` that must be used for modeling. The output of this function is the P-value of the model. I need to set the value of input parameter `Dose`in order to achieve the lowest P-value as possible and so I decided to create a function that finds the value of the `Dose` resulting in the lowest value of P-value achieved by `glm.Vdose`.                   fit.Vdose <- function(dvh.matrix, outcome, start.Vdose=20) {           out<-nlminb(start=start.Vdose, objective=glm.Vdose, dvh.matrix=dvh.matrix, outcome=outcome)           return(list("Vdose"=out$par, "P"=out$objective*1e-12))         }      Now I would like to introduce also different parameters in the `glm.Vdose` function, so allowing to get a 'multivariate' approach in giving the final optimized model regard to the value of `Dose` parameter. So I changed the two functions as follow:                   glm.Vdose <- function(dvh.matrix, InputDose, formula) {           D<-Vdose(dvh.matrix=dvh.matrix, Dose=InputDose)           f <- update.formula(formula, ~ . + D)             model <- glm(formula=f, family=binomial(link="logit"))           P<-summary(model)$coefficients[, "Pr(>|z|)"][2]           return(P*1e12)         }              fit.Vdose <- function(dvh.matrix, formula, start.Vdose=20) {           formula <- as.formula(formula)           out<-nlminb(start=start.Vdose, objective=glm.Vdose, dvh.matrix=dvh.matrix, formula=formula)           return(list("Vdose"=out$par, "P"=out$objective*1e-12))         }      I'm trying to pass the value of `formula` in order to allow to introduce more parameters to be fitted. The first two functions run correctly, the second two give an error that is:                   Error in model.frame.default(formula = f, drop.unused.levels = TRUE) :            invalid type (closure) for variable 'D'      I don't know how can I correct this bug. Anybody knows?


title: Modeling a 1D random walk with nonconstant probability in a point

text: There is a 1D discrete random walk system where the probabilities at all points are $\frac{1}{2}$ (probability of going forward and backward) except one point at on $l$ (which is an integer). The probability of going right is $p$ at this point ($p$ isn't necessarily $\frac{1}{2}$). There is just one number where the probability of going right and left is different from the other numbers.   * How can I model this process?   * For fixed $N$ (number of steps) how can I find the probability of reaching a specified point by starting at $0$? For example, with $N=10$, $l=3$, and $p=1/3$, what is the probability of ending at (say) $4$ after $10$ steps?


title: K-Fold Cross Validation for mixed-effect models: how to score them?

text: I'm trying to use k-fold cross validation for model selection for a mixed- effect model (fitted with the `lme` function). But, what exactly do I use as the score for each fold? Presumably I don't just fit each candidate model to the validation subset, calculating new coefficients based on the new data. If I understand correctly, I'm supposed to score the models according to how well a model with coefficients calculated using the training data fits the validation data. But how does one calculate AIC, BIC, logLik, adjR^2, etc on an artificial model that gets its coefficients from one source and its data from another? With so many people advocating cross-validation, I thought there would be more information and code available for calculating the scores by which models will be compared. I can't be the first one trying to cross-validate `lme` fits in R, yet I see absolutely nothing about what to use as the score... how does everyone else do this? What am I overlooking?


title: Show power divergence statistic converges to LRT statistic

text: This is a homework problem out of Agresti's Analysis of Categorical Data (1.34b) For counts $\\{n_j\\}$ the power divergence statistic for testing goodness of fit is $$ \dfrac{2}{\lambda(\lambda+1)}\sum_{i=1}^m n_i[(n_i/\hat{\mu}_i)^\lambda-1] \quad \text{for }-\infty < \lambda < \infty$$ As $\lambda$ approaches 0, show that it converges to the LRT statistic $G^2 = 2\sum_{i=1}^m n_i\log{\left(\dfrac{n_i}{\hat{\mu}_i}\right)}$. I have no idea where to start. One of the original papers that discussed the Power Divergence Statistic simply states this as a fact (http://pages.stern.nyu.edu/~dbackus/Exotic/1Robustness/CressieRead%20JRSSB%2084.pdf) Any clues?


title: Nonparametric equivalent of ANCOVA for continuous dependent variables

text: I have an independent categorical variable ($X$ with two categories, $x_{1}$ and $x_{2}$) and two continuous dependent variables ($y$ and $z$). Using a Mann Whitney test, I know that $y$ is significantly associated with $x_{1}$ and $z$ is likewise significantly associated with $x_{1}$. However, it could be that either $y$ confounds the relationship seen between $x_{1}$ and $z$, or _vice versa_ , i.e. $z$ confounds the relationship seen between $x_{1}$ and $y$. What distribution-free tests can I use to try account for each factor in tests of $y$ _versus_ $X$ and $z$ _versus_ $X$? How can I achieve this in R and SPSS?


title: Can anyone explain “cond_indep_fisher_z.m” implementaion and application in Bayesian Toolbox?

text: I've just drawn some samples from the ASIA Bayesian network using BNToolbox. Now I need to test the conditional independence between two arbitrary random variables given the third one but there is no guide how to do it according to available samples. I noticed a function exists in BNToolbox named cond_indep_fisher_z which do my desired task. There is its implementation. I would be thankful to indicate me how to use it. I just have five random variables and 1000 samples at hand.               function [CI, r, p] = cond_indep_fisher_z(X, Y, S, C, N, alpha)     % COND_INDEP_FISHER_Z Test if X indep Y given Z using Fisher's Z test'     % CI = cond_indep_fisher_z(X, Y, S, C, N, alpha)     %     % C is the covariance (or correlation) matrix     % N is the sample size     % alpha is the significance level (default: 0.05)     %     % See p133 of T. Anderson, "An Intro. to Multivariate Statistical Analysis", 1984          if nargin < 6, alpha = 0.05; end          r = partial_corr_coef(C, X, Y, S);     z = 0.5*log( (1+r)/(1-r) );     z0 = 0;     W = sqrt(N - length(S) - 3)*(z-z0); % W ~ N(0,1)     cutoff = norminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95     %cutoff = mynorminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95     if abs(W) < cutoff       CI = 1;     else % reject the null hypothesis that rho = 0       CI = 0;     end       p = normcdf(W);


title: Taking logarithm transformation of variable for a regression

text: I have read other links on here about when to take the logarithm of a dependent (or independent) variable. I understand that the log transformation gives a % interpretation rather than a levels interpretation, which might be more desirable because a 1 unit increase means different things depending on the initial starting point. I also understand that if a variable is log- normally distributed like income, then taking the log makes sense. Also, if the theory says log is appropriate, sure. But, what about variables like population density in a region or the child- teacher ratio for each school district or the number of homicides per 1000 in the population. I have seen professors take the log of these variables. It is not clear to me why. For example, isn't the homicide rate already a percentage? The log would the the percentage change of the rate? Why would the log of child-teacher ratio be preferred? Suppose there are no theoretical models supporting that the relationships should be log. Is it the case that we view a class-size of say 20 students so differently from say 25 students that a 1-student increase in class sizes means entirely different things depending on if we are 20 or 25? The same with the homicide rate? A homicide rate of 1 murders per 1000 means something so different from 10 per 1000 that an increase of 1 has very different meanings?


title: How do you determine whether a time series is stationary with reasonable confidence?

text: There are many time series which are obviously not stationary (I mean, by definition), some of them are even determinstic, yet the unity root tests like ADF/PP test, etc., happily reject the existence of an unity root in these processes. So my question is what kind of methods you used to determine whether a process is a truly stationly stochastic process with reasonable confidence?


title: Generating random variable from density function

text: How can I generate a random variable of size n= 2914 if I have the density function?. So the problem is that I have the density f(x) (function well defined)               P<-function(a,e) { ( (1/6)*(1^3) )-((a/2)*(1^2)) +(((((a)^2)/2)+e)*1)}          D<-function(u,mu,sigma) {dlogis(u,mu,sigma)}          K<- function(u,a,e) {(((1/2)*(u^2))- (a*u) +(((a^2)/2)+e))}          H<-function(u,mu,sigma){ plogis(u,mu,sigma, lower.tail = TRUE)}          Fprim<- function(u,a,e,mu,sigma) (1/P(a,e))*(D(u,mu,sigma))*(K(H(u,mu,sigma),a,e))          Fprim(1,a,e,mu,sigma)           df<- function(u)  Fprim(u,a,e,mu,sigma)          #### Parameter n,a,e,mu,sigma      n<-2914     mu<- -0.42155226     sigma<- 0.60665552     a<- 0.43218138     e<- 0.02149706      I think I need to reverse and to use Monte Carlo, I don't know how to do?


title: Uniqueness of $x'\beta$ even when $\mathbb{E}(x^Tx)$ is not invertible

text: As discussed in user25658's answer to this question, when one wants to compute $$ \beta = \mathbb{E}(x^Tx)^{-1} \mathbb{E}(x^TY) $$ but $\mathbb{E}(x^Tx)$ is not invertible, $\beta$ is not uniquely identified. A characterization of all possible values of $\beta$ is given by $$\beta= \mathbb{E}(x^Tx)^g\mathbb{E}(x^Ty)$$ where $\mathbb{E}(x^Tx)^g$ is the generalized inverse of $\mathbb{E}(x^Tx)$. In Hansen, Econometrics, pg 34, section 2.18 (as of today's version on Hansen's webiste) it is written that $x^T\beta = x^T \mathbb{E}(x^Tx)^g\mathbb{E}(x^Ty) $ is however uniquely identified. This result is not completely counterintuitive to me. I understand that that this is matrix multiplication and that the fact that for every $x$ there may be multiple $\beta(x)$, say $\beta(x)_1 \neq \beta(x)_2$, does not imply $x^T\beta(x)_1 \neq x^T\beta(x)_2$. Nevertheless, I have been unsuccessful in trying to prove the claim. Could someone help me with that?


title: Proof for bivariate conditional mean of Gaussian dist

text: I see that a lot of questions are answered here for multivariate and bivariate conditional distributions. But I did not find the proof of these equations (I need just for bivariate case). ![enter image description here](http://i.stack.imgur.com/6SCb1.png) ![enter image description here](http://i.stack.imgur.com/sxyJi.png) to get this ![enter image description here](http://i.stack.imgur.com/TqUUf.png) I appreciate your help


title: Logistic regression vs. estimating $\beta$ using linear regression and applying the inverse-logit function

text: If I want to perform a logistic regression analysis, can I use the linear regression model to compute the values of the parameters and then apply the formula: $$\frac{1} {1 + \exp-(\beta_0+x_1\beta_1+\ldots +x_i\beta_i)}?$$


title: Degrees of freedom for standard deviation of sample

text: would someone please explain why the degrees of freedom for a random sample is _n-1_ instead of _n_ ? I'm looking for an explanation that is intuitive and easily understood by a high school student. ![http://www.statsdirect.com/help/image/stat0019_wmf.gif](http://i.stack.imgur.com/xb1sb.gif)


title: Standard deviation and n-1 : degrees of freedom or unbiased?

text: I am trying to understand if the division for n-1 in the calculation of standard deviation related to degrees of freedom, unbiasing of the estimator, or both. I understand what is a degree of freedom when we talk about t, chi, f distribution and so on. I was also trying to understand what is a degree of freedom in respect of the standard deviation. I became aware of another explanation, which relates to the bias of the estimation of $\sigma$. Now I am confused: why do we divide by $n-1$? Is it because of the degrees of freedom (and in this case, what does it actually mean?) or is it because we want to have an unbiased estimator? Or are the two explanations connected in some way? (Is it normal that the concept of DOF is so slippery and difficult to understand? Or is that just me?) \--I have read the suggested topic, and I think that it only apparently answered my question. The OP of that topic asked about degrees of freedom, and the answer talked about biased estimator. I don't understand how they relate to each other.


title: Non-random small sample - What conclusions can I get?

text: For my thesis I developed a psychotherapeutic method to address social competence in aggressive/withdrawn children. I have obtained a sample of 10 subjects (not random, the professors chose the children by their behavior in the classroom) which I divided in two groups of 5 each (again not random, I couldn't choose, the children were from different schools and I couldn't, by any means, do a random selection).One of the groups were designed to be the experimental and the other control. I have applied to each children, parent and teacher a questionnaire to assess social problems before the application of the method and in the end. Now my question is: My objective is not to generalize the results because this is an exploration study but to compare the results of each group in the two moments. Do the sampling problem and the little size of the sample compromise any results or can I do any tests? Is it possible to use parametric tests? Thanks in advance


title: Incremental IDF (Inverse Document Frequency)

text: In a text mining application, one simple approach is to use the $tf-idf$ heuristic to create vectors as compact sparse representations of the documents. This is fine for the batch setting, where the whole corpus is known a-priori, as the $idf$ requires the whole corpus $$ \mathrm{idf}(t) = \log \frac{|D|}{|\\{d: t \in d\\}|} $$ where $t$ is a term, $d$ is a document, $D$ is the document corpus, and $T$ (not shown) is the dictionary. However typically new documents are received over time. One option is to keep using the existing $idf$ until a certain number of new documents have been received, and the recalculate it. However this seems rather inefficient. Does anyone know of an incremental update scheme that (possibly approximately) converges to the value if all the data were seen in advance? Or alternatively is there another measure which captures the same notion but can be computed in an incremental fashion? There is also a related question of whether the $idf$ remains a good measure over time. Since the idf captures the notion of the corpus word frequency, it is conceivable that older documents in the corpus (say for example, that my corpus includes over 100 years of journal articles), as the frequencies of different words change over time. In this case it might actually be sensible to throw out older documents when new ones come in, in effect using a sliding window $idf$. Conceivably, one could also store all previous $idf$ vectors as new ones are calculated, and then if we wanted to retrieve documents from say 1920-1930, we could use the $idf$ calculated from documents in that date range. Does this approach make sense? Edit: There is a separate but related issue about the dictionary $T$. As time evolves, there will be new dictionary terms that didn't appear before, so $|T|$ will need to grow, and hence the length of the $idf$ vector. It seems like this wouldn't be a problem, as zeros could be appended to old $idf$ vectors.


title: Do these residual plots indicate that my least squares regression coefficient estimates may be biased?

text: Lets say I have a linear regression: $$y \sim 1 + x_1+x_2$$ where the range of $x_2$ is $[0,10]$. I fit this model using `lm` or `rlm` with regression weights in `R`. I collect the residuals and plot them against $x_1$, I found that the residuals show a pattern with respect to the variable $x_1$. The $R^2$ of regressing the residuals onto $x_1$ is $20\%$. Is that possible? What could be the causes? * * * After the same linear regression as above, if I take a smaller portion of the data, say all the data with $x_2<6$. Then I collect the residuals and $x_1$ of this subset and plot the subsetted residuals against the subsetted $x_1$. I found that the residuals still show a pattern with respect to the variable $x_1$. (The two 20% above are just for example... they are not related ... and maybe there is a theory saying that one should be definitely larger than the other, etc. ) Is that possible? What could be the causes? * * * **Edit:** Let me try to describe the shape of the pattern. Lets say the range of $x_1$ is $[0, 100]$.   * At around $x_1=1$, the residuals are in a vertical band of $[-0.1, 0.1]$.   * At around $x_1=10$, the residuals are in a vertical band of $[-1, 1]$.   * ...   * At around $x_1=100$, the residuals are in a vertical band of $[-10, 10]$. I intentionally put these numbers so you see the upper-band and lower-band are growing somewhat linearly as $x_1$ increases. I know this is heteroskedasiticity. But I guess since I am concerned about "bias", not inference... so I don't worry about the heteroskedasiticity...


title: Fitting log-normal distribution in R vs. SciPy

text: I've fitted a lognormal model using R with a set of data. The resulting parameters were:               meanlog = 4.2991610      sdlog = 0.5511349      I'd like to transfer this model to Scipy, which I've never used before. Using Scipy, I was able to get a shape and scale of 1 and 3.1626716539637488e+90 -- very different numbers. I've also tried to use the exp of the meanlog and sdlog but continue to get bizarre graph. I've read every doc I can on scipy and am still confused about what the shape and scale parameters mean in this instance. Would it just make sense to code the function myself? That seems prone to errors though, as I'm new to scipy. SCIPY Lognormal (BLUE) vs. R Lognormal (RED): ![Scipy Lognormal \(BLUE\) vs. R Lognormal \(RED\)](http://i.stack.imgur.com/ORL1y.png) Any thoughts on what direction to take? The data is fit very well with the R model, by the way, so if it looks like something else in Python, feel free to share. Thank you! Update: I'm running Scipy 0.11 Here's a subset of the data. The actual sample is 38k+, with a mean of 81.53627: Subset: > > > x   >  [60, 170, 137, 138, 81, 140, 78, 46, 1, 168, 138, 148, 145, 35, 82, 126, > 66, 147, 88, 106, 80, 54, 83, 13, 102, 54, 134, 34]   >  numpy.mean(x)   >  99.071428571428569   > Alternatively: I am working on a function to capture the pdf:               def lognoral(x, mu, sigma):         a = 1 / (x * sigma * numpy.sqrt(2 * numpy.pi) )         b = - (numpy.log(x) - mu) ^ 2 / (2 * sigma ^ 2)         p = a * numpy.exp(b)         return p      However, this give me numbers the following (I tried several in case I was getting the meaning of sdlog and meanlog mixed up):               >>> lognormal(54,4.2991610, 0.5511349)     0.6994656085799437      >>> lognormal(54,numpy.exp(4.2991610), 0.5511349)     0.9846125119455129     >>> lognormal(54,numpy.exp(4.2991610), numpy.exp(0.5511349))     0.9302407837304372      Any thoughts? Update: rerunning with "UPQuark's" suggestion: > > > shape, loc, scale (1.0, 50.03445923295007, 19.074457156766517) The shape of the graph is very similar, however, with the peak happening around 21.


title: In calculating the F-measure with precision and recall, why is the harmonic mean used?

text: The article for F-measure in Wikipedia says: > The traditional F-measure or balanced F-score (F1 score) is the harmonic > mean of precision and recall: $F_1=2\times\frac{precision \times > recall}{precision+recall}$ Why is the _harmonic mean_ used in particular, and not the arithmetic mean or geometric mean or any other type of averages? What exactly does it mean, to calculate an harmonic mean?


title: Calculating weight as 1/(stnd error) for weighted regression if stnd error = 0

text: I'm assisting a colleague with a weighted regression of average length of stay (LOS) measured in days vs. inpatient admit rate in a dataset consisting of inpatient records from 30 hospitals. We've calculated the weight at each data point as the inverse standard error of patient's LOS at the hospital. Is there a standard procedure for estimating the inverse stnd error weight where the stnd error = 0? In several hospitals with low patient volume the LOS is identical for all patients during our analysis time period, and therefore the stnd error = 0, producing a weight equal to infinity. We could drop these data points from the regression (or avoid a weighted regression entirely), but in principle it seems there ought to be an accepted technique for calculating weights in special cases where the variance = 0. I haven't had any luck checking my stats textbooks. * * * Thank you for the advice, Michael and whuber. Most of the total error sum of squares can be attributed to measurement error in my case when I run a simple unweighted least squares regression (RSS=44.5, ESS=168.9, TSS=213.4). So if I were to construct my own weighting scheme it might entail, at one extreme if variance=0 then the weight = # obs in that hospital, and at the other extreme, if variance=infinity, the weight=0. Perhaps a handy formula could be weight_i = N_i/(N_i^CV_i), where weight_i = weight for hospital i, N_i = # obs for hospital i, and CV_i = Coefficient of Variation of observed LOS for hospital i? ![ALOS vs. Admit Rate](http://i.stack.imgur.com/N0yKb.png)


title: How to treat holidays when working with time series data?

text: I am analysing hourly demand data for electricity. To make my forecasts more accurate, how should I treat the national holidays in the data? In particular, how should I treat them in R?


title: Increasing number of features results in accuracy drop but prec/recall increase

text: I am new to Machine Learning. At the moment I am using a Naive Bayes (NB) classifier to classify small texts in 3 classes as positive, negative or neutral, using NLTK and python. After conducting some tests, with a dataset composed of 300,000 instances (16,924 positives 7,477 negatives and 275,599 neutrals) I found that when I increase the number of features, the accuracy goes down but the precision/recall for positive and negative classes goes up. is this a normal behavior for a NB classifier? Can we say that it would be better to use more features? Some data:               Features: 50         Accuracy: 0.88199     F_Measure Class Neutral 0.938299     F_Measure Class Positive 0.195742     F_Measure Class Negative 0.065596          Features: 500        Accuracy: 0.822573     F_Measure Class Neutral 0.904684     F_Measure Class Positive 0.223353     F_Measure Class Negative 0.134942      Thanks in advance... **Edit 2011/11/26** I have tested 3 different feature selection strategies (MAXFREQ, FREQENT, MAXINFOGAIN) with the Naive Bayes classifier. First here are the Accuracy, and F1 Measures per class: ![enter image description here](http://i.stack.imgur.com/VVFgg.png) Then I have plotted the train error and test error with an incremental training set, when using MAXINFOGAIN with the top 100 and the top 1000 features: ![enter image description here](http://i.stack.imgur.com/Nto1W.png) So, it seems to me that **although the highest accuracy is achieved with FREQENT, the best classifier is the one using MAXINFOGAIN, is this right** ? When using the top 100 features we have bias (test error is close to train error) and adding more training examples will not help. To improve this we will need more features. With 1000 features, the bias gets reduced but the error increases...Is this ok ? Should I need to add more features ? I don't really know how to interpret this... Thanks again...


title: testing significance of difference between dependent correlations

text: When testing the correlation between xy and xz, you also need to include the yz in the formula, plus sample size. Well what sample size do you use if these variables are items asked on a questionnaire that only some participants in the sample answered so there is a different N for each correlation?


title: Suggestions on how to visualize some data

text: Basically, I have several design options with 22 variables with a 0-1 option. For each design I have an output which is a real number. So...I would like to know a way to express as clear as possible the effect of all the variables in the output, expressing how well each of the designs did. Any suggestions? Clarification: I have (potentially) 2^22 possible designs correct? And I don't know how the variables correlate. So to get a starting general idea I want to display the results graphically in a way that captures information about each design as whole. Right now I am displaying simply 22 bars showing the variation of the dependent variable when you switch each variable and hold all the rest in a constant "nominal" configuration. This however makes me lose information about how mutually changing variables play a role in the process I am aware this is an oddly open-question, and I was just looking for suggestions that more experienced people in the field might have


title: Can we use CFA (confirmatory factor analysis) to ascertain validity coefficient?

text: Validity coefficient is an important measure for test development especially in psychology and sociology wherein we have too indirect measurements. Perhaps we have different ways to know its size.