title: problems to build the linear multiple regression in R and to explain the output

text: I have a problem to build and to explain the linear multiple regression. I have a data set called `Cars93` with 26 variables (numeric and not numeric) and 93 observations. This data set you can find in the `MASS` R package. I want to build a linear regression model for predicting the price of a car. Then I have to do a variable selection (forward and backward stepwise) using AIC and BIC in R. My knowledge in R is too little that´s why I have some problems. I really hope you can help me! 1) The data set has some missing values I just solved this problem like this:                Cars93 [! Complete.cases (Cars93)]       Cars93new <- na.omit (Cars93)       Cars93 = Cars93new       I think some informations are going lost. Is there another solution to eliminate the missing values? 2) Some variables from the dataset are not numeric I tried to convert these values into numerical values like this:               Cars93 $ airbags = factor (Cars93 $ airbags, labels = c (2,1,0))      Cars93 $ airbags        [1] 0 2 1 2 1 1 1 1 1 1 2 0 1 2 0 1 2 2 1 0 1 1 1 1 0 2 0 0 0 1 1 1 1 0 1 1 2 2      [39] 0 0 0 0 1 1 2 2 2 0 0 1 1 2 1 0 0 1 1 1 1 0 1 1 0 0 0 2 0 2 1 1 0 0 1 0 1 1      [77] 1 0 0 0 1 2      Levels: 2 1 0       I did the same with other not numeric variables. Afterwards I tried to build a linear model regression with all variables:               Modell=lm(Price~Horsepower+EngineSize+MPG.city+MPG.highway+Rev.per.mile+Man.trans.avail+Fuel.tank.capacity+Passengers+Length+Wheelbase+Width+Turn.circle+Weight+Rear.seat.room+Luggage.room+Origin+AirBags+Type+Cylinders+Weight+PRM)     summary(Modell)      But the output does make any sense:               Call:     lm(formula = Price ~ Horsepower + EngineSize + MPG.city + MPG.highway +          Rev.per.mile + Man.trans.avail + Fuel.tank.capacity + Passengers +          Length + Wheelbase + Width + Turn.circle + Weight + Rear.seat.room +          Luggage.room + Origin + AirBags + Type + Cylinders + Weight +          RPM)          Residuals:         Min      1Q  Median      3Q     Max      -9.4893 -2.3664 -0.0062  2.1180 18.1112           Coefficients:                         Estimate Std. Error t value Pr(>|t|)         (Intercept)        81.335018  37.993697   2.141 0.036826 *       Horsepower          0.123535   0.049355   2.503 0.015372 *       EngineSize         -0.615828   3.047223  -0.202 0.840602         MPG.city           -0.392888   0.470385  -0.835 0.407259         MPG.highway         0.013646   0.428978   0.032 0.974740         Rev.per.mile        0.001498   0.002511   0.597 0.553206         Man.trans.availYes -1.600967   2.480497  -0.645 0.521387         Fuel.tank.capacity  0.462731   0.572169   0.809 0.422219         Passengers          0.615593   1.823089   0.338 0.736925         Length              0.074875   0.130511   0.574 0.568547         Wheelbase           0.740146   0.343760   2.153 0.035796 *       Width              -1.745792   0.571082  -3.057 0.003473 **      Turn.circle        -0.695287   0.415708  -1.673 0.100203         Weight             -0.004068   0.006255  -0.650 0.518171         Rear.seat.room      0.101150   0.420050   0.241 0.810619         Luggage.room        0.176183   0.367199   0.480 0.633306         Originnon-USA       1.881047   1.762845   1.067 0.290696         AirBagsDriver only -3.294049   1.888346  -1.744 0.086777 .       AirBagsNone        -8.535307   2.289737  -3.728 0.000464 ***     TypeLarge          -1.692122   3.999146  -0.423 0.673887         TypeMidsize         2.684947   2.639047   1.017 0.313504         TypeSmall           1.913341   2.896592   0.661 0.511710         TypeSporty          4.686129   3.268426   1.434 0.157407         Cylinders4         -3.126727   4.554852  -0.686 0.495360         Cylinders5         -4.732933   7.498898  -0.631 0.530605         Cylinders6          0.224795   5.695793   0.039 0.968664         Cylinders8          4.020677   7.255406   0.554 0.581755         RPM                -0.002778   0.002450  -1.134 0.261805         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 5.009 on 54 degrees of freedom       (11 observations deleted due to missingness)     Multiple R-squared:  0.8313,    Adjusted R-squared:  0.747      F-statistic: 9.859 on 27 and 54 DF,  p-value: 1.014e-12      I have 4 times "Cylinders" and "Type" and twice "AirBag" in the summary. I dont know why... And only 4 variables are significant in the model. Can somebody tell me, where is the mistake in my model? I also would like to know how to test other assumptions in R for multiple linear model.


title: Dealing with missing data due to variable not being measured over initial period of a study

text: I was recently consulting a researcher in the following situation. **Context:**   * data were collected over four years at around 50 participants per year (participants had a specific diagnosed clinical psychology disorder and were difficult to obtain in large numbers); participants were only measured once (i.e., it's not a longitudinal study)   * all participants had the same disorder   * the study involved participants completing a set of 10 psychological scales    * the 10 scales measured various things like symptoms, theorised precursors, and related psychopathology: the measures tended to intercorrelate around $r = .3$ to $.7$.   * in the first year one of the scales was not included   * the researcher wanted to run structural equation modelling on all 10 scales on the entire sample. Thus, there was an issue that around a quarter of the sample had missing data on one scale. **The researcher wanted to know:**   * What is a good strategy for dealing with missing data like this? What tips, references to applied examples, or references to advice regarding best practice would you suggest? I had a few thoughts, but I was keen to hear your suggestions.


title: Getting seRious about time series with R

text: If you think back, to when you first started with time series analysis. What tools, R packages and internet resources do you wish you had known about? What I'm trying to ask is, where should one start? Specifically, are there any resources for R that really boil it down for one who is "new" to time series analysis with R.


title: ROC for more than 2 outcome categories

text: How do you construct ROC Curves when there are more than two outcome categories (in my case, I have four)? I've heard you should do this for the most popular group. Are there any other ideas? Are there functions in R to help with this?


title: Graphical representation of variance

text: I'm learning statistics and I can imagine pretty well what the standard deviation looks like (image here). But, knowing that the standard deviation is the square root of the variance, I just can't figure out what that looks like. Can anybody provide me with an illustration or a plot to help me understand that?


title: Alternative denominator for F test in nested model comparison?

text: In comparing two nested multivariate models, the usual F test statistic is in this form: $\frac{(SSR_2 - SSR_1)/(p_2 - p_1)}{MSE_2} \sim F_{p_2-p_1, N - p_2}$ assuming model 1 is nested within model 2, where $SSR_i$: Residual sum of squares,   $SST$: Total sum of squares,   $p_i$: model dimension;   $MSE_i$: mean square error, calculated as $\frac{SST - SSR_i}{N - p_i}$. I'm wondering if it is valid to change denominator of this statistic to $MSE_1$ (i.e. $MSE$ from reduced model) and use $F_{p_2-p_1, N - p_1}$ instead?


title: Mean of Posterior distribution

text: If $X_1,\ldots,X_n$ be iid $\sim N(\theta,\sigma^2)$, and let $\theta$ has double exponential distribution, $\pi(\theta) =\frac{e^{-|\theta|/a}}{2a}$, with $a>0$ known. Find mean of the posterior distribution. **My Work** We need to find joint distribution of $f(\theta|x)$ , where $x$ is vector $f(\theta|x) = \frac{e^{-|\theta|/a}}{2a}\frac{1}{(2\pi\sigma^2)^{n/2}}. e^{-n(x-\theta)^2/{2\sigma^2}}$ $f(\theta|x) \propto e^{\frac{-n(x-\theta)^2}{2\sigma^2}-\frac{|\theta|}{a}}$ I understand that I need to consider two cases, one for positive $\theta$ and one for negative. But I am finding hard time to go further. Please suggest some feedback or help.


title: Unit root tests: how to decide if to include a trend and/or a constant

text: Applying a test to univariate time series data for checking if the series has a unit root or not, one is faced with a decision if one would like to test if the series is stationary around a constant or a trend. But I am not pretty sure, based on what sign should I infer which test is the reasonable one? Is it correct to look at the plot of the data and see if there is a trend visible, and if yes, then opt for a test of stationarity around a trend? And with the constant the same way: if there is no trend visible then see if the series might have a mean different from zero?


title: Average difference in averages and t-test

text: I'm currently working with some individual level data and have run up against a wall in figuring out the best way to analyse it. Let me describe: We have a sample of N=9000 Physicians. For each Physician, I've counted the number of nights they worked (N-Nights) and the number of Patients they saw (N-Night-Patients). I've arrived at an average Patients-per- Night metric by dividing the N-Night-Patients by N-Nights. For each Night a Physician works, they could have also worked directly into the next day. I've arrived at similar counts and a similar average: N-Day- Patients / N-Days. Note that N-Days will be less than or equal to N-Nights because I'm only interested in the days that followed directly from a Night. My goal is to determine if there is a difference in the average number of patients treated on a night is different from the average number of patients treated on a day. My analysis is at the physician level, where I have the two averages, two numerators and two denimonators in each record. My intuition tells me that a paired t-test might be suitable. My independent variable is Type, where type can be Night or Day.My question concerns whether it's legitimate to use averages as the dependent variable. The results would then be the "average difference in averages". Is this appropriate? If not, are there any ways you could recommend to analyze the data? Is there possibility for regression with weighting? Thanks in advance.


title: Proof for mean and variance of a multidimensional half normal distribution

text: Given a $p$-dimensional random variable $\xi \sim \mathcal{N}\left( 0, I_{p} \right)$ and a $p\times p$ dimensional diagonal matrix $\Lambda$ with diagonal entries $\lambda = \left( \lambda_{1},\cdots,\lambda_{p} \right)$. Can anybody point me to a paper or help me with the proof that: mean: $\mathbb{E} \left(\Lambda \left| \xi \right|\right)$ = $\sqrt{\frac{2}{\pi}}\lambda$ and $\mathbb{V}ar\left( \Lambda \left| \xi \right| \right)$ = $\left( 1 - \frac{2}{\pi}\right) \Lambda^{2}$. On Wikipedia, one can find the one-dimensional results which are $\sqrt{\frac{2}{\pi}}\lambda$ and $\left( 1 - \frac{2}{\pi}\right)\lambda^{2}$, respectively.


title: Principal Components Regression: Transform 95% CIs back to original space

text: I have a set of predictors that clearly suffer from some amount of multicollinearity, so I am using PCA to make the columns of X orthogonal. I am also using this as a way to regularize the subsequent regression by removing components that account for ~0% of the variance. For example, if I run OLS regression on PCA-transformed data that has 8 predictors, I am then able to use the eigenvectors from the original PCA transformation to get the beta weights for the original 12 predictors. So far, so good. However, to be able to evaluate the contributions of these predictors to the model fit, I'd like to transform the 95% confidence intervals back into the original space of the 12 predictors. That way, I can use the overall R^2 and associated p-values for the full model to find regressions that are significant, where specific predictors have non-zero contributions. It is unclear to me how to transform the 95% confidence intervals. If that's not possible, is there another way to evaluate the significance of specific predictors in the original space? Thanks


title: Efficient data flow for R?

text: There are already some good answers on project management (eg. How to efficiently manage a statistical analysis project?). These are great, and make life a lot easier. In particular, the workflow that runs along the lines of   1. Load raw data.   2. Manipulate raw data into useful forms. (Unload raw data to save space).   3. Perform analysis and store results.   4. Make and store figures. One thing that's still hard is to get from raw data to final results and figures, without having to run the whole project from scratch every time something changes. What I would like to be able to do is something like this:   * Attempt to create figure. Is manipulated data is available in memory?      * Yes: create figure.     * No: Does manipulated data exist on disk?        * Yes: Load manipulated data.        * No: Load raw data, manipulate, and save, then unload raw data objects. Also, in each step it would be nice to have some trigger to force a full re- load, if the raw data has been updated. Is there an existing framework in R for doing something like this? Or are there any recommended ways of doing this? As it stands, I often have to run everything from scratch, which can take ages (large data files, complex manipulations), and is a waste of resources


title: Measuring statistical significance of machine learning algorithms comparison

text: Let us consider a comparison of two machine learning algorithms (A and B) on some dataset. Results (root mean squared error) of both algorithms depend on randomly generated initial approximation (parameters). Questions:   1. When I use the same parameters for both algorithms, "usually" A slightly outperforms B. How many different experiments ( _with different parameters /_ updated _/_ ) have I to perform to make "sure" that A is better than B?   2. How to measure significance of my results? (To what extent I am "sure"?) Relevant links are welcome! PS. I've seen papers in which authors use t-test and p-value; but i'm not sure if it is ok to use them in a such situation. **UPDATE.** The problem is that A (almost) always outperforms B if initial params and learning/validation/testing sets are the same; but it doesn't neccessarily hold if they differ. I see the following approaches here:   * split data into disjoint sets D_1, D_2, ...; generate parameters params_1; compare A(params_1, D_2, ...,) and B(params_1, D_2, ...,) on D_1; generate params_2; compare A(params_2, D_1, D_3,...) and B(params_2, D_1, D_3,...) on D_2 and so on. Remember how often A outperforms B.   * split data into disjoint sets D_1, D_2, ...; generate parameters params_1a and params_1b; compare A(params_1a, D_2, ...,) and B(params_1b, D_2, ...,) on D_1; .... Remember how often A outperforms B.   * first, do cross-validation for A. Then, independently, for B. Compare results. Which approach is better? How to find significance of the result in this best case?


title: Histogram question: How do we choose a perfect histogram?

text: 20 students from a school were randomly selected to study the dollar amount that they spend on purchasing soft drink in a typical week. The data collected (in dollars):               $24, $35, $17, $21, $24, $37, $26, $46, $58, $30,     $32, $13, $12, $38, $41, $43, $44, $27, $53, $27      Some graphical representations are created based on the above data Figures 1 to 4 are shown below. Which do you consider an appropriate histogram? You can choose more then 1. ![Histogram.](http://i.stack.imgur.com/4WWWH.png) This is an actual quiz question but I find it to be so poorly worded I cannot understand a single bit. Can anyone help?


title: Set a threshold for a sequence

text: I am now having a sequence of numbers, e.g. $x_1, x_2,\dots , x_N$, ($0\leq x_i\leq 1$), which may represent the correlation between a couple of measurements. I want to set a threshold for this sequence, that for $x_i$ which are greater than the threshold may have a significant difference with $x_j$ which are smaller than the threshold. Or in another form, assuming these values are obtained from a distribution which is unknown, I want to determine a threshold, that $\Pr(x_i > \text{threshold}) = 0.1$, any non-parametric method can determine this threshold?


title: Translate R to C++ (eventually with Rcpp)

text: I would like to learn to use **Rcpp**. I went through the docs on the package's CRAN website, but i feel working on a practical example (second practical, considering convolve3) would be more helpfull . I propose the following code, from the the **robustbase** package because it's neither too long nor too short, uses a combination of R types and R functions and has one of those small arithmetic iterations that are too slow in R. How would you go about Rcpp-ing it ?               scaleTau2<-function (x, c1 = 4.5, c2 = 3, consistency = TRUE, mu.too = FALSE){     n <- length(x)     medx <- median(x)     x. <- abs(x - medx)     sigma0 <- median(x.)     mu <- if (c1 > 0) {         x. <- x./(sigma0 * c1)         w <- 1 - x. * x.         w <- ((abs(w) + w)/2)^2         sum(x * w)/sum(w)     }     else medx     x <- (x - mu)/sigma0     rho <- x^2     rho[rho > c2^2] <- c2^2     if (!identical(consistency, FALSE)) {         Erho <- function(b) 2*((1-b^2)*pnorm(b)-b*dnorm(b)+b^2)-1         Es2 <- function(c2) Erho(c2*qnorm(3/4))         nEs2 <-ifelse(consistency == "finiteSample",n-2,n)*Es2(c2)     }     else nEs2 <- n     c(if (mu.too) mu, sigma0 * sqrt(sum(rho)/nEs2))     }      Please explain as much as you can. **EDIT** It's really the idea of a step by step explanation of how you would go about converting a well written (and documented) R code (so at least the foundations are okay) unto an efficient implementation. The choice of the code is arguable a bit random but i think it reflects the arch-typical script on our blueprints (calls R functions that one doesn't want to translate, uses arithmetic loops....). **EDIT2** from the comments i realize this may actually be a big work to do in C++ (i didn't realize it when posting the code). In regard of this, using individual pieces as pedagogical devices is ok. I'll eventually parse the pieces together by editing the question.


title: Derivation of the formula for partial correlation coefficient of second order

text: I came across this formula in some online resources. $$r_{12.34} = \frac {r_{12.3} - r_{14.3}r_{24.3}}{ \sqrt {(1- r_{14.3}^2 )(1-r_{24.3}^2 }}$$ I can use this but I wanted a proof of the formula. Can you prove the formula or please direct me to any link with a proof?


title: How to handle the difference between the distribution of the test set and the training set?

text: I think one basic assumption of machine learning or parameter estimation is that the unseen data come from the same distribution as the training set. However, in some practical cases, the distribution of the test set will almost be different from the training set. Say for a large scale multi-classification problem that tries to classify product descriptions into about 17,000 classes. The training set will have highly skewed class priors, such that some class might have many training examples, but some might have just a few. Suppose we are given a test set with unknown class labels from a client. We try to classify each product in the test set into one of the 17,000 classes, using the classifier trained on the training set. The test set would probably have skewed class distributions but probably very different from that of the training set, since they might be related to different business areas. If the two class distributions are very different, the trained classifier might not work well in the test set. This seems especially obvious with the Naive Bayes classifier. Is there any principled way to handle the difference between the training set and a particular given test set for probabilistic classifiers? I heard about that "transductive SVM" does similar thing in SVM. Are there similar techniques to learn a classifier that performs best on a particular given test set? Then we can retrain the classifier for different given test sets, as is allowed in this practical scenario.


title: About P-value for ARIMA in R

text: I know I could obtain $p$-value by calculating coef/standard error using R command `pnorm()`. When the $p$-value of each coefficient is less than 0.05 (confidence interval 95%), can I reject those coefficients? For example, given these $p$-value:               ar1: 0.003     ar2: 0.432     ar3: 0.04      do I have to reject `ar2` and run my analysis again? Or, since the last $p$-value of `ar3` is less than 0.05, can I accept all?


title: Estimating joint distributions using copula package in R

text: I am trying to estimate the joint distribution of stock returns using the `copula` package. I have read a couple of papers on copulae, but alas my lack of math knowledge prevents me from understanding beyond the basics if that! I am thinking of estimating the joint distribution using the following steps and have a couple of questions about the process in general:   1. Download stock prices for two companies   2. Convert prices to continuous returns   3. Transform returns to marginal uniform distribution.      1. How do I do this??   4. Create the copula object       1. When creating the Copula object using `ellipCopula()` or `archmCopula()`, I need to specifiy the type of copula and the parameters. How do I know which copula to use? And how do I obtain the parameters... I thought the point of fitting a copula to the data was to obtain the estimated parameters, so why do I need to specify here? Am I just providing initial estimates of the parameters which are then corrected in the actual fitting?    5. Specifiy the bivariate distribution      1. When specifying the bivariate distribution, I need to specify the margins in `mvdc()`, since, according to a number of sources, I should convert the data (returns) to uniform distributions, do I specify the margins as 'unif' in the function?   6. Define the data to be fitted      1. I am assuming the will be a matrix with 2 columns in this case, each column being the uniforms   7. Fit the data As you can see I am very new at this and any help would be greatly appreciated!


title: Cross validation with nonparametric smoothing regressions

text: When I use regression models I like to explore functional relationships using nonparametric smoothing regression (e.g. generalized additive models, lowess, running line smoothers, etc.) before estimating a parametric model using, as appropriate, nonlinear least squares regression to estimate parameters for functions suggested by the nonparametric model. What is a good way to think about performing cross validation in the nonparametric smoothing regression phase of such an approach? I wonder if I might encounter a situation where in random holdout sample A a relationship approximated by a "broken stick" linear hinge function might be evident, while holdout sample B suggests a relationship that would be better approximated by a parabolic threshold hinge function. Would one take a non-exhaustive approach hold back some randomly selected portion of the data, perform the nonparametric regression, interpret plausible functional forms for the result, and repeat this a few (human-manageable) number of times and mentally tally plausible functional forms? Or would one take an exhaustive approach (e.g. LOOCV), and use some algorithm to 'smooth all the smooths' and used that smoothest of smooths to inform plausible functional forms? My applications will typically entail human-manageable numbers of predictor variables (a handful to a few dozen, say), but my sample sizes will range from from a few hundreds to a few hundred thousand. My aim is to produce an intuitively communicated and easily translated model that might be used to make predictions by people with data sets other than mine, and which do not include the outcome variables. References in answers very welcome.


title: Repeated Measures ANOVA pairing counterbalance issue

text: I have been conducting all of my analyses on one DV (Quality) using a repeated measures ANOVA, with two variable (Time and Gender) at two levels each. However, it occurred to me today, that this might be incorrect, and here's the problem: Though Gender is one of my variables, the truth is that participants don't find out gender until Time two. The setup is such that they read two essays (A and B) and evaluate them both (Quality). Then they find out that each of these essays is Male (M) or Female (F) and rate them again for Quality, however, I had counterbalanced the essays so that sometimes they would see AM or BM (or AF / BF). Thus it seems that using a normal repeated measures ANOVA, the pairing would be wrong so the model assumes that A always goes with F and B always goes with M, so it's taking the difference between A and F and analyzing it agains the diff between B and M, correct? MY ROUGH SOLUTIONS: Since it's nearly impossible for me to tell which essay was paired with which gender, my solution was to take the average of essay A and B and make a new variable (AB), and compare F and M to AB. However, it doesn't seem like I can do that using the normal model (at least not in SPSS). Manually, I can do separate paired t-tests, or I could subtract the Grand Mean of (AB) from both F and M and then do a paired T-test on those, but again, that doesn't seem right. Unless, maybe it is valid? Or, would I just put in Time as a 3 level variable, and look at the pairwise contrasts after do the repeated measures without crossing it with Gender? VALID SOLUTION? Is there any way to do a repeated measures ANOVA with all of these variables without having to do manual calculations or series of t-tests in which I would be using the wrong dfs? I appreciate your help! PS: I use SPSS Thanks!


title: Help needed with interpreting mixed-model and factorial ANOVA

text: I really need help with an analysis I am doing. The study involves a comparison of two different treatments for bulimia (weight concern treatment and perfectionism treatment) and a combination of the 2 treatments, as well as a control group over time. In other words, participants were randomly allocated to one of four groups: (a) an active control group which involved group counselling sessions but no targeting of perfectionism or weight concern; (b) a perfectionism intervention group who received a treatment aimed at reducing perfectionism; (c) a weight concern treatment group who received a treatment aimed at reducing perceived weight concern; and (d) a combined treatment program where both perfectionism and weight concern were targeted. I therefore have a factorial design with two independent variables: perfectionism treatment (yes, no) and weight concern treatment (yes, no). The within-subjects factor is Time (i.e. T1 before treatment and T2 post treatment). The dependent variable is disordered eating symptoms (continuous). I initially ran a 3-way mixed-model ANOVA, and found a significant 3-way interaction. ![Mixed model ANOVA multivariate tests](http://i.stack.imgur.com/mIBZT.png) ![Mixed model ANOVA tests of within subjects effects](http://i.stack.imgur.com/nvvn5.png) ![Mixed model ANOVA tests of within subjects contrasts and tests of between subjects effects](http://i.stack.imgur.com/hbLlU.png) To further test this interaction, I then ran a factorial ANOVA. To do this, I computed a new dependent variable which was the change in disordered eating symptoms between T1 and T2 (I did T1 – T2, so a higher positive change score indicates greater improvement in disordered eating symptoms). Here are the results from my factorial ANOVA (using change in disordered eating scores as my DV): ![Factorial ANOVA Descriptive statistics](http://i.stack.imgur.com/dym24.png) ![Factorial ANOVA tests of between subjects effects](http://i.stack.imgur.com/3eSxS.png) ![Interaction line graph](http://i.stack.imgur.com/xTdlI.png) My problems are as follows:   1. I’m not sure how best to present the results (i.e. graphically, table??). I think it would make sense to present the change over time in a graph, a bit like the one I have created, above. I’m not sure if my graph is correct, and if it shows the interaction? The difficulty I’m finding is dealing with the fact that the treatment conditions are not part of one single independent variable with several levels, but rather are separate independent variables. I am used ot seeing them as IV of “treatment group” (perfectionism treatment, weigh concern treatment, combination treatment, control), but instead here we have “perfectionism treatment” (yes, no) and “weight concern treatment” (yes, no)   2. I don’t know whether I need to do further tests, i.e. simple effects analysis following this factorial ANOVA or if that’s it? I tried running 2 t-test but just got more confused    3. I'm also confused by the fact that in the initial mixed-model ANOVA I got a barely significant interaction between the 2 treatments in the Test of Between Subjects Effects box (p = 0.55), however the effect was significant in the factorial ANOVA. Also, the main effect of weight concern treatment was not significant in the Test of Between Subjects Effects box (p = 0.085) in the mixed-model ANOVA, however it was significant in the factorial ANOVA (p < 0.001). Can I conclude that the weight concern treatment IV has a main effect as p is significant in the factorial?? What about the fact that it showed as non-significant in the 3-way?   4. Am I correct in thinking that there wasn’t much difference between the effect of the perfectionism treatment and the combination treatment and the difference in means is small (and the lines on the graph are nearly on top of each other and parallel)? Where can I get the significance measure for this and for the other treatment groups? Any help would be greatly appreciated!


title: What is the equation for combining standardized mean differences across multiple studies?

text: I am trying to come up with a combined/pooled standardized mean difference for meta-analysis of a dozen or so studies that have the same binary predictive groups and the same continuous outcome measure. I am able to retrieve the sample size, the mean, and the standard deviation for the control and experimental group in each study, which several sources such as Chochrane Reviews (http://130.226.106.152/openlearning/html/modA1-3.htm) have said I need -- but I can't find any resources that explain what to do with those values. Here is some fake sample data similar to what I have, in case it is useful for example purposes:               Study, Treatment N, Treatment Mean, Treatment SD, Control N, Control Mean, Control SD     A, 30, 15.5, 3.2, 31, 17.9, 3.0     B, 60, 17.3, 1.1, 60, 18.2, 0.9     C, 40, 15.1, 3.5, 44, 18.0, 2.8


title: Estimating beta-binomial distribution

text: Suppose that I culture cancer cells in _n_ different dishes _g₁_ , _g₂_ , … , _g n_ and observe the number of cells _n i_ in each dish that look different than normal. The total number of cells in dish _g i_ is _t i_. There is individual differences between individual cells, but also differences between the populations in different dishes because each dish has a slightly different temperature, amount of liquid, and so on. I model this as a beta-binomial distribution: _n i_ ~ Binomial( _p i_, _t i_) where _p i_ ~ Beta( _α_ , _β_ ). Given a number of observations of _n i_ and _t i_, how can I estimate _α_ and _β_?


title: Why are small sample sized avoided in meta-analyses?

text: In most meta-analyses that I've read, the authors choose to exclude studies with really small sample sizes (e.g n=10). Why is that? ## Speculation I could speculate that one reason is that since the effect size from a study depends on the means of the groups as well as the pooled standard deviation, like so $$\frac{\bar{X}_1-\bar{X}_2}{SD_{pooled}}$$ and small samples cannot be assumed to achieve enough of a normal distribution among their results, the pooled standard deviation will be terribly off, thereby possibly messing up the overall effect size. Further, some authors argue that studies with large sample sizes, among other things, tend to be of higher quality than studies with smaller sample sizes, and that there might be a publication bias going on, with smaller studies having to have positive results to get published. For example, this is the argument laid forward by Nüesch et al., 2010, although small studies are here defined as $n<100$. However, I guess one could extrapolate from results like these to what to do with even smaller studies. ## References   * Nüesch et al., 2010. Small study effects in meta-analyses of osteoarthritis trials: meta-epidemiological study. _BMJ, 2010 July 16, 10.1136/bmj.c3515._


title: What is a "saturated" model?

text: What is meant when we say we have a saturated model?


title: What is the difference/relationship between method of moments and GMM?

text: Can someone explain to me the difference between method of moments and GMM (general method of moments), their relationship, and when should one or the other be used?


title: can PCA be applied on grid datasets (variables) without observations

text: hopefully there is someone with a GIS bacground that could provide some clairty! i have 15 raster (grid) data layers derived from a DEM using the terrain analysis function in SAGA. naturally these derivations are strongly correleted with the DEM as they were derived from the dem, however i am trying to ascertain which are the most influential variables for optimising a sampling approach for soils surveys. i have read some papers wher guys have used PCA or PFA for assessing the covariance between variables and then reducing the variables accordingly, however, i do not have any observations to start with (as would be the case if you were going to a new study site and wanted to decide on the optimall sampling locations) as i am trying to model optimal the location of siting a survey point using the PCA (ultimately using Conditioned latin hypercube sampling. does this make sense?


title: Reason for not shrink the bias term

text: For linear model, $y=\beta_0+x*\beta+\varepsilon$, the shrinkage term is always like $P(\beta) $. What's the reason we do not shrink the bias term $\beta_0$? Comparatively, should we shrink the bias term in the Neural network model?


title: two inversly correlated indepent variables...how to combine?

text: I have a few instances where two independent variables are almost perfectly correlated inversely usually in such an instance i would remove one in this instance the two variables each comprise:               Variable A: 50% zero values, rest a normalish distribution separate to the zeros;     Variable B: same as A      the zeros in A correspond to non zeros in B and vice versa originally i planned to remove one but i would loose the descriptive power of the non zero values of one of the variables.. if i dropped Var B i would effectively be putting its **distributed** non zero values into the 0 value of Var A. What is the correct way to proceed? I could of course run models on each but i have more than one instance of this


title: How to find an unbiased estimator?

text: Suppose $X_1, X_2, ...,X_n$ are samples from a uniform discrete distribution with probability 1/3 on each of the points $\theta-1, \theta, \theta+1$, where $\theta\in\mathbb{Z}.$ From "Theory of point Estimation" by Lehman and Casella. Page130 Problem 1.9 I want to find all unbiased estimators of zero. How should I do that? My understanding:   The definition of unbiased estimator says, $E[\delta(x)]=g(\theta)$. which means, I want anything that has $bias=E[\delta(x)]-g(\theta)=0$ But how exactly do I find this guy, $\delta(x)$? And what is the parameter of interest here? My work:   $f(x;3)= \frac{1}{3}, \text{for } \theta-1\le x \le \theta+1, \text{where }\theta\in\mathbb{Z}$ Could someone point out what I am missing here? I feel that I am pretty close to the answer.   Thanks in advance!


title: Properties of conditional probability distributions

text: This is a problem from a machine learning pset that I'm self-learning from http://www.seas.harvard.edu/courses/cs281/assignment-1.pdf. > Suppose we are provided with a hierarchy of three distributions: > $p(\alpha)$, $p(\theta | \alpha)$ and $p(y | \theta, \alpha)$. Write > expressions to compute the following related distributions: > > $p(y|\theta)$, $p(y|\alpha)$, $p(y)$, $p(\theta|y,\alpha)$, > $p(\alpha|y,\theta)$, $p(\theta)$, $p(\theta|y)$, $p(\alpha|y)$ I see that I integrate over the joint distribution to get the marginal probability $$p(y|\theta) = \int p(y|\theta,\alpha)p(\alpha)d\alpha$$ This gives me the first expression in terms of the givens. However, I don't know how to get the others. I know I can set up various equations using Bayes' Theorem $$p(\theta|\alpha) = p(\alpha|\theta) * p(\theta)/p(\alpha)$$ But the above equation has two unknowns ($p(\theta)$ and $p(\alpha|\theta)$). Is there another property that I can exploit besides marginalizing the joint distribution and utilizing Bayes' Theorem?


title: Validating a logistic regression for a specific $x$

text: I have a logistic regression model for 0/1 binary response data that is built from samples $(x_1,Y_1),\dots,(x_m,Y_m)$, where $x_1,\dots,x_m$ are, fixed, nonrandom, real values and $Y_1,\dots,Y_m$ are the corresponding binary response values. So, my model is: $$\log\bigg(\frac{\text{P}(Y=1|x)}{1-\text{P}(Y=1|x)}\bigg) = \alpha + \beta x.$$ Using our data, we estimate $\alpha$ and $\beta$ to be $\hat{\alpha}$ and $\hat{\beta}$. Now, someone wants to test my model. They bring in $n$ new values of $Y$, say, $\hat Y'_1,...,\hat Y'_n$, all corresponding to a particular value of $x$, say $x = x_0$. They claim that if my model is correct, the proportion of $\hat Y'_i$'s that are equal to 1 must be close to $$\frac{1}{1+\exp(\hat{\alpha} + \hat{\beta} x_0)},$$ (say within a 95% confidence interval for this quantity, which we can obtain using the asymptotic normality of $\hat{\alpha}$ and $\hat{\beta}$ and their asymptotic covarriance matrix). But we know that even if the model is perfectly correct and $$\text{P}(Y = 1 | x_0) = \frac{1}{\exp(\hat{\alpha} + \hat{\beta} x_0)},$$ the rate of convergence of $(\hat Y'_1+\ldots+\hat Y'_n)/n$ to $\text{P}(Y = 1|x_0)$ is of order $1/\sqrt{n}$ (please correct me if I'm wrong on this). Therefore, even a sample of size 100 (which would be considered quite large) isn't really large enough for the sample proportion to be close to its presumed expected value that the logistic regression model claims. How do we really test the claim that the model is correct using this newly obtained data? I've looked around quite a bit, and the Bayesian approach seems to be to use simulation to create a prediction interval for the sample proportion of a sample of size $n$ (see for instance chapter 7 in Gelman and Hill's regression book), but I haven't been able to find anything in the frequentist world to address this question. I'd very much appreciate any comments, suggestions, or pointers to the literature.


title: Linear regression: explanation of slope constraints

text: I am given regression equations, one showing linear regression of x on y, and the other y on x. Both refer to the **same set of data-points**.               4X - 5Y + 33 = 0     20X - 9Y - 107 = 0        * Taking the first to be X (dependent) on Y (independent) => of the form `X = a + bY`. I end up with slope value of 1.25    * Taking the second to be Y (dependent) on X (independent) => of the form `Y = a + bX`. I end up with slope value of 2.22 Now, these values **somehow imply** that the system of regression lines is _invalid_? How is that? I am unable to visualize this. The book states the "rule" being that:   1. Both coefficients (slopes) must be less than 1   2. Both coefficients (slopes) must be of the same sign Are these correct? What other constraints exist for a set of regression equations to be valid?


title: Calculation of a time until the absorbing state given chemicals are initially separate

text: Biochemist from Siberia has just patented a new super genome formed by amalgamating genome X and genome Y. These genomes are harmless when separate, but when mixed there is a 45% chance of an unwanted side-reaction and a 10% chance of no reaction at all. Calculate the time until the absorbing state given that the chemicals are initially seperate. The sketch is as follows: ![image](http://s9.postimg.org/opq36p45r/1task.jpg) The time before reaching the absorbing state is also called the mean time to failure (MTTF). This is the “Laplace transform” of the rate of the disaster. I would like to understand this, not to be solved, so please give me only directions where to go.


title: How to test group differences when neither parametric nor nonparametric assumptions are met?

text: Consider following example. We asked 50 men and 50 women one question and they answered on 5 point likert scale (1 is completely disagree and 5 completely agree). Suppose that 50 men responded 20 times 1, 3 times 2, 4 times 3, 20 times 4 and 3 times 5 and 50 women responded 15 times 1, 15 times 2, 15 times 3, 2 times 4 and 3 times 5. I would like to test whether mean (or median) answers in both groups are significatly different. The problem is that I can not use parametric (t-test) neither nonparametric (Mann-Whitney) test since data are nonnormal (thus t-test is not appropriate), ordinal and thus not continuous, multimodal, with different variances and with different shapes (thus Mann-Whitney test is not appropriate). How can I test my hypothesis? Thank you for any answer.


title: How to interpret Dickey Fuller (DF) test results in R (for unit test)

text: I read 1) Intuitive explanation of unit root and 2) http://www.r-bloggers.com/unit-root-tests/ for doing unit root test. I have basic questions: 1)should I check for unit root on both 'x' and 'y' variables of my lm equation or just one of them? For AR process, it doesnt matter (I suppose), but my "lm" equation has 1) few variables(x's) on right hand side, and also 2)one lag of y. Should I check ur.lm for all x's as well? 2)I am confused on what parameters should I pass to ur.df? what does type parameter mean(none/trend/drift) and how do I decide which of none/trend/drift is the right one 3)How would I interpret the following df test result(what is z.diff.lag? Is it correct to say z.diff.lag does NOT exhibit unit root given its high statistic, hence z.diff.lag is stationary; but what is z.diff.lag anyway)                   > summary(ur.df(temp3))              ###############################################           Augmented Dickey-Fuller Test Unit Root Test #           ###############################################               Test regression none          Call:               lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)         Residuals:               Min      1Q  Median      3Q     Max                 -5.7048 -0.7167 -0.1592  0.4686  7.4964           Coefficients:               Estimate Std. Error t value Pr(>|t|)                    z.lag.1    -0.0002636  0.0001407  -1.874   0.0609 .                  z.diff.lag -0.3607210  0.0048720 -74.040   <2e-16 ***


title: Strategy for finding cluster solution with positive values for Cubic Clustering Criterion

text: I am trying to cluster a data set with approximately 3500 entries, consisting primarily of likert-10 scale survey answers. Most of the survey answers are skewed - a typical example is below: ![distribution](http://i.imgur.com/pDo18EM.jpg) I've done the following:   * Range-standardized the data   * Removed rows where the survey respondent answered all questions with a 5/neutral or 10/strongly agree   * Used SAS's PROC VARCLUS and PROC FACTOR to reduce the variables used in the clustering solution   * Eliminated variables with the greatest skew   * Use the TRIM= option to eliminate cluster outliers. Despite all of this, all of my cluster solutions are consistently showing negative Cubic Clustering Criterion values. Here is a typical example: ![Negative CCC](http://i.imgur.com/dzrBv09.jpg) Can anyone suggest a strategy for coming up with a cluster solution with positive CCC? I'm not sure what to try next.


title: Frequency based discounting of values

text: Let's say there are two products `A` and `B` both of which have some value. Consider a population of `10.000` people choosing between these two products. I would like to come up with a function that adjusts the price of `A` and `B` based on the frequency of people choosing the products. The function should satisfy the following criteria: `1.` The value of A and B should always sum to 1 and be in the range between 0 and 1. `2.` The more people are choosing a product the lower should its value be and similarly the fewer people are choosing a product the higher should its value be. To illustrate let's say that `on day 1:` `6000` people are choosing product `A` and `4000` are choosing product `B` and their value is `0.6` and `0.4` respectively. `on day 2:` `7000` people are choosing product `A` and `3000` are choosing product `B` the value of `A` should decrease and the value of `B` should increase. `3.` The adjustment of the value should not be linearly increasing or decreasing based on the frequency of people using the product but the increase or decrease should be steeper when the proportion of people using the products is close to `50%` but should increase or decrease slower when the people choosing the products is skewed (e.g. `90` vs. `10%`).


title: How to test the ARIMA coefficients?

text: Which test is required to test whether coefficients estimated as part of ARIMA procedure is different from 0? And how does one compute this test? I am reading some procedures regarding the inversion of Hessian matrix during MLE. Is this some direction I can pursue? I really do not want to get into bootstrapping or some other resampling as my dataset is large and takes forever to compute.


title: Regression Analysis

text: hello can someone help me with this problem? I am using R, so what functions I should use? Thank you![enter image description here](http://i.stack.imgur.com/IurEa.png)


title: Sample autocovariance of Durbin–Watson test

text: I understand Durbin–Watson test, but I can't understand this sentence. I cannot prove it. > The Durbin-Watson test statistics is asymptotically equivalent to (rootT*C), > where C is the sample first-order autocovariance of the least squares > residuals. In my opinion, we usually identify the sample autocorrelation, not autocovariance. So, I cannot take out the sample autocovariance from Durbin–Watson equation. Can anyone explain this sentence?


title: Dropping data from people who have "perfect" scores

text: OK, so I have data from a class that had a preparatory self-test to see how prepared they were for the class, and the final results for the class. The preparatory self-test had a range from 0..13 and the final had a range from 0..100. (don't mind the normal, none, distinction legend, that's just the class (Coursera) had a threshold for certificates at 50%, and you'd get a "with distinction" above 80%) But, here, look at these two graphs: ![Before culling the 13s](http://cablemodem.hex21.com/~binesh/coursera/useful- genetics-1/scatterplot-full.png) ![After culling the 13s](http://cablemodem.hex21.com/~binesh/coursera/useful- genetics-1/scatterplot.png) The correlation coefficient for the full data set is r=0.29782 and if I remove the 13s from the preparatory self-test becomes r=0.294383, so I guess the correlation doesn't change all that much. But, I guess my question was more sort of, is it valid/a good idea to remove the 13's from the preparatory test? Here's what I was thinking: The final seems to be valid for the set of students being tested. But, you really can't distinguish between one student who got a 13 and another who got a 13, because they "maxed" out the test. (Like if you had a scale that measured up to 100 grams, and you put an elephant, a whale and a chicken on it, and declared them all to be "100 grams") What should one do if the test is being maxed out by certain people? Is it better to leave it in, to remove it, or to do something else altogether with it? Thanks!


title: Classification of 10 or more ordinal classes

text: We have to classify new data in 11 ordinal classes based on a training set. Do I use an ordinal approach such as ordered logit; or methods for continuous data? But how do I cut the new data then into classes? Do I just round them to the nearest class? For example, we have to guess years of education. I could do simple OLS and round to the nearest value. But what if my model predicts -3? or +20? My intuition is that in this case it is also better to work with a function in the range of 0 to 1 and estimate probabilities rather than absolute values. Is that right?


title: Joint distribution calculation of 2 variables in R

text: I'm trying my hands on doing Bayesian analysis. So need some guidance. I have a dataset in the following form: ![enter image description here](http://i.stack.imgur.com/ykQt0.jpg) Will each row have a joint distribution value. How do we calculate it in R, what packages should I use?


title: Relationship between entropy and information gain

text: Based on papers :1. Deniz Erdogmus, Member, IEEE, and Jose C. Principe, An Error-Entropy Minimization Algorithm for Supervised Training of Nonlinear Adaptive Systems   1. J. Principe, D. Xu, and J. Fisher, “Information theoretic learning,” in Unsupervised Adaptive Filtering, S. Haykin, Ed. New York: Wiley, 2000, vol. I, pp. 265–319. Entropy (Shannon and Renyis) has been used in learning by minimizing the entropy of the error as an objective function instead of the Mean Square error. Rationale is that minimizing entropy = maximizing mutual information. Now, entropy = disorder = uncertainty. Higher the uncertainty, more is the entropy. Also, higher entropy = high information content (used in compression), hence we cannot compress a signal with high entropy. So, in view of the above I have questions which is Is mutual information another name for information gain? Next, if max entropy implies high information content then why do we minimize entropy of the error between output signal and the desired signal. Is there any proof which shows that minimizing entropy of error when used as a fitness function would mean that we are getting close to the true estimate of the unknown parameter?


title: Replication in multivariate logit in R

text: I am new to both logit models and GLMs, but I think one of those two model classes might be the correct analysis for my data set: I am interested in comparing the composition of the diet of fish collected in 3 habitat types - thus "habitat" is a categorical predictor with 3 levels. I want to test 2 types of null hypotheses:   1. That the relative proportions of food types A, B and C are the same among the 3 habitat types    2. Within each habitat, that the biomass of food is equally distributed between prey types A, B, and C. In other words, to be able to make evaluations such as "The proportional composition of fish diet varied significantly among habitats.... The relative contribution of food A to fish diets was significantly higher in Habitat 1 than in Habitat 3... Food type A comprised a larger proportion if the diet of fish in Habitat 1 than did food type B." The response data are proportions of food in categories A, B or C. Of course the proportions are bounded between 0 and 1 and within an individual fish the proportions sum to 1. The raw data are actually proportions, not counts of food items. I've taken the raw food item counts and used length/width relationships for the food types to yield biomass estimates for each food type - so the data in the "cells" for columns A, B and C are proportions of the total biomass in the gut of an individual fish. I have analyzed multiple fish from each habitat type (22 from Habitat 1, 20 from Habitat 2, 8 from Habitat 3). In reading on log-linear functions and multinomial logistic regressions, I now have myself thoroughly confused. The elements of this application:   1. I have a definite predictor variable (habitat)    2. At the level of the individual fish (the replicates), the 3 food categories (A,B, and C) are not mutually exclusive.    3. Because of the biomass conversion, the raw data are proportions, not counts. What it the best way to evaluate the null hypotheses for these data? What post-hoc comparison approaches are appropriate? General suggestions are welcome, suggestions tailored to R are even more welcome.


title: Use of text mining document classification via decision rules?

text: Edit: could someone please explain decision rules in text mining and what they are useful for and how to obtain them from text?


title: Canonical Correlation

text: I am working on research involving the collection of data from employees and customers. Specifically, my research question deals with the influence of employee satisfaction on the quality of customer service. I will collect data from both employees and customers, and plan on running a correlation analysis on the results. I know i will run **canonical correlation** , but my question is, if i have collected 200 questionnaires from employees, and I have 400 from customers, in a single SPSS file, there is an issue of missing data, since employees data is 200 less than the customers data. Does anyone have a suggestion for how I should handle this?


title: Clustering large scale data in fine-grained clusters

text: I've got large amount of data (e.g. 100K) and I want to cluster them in very fine-grained clusters (e.g. 10K). I look for an appropriate algorithm that uses the similarity function instead of whole data and also I can't compute whole similarity matrix because of data scale. Python implementation is appreciated.


title: Chi Square Test for QC Samples

text: An analytical procedure is run for 112 consecutive days. Each day, the analyses are run in four batches in order -- Batch-1, Batch-2, etc. In each batch is a QC sample, labeled QC1, QC2, etc. The following R code generates a contingency table for the data, showing the number of times each QC sample was a particular rank. For example, QC1 had the highest result 43 days out of 112 and QC4 had the lowest result 28 days out of 112. The analyst has the feeling that the first batch of the day tends to have the highest results. What statistical test should I use to check out whether or not the ranks of the QC samples are random from day to day? I understand that I should not use a chi- square test of independence between the QC sample and the ranks, but why not?               count <- c(43,26,26,17,24,33,28,27,18,27,27,40,27,26,31,28)     mat.obs <- matrix(count,nrow=4,byrow=TRUE)      rownames(mat.obs) <- c("QC1","QC2","QC3","QC4")     colnames(mat.obs) <- c("1","2","3","4")     mat.obs


title: Nesting, sort-of: Multiple observations mapping to a single DV measure?

text: This is a made up example, but it gets to a point I am trying to figure out. I know that you can use a nested model if your IVs are collected in clusters -- for example, if you want to model the effect of "the # of hours a student spends working" on his "grades", you can cluster by school and maybe teacher. But, what if you want to test the effect of "the # of hours students spend working" on "the ranking of the school"? Then, you have all observations within a school mapping to a single DV measure. This seems pretty different. How should it be modeled?


title: Nominal/ordinal Variable Measures of Association

text: I have an ordinal variable (employee performance ratings on 5-point scale) and I want to do a crosstab with a nominal variable (4 levels of race/ethnicity). What is the best measure of association for these two variables in SPSS?


title: Cox regression to see if distribution of disease predicts survival

text: I'm working with clinical data from a few hundred patients. Every patient was measured at 6 different timepoints with regular intervals. Each patient is assessed on 3 different domains and given a score for each domain ranging from 0 - 11. On this scale 11 is perfect and 0 is the worst score. A very focal distribution of scores is 11,11,0 and an example of a very diffuse distribution is 7,7,7. The research question is whether the distribution of the scores is related to survival length of the patient. My current approach is the following: I subtract the lowest score from the two other scores and square each. For example: a score of 0,2,11 yields a value of (11-0)2 + (2-0)2 = 125. Subsequently, I perform a cox regression to look at the correlation with survival. Furthermore, I want to match the pattern of progression of the patient to a very focal progression of disease, i.e. someone who progresses to 0,11,11 rather than 8,7,7 for a sum of 22. This would allow to incorporate the data of multiple timepoints in the analysis. Do you think the cox regression is a fitting approach? The problem is that there is a different outcome possible for each sum of the three domains. Secondly, what analysis would I have to do to relate the degree to which the pattern of progress matches with survival?


title: Statistical test for positive and negative predictive value

text: I was reading through a paper and I saw a table with a comparison between PPV (Positive Predictive Value) and NPV (Negative Predictive Value). They did some kind of statistical test for them, this is a sketch of the table:               PPV    NPV    p-value     65.9   100    < 0.00001     ...      Every rows refers to a particular contingency table. What kind of hypothesis test did they do? Thanks!


title: Article about misuse of statistical method in NYTimes

text: I am referring to this article: http://www.nytimes.com/2011/01/11/science/11esp.html > Consider the following experiment. Suppose there was reason to believe that > a coin was slightly weighted toward heads. In a test, the coin comes up > heads 527 times out of 1,000. > > Is this significant evidence that the coin is weighted? > > Classical analysis says yes. With a fair coin, the chances of getting 527 or > more heads in 1,000 flips is less than 1 in 20, or 5 percent, the > conventional cutoff. To put it another way: the experiment finds evidence of > a weighted coin “with 95 percent confidence.” > > Yet many statisticians do not buy it. One in 20 is the probability of > getting any number of heads above 526 in 1,000 throws. That is, it is the > sum of the probability of flipping 527, the probability of flipping 528, 529 > and so on. > > But the experiment did not find all of the numbers in that range; it found > just one — 527. It is thus more accurate, these experts say, to calculate > the probability of getting that one number — 527 — if the coin is weighted, > and compare it with the probability of getting the same number if the coin > is fair. > > Statisticians can show that this ratio cannot be higher than about 4 to 1, > according to Paul Speckman, a statistician, who, with Jeff Rouder, a > psychologist, provided the example. **First question:** This is new to me. Has anybody a reference where I can find the exact calculation and/or can YOU help me by giving me the exact calculation yourself and/or can you point me to some material where I can find similar examples? > Bayes devised a way to update the probability for a hypothesis as new > evidence comes in. > > So in evaluating the strength of a given finding, Bayesian (pronounced BAYZ- > ee-un) analysis incorporates known probabilities, if available, from outside > the study. > > It might be called the “Yeah, right” effect. If a study finds that kumquats > reduce the risk of heart disease by 90 percent, that a treatment cures > alcohol addiction in a week, that sensitive parents are twice as likely to > give birth to a girl as to a boy, the Bayesian response matches that of the > native skeptic: Yeah, right. The study findings are weighed against what is > observable out in the world. > > In at least one area of medicine — diagnostic screening tests — researchers > already use known probabilities to evaluate new findings. For instance, a > new lie-detection test may be 90 percent accurate, correctly flagging 9 out > of 10 liars. But if it is given to a population of 100 people already known > to include 10 liars, the test is a lot less impressive. > > It correctly identifies 9 of the 10 liars and misses one; but it incorrectly > identifies 9 of the other 90 as lying. Dividing the so-called true positives > (9) by the total number of people the test flagged (18) gives an accuracy > rate of 50 percent. The “false positives” and “false negatives” depend on > the known rates in the population. **Second question:** How do you exactly judge if a new finding is "real" or not with this method? And: Isn't this as arbitrary as the 5%-barrier because of the use of some preset prior probability?


title: What is the sum of squared t variates?

text: Let $t_i$ be drawn i.i.d from a Student t distribution with $n$ degrees of freedom, for moderately sized $n$ (say less than 100). Define $$T = \sum_{1\le i \le k} t_i^2$$ Is $T$ distributed nearly as a chi-square with $k$ degrees of freedom? Is there something like the Central Limit Theorem for the sum of squared random variables?


title: Are there statistical lessons from the "Bible Code" episode

text: Although this question is somewhat subjective, I hope it qualifies as a good subjective question according to the faq guidelines. It is based on a question that Olle Häggström asked me a year ago and although I have some thoughts about it I do not have a definite answer and I would appreciate some help from others. ## Background: A paper entitled "Equidistant letter sequences in the book of Genesis," by D. Witztum, E. Rips and Y. Rosenberg made the extraordinary claim that the Hebrew text of the Book of Genesis encodes events which did not occur until millennia after the text was written. The paper was published by "Statistical Science" in 1994 (Vol. 9 429-438), and was offered as a "challenging puzzle" whose solution may contribute to the field of statistics. In reply, another paper entitled "Solving the Bible code puzzle" by B. McKay, D. Bar-Natan, M. Bar-Hillel and G. Kalai appeared in Statistical science in 1999 (Vol. 14 (1999) 150-173). The new paper argues that Witztum, Rips and Rosenberg's case is fatally defective, indeed that their result merely reflects on the choices made in designing their experiment and collecting the data for it. The paper presents extensive evidence in support of that conclusion. (My own interests which are summarized in Section 8 of our paper are detailed in another technical report with Bar Hillel and Mckay entitled "The two famous rabbis experiments: how similar is too similar?" See also this site.) ## The questions: Olle Häggström's specific question was: > "I once suggested that your paper might be useful in a statistics course on > advanced undergraduate level, for the purpose of illustrating the pitfalls > of data mining and related techniques. Would you agree?" In addition to Olle's question let me ask a more general question. > Is there something related to statistics that we have learned, (including > perhaps some interesting questions to ask) from the Bible Code episode. Just to make it clear, my question is restricted to insights **related to statistics** and not to any other aspect of this episode.


title: Use Empirical CDF vs Distribution CDF?

text: _NOTE: I purposely did not label the axis due to pending publications. The line colors represent the same data in all three plots._   I fitted my data using a negative binomial distribution to generate a pdf. I am happy with the pdf and meets my research needs. PDF plot: ![alt text](http://i.stack.imgur.com/ttDnS.png) * * * **For when reporting the CDF, should I use the empirical or fitted CDF?** There are slight differences between the empirical and fitted CDF, specifically at x = 40, the yellow and cyan lines intersect in the empirical distribution, but not the fitted. Empirical:   ![alt text](http://i.stack.imgur.com/j0ELI.png) Negative Binomial CDF: ![alt text](http://i.stack.imgur.com/13yyc.png)


title: What is the difference between "multinomial" and "polytomous"?

text: I would like to know the difference a multinomial model and a polytomous model. Both a multinomial model, for example multinomial logit / probit model, and a polytomous model are used for estimation of quantitative data. But, a polytomous model is seldom used. I think that the difference is whether to estimate parameters for every category. So, in polytomous model, the parameters of each category is the same. Is it correct? Moreover, I would be interested in understanding the estimation method for the polytomous model. If someone understands it, I would appreciate it someone could point that out to me.


title: Is binomial dispersion the same as binomial distribution?

text: I have been asked by a client to write a program in R which calculates the "binomial dispersion test." I have not found anything for this, but I have found that you can calculate the binomial distribution. Are they the same thing? I suspect they are.


title: Power analysis for moderator effect in regression with two continuous predictors

text: Related to an earlier question on power analysis for multiple regression, a social science researcher asked me about power analysis for moderator regression (i.e., an interaction effect). The researcher asked me: > I seem to recall that power of tests for moderation with two continuous > predictor variables is low - do you know the minimum sample size requirement > in this context? From the context, it can further be assumed that this is an observational study (not an experimental study) and that the dependent variable is continuous. ### Question   * What advice would you give regarding calculating the minimum sample size required?   * Are there any caveats that you would present?


title: How do I calculate a confidence interval for the mean of a log-normal data set?

text: I've heard/seen in several places that you can transform the data set into something that is normal-distributed by taking the logarithm of each sample, calculate the confidence interval for the transformed data, and transform the confidence interval back using the inverse operation (e.g. raise 10 to the power of the lower and upper bounds, respectively, for $\log_{10}$). However, I'm a bit suspicious of this method, simply because it doesn't work for the mean itself: $10^{\operatorname{mean}(\log_{10}(X))} \ne \operatorname{mean}(X)$ What is the correct way to do this? If it doesn't work for the mean itself, how can it possibly work for the confidence interval for the mean?


title: How do i compare NDCG values?

text: I want to compare different ranking models using the NDCG metric. I average the NDCG of the models, then use the difference between the avg. NDCGs to compare the models. My questions:   1. What does it mean for a ranker to have 0.1 (avg.) NDCG more than another ranker? Can i say more than: "It ranks a little better"? Especially, how big a difference in NDCG would be "significant" in absolute terms?   2. Can i say something about statistical significance of the difference of the average NDCG w.r.t the 0-Hypothesis? I'm quite unsure about whether normality assumption holds here... Especially since i'm not looking at the average of the differences, but the difference of the averages.


title: Data mining approaches for analysis of sequential data with nominal attributes

text: Question for the experienced data miners out there: ### Given this scenario:   * There are N shopping carts   * Each shopping cart is filled with an arbitrary number of M items from an infinitely large set (with the current amount of data I have, that arbitrary number can hit numbers around 1500)   * The order in which each cart is filled is significant   * There are other attributes such as geolocation of shopper, but these can be (and currently are) thrown out in favor of making the algorithm simpler ### I need to:   * At a particular point in time, given only the ordered sets of items in each cart, identify 'similar' carts without prior knowledge of class labels   * After a certain amount of data has been collected and a drudge works through the data and assigns labels, create a classifier that can work quickly with future unseen data ### Initial approach:   * So far, my approach has been focused on the first point. My method uses k-means clustering and handles the sequential nature of the data by using a distance matrix generated by calculating the Hamming distance between carts. In this way, [apple, banana, pear] is different from [pear, apple, banana], but [apple, banana, pear] is less different from [apple, banana, antelope]. The appropriate value of k is determined through investigation of the silhouette coefficient. The clusters generated from this seem to make sense, but the runtime of my method will definitely be prohibitive as my dataset scales. ### Question:   * Would anyone happen to have any suggestions for a novice data miner for this problem? ### Edits with more info:   * I've found suggestions that consider using n-gram features and comparing them pair-wise. A concern I have about this is order: will the order of the sequences be maintained if n-gram models are used? Also, I see performance issues being a larger possibility with this method.


title: normalized cointegrating equations in johansen cointegration test

text: I have 3 cointegrating equations. When looking at the long term relationship among the variables, I guess I have to look at "Normalized cointegrating coefficients (standard error in parentheses)" for "3 Cointegrating Equation(s)". How can I then formulate the formula? because most variables are normalized to 1.000 Thank you!


title: Deep Belief Network (number of layers)

text: So we have "several RBMs" > Deep Belief Network >   >   >     A deep belief network is obtained by stacking several RBMs on top of > each >     other. The hidden layer of the RBM at layer `i` becomes the input of the >     RBM at layer `i+1`. The first layer RBM gets as input the input of the >     network, and the hidden layer of the last RBM represents the output. > When >     used for classification, the DBN is treated as a MLP, by adding a > logistic >     regression layer on top. >   but why more RBM's don't means better? How to determine best number of layers? Seems related: Deep belief network performs worse than a simple MLP


title: Permutations exercises with alike objects; check my homework

text: could you please review my answers and confirm my reasoning. This is self- teaching work. In how many ways can 8 people be seated in a row if: (a) there are no restrictions: ### Answer: 8! (b) persons A and B must sit next to each other: ### Answer: 8!/2! Now, the demoninator is 2! because the couple are considered 'a like', right? (c) there are 4 men and 4 women and no 2 men and 2 women can sit next to each other ### Answer: 8!/(4!*4!) I think this is correct, can you explain why? (d) there are 5 men and they must sit next to each other. ### Answer: 8!/5! I can't explain why this is correct (e) there are 4 married couples and each couple must sit together ### 8!/(2!X4) the 4 couples are alike.


title: What is the interpretation of "generalized" partial correlations?

text: The usual partial correlation between X and Y given the set of variables Z is the Pearson correlation between residuals resulting from the liner regression of X on Z and Y on Z. It can be computed using the recursive formula (see wikipedia) from the correlation matrix. My question is what is the interpretation of the partial correlation if the correlation matrix contains other correlation coefficients, like Kendall's tau or Spearman's rho, and not Pearson correlation coefficients.


title: Independent variable misclassification and statistical tests

text: Suppose we have our standard DGP, $y=\alpha+\beta x+\varepsilon,$ where $x$ is binary. Let's say the observed $x$ is actually measured with error, so that the explanatory variable is misclassified for some fraction of the population. This is very non-classical measurement error. I am trying to summarize the effects of this type of measurement error. From Bound, Brown, and Mathiowetz (2001), I've learned that as long as the misclassification is non-differential (and none of the other classical measurement error assumptions are violated), the coefficient on $x$ will be attenuated. Here non-differential means the classifications errors are independent of the dependent variable $y$. Two questions remain:   1. My simulation-based intuition is that size goes up and the power of the test goes down as the degree of misclassification rises. I am somewhat more confident of the second.   2. If $x$ is vector (but we only have measurement error in one variable), is the bias still unambiguously toward zero?  Bound, John, Charles Brown, and Nancy Mathiowetz, 2001. “Measurement Error in Survey Data.” In James J. Heckman and Edward E. Leamer (Eds.), _Handbook of Econometrics_ , Volume 5. Amsterdam: Elsevier Science, 2001. Link to working paper version.


title: Merge text collection subsamples for cross-validation

text: I'm doing a research and to test our methods we're using 5-fold cross- validation. Our collections consist of some pre-classified text. To use the data in Weka, I'm pre-processing it with the `StringToWordVector` filter. So far, so good... The problem is: both me and my supervisor are running the algorithms. I use Weka, and he uses his algorithms implemented in Matlab. In order to use the same subsamples for the 5-fold cross-validation in both Weka and Matlab, I split each of our collections into 5 subsamples. After doing that, for each subsample I created the `arff` file and applied the `StringToWordVector`. **The question is** : how can I merge these subsamples in order to get the train and test files for each fold? I can't merge the text files manually and then apply the `StringToWordVector`, because the words on the training file would be different than the words on the test file. I could also apply the `StringToWordVector` filter in the collection, without splitting it into subsamples, and then use the `StratifiedRemoveFolds` filter to generate the train and test files, but I want to use the same subsets on both Weka and Matlab and, if possible, I want to use the subsamples that I've already created.


title: Extending 2-class models to multi-class problems

text: This paper on Adaboost gives some suggestions and code (page 17) for extending 2-class models to K-class problems. I would like to generalize this code, such that I can easily plug in different 2-class models and compare the results. Because most classification models have a formula interface and a `predict` method, some of this should be relatively easy. Unfortunately, I haven't found a standard way of extracting class probabilities from 2-class models, so each model will require some custom code. Here's a function I wrote to break up a K-class problem into 2-class problems, and return K models:               oneVsAll <- function(X,Y,FUN,...) {         models <- lapply(unique(Y), function(x) {             name <- as.character(x)             .Target <- factor(ifelse(Y==name,name,'other'), levels=c(name, 'other'))             dat <- data.frame(.Target, X)             model <- FUN(.Target~., data=dat, ...)             return(model)         })         names(models) <- unique(Y)         info <- list(X=X, Y=Y, classes=unique(Y))         out <- list(models=models, info=info)         class(out) <- 'oneVsAll'         return(out)     }      Here's a prediction method I wrote to iterate over each model and make predictions:               predict.oneVsAll <- function(object, newX=object$info$X, ...) {         stopifnot(class(object)=='oneVsAll')         lapply(object$models, function(x) {             predict(x, newX, ...)         })     }      And finally, here's a function to turn normalize a `data.frame` of predicted probabilities and classify the cases. Note that it is up to you to construct the K-column `data.frame` of probabilities from each model, as there is not a unified way to extract class probabilities from a 2-class model:               classify <- function(dat) {         out <- dat/rowSums(dat)         out$Class <- apply(dat, 1, function(x) names(dat)[which.max(x)])         out     }      Here's an example using `adaboost`:               library(ada)     library(caret)      X <- iris[,-5]     Y <- iris[,5]     myModels <- oneVsAll(X, Y, ada)     preds <- predict(myModels, X, type='probs')     preds <- data.frame(lapply(preds, function(x) x[,2])) #Make a data.frame of probs     preds <- classify(preds)     >confusionMatrix(preds$Class, Y)     Confusion Matrix and Statistics                      Reference     Prediction   setosa versicolor virginica       setosa         50          0         0       versicolor      0         47         2       virginica       0          3        48      Here is an example using `lda` (I know lda can handle multiple classes, but this is just an example):               library(MASS)     myModels <- oneVsAll(X, Y, lda)     preds <- predict(myModels, X)     preds <- data.frame(lapply(preds, function(x) x[[2]][,1])) #Make a data.frame of probs     preds <- classify(preds)     >confusionMatrix(preds$Class, Y)     Confusion Matrix and Statistics                      Reference     Prediction   setosa versicolor virginica       setosa         50          0         0       versicolor      0         39         5       virginica       0         11        45      These functions should work for any 2-class model with a formula interface and a `predict` method. Note that you have to manually split up the X and Y components, which is a little ugly, but writing a formula interface is beyond me at the moment. Does this approach make sense to everyone? Is there any way I can improve it, or is there an existing package to solve this issue?


title: Name of mean absolute error analogue to Brier score?

text: Yesterday's question Determine accuracy of model which estimates probability of event got me curious about probability scoring. The Brier score $$\frac{1}{N}\sum\limits _{i=1}^{N}(prediction_i - reference_i)^2$$ is a mean squared error measure. Does the analogous mean absolute error performance measure   $$\frac{1}{N}\sum\limits _{i=1}^{N}|prediction_i - reference_i|$$ have a name, too?


title: Probability that a draw from one sample is greater than a draw from a another sample

text: I have a sample of `b` for id=1 and I want to compare this sample with the sample of id={2,3} of `a` i.e. I want to compute the probability that a draw from sample `b` (id=1) is greater than a draw from sample `a` (id={2,3}).                   library(data.table)         df <- data.table(id=c(rep(1,500),rep(2,500),rep(3,500)),a=rnorm(1500),b=rnorm(1500,1))         ot <- df[id!=1,a]         ow <- df[id==1,b]         sa <- df[id!=1,sample(a,500)]         #all permutations of eig and sam         fall <- expand.grid(ow,sa)         #all permutations of eig and other         all  <- expand.grid(ow,ot)              #1. best possible empirical measure         sum(all[,1]>all[,2])/nrow(all)         #2.          sum(fall[,1]>fall[,2])/nrow(fall)         #3.           sum(ow>sa)/length(sa)      Nb. 1 would obviously be the best solution as in my real data set the distributional properties of a and b are unknown. Nb. 1 is however not feasible because I would have to repeat this procedure for every `ID`. My data set is huge. The solution I used until now was Nb. 3: I a draw a random sample from `a` (id={2,3}) in the size of the sample for `b` (id=1). But this is of course imprecise. Solution Nb. 2 would be a compromise but is also not feasible computation-wise. I would like to find a data table (this is a faster implementation of data frames in the form of the R package "data.table") solution because those are normally much faster. I didn't post this question on Stackoverflow, because I also hoped that one of you could maybe hint me to a mathematical trick to overcome this problem, but of course also coding advises are welcome. Thanks in advance!


title: Modeling paired data: Philosophical difference between using mixed effects model and using grouping-variable-level outcome

text: I am facing a problem of detecting change in infection rate after surgery in about 400 eye surgery patients. Crucial medical background includes the following-- About 30% patients have both eyes undergoing surgery, 70% have one eye undergoing surgery; Infection status (yes/no) is collected on day 14 after surgery. All patients are treated with antibiotics on day 1 after surgery to prevent infection. There is no difference in dosage between those undergoing surgery on one or two eyes. The antibiotics' effect is systemic (whole-body) not local (only effective at the eye having operation). What is the philosophical / interpretive difference between-- Using infection rate per eye as outcome and mixed-effects logistic regression model with patient ID as random-effects variable; and Collapsing the paired data into per-person infection status, create a categorical variable "number of eyes undergoing surgery", and use infection rate per patient as outcome and "number of eyes undergoing surgery" as one of the covariates in a logistic regression? An obvious observation is that if collapse the paired data, the infection per- patient prevalence will look higher than the per-eye prevalence. Which model makes more sense (and why) please? Thank you!


title: Using the calibrate function in rms package for Cox model: use 'hare' method or 'KM'?

text: I'm creating a calibration plot for my breast cancer prognostic Cox model, which doesn't include any fancy transformations, using the `calibrate()` function in the `rms` package for `R`. Should I use the `hare` method or `KM`? Similar literature to what I'm doing always uses the `KM` method, but I'm curious if `hare` is more advantageous. There are 3 options I have: use the `KM` method, use the `hare` method, or plot both `hare` and `KM` on one figure. Here's sample code:               library(rms)     library(survival)     library(mice)          remove(veteran)     data(veteran)     veteran$trt=factor(veteran$trt,levels=c(1,2))     veteran$prior=factor(veteran$prior,levels=c(0,10))     survmod=with(veteran,Surv(time,status))     fit=cph(survmod~celltype+karno,data=veteran,surv=TRUE,u=100,x=T,y=T)     cal.KM <- calibrate(fit, u=100, cmethod='KM', m=20, B=20)     cal.hare=calibrate(fit,u=100,cmethod='hare',B=20)     plot(cal.hare)     plot(cal.KM,add=TRUE)      Thank you for your help!


title: Goodman Gamma for nominal variables

text: I need to give a direction to the association between an ordinal (independent) and a nominal variable (dependent). This, after having performed a Chi-square test. The ordinal variable is like a blood test score with values (high, low, average) and the nominal variable is yes/no (reflecting the possibility of being sick). Is it OK to use the Goodman gamma to do this? If not because one variable is nominal, which is a good alternative?


title: If the t-test and the ANOVA for two groups are equivalent, why aren't their assumptions equivalent?

text: I'm sure I've got this completely wrapped round my head, but I just can't figure it out. The t-test compares two normal distributions using the Z distribution. That's why there's an assumption of normality in the DATA. ANOVA is equivalent to linear regression with dummy variables, and uses sums of squares, just like OLS. That's why there's an assumption of normality of RESIDUALS. It's taken me several years, but I think I've finally grasped those basic facts. So why is it that the t-test is equivalent to ANOVA with two groups? How can they be equivalent if they don't even assume the same things about the data?


title: what can be inferred from the error plot of my classifier?

text: I am working on a classification task where I use 12 features. My training set has 400 samples of positive data and 2000+ samples of unlabelled data. While testing has 34 positive samples and 999 samples of unlabelled data which is mostly negative ( but not known for sure how many ). Before using it for training I didn't standardised them for two reasons... I read somewhere that the variance among the samples has information ( here it is a biological data ) and standardisation kills information. Secondly when I tried to train and test classifier after standardisation it remained poorer in performance. I am using semi supervised learning technique and tool i use for the same is http://vikas.sindhwani.org/svmlin.html ( parameters chosen A=3 and W=0.3) I prepared an input file with positive and unlabelled datasets initially. After training I calculated the error and plotted the same for training sets and test sets. Here is the plot i get https://www.dropbox.com/s/uqma194uiipm9n8/sample_plot.png I have three issues: 1) What does this error plot infer ? Red is testing error and blue is the training error. x-axis is data size, y is error . Why I am having such high fluctuations in error while a standard classifier error plot has less ups and downs and both come closer with training 2) It is known that my data has non-linear characteristics while the classifier I am using has linear kernel. Will this effect the classification results drastically ? When I say drastically , it means that my classifier performs fair on my independent test sets, but fails miserably on the same kind of but different independent dataset from other sources. 3) How to chose threshold in semi-supervised learning technique and in general for any classifier ? In co-training method we need to separate out negatives from unlabeled samples. But I couldn't find how one can separate out and on what basis see Table 1 from this paper \--->> www.dropbox.com/s/7rhis2gl5ztm3wn/BlumMitchell98.pdf 4) Is it necessary to standardise the data in a classifier, knowing that biological data has non-linear form and standardisation may cause to loose information ?


title: Test for equal variance

text: I need to test for equal variance. I tried to read up and it mentioned there are a few tests I can use to test for equality of variance such as F Test (F ratio?), Levene Test, Bartlett Test etc. I was wondering if there is a qualitative method of telling whether there is equal variance by using the SD and the mean of the two data sets (both normally distributed) My lecturer was saying I can assume equal variance by just looking the mean and the SD? I do not know how to tell... can anyone help?


title: Which vocabulary to use for a special version of the Rocchio classification algorithm?

text: My issue is a bit hard to explain in this question's title, so hopefully I can make clear what my problem is about in this text. I'm dealing with partially supervised text classification. I have a set of positive documents and a set of unlabeled documents (which contains both positive and negative documents). My goal is to identify documents in this unlabeled set which are most probably negative documents. As soon as I've identified them, I use those and the positive set to classify the rest of the unlabeled documents. In order to identify the set of reliable negative documents, I use a special version of the Rocchio classification algorithm which is explained in this paper: > Xiao-Li Li, Bing Liu, See-Kiong Ng (2010) _Negative Training Data can be > Harmful to Text Classification_. In: Proceedings of the Conference on > Empirical Methods in Natural Language Processing (EMNLP-2010). In the upper left corner of page 6 (figure 4), there is pseudocode explaining the algorithm. Below, are the two relevant lines that my question is about. $$ \begin{aligned} 3.~~~~~p &= \alpha\frac{1}{|P|}\sum_{d_j\in P}\frac{d_j}{||d_j||}-\beta\frac{1}{|PN|}\sum_{d_i\in PN}\frac{d_i}{||d_i||} ; \\\ 4.~~~~~n &= \alpha\frac{1}{|PN|}\sum_{d_i\in PN}\frac{d_i}{||d_i||}-\beta\frac{1}{|P|}\sum_{d_j\in P}\frac{d_j}{||d_j||} ; \end{aligned} $$ At this stage of the algorithm, I have a set _P_ of positive documents and a set _PN_ of potential negative documents which were identified in a previous step. Each document in these two sets is represented as a vector (in bold letters) of TF-IDF values of the word vocabulary from the respective set. In the first line of the code above, I subtract the _PN_ -vector from the _P_ -vector. In the second line, I do it vice versa. The goal is to create a positive prototype vector **p** and and a negative prototype vector **n**. **My question is the following:** Which vocabulary do I have to take into account for each of these two subtractions? Do I have to create all feature vectors from the entire vocabulary of both the positive and the potential negative set? Or do I have to use only the vocabulary from the positive set in the first line and the vocabulary from the potential negative set in the second line? Or anything completely different from that? Unfortunately, this isn't explained anywhere.


title: Regressing on only the positive part of a vector

text: I have an interesting question, with its original application in finance. Suppose I have a stock return $Y$, and a set of independent variables (other tradable assets) $X$. Typically, one _hedges_ Y with X by finding the least squares solution to regressing $Y$ on $X$. Thus $Y - X\beta$ is approximately zero. But what if we instead decompose $Y = Y^{+} + Y^{-}$, where $Y^{+}$ is the vector of _positive_ components of $Y$ and zero otherwise, likewise for $Y^{-}$. If we then regress $Y^{-}$ onto $X$, then $Y - X\beta_{neg}$ only hedges the negative values, and exposes us to gains from positive values of $Y$ as $$Y - X\beta_{neg}$$ $$ = Y^{+} + ( Y^{-} - X\beta_{neg} )$$ $$ \approx Y^{+} $$ This sounds great to hedge downside risk, but it doesn't really work in practice unfortunately. Any advice?


title: Model simplification results in non normal residuals

text: I have a data set that includes a response variable (counts), two categorical explanatory variables, and one covariate. The response variable being counts, I ran a GLM with poisson distribution, and simplified the model. While the saturated model showed normal and homoskedastic residuals, as I reach the minimum adequate model, they fail normality and homoskedasticity. I'm a bit confused as to what the best procedure would be, and/or what the most important assumption. If I'm just looking to assess whether the interactions have a significant effect, would it be appropriate to simply transform the data in order to do an ANCOVA? Thanks for your help.


title: Restricted Boltzmann Machines for regression?

text: I'm following up on the question I'd asked earlier on RBMs. I see a lot of literature describing them but none that actually talks of regression (not even classification with labelled data). I get a feeling that it is used for unlabelled data only. Are there any resources for handling regression? Or is it as simple as adding another layer on top of the hidden layer and run the CD algorithm up and down? Thanks much in advance.


title: Can I use t-test to compare my control and sample mean to confirm significant difference?

text: My control is a sample without any drug treatment. My samples are treated with drugs. They were both randomly taken from a normal distribution. The hypothesis is $\text{H}_0: \mu = \mu_0$ and $\text{H}_\text{a}: \mu \neq \mu_0$ The calculation would be $t = (X-\mu_0)/(S/\sqrt n)$ and then confirm whether or not the control and the sample mean are significantly different. $\mu_0$ is based on the mean taken from the control, and X comes from sample mean Does this make sense?


title: How to plot an interaction term, using a model's coefficients, of three-factorial GEE model with full-order interactions (geepack package)?

text: Using function `geeglm` from package `geepack` (Generalized Estimating Equation), I have modeled counts as being dependent on two nominal (factor) variables, one continuous variable with 3rd-order interactions and one grouping variable. This is the model:               m1<-geeglm(formula = dependent.var  ~  cat.var1 *  cat.var2 * contin.var,                  family = poisson, id = group, corstr = "exchangeable")      Factor variable _cat.var1_ has two levels (CD and WL), factor variable _cat.var2_ has two levels (SAH and SKA). _Group_ is a grouping variable to account for autocorrelation among somehow related subjects. Using `anova` function, I found that the only significant model’s interaction terms are **cat.var1 : contin.var** and **cat.var2 : contin.var**. I would like to plot these two interactions: I have used the codes below, but I am not happy with the fitted lines. **Line for CD-level of _cat.var1_ in the first figure is the same as line for SAH-level of _cat.var2_ in the second figure (see figs) – therefore I think that they are not correct and I do not know how to solve this. Which model's coefficients I have to use to construct that lines?** The model’s coefficients used for producing of the figures are thereinafter. ![cat.var1 : contin.var](http://i.stack.imgur.com/5omoS.jpg) ![cat.var2 : contin.var](http://i.stack.imgur.com/e414Y.jpg) Codes:               # 1st plot: cat.var1 : contin.var     par(mfrow=c(1,1))     plot(contin.var, dependent.var, type="n")     points(contin.var[cat.var1=="WL"], jitter(dependent.var[cat.var1=="WL"]))     points(contin.var[cat.var1=="CD"], jitter(dependent.var[cat.var1=="CD"]),            pch=16)     x <- seq(3,7,0.1)     y1 <-exp(-0.1584 + 0.3474*x)     lines(x,y1, lty=2)                                 # CD     y2 <-exp(-0.1584 + 3.7293 +(-0.6685 + 0.3474)*x)     lines(x,y2, lty=1)                                 # WL          # 2nd plot: cat.var2 : contin.var     par(mfrow=c(1,1))     plot(contin.var, dependent.var, type="n")     points(contin.var[cat.var2=="SKA"], jitter(dependent.var[cat.var2=="SKA"]))     points(contin.var[cat.var2=="SAH"], jitter(dependent.var[cat.var2=="SAH"]),             pch=16)     x <- seq(3,7,0.1)     y1 <-exp(-0.1584 + 0.3474*x)     lines(x,y1, lty=2)                                 # SAH     y2 <-exp(-0.1584 + -6.9490 + (1.2249 + 0.3474)*x)     lines(x,y2)                                        # SKA      Here are my model's coefficients, used for producing of figures, outputs from `summary(m1)` (statistically significant terms - outputs from `anova(m1)` \- are denoted with asteriks **):               Coefficients:                                              Estimate          (Intercept)                                  -0.1584        cat.var1WL                                    3.7293        cat.var2SKA                                  -6.9490       contin.var                                    0.3474   *     cat.var1WL: cat.var2SKA                       3.9970        cat.var1WL: contin.var                       -0.6685   *     cat.var2SKA: contin.var                       1.2249   **     cat.var1WL: cat.var2SKA: contin.var          -0.6860             Factor variable *cat.var1* has two levels (CD and WL),     factor variable *cat.var2* has two levels (SAH and SKA).      Here are the data: https://docs.google.com/file/d/0Bz8ojhHeiNclVi1oT0ZwTEtEN2s/edit


title: ANCOVA/ANOVA: Testing different models against each other

text: I have a dependent variable $Y$, an independent variable $X$, and a categorical independent variable T (which can take 2 levels). Now I have 4 models: **MI:** different slopes, different intercepts $$y_{i1}=\alpha_1+\beta_1x_{i1}, y_{i2}=\alpha_2+\beta_2x_{i2}; \text{or}\; y_{ij}=\beta_0+\beta_1X_{ij}+\beta_2T_i+\beta_3X_{ij}T_j$$ **MII:** same slope, different intercepts $$y_{i1}=\alpha_1+\beta x_{i1}, y_{i2}=\alpha_2+\beta x_{i2}; \text{or}\; y_{ij}=\beta_0+\beta_1X_{ij}+\beta_2T_j$$ **MIII:** same intercept, different slopes $$y_{i1}=\alpha+\beta_1x_{i1}, y_{i2}=\alpha+\beta_2x_{i2}$$ **MIV:** one regression line for both $$y_{ij}=\alpha+\beta x_{ij}; \text{or}\; y_{ij}=\beta_0+\beta_1X_i$$ Now I am trying to compare these models one against the other. I am not sure, if I am right with what I think:   Comparing MI to MII tests $H_0:\, \beta_1=\beta_2$ (this I do know), which means that $H_0$: same slope against $H_a$: different slopes (which I am not sure of). If I compare these in R, I get a p-value. Now what I think is that, if the p-value < alpha level, I can reject $H_0$, which would mean I should stick to the first model. If it is a large value, I cannot reject the null hypotheses, which means I should stick with the second model (as there I have the same slope). I am not sure if I am right with that.   Plus I cannot write down the second equation for the model third model, can someone help out?


title: Good GLIM Books

text: I finished my B.S in Statistics in the early 80's. Recently, for a variety of reasons, I've got interested in Generalized Linear Models. Would you guys please point me some good books on the theory and applications of GLIM's? Any references with lots of worked examples, proofs of the main results, and many exercises would be great. A solutions manual for self study would be very useful. Thanks!


title: How to decide between pearson, distance correlation or maximal information coefficient for currency data?

text: I have one dataset with currency exchange ratios for several years. The exchange ratios are for two Asiatic currencies and the Canadian dollar, all against American dollars. My (perhaps naïve) thinking is that I can observe which economies are more strongly tied by looking to the association strength between the ratios. So, I have got an Excel add-in that can calculate distance correlation and maximal information coefficient, plus with some tweaking it can even give p-values of Pearson. However, when I calculate the three measures for all pairs of variables and their p-values against the null hypothesis that the variables are independent, they give utterly discordant results. Which measure should I believe better? Is there anything wrong with my analysis? ![Screenshot with the problem](http://dl.dropbox.com/u/5363697/currency_screeshot.png) It is possible to check both the Excel worksheet and a PDF version of it. Thanks in advance for any help.


title: Entering values into a linear regression / gradient descent algorithm

text: The instructor of the machine learning course I've been taking whipped this algorithm at my head without explaining how to apply it to a training set. ![enter image description here](http://i.stack.imgur.com/k8jlS.png)               m = number of training sets (3)     j = subscript 0 or 1     x(i) = x from the ith training set     y(i) = y from the ith training set     alpha = ?     h theta = ?      If I had a training set of (1,1), (2,2), (3,3), how would I go about entering it into this algorithm?