title: Why does bootstrap gives different p-values for CC and my regression coefficient in simple linear regression?

text: I'm using R and the package boot. my boot function returns the coefficients and the r.square and sqrt of the r.square                rsq2 <- function(formula, data, indices) {          d <- data[indices,] # allows boot to select sample           fit <- lm(formula, data=d)          return(c( coef(fit),summary(fit)$r.square,sqrt(summary(fit)$r.square)))      }          boot(data=mtcars,rsq2,1000,formula="mpg~wt")          ORDINARY NONPARAMETRIC BOOTSTRAP               Call:     boot(data = mtcars, statistic = rsq2, R = 1000, formula = "mpg~wt")               Bootstrap Statistics :           original        bias    std. error     t1* 37.2851262  0.1170681005  2.32470420     t2* -5.3444716 -0.0557873180  0.70523904     t3*  0.7528328  0.0019940314  0.05758118     t4*  0.8676594  0.0004980607  0.03362469          I determine a pvalue using the original value and the std. error determined by the boot program          pvalue of the regression coefficient     2*pnorm(-abs(-5.3444716/0.70523904))     [1] 3.502706e-14     pvalue of the CC      2*pnorm(-abs(0.8676594/0.03362469))     [1] 7.947938e-147      Why are these not the same value? In a simple univariable the linear regression of the p-value of the coefficient and Pearson's CC are the same value?


title: Risk assesment and non-statistician's perception of percentages

text: I'm planning on doing quantile regression for estimating pain relief after hip surgery. My aim is to construct for my patients information in this format: > 4 out of 5 patients experience with your background experience at least a 50 > % decrease in pain intensity one year after surgery My problem is: Is 80 % perceived as very probable? Is 90 % much better? What I'm looking for is some kind of **reference article** on _how people perceive percentages_ to determine a good quantile. I've been looking on Google scholar for something useful but strangely I haven't found anything that I can use. There is plenty of articles on how to communicate risk, a good summary can be found here, but none of them discuss how percentages in this fashion. I guess what I'm looking for is something like the p-value graph (originally posted in this question) but with "John Doe's" instead of "Scientist's". ![Magical p-values](http://i.stack.imgur.com/PNBph.gif)


title: Skills & coursework needed to be a data analyst

text: I am graduating with a bachelor's degree in applied math and I will pursue a master's degree in statistics this fall. There are many specialized fields in applied statistics. I realize I may be more interested in becoming a data analyst, working with big data in market research or IT, but I'm not sure what kinds of skills are needed and what are the core courses I should take to become a data analyst. I've learned some MatLab, R, SPSS, C, ruby, can't say I'm an expert on these. I'm learning SAS, because SAS seems to be what every stats master should know. (Is this really true?) I'm a little confused right now and need some specific advice about what skills I should intensively learn and what courses I should take. TIA.


title: Linear factor analysis to find stationary factors?

text: I have a multivariate time series. My model assumption is that it was generated from underlying factors, some of which are stationary, using a linear mixing. How can I access the stationary factors? I have tried independent component analysis (ICA), but the independence assumption is too strong.


title: Hausman test in Stata : rank problem

text: When I run `hausman` in Stata, there is the following message:               Note: the rank of the differenced variance matrix (8) does not equal the number of      coefficients being tested (9); be sure this is what you expect, or there may be problems      computing the test.  Examine the output of your estimators for anything unexpected and      possibly consider scaling your variables so that the coefficients are on a similar scale.        I searched on the internet and they suggested me to use `hausman, sigmamore` or `hausman, sigmaless` but there is the same message. What should I do to fix this?


title: positive samples in dataset for classifier training - how should environment figure in?

text: I have a question about preparing the dataset of positive samples for a cascaded classifier that will be used for object detection. As positive samples, I have been given 3 sets of images:   1. a set of **colored** images in full size (about 1200x600) with a **white background** and with the object displayed at a different angles in each image   2. another set with the same images in **grayscale** and with a **white background** , scaled down to the detection window size (60x60)   3. another set with the same images in **grayscale** and with a **black background** , scaled down to the detection window size (60x60) My question is that in Set 1, should the background really be white? Should it not instead be an **environment** that the object is likely to be found in in the testing dataset? Or should I have a fourth set where the images are in their natural environments? How does environment figure into the training samples?


title: Please help to understand Agresti's book

text: I need some help to understand some topics in Agresti's _Categorical Data Analysis_. In section 6.3.1 (p 231), he provided a model like:   $$ \text{logit}(\pi_{ik})=\alpha+\beta x_i+\beta_k^Z $$ where $i$=1,2, $k$ = 1, $\ldots$, $K$ Basically we have 3 categorical variables $X$, $Y$, & $Z$ ($Z$ is a confounding variable), and we are trying to test the conditional independence of $X$ & $Y$ for each level of $Z$. He says that "this model assumes that [the] $XY$ conditional odds ratio [is the] same at each category of $Z$ namely exp($\beta$)" I cannot understand this statement. How is the odds ratio exp($\beta$)? What about $\beta_k^Z$? isn't this also varying for each level of $Z$? Then the odds ratio must be different $Z$-wise? What I am missing here? Thanks for your time


title: Subspace clustering in R using package orclus

text: Currently I am working on some subspace clustering issues. I found one useful package in R called orclus, which implemented one subspace clustering algorithm called orclus. As stated in the package description, there are two key parameters to be determined. One is the subspace dimensionality and the other one is the cluster number. It is stated that to determine the optimal value of subspace dimensionality, one statistic, the cluster sparsity coefficient can be used. The closer the statistic to zero, the better the performance. However, when actually trying this implementation, I found that the statistic is minimal when subspace dimensionality is 1; and the larger the subspace dimensionality, the larger the statistic. Does it make sense? I was not expecting such monotonic trend.


title: Regarding kernel-based naive Bayesian classifier

text: Are there any good references for kernel-based Naive Bayesian classifier?


title: Factor analysis for ordinal variables that have different categories

text: I have a data set that contains about 40 categorical variables that are taken as independent variables (and believed to be related to some unobservable human resource factors) and 4 categorical variables (like a company's turnover, competition for job etc.) that are taken to be the dependent variables. All these variables categories can be ordered although the variables have different number of categories. I want to see how the factors actually affect the dependent variables. I know how to do factor analysis with continuous variables. How do I do it with these kinds of categorical variables with different categories? Besides there is another categorical variable "size" representing the number of full time staffs ($\le 25$ is given 1, $26-99$ is given 2 and $>100$ is given 3) which is planned to use as a mediating or moderating variable. I want to see how the human resource factors affects the dependent variables. I mean, to make comments like “For small companies (smaller than 25 employees), human resource factor-1 impacts on profitability quite dramatically, however, there is no relationship between voice and profitability in medium and large hospitals...bla bla bla” (for example) How do I relate the factors with the categorical dependent variables? Can I do whatever you suggest in R or in SPSS?


I'm trying to do a logistic regression on some data. Here's a simplified version of the situation: I'm trying to predict student success based on their history, etc. One of my predictors is the percentage of the courses they've passed in the past. If they haven't taken any courses, I don't want to set that to zero, because that's obviously different from having failed all their courses. Right now, these cases are set as NaN, but when I use the glm function in R, I get the following error: > Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, : > NA/NaN/Inf in foreign function call (arg 1) **How do I predict performance for individuals who haven't taken any courses yet?**


title: PLINK meta-analysis method

text: I am trying to determine exactly what method PLINK uses for meta-analysis. I can't find it in the program documentation/FAQ anywhere. Inverse variance? (given the meta command only requires effect size and SE)? ALso, what method its using to generate a P-value?


title: At least, how many times an experiment should be replicated?

text: Is there a recommendation on the number of times that an experiment should be replicated? As many of you know, is not always possible to make many replicas. What would be the recommended minimum? Is there some references to support it? In my particular case (animal reproduction), for reasons of seasonality, I can only replicate experiments 3 times and I have sometimes been criticized for the low number of replicates performed. Could be considered appropriate to assess the effect that a parameter measured 3 times in the same individuals have on the performance of these individuals?


title: Whether to use multinomial logistic regression with a three-category outcome variable and many binary predictors?

text: I have a dependent variable made up of 3 categories and 14 binary predictor variables. I have tried using `mlogit` and `nnet/multinom` packages in R. **Is there a better approach than multinomial logistic regression for this particular scenario?**


title: Is it possible to apply Bayes Theorem with only samples from the prior?

text: I was just wondering-- is it possible/practical to apply bayes' theorem without an analytical expression for the prior, only samples? For example, say you have sufficient draws from a posterior distribution from a previous experiment via MCMC methods, and you'd like to use that posterior as the prior for a new one. You have an analytical expression for the likelihood as before, but now only samples from the (new) prior. How would you proceed?


title: Confidence and prediction intervals of linear regression model

text: Okay, so I am trying to understand linear regression. I've got a data set and it looks all quite alright, but I am confused. This is my linear model- summary:               Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) 0.2068621  0.0247002   8.375 4.13e-09 ***     temp        0.0031074  0.0004779   6.502 4.79e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.04226 on 28 degrees of freedom     Multiple R-squared: 0.6016, Adjusted R-squared: 0.5874      F-statistic: 42.28 on 1 and 28 DF,  p-value: 4.789e-07       so, the p-value is really low, which means it is very unlikely to get the correlation between x,y just by chance. If I plot it and then draw the regression line it looks like this: http://s14.directupload.net/images/120923/l83eellv.png (Had it in as a picture but I am - as a new user - currently not allowed to post it) Blue lines = confidence interval Green lines = prediction interval Now, a lot of the points do not fall into the confidence interval, why would that happen? I think none of the datapoints falls on the regression line b/c they are just quite far away from each other, but what I am not sure of: Is this a real problem? They still are around the regression line and you can totally see a pattern. But is that enough? I'm trying to figure it out, but I just keep asking myself the same questions over and over again. What I thought of so far: The confidence interval says that if you calculate CI's over and over again, in 95% of the times the true mean falls into the CI. So: It it is not a problem that the dp do not fall into it, as these are not the means really. The prediction interval on the other hand says, that if you calculate PI's over and over again, in 95% of the times the true VALUE falls into the interval. So, it is quite important to have the points in it (which I do have). Then I've read the PI always has to have a wider range than the CI. Why is that? This is what I have done:               conf<-predict(fm, interval=c("confidence"))     prd<-predict(fm, interval=c("prediction"))      and then I plotted it by:               matlines(temp,conf[,c("lwr","upr")], col="red")     matlines(temp,prd[,c("lwr","upr")], col="red")      Now, if I calculate CI and PI for additional data, it does not matter how wide I choose the range, I get the exact same lines as above. I cannot understand. What does that mean? This would then be:               conf<-predict(fm,newdata=data.frame(x=newx), interval=c("confidence"))     prd<-predict(fm,newdata=data.frame(x=newx), interval=c("prediction"))      for new x I chose different sequences. If the sequence has a different # of observations than the variables in my regression, I am getting a warning. Why would that be?


title: What does $E_yT_y$ mean?

text: I am confused about the notation in my Stochastic Processes class, and I can't find a place in the textbook that explicitly defines this notation. He uses $E$ to mean simple expectation i.e. for tail sum formula for expectation it is written $$EX = \sum_{k = 1}^\infty P(X \ge k)$$ Then in the next paragraph, he writes that the probability of returning at least $k$ times $\\{N(y) \ge k\\}$ is $$E_xN(y) = \sum_{k = 1}P(N(y) \ge k) = \cdots = {\rho_{xy} \over 1 - \rho_{yy}}$$ where I added the `\cdots` because the calculations aren't relevant. As a direct consequence of this, I am confused about what the "limiting fraction of time we spend in each state" theorem is saying. It goes as follows: > **Theorem 1.21 (Asymptotic Frequency)**. Suppose our Markov Chains is > irreducible and all states are recurrent. If $N_n(y)$ is the number of > visits to $y$ up to time $n$, then $${N_n(y) \over n} \to {1 \over E_yT_y}$$ where $T_y$ is defined to be $T_y^1$ for $$T_y^k = \min\\{n > T_y^{k-1}: X_n = y\\}$$ i.e. the time of first return to $y$.


title: Model selection: where to start?

text: For a general modeling problem, there are literally at least a dozen choices of statistical and algorithmic models to choose from. Off the top of my head, choices could be: regression (and its variants), neural networks (and its variants), SVM (and its variants), random forests (and its variants), etc. So what I'm asking about isn't much about what is referred to 'model selection'. I'm looking for a guide that can help me (intelligently) choose candidate models. Or is this process more of an art? I'm not looking for models specialized for certain data types. I'm looking for models that generically model tabular data. I tried consulting with Google but it wasn't of much help. Surely, there must be a guide on how to begin the model selection process.


title: How is it possible that the normal independent variable is not significant anymore after mean centering in a curvilinear relationship?

text: > **Possible Duplicate:**   >  Two negative beta's in a curvilinear regression when mean centered or using > standardized values The problem I encounter is the following: Imagine a (perfect) inverted U-shaped relation between an independent variable and a dependent variable. When you look at the curve estimation there is indeed a perfect inverted U-shaped relation. If you conduct a regression analysis you have the normal independent variable and the squared term of the independent variable. If these are both significant, this indicates that there is an inverted U-shaped relationship. However, it is well-known that in this way you have a very high VIF value (multicolliniearity), because there is an almost perfectly positive correlation between the normal independent variable and independent variable squared (because you square it of course). The solution for reducing the multicolliniearity is to mean center the independent variable, and after then taking the squared term of the mean centered independent variable. If I perform a regression now (with both of the two predictors mean centered) I do not have problems anymore with multicollinearity, but the normal independent variable is not significant anymore (the independent variable squared is still significant). How is this possible? What is wrong? Is it still possible to prove the inverted U shaped relation when only the independent variable squared is significant? Or is it possible to ignore the multicollinearity issues?


title: How can I choose the 10 variables that explain the most variation in a wealth index?

text: I have household survey data with 32 questions about assets the household has or doesn't. I assume that taken together the answers to these asset questions (e.g. how many televisions does the household own) are an indication of wealth, and could be used to make a good index of wealth, e.g. using the first component in a principal components analysis. What I want to do, however, is to choose 10 of these variables that jointly explain the largest possible proportion of the variation in wealth and use those as the questions in a shorter questionnaire that I am developing. What is the best way of doing this? One possibility that has occurred to me is to calculate the wealth index using PCA then regress this on every possible combination (60 million or so I think) of 10 variables from the 32, and see which gets the highest R-squared. I'm hoping there's an easier way. Ideally I'm looking to implement this in Stata.


title: Dealing with a variable amount of entries for a given attribute

text: I want to analyze a .csv file with the results of a customer survey. One of the questions is: "Tell us the top 3 reasons you chose us.". The answers are short paragraphs of text listing a ranging number of reasons (not always 3). An example answer could be "1. Speed, 2. Customer service, 3. Price". Another might be: "I liked your customer service, and it wasn't that expensive!" I would like to be able to turn this data into an overview like this:               REASONS     Customer Service: 2     Price: 2     Speed: 1      I'm facing two problems:   1. I'm interpreting and splitting up the reasons manually in a spreadsheet. Any tips on how to do this more efficiently are very welcome.   2. I have a field "Reasons" which I want to contain a variable amount of values per record. I don't think spreadsheet software such as Libre Calc supports that, and I suspect I might need to export this to SQL. Am I right in suspecting that? And how would I achieve the desired multiple-fields-per-attribute requirement in SQL?


title: Deriving ${\rm var}(\overline{X})$ from expected value definition

text: I want to show ${\rm var}(\overline{X})$ is $\sigma^2/ n$ for $n$ independent random variables each with mean $\mu$ and standard deviation $\sigma$. I would like to show this using: $$ E(\overline{X} - \mu) ^ 2 $$ I'm doing something wrong in my own calculations because it's not reducing to $\sigma ^ 2/ n$. Any thoughts/comments? Thanks! Diane


title: IV (2SLS) with ordered/binary dependent and endogenous variable

text: thanks for taking time to read this first of all. I am looking for an IV procedure if the endogenous variable is binary or ordinal (0 to 5). I see the problem of using 2SLS if the endogenous variable is ordinal, but I fail to find a solution. I am convinced there is a solution and many people have faced this problem before, so sorry if this is too obvious. Thanks everyone!


title: Numbers 1-100, calculating how many I must pick to get x% of reliability

text: I would like to know the formula or method of calculating this. For example, I want to know how many numbers I must pick in order to have 95% confidence that I will be within 2% of the average. I'm also interested in doing a similiar analysis for the median as well. Update: Numbers are not put back in the pool. I think this has to do with combinations rather than trials? For example, if I pick 1 number or 99 numbers there are 100 combinations, if I pick 100 numbers there is 1 combination etc. I think you could figure out all of the combinations when you pick 1 ball, 2 balls, 3 balls etc and see how many are within 2% of the median. Im interested in seeing just how reliable the average price statistic is for the housing market. For example, in a single month 5% of the houses may sell, or maybe 20% will. Housing prices are all over the map but I wanted to do something simplistic like this. And now I have thought of something else...let's say I have 1-500 instead of 1-100. If I pick from 1-100 and pick 10 balls(10% of them), Im guessing that the answer would differ from taking 50 balls out of 1-500(10% of them as well)?


title: Combining exponentially smoothed variances

text: I'm currently tracking EWMA and EWMV with $M_k=(1-\alpha)*M_{k-1}+\alpha*x$ $S_k=(1-\alpha)*(S_{k-1}+\alpha*(x-M_{k-1})^2)$ Now, I'm tracking these values in data set with 4 different dimensions and would like to "fold" these dimensions to get aggregated views, but I'm not sure how I can combine mean and variance of multiple data sets together when they are exponentially smoothed, as I don't have population size. As dimensions have different scales, surely I just can't assume population of 1 for each partial set?


title: Calculating AIC “by hand” in R

text: I have tried calculating the AIC of a linear regression in R but without using the `AIC` function, like this:               lm_mtcars <- lm(mpg ~ drat, mtcars)          nrow(mtcars)*(log((sum(lm_mtcars$residuals^2)/nrow(mtcars))))+(length(lm_mtcars$coefficients)*2)     [1] 97.98786      However, `AIC` gives a different value:               AIC(lm_mtcars)     [1] 190.7999      Could somebody tell me what I'm doing wrong?


title: How to handle Missing Data

text: Am conducting a study where I teach my participants 10 lessons. After each lesson they have to answer questions to check understanding of the lesson. Some participants missed between 2-4 lessons. How do I handle the missing data? Do I drop them?


title: Pros and Cons of Fitting a Spatial Regression to Cumulative Data

text: When using regression to estimate a relation between a variable of interest and distance from a certain point (eg a distance-decay curve), what are the advantages and disadvantages of fitting the regression line to cumulative data rather than data by band, and what is the balance of advantage? The following is a slightly contrived example designed to illustrate the question while avoiding extraneous issues. Suppose a researcher tags in some way a random sample of 100 nuts growing on a nut tree. Later, when the nuts have fallen, the researcher identifies where they are, measures their distances from the tree, and counts the numbers by 5 metre distance bands, obtaining the following (hypothetical) data: ![enter image description here](http://i.stack.imgur.com/7Oe5c.png) The researcher wants to use the sample to estimate the distance-number relation for the whole population of nuts on the tree. One approach would be to estimate a regression of N on MD. Another would be to estimate a regression of CN on D. The results of either approach can easily be converted to the other form by summing or taking differences. Is there any advantage in terms of accuracy in the latter approach? The advantages that occur to me are these:   1. The mean distance of the nuts within any band might not be the mid-point of the band. If, say, the mean distance is generally less than the mid-point, regression of N on MD will result in bias. By contrast, regression of CN on D is unaffected by the distribution of distances within bands.   2. Both approaches depend upon banding, raising the risk that their results will depend on the particular banding structure chosen (see Openshaw, The Modifiable Areal Unit Problem). However, the cumulative approach seems less dependent: considering N and CN at MD = 17.5 and D = 20, for example, suppose 1 nut at distance 19.99 metres had been excluded by a marginally smaller band. In that case, N would be reduced from 11 to 10, a difference of 9%, while CN would be reduced from 90 to 89, a difference of only 1.1%   3. Using CN avoids the complication of what to do with the zero values of N if a logged functional form appears appropriate. An important disadvantage, as whuber has pointed out in his comment, is interdependency between the CN values. Furthermore, heteroscedasticity in the N values will affect the CN values in a complex way. I am not aware of an estimation method that can handle these features - any suggestions would be appreciated.


title: forecasts ahead of t+1 in time-varying linear state space models

text: When the matrices in the model are constant, then performing forecasts is straight-forward. However, when using a time-varying model like dynamic regression I'm not sure how to proceed since we don't observe the predictors either (are time series too) so my intuition is that the forecasts would be those of the constant part of the model (ignoring the dynamic part) for steps ahead further than t+1. Any thoughs would be appreciated.


title: Overdispersion in logistic regression

text: I'm trying to get a handle on the concept of overdispersion in logistic regression. I've read that overdispersion is when observed variance of a response variable is greater than would be expected from the binomial distribution. But if a binomial variable can only have two values (1/0), how can it have a mean and variance? I'm fine with calculating the mean and variance of successes from x number of Bernoulli trials. But I cannot wrap my head around the concept of a mean and variance of a variable that can only have two values. Can anyone provide an intuitive overview of:   1. The concept of a mean and variance in a variable that can only have two values    2. The concept of overdispersion in a variable that can only have two values


title: Multinomial model on victim reporting behaviour

text: I am trying to fit a multinomial logit model on the data below. I want to see whether gender, age, population and the province have an influence on the reporting behavior of a victim of crime. This is my model:               library("mlogit")     #source("~/Forge/mlogit/chargement.R")          R <- mlogit.data(data, shape="wide", choice = "choice_of_reporting")     m <- mlogit(choice_of_reporting ~ gender_household_head + age_household_head                   + pop + prov_name, R)     summary(m)               mlogit(formula = choice_of_reporting ~ gender_household_head + age_household_head                + pop + prov_name, data = R, method = "nr",               print.level = 0)      But I get the following error:               Error in solve.default(H, g[!fixed]) :        Lapack routine dgesv: system is exactly singular: U[4,4] = 0      Can someone tell me what this error mean and how I can solve it? Thank you for your assistance! An extract from my data:                         DEPENDENT VARIABLE               GENDER             AGE     POP   PROVINCE                Did not report                     1                 56    black EASTERN CAPE     5566       SAPS and other                     1                 42    black EASTERN CAPE     5567   Other but not SAPS                     0                 66    black EASTERN CAPE     5568            SAPS only                     0                 68    black EASTERN CAPE     5576            SAPS only                     0                 60    black EASTERN CAPE     5584       SAPS and other                     0                 76    black EASTERN CAPE     5600            SAPS only                     0                 35    black EASTERN CAPE     5611            SAPS only                     1                 68    black EASTERN CAPE     5612            SAPS only                     0                 52    black EASTERN CAPE     5615       Did not report                     0                 46    black EASTERN CAPE     5667   Other but not SAPS                     1                 59    white EASTERN CAPE     5668   Other but not SAPS                     1                 59    white EASTERN CAPE     5685       Did not report                     0                 24    black KWAZULU-NATA     5689       Did not report                     0                 62    black KWAZULU-NATA     5697            SAPS only                     0                 71    black KWAZULU-NATA     5698       Did not report                     0                 63    black KWAZULU-NATA     5703       Did not report                     0                 85    black KWAZULU-NATA     5704       Did not report                     1                 59    black KWAZULU-NATA     5705       Did not report                     1                 81    black KWAZULU-NATA     5708       Did not report                     1                 42 coloured KWAZULU-NATA     5709       Did not report                     0                 33    black KWAZULU-NATA     5721       SAPS and other                     1                 46    black EASTERN CAPE


title: Augmented split plot experimental design

text: I planned a field experiment with different genotypes according to an augmented design (check repeated entries plus unrepeated entries). Since I am testing a segregating population, each plot split into different plant colours. Therefore, I would like to analyze singularly each colour group within each plot. This looks like a kind of augmented split plot design. How can analyze such an experiment in R?


title: Calculation of confidence interval for estimation of demand curve - how many degrees freedom?

text: I'm trying to understand how uncertain an estimation of a demand distribution of a good is. If I take a random sample of n people from the representative consumer population who give me their personal cutoff/choke prices (assuming they accurately report them), I can construct a demand curve with n points representing an estimation of the proportion of the population who would buy the good at that price. Given that I used n independent pieces of information to calculate n points, should I be using a t-distribution with 1 degree of freedom to calculate the confidence interval for each point? And, if I want to reduce this interval by increasing the degrees of freedom, I either need to take more samples, or increase the granularity of the demand curve to get more d.f. at each point?


title: Cost Benefit Analysis of Pre-screening Widgets for Faults before they Fail

text: I want to build a model that determines whether to pre-screen my widgets for defects. If I do pre-screen, it costs a fixed amount per check and I resolve the problem 100% of the time. If I don't pre-screen, its cost will be the cost of fixing the problem, plus the cost of materials for expired inventory. I have a whole bunch of variables that impact whether a defect will occur in each particular widget, but if it does occur and I didn't pre-screen it, it will definitely result in a fixed repair cost and a variable expired inventory cost based on cost of materials. For each widget, I need to know if the expected cost of the defect exceeds or is less than the cost of pre-screening. Things I know:   * A - Cost of pre-screening   * B - Model to Predict Whether it will fail - based on test data - of ones it says it won't have a defect, it is right 98.5% of the time (false negatives).   * C -Model to Predict expected cost of Raw Materials for each widget produced   * D - % of repairs that will be successful and result in no expired inventory cost   * E - Fixed Cost of Repairs Therefore, I should screen if: A - B*(E+C*D) > 0 I have all of these pieces, but I don't know how to build this into a model. Anyone with R or Stata skills, I would appreciate being pointed in the direction of a good package to use for building this model.


title: Using residualized predictors outside the linear model context

text: Can anyone point me towards a good explanation of when a residualized variable in a regression will give you the same answer as using a non-residualized variable with controls? For instance, say I want to know the effect of a variable x on y and need to control for a and b. In a classic linear model framework I can either add a and b as covariates (i.e., control variables) to the model of y on x, or I can first regress x on a and b, and then use the residuals from this regression (the residualized x) to predict y. Both will give the same coefficient for x. This works in the linear model case, but does a residualized x give the same coefficient as x with controls for other types of models, e.g., logit models or poisson models? My own simple simulations suggest they do not (see R code below), but I am trying to understand why, and if residualization can ever be used in place of adding controls outside of the linear model framework. Can anyone point me towards a good explanation?               #generate the data     n=10000     set.seed(3345)     a=rnorm(n); b=rnorm(n)     x = .4*a + .4*b*b + rnorm(n)     y = .5*x + .3*a + .3*b*b + rnorm(n)          ## LINEAR MODEL ####     #a model with controls gets the right coefficient     summary(lm(y ~ x + a + I(b^2)))     residmod=lm(x ~ a + I(b^2))     x.resid=resid(residmod)     #using a residualized variable gets the same coefficient     summary(lm(y ~ x.resid))          ## LOGIT MODEL ####     y=.5*x + .3*a + .3*b*b + rlogis(n)     ydichot=ifelse(y >0, 1, 0)     #a model with controls gets the right coefficient     summary(glm(ydichot ~ x + a + I(b^2), family=binomial))     #using a residualized variable does NOT get the same coefficient     summary(glm(ydichot ~ x.resid, family=binomial))          ## POISSON MODEL ####     mu=exp(.5*x + .3*a + .3*b*b)     ycount=rpois(n, mu)     summary(glm(ycount ~ x + a + I(b^2), family=poisson))     #using a residualized variable does NOT get the same coefficient     summary(glm(ycount ~ x.resid, family=poisson))


title: Open source libraries in Science

text: Is it a good or a bad practice to use R packages from CRAN for research? I'm talking about the common packages like: simple models for regression, estimation, econometrics Most of them use function that can be written easily on your own. My yes arguments:   * Time to focus on the main part of the research   * The community is a good quality control   * R offers a lot of modern methods   * Some models are too complicated to write them on your own My con arguments:   * Not every package of R has been created in an academic environment   * There could be bugs that influence the outcome and I do not know it   * Most models can be written in a short time period from scratch How can we use open source without risking failures in the outcome? Are there certain quality indicators for packages in general and for R?


title: Is it possible to calculate AIC and BIC for lasso regression models?

text: Is is possible to calculate an AIC or BIC values for lasso regression models and other regularized models where parameters are only partially entering the equation. How does one determine the degrees of freedom? I'm using R to fit lasso regression models with the glmnet function in the glmnet package, and I'd like to know how to calculate AIC and BIC values for a model. In this way I might compare the values with models fit without regularization. Are there any functions that can do this?


title: Difference between post-hoc test after Kruskal-Wallis and Bonferroni corrected Mann-Whitney tests

text: I have some non Gaussian distributed variable and I need to check if there are significant differences between the values of this variable in 5 different groups. I have performed Kruskal-Wallis one way analysis of variance (which came up significant) and after that I had to check which groups were significantly different. Since the groups are kind of sorted (the values of the variable in the first group are supposed to be lower than the values of the variable in the second group which are supposed to be lower than the values of the variable in the third group and so on) I only performed 4 test:               Group 1 vs Group 2     Group 2 vs Group 3     Group 3 vs Group 4     Group 4 vs Group 5      I have performed this analysis with two different methods. I started by using the Dunn's Multiple Comparison Test but nothing came up significant. On the other hand if I use a Mann-Whitney test and correct for the number of tests (4) using Bonferroni, 3 tests come up significant. What does it mean? Which results should I trust? Any help is greatly appreciated.


title: Licensing of analytical code and data - which license to use, and what issues?

text: Not sure if this belongs to CrossValidated, Computational Science, or something else entirely. I've got some code and data that's going into a publication, but I'd like to make said code and data publicly available so that the analysis can be replicated. If I was having it hosted by the publisher, I'd ignore this and let them deal with it, but for a number of reasons to code and data is being hosted externally on my own site. I presume I should put some sort of license on said code and data. Is there a Creative Commons or other type of license that lends itself particular well to open data? Something else I should consider? Just not bother?


title: Should I remove a variable from a PCA analysis if it is the only variable to load on a component?

text: I have four continuous variables drawn from secondary data of companies. The combination of these four variable represent my dependent variable. These four variables are real values. They are just only four values with no items. I am using PCA to combine them together. The results I have got are that three variables are categorized under one component and fourth is categoarized alone under the second component. **Should I consider these two componenents as two dependent variables, one with three variables and the second one with one variable or should I just delete the second component?**


title: Comparing employee sales monthly data despite regional differences

text: I'm not exactly sure how to title this thread, or what tags to use since advanced statistics isn't really my thing right now. Please change the title and tags as needed to improve this question. Say I have 10 different employees who sell newspaper subscriptions, and each employee covers a different geographical region of the state. I have data on the total number of newspapers sold for each employee at the end of each month for a 12 month period, January-December. What I want to know, is how do I normalize the numbers and get a performance evaluation of the employee with complete understanding that some geographical regions naturally perform better than others based on various demographical information. One employee might be doing very well in an area that traditionally has lower volume due to a lower population density, whereas another employee might have sold a ton of newspaper subscriptions but they operate in a region that has a much higher population density than the aforementioned employee. Where do I start? What information do I need to make this happen? How much historical data do I need?


title: Interpreting main effects in the presence of an interaction between other variables

text: I could really use some advice with an interaction interpretation. I have conducted a 3 way repeated measures anova with terms A, B and C. A has 2 levels, B has 3 levels and C has 4 levels. My results show a significant main effect of A and a significant interaction effect of BxC. I am aware that if there is an interaction effect then one should not interpret the main effect. Does that hold true in this situation when the variable A with a main effect does not experience any significant interactions?


title: How is the working correlation matrix estimated for GEE?

text: From a slides > How does GEE work? > >   * First, a naive linear regression analysis is carried out, assuming the > observations within subjects are independent. >   * Then, **residuals are calculated from the naive model (observed- > predicted) and a working correlation matrix is estimated from these > residuals**. >   * Then the regression coefficients are refit, correcting for the > correlation. (Iterative process) >   * The within-subject correlation structure is treated as a nuisance > variable (i.e. as a covariate) > I was wondering in step 2, how is the working correlation matrix estimated for GEE? What is "the naive model"? Thanks!


title: Predicting count data (clicks) based on the words used in a title

text: So I have a data set of articles like so.               ‘Difference of Python and Ruby’    537     ‘Advanced Tutorial on Python’     1438     ‘HTTP library’                     134      The numbers next to them are click data. I am looking to do regression analysis between the words in each of the titles and the number of resulting view the article received, i.e. to draw conclusion on which words are best to use when titling an article for a specific community or group. Also, I would love to hear other ideas as well as to what other cool conclusion I could draw from such data. **Edit** Requested Data I also added unique visits data as well.               [('Travis CI: Announcing Python and Perl support on Travis CI',       {u'clicks': 81, u'unique': 74}),      ('Software Maniacs blog \xc2\xbb Django as a micro-framework',       {u'clicks': 333, u'unique': 318}),      ('Interpreted Languages: Perl, PHP, Python, Ruby (Sheet One) - Hyperpolyglot',       {u'clicks': 236, u'unique': 221}),      ('PyCon 2012 has gone mobile on Guidebook!', {u'clicks': 1, u'unique': 1}),      ('Twitter', {u'clicks': 17, u'unique': 17}),      ('Static Modification of Python with Python: the AST Module - Blueprint Forge.',       {u'clicks': 122, u'unique': 109}),      ('Python for iOS - A Python interpreter and development environment for the iPod Touch, iPhone, and iPad',       {u'clicks': 151, u'unique': 142}),      ("Vivek's blog - Roll your own autocomplete solution using Tries.",       {u'clicks': 211, u'unique': 203}),      ('Twitter', {u'clicks': 4, u'unique': 4}),      ('Ian Bicking: a blog :: Python Application Package',       {u'clicks': 114, u'unique': 111}),      ('PEP 414 -- Explicit Unicode Literal for Python 3.3',       {u'clicks': 50, u'unique': 50}),      ('Static Modification of Python with Python: the AST Module - Blueprint Forge.',       {u'clicks': 0, u'unique': 0}),      ('Requests: HTTP for Humans \xe2\x80\x94 Requests 0.10.6 documentation',       {u'clicks': 18, u'unique': 17}),      ('Haskell, Python and readability : programming',       {u'clicks': 246, u'unique': 236}),      ('PyCon 2012 has gone mobile on Guidebook!', {u'clicks': 40, u'unique': 38}),      ('Pycoders Weekly (pycoders) on Twitter', {u'clicks': 12, u'unique': 12}),      ('SWIX social media optimization, and social network marketing ROI',       {u'clicks': 1, u'unique': 1}),      ('PEP 414 -- Explicit Unicode Literal for Python 3.3',       {u'clicks': 0, u'unique': 0}),      ('[Python-Dev] PEP 414 - Unicode Literals for Python 3',       {u'clicks': 82, u'unique': 82}),      ('Google App Engine Blog: Announcing the General Availability of the Python 2.7 Runtime for App Engine',       {u'clicks': 57, u'unique': 56}),      ('What\xe2\x80\x99s new in Celery 2.5 \xe2\x80\x94 Celery 3.0.21 documentation',       {u'clicks': 237, u'unique': 231}),      ('Python Closures and Decorators (Pt. 2)', {u'clicks': 59, u'unique': 57}),      ('Python for iOS - A Python interpreter and development environment for the iPod Touch, iPhone, and iPad',       {u'clicks': 0, u'unique': 0}),      ('Travis CI: Announcing Python and Perl support on Travis CI',       {u'clicks': 1, u'unique': 1}),      ('Quiz & Learn Python - Learn Python on Your Mobile Device',       {u'clicks': 259, u'unique': 244}),      ('Python Closures and Decorators (Pt. 1)', {u'clicks': 1, u'unique': 1}),      ('Python Closures and Decorators (Pt. 1)', {u'clicks': 219, u'unique': 207}),      ('Trie - Wikipedia, the free encyclopedia', {u'clicks': 18, u'unique': 15}),      ('Ian Bicking: a blog :: Python Application Package',       {u'clicks': 0, u'unique': 0}),      ('Software Maniacs blog', {u'clicks': 42, u'unique': 41}),      ('Quiz & Learn Python - Learn Python on Your Mobile Device',       {u'clicks': 1, u'unique': 1}),      ('Python Closures and Decorators (Pt. 2)', {u'clicks': 0, u'unique': 0}),      ('SWIX social media optimization, and social network marketing ROI',       {u'clicks': 71, u'unique': 69}),      ('Twitter', {u'clicks': 6, u'unique': 6}),      ('Software Maniacs blog \xc2\xbb Django as a micro-framework',       {u'clicks': 4, u'unique': 4}),      ('Haskell, Python and readability : programming',       {u'clicks': 0, u'unique': 0})]


title: Confusion related to derivation of posterior of gaussian processes

text: I was reading this paper related to gaussian processes http://ipg.epfl.ch/~seeger/lapmalmainweb/papers/bayesgp-tut.pdf. I didn't get how they derived the posterior of the gaussian process given some observed data I understand that y follows normal having mean u and std $\epsilon$. Both y and u have normal distribution with mean 0. Only the covariance is different in the case of y an u.Can anyone please explain how the posterior was derived. I didn't get how the $K$ part came in the mean of the posterior and the covariance of the posterior Well u ~ 0 K K_uy y ~ 0 K_uy K+$\sigma^2$ I didn't get how the part K_uy is calculated. They have placed it to be K. I am not sure why. Clarificiations? Also can you explain how the prediction process is derived. I know the integration makes sense. But how did they get the final expression? ![enter image description here](http://i.stack.imgur.com/qpAlE.png)


title: Cross correlation vs mutual information

text: What is the difference between Cross Correlation and Mutual Information. What kind of problems can be solved using these measures and when is it appropriate to use one over the other. Thanks for the comments. To clarify, the question is prompted by an interest in iomage analysis rather than time series analysis although any enlightenment in that area would also be appreciated


title: Naive Bayes with invalid independence assumption

text: I'm trying to understand the effects of adding non-conditionally independent features to a naive Bayes classifier. Let's say I have the features vector $X = [x_1,x_2,x_3,x_4]$ and that for each value of $x_3$ I get the same value for $x_4$: For all $i \in \\{samples\\}$, $x_{3}^{i} = x_{4}^{i}$ I could say that the conditionally independent assumption of $x_n$ given the class $Y = y_k$ does not hold anymore since the value of $x_{3}^{i}$ foresee $x_{4}^{i}$, and that naive Bayes classifier may not produce the expected results. I'm not really sure about that explanation and I would appreciate your point of view about it.


title: What is the point of measuring statistical distance?

text: On pg. 378 of "Cryptography with Tamperable and Leaky Memory", Kalai et al. claim two probability distributions are $e(k)$ close if the distance between them is at most $e(k)$. What is significance of two distributions X and Y being "close to" or "far from" each other? Why would anybody care, especially in cryptography?


title: I can't tell whether my data are seasonal or just affected by calendar effects

text: I used to think my time-serie was seasonal, but then I realized it simply needed some calendar-effect adjustment. I tried that and I'm now in doubt there might still be some seasonality left. I tried adding a 12 month differencing, and it seemed to work, but then I thought: why not just try 12 month differencing in the first place? Shouldn't this account for calendar- effects as well, all in one "bundle"? There you go, runplot + monthly subseries plot + ACF plot of   * raw (top),    * adjusted (top middle),    * 12mths differenced (bottom middle) and    * adjusted+differenced (bottom)  data. NOTE: the order of the images does not follow the order of my reasoning, I put seasonal differencing _before_ calendar adjustment + seasonal differencing: ![enter image description here](http://i.stack.imgur.com/ozT7M.png) ![enter image description here](http://i.stack.imgur.com/h9COn.png) ![enter image description here](http://i.stack.imgur.com/xqWI5.png) ![enter image description here](http://i.stack.imgur.com/bavIy.png) ![enter image description here](http://i.stack.imgur.com/HUJuC.png) ![enter image description here](http://i.stack.imgur.com/RqOWw.png) How do I know which method gave the correct results (if any)? Is this "science" or just "messing with numbers"? Would some further cross-validation help (e.g. stl decomposition, seasonal dummies)? Thank you, as always, for whatever help you can grant me!


title: MCMC for infinite variance posteriors

text: My question originated from Xi'an's suggestion to check integrability against the posterior in my nonlinear hierarchical model. I did not check it, but had possible infinity in mind and found out that one of my conditionals (which is inverse gamma distribution) have infinite variance. So that sampling resulted in chains often being far away at the tails. So, this raised a question for me: How to properly sample from distributions with infinite variances like inverse gamma distribution (with $ \alpha = 2 $) or Levy distribution or any other distribution with infinite variance? What does MCMC offer? I tried to search for papers dealing with such issues, but still no luck.


title: Can I use generalised least squares with a binomial distribution and a nested structure?

text: I'm trying to fit linear models to my data in R. I need to use a generalised least squares method as I have heterogeneity of variance in one of my variables. I was planning to use **varIdent** , as the variable is nominal. But - I also want to have a random term as my data have a nested structure, and my data need to be modelled using a binomial distribution. I can't find any information on whether this is possible using the **gls** command in package **nlme** , does anyone have any information that could help me, please?


title: Why is ezANOVA not giving me Mauchly´s test of sphericity, even though I belive I have more than 2 levels of within-Ss variables?

text: I am trying to do learn how to do a repeated measure in R using `ezANOVA` and I have used the following command:               > ezANOVA(data=NL2, dv=.(ln.conc.), wid=.(sample), within=.(time),     > between=.(pH), type=1, detailed=TRUE)      This is my output and my question is why Mauchuly's test is not appering? I believe I have 9 levels in my within-Ss (`time=A-I`, see dataset below). The p values I get is the same as in SPSS so the test must be sort of correct...:               $ANOVA     Effect   DFn DFd    SSn         SSd         F        p           p<.05    ges     pH        2   6    0.09673652 0.9441442 0.3073784 7.462977e-01           0.04724473     time      8  48    1.08726236 1.0066814 6.4802768 1.065187e-05     *     0.35787718     pH:time  16  48    0.26690606 1.0066814 0.7954038 6.831493e-01           0.12035093      Dataset used (I have converted `time` and `pH` into letters so that R understand that it is a factor and not numericals). `ln(conc)` is measured at 3 `pH`s (`A-C`), 3 different sample sites per `pH` (`A-C`), at 9 different times (`A-I`):               time pH sample  ln(conc)     A   A   A   -0.862749964946125     B   A   A   -0.725670372265505     C   A   A   -0.636766847123838     D   A   A   -0.544727175441672     E   A   A   -0.681218609694672     F   A   A   -0.650087691099498     G   A   A   -0.811930716549912     H   A   A   -0.744440474947496     I   A   A   -0.881889305156823     A   A   B   -0.603306476560156     B   A   B   -0.592397277459802     C   A   B   -0.555125882662571     D   A   B   -0.534435489405125     E   A   B   -0.834710744881732     F   A   B   -0.869884359059999     G   A   B   -0.923818998294947     H   A   B   -1.14256417619729     I   A   B   -1.08767234862978     A   A   C   -0.818710403535291     B   A   C   -0.7339691750802     C   A   C   -0.755022584278033     D   A   C   -0.76142602131324     E   A   C   -0.738144546490681     F   A   C   -0.825536368605691     G   A   C   -0.936493439191674     H   A   C   -1.04412410338404     I   A   C   -1.1679623668029     A   B   D   -0.681218609694672     B   B   D   -0.449416995637347     C   B   D   -0.903868211875598     D   B   D   -0.659712404473708     E   B   D   -0.959720289801491     F   B   D   -0.657780036722654     G   B   D   -0.681218609694672     H   B   D   -0.460449416440924     I   B   D   -0.99695863494161     A   B   E   -0.358104536748327     B   B   E   -0.187535123846842     C   B   E   -0.379797361359587     D   B   E   -0.454130280089445     E   B   E   -0.574475650842447     F   B   E   -0.707246104939447     G   B   E   -0.802962046567152     H   B   E   -0.45886588483528     I   B   E   -1.14570389620196     A   B   F   -0.951917909517306     B   B   F   -0.889162064485903     C   B   F   -0.685179010910768     D   B   F   -0.685179010910768     E   B   F   -0.860383099935859     F   B   F   -0.913793851675568     G   B   F   -0.95711272639441     H   B   F   -0.944175935363691     I   B   F   -0.881889305156823     A   C   G   -0.532730459154041     B   C   G   -0.579818495252942     C   C   G   -0.603306476560156     D   C   G   -0.562118918153541     E   C   G   -0.765717873394781     F   C   G   -0.755022584278033     G   C   G   -1.03282454813011     H   C   G   -0.877070018720874     I   C   G   -1.02443289049386     A   C   H   -0.53102833108351     B   C   H   -0.563874844855806     C   C   H   -0.3147107448397     D   C   H   -0.631111789640493     E   C   H   -0.558616287602339     F   C   H   -0.703197516413447     G   C   H   -0.655851395816248     H   C   H   -0.894040122939335     I   C   H   -0.860383099935859     A   C   I   -1.02722229258144     B   C   I   -0.881889305156823     C   C   I   -0.959720289801491     D   C   I   -0.972861083362549     E   C   I   -0.881889305156823     F   C   I   -0.954511944694353     G   C   I   -0.936493439191674     H   C   I   -1.00785792539965     I   C   I   -0.911303190363116


title: How to test (and accept) that a coefficient in a linear regression model equals zero

text: I understand that in a linear regression model like: $y_i = b_0 + b_1 x_i + \epsilon_i$ I can have a null and an alternative hypothesis: $H_0: b_1 = 0$ and $H_1: b_1 \neq 0$. And then I can reject $H_0$ or fail to reject $H_0$. But what if I want to accept that $b_1 = 0$?


title: Expected value of an indicator function

text: We have $$ I_i (x) = \begin{cases} 1 \quad x-h<X_i<x+h \\\ 0 \quad \text{otherwise} \end{cases} $$ Does that indicator function follow the Bernoulli distribution? If so, the expected value can be computed easily.


title: Is a one sample t test the right approach for this question?

text: This is an example from my book that I'm trying to figure out but I'm totally lost. Can someone complete this step by step so I understand? > A new meditation technique based on making barnyard animal noises promises > to calm people, but you might be skeptical. You can imagine that it could > cause stress to sit around oinking and mooing. To test the therapy, you ask > 16 people to make animal noises, then ask them to rate emotion on a scale > from 1-7, with 1 representing “Very Calm,” and 7 representing “Very > Stressed.” You happen to know that people who don’t oink and moo have an > average emotion rating of 4, with a standard deviation of 2.57. Below are > the ratings from the 16 people who made animal noises.                4  5   6   5   4   3   7   4              3  4   4   3   5   6   5   4      > Write the null and alternative hypotheses for this example.   >  Write the specific type of analysis you will use.   >  Zobt? Z crit?   >  Provide a statement of a significant or nonsignificant effect.   >  Write the obtained value. I believe my null and alternate hypotheses are: $\text{H}_0: \overline{X}= μ\quad\text{vs}\quad \text{H}_1:\overline{X} \neq μ$ for $t$ I did $\frac{4.5-4}{2.57/\sqrt{16}}$ and got $.778$? Is this correct? My $z_\text{crit}$ would be $2.131$?


title: Is there an implementation of Item-Item Collaborative Filtering in R?

text: I was reading about item-item collaborative filtering at: http://guidetodatamining.com/guide/ch3/DataMining-ch3.pdf But I did't find an implementation in R of this algorithm. Is there one?


title: Interpreting fixed effects coefficients with categorical independent variables

text: I am trying to find out if change in marital status have effect on person's health. I am using 2 data points panel data. ![Random and Fixed Effects](http://i.stack.imgur.com/SLM2H.jpg) My question is, how to interpret fixed coefficients if independent variable is categorical? Does coefficient in fixed effect mean that if person change marital status from married (reference category) to divorced, his health get better by .25? Reason why I ask (other than this within change disproves my theory), if I make for example "widowed" category reference, I still get fixed effects. It's really unlikely for person to get from "widowed" to "divorced" in 4 years (50+ population). So, am I interpreting the fixed effects with categorical independent variables wrong? Thanks! * * *


title: Best statistical solution for text mining

text: I'd like to choose a statistical solution for the next years which can manage the major forms of text mining analysis, like mood detection. Is there a solution which stands out from others (R, SAS, SPSS,...) ?


title: How can I statistically formulate football player attributes similar to Football Manager games or footballdatabase.eu?

text: I have collected player match history statistics like runs, shoots, passes and such. I would like to come up similar attribute values (psychical, mental, technical) to Football Manager game or footballdatabase.eu . What would you suggest to come up such a statistical formulation for each attribute based on known numerical statistics? What are the resources and headings that might be helpful to me ? Any recommendation, suggestion or comment is very appreciate :) Here is a example page from football manager. ![enter image description here](http://i.stack.imgur.com/dO2SC.png)


title: Transformation of the data

text: Which Transformation should i apply for the data which having 0's and large values? min value is - 0 and max value is - 2254 here we should not delete any cases from the data.


title: Clustering with K-Means and EM: how are they related?

text: I have studied algorithms for clustering data (unsupervised learning): EM, and k-means. I keep reading the following : > k-means is a variant of EM, with the assumptions that clusters are > spherical. Can somebody explain the above sentence? I do not understand what spherical means, and how kmeans and EM are related, since one does probabilistic assignment and the other does it in a deterministic way. Also, in which situation is it better to use k-means clustering? or use EM clustering?


title: Random forest importance differs between rf$importance and importance()

text: My model is working ok (the AUC is 0.7) but the importances from a `randomForest` run for my binary classification problem differ depending on how I retrieve them. Is this normal? They seem to be scaled up when I call `importance(rf)` (by order 105). More importantly, the first three columns don't match at all, even discounting the magnitude! I'm unsure how to interpret this at the moment having read `?importance` and `?randomForest`. Here's my randomForest call:               rf <- randomForest(trnX, as.factor(trnY), ntree=1000, importance=TRUE)      and here's a portion of the outputs: `rf$importance`:                                  0             1 MeanDecreaseAccuracy MeanDecreaseGini     var1   -2.308793e-05 -6.124117e-05        -3.642557e-05         30.47050     var2    6.169346e-04 -3.947637e-04         2.563570e-04         48.04932     var3    3.621287e-03  1.834355e-03         2.981152e-03         68.28302     var4    3.655234e-03  1.981978e-03         3.057267e-03         79.12254     var5    5.350649e-03  1.041555e-03         3.812376e-03         84.88832     var6    3.366199e-03  1.707144e-03         2.773894e-03         78.02293     ...      `importance(rf)`:                               0          1 MeanDecreaseAccuracy MeanDecreaseGini     var1   -0.4739717 -0.8497967           -0.8952541         30.47050     var2    7.9028043 -3.6680758            4.1625329         48.04932     var3   17.1279919  5.9733456           20.8354079         68.28302     var4   17.3895593  7.2370419           22.4948931         79.12254     var5   20.7470555  3.2702612           24.3069376         84.88832     var6   16.9071520  6.2532409           21.1648686         78.02293     ...      My data has 8000 data points and 60 variables. There is definitely correlation between variables and they are of different types: some are binary flags, other are integers, others continuous. Some variables have missing values (which I'm fixing with `na.roughfix()`). Additionally, a run from gbm:               gbmBernModel <- gbm(formula = var1 ~ ., distribution="bernoulli",                          data = trn, verbose="CV", n.trees=10000, interaction.depth=2)      produced relative importances (viewed by `summary(gbmBernModel)`) that made logical sense to me. Any pointers would be greatly appreciated. Thank you.


title: Truncated Explanatory Variable

text: I am running a model where one of the explanatory variables is truncated. In particular, the variable measures the duration of unemployment (retrospectively) in months and it is truncated at 2 years, which happens for around 10% of the sample. In my model I want to treat the variable as continuous (as I already have quite a few categorical variables), but I am not sure how to deal with the truncated values - it seems that I have a continuous variable, which is also categorical at the same time (i.e. due to truncation). The model will be something like: Y = A + B*X + C*Z + E, where Z<24 Y = A + B*X + C*24 +E, where Z>=24 where X includes some covariates and Z is the truncated variable of interest (continuous) Issue: How to deal with Z in the regression?


title: Matlab - Single Layer Perceptron - output?

text: I'm trying to train single layer perceptron. I don't understand one thing. How to calculate output for vector input x? For example let say we have this data to use for training:               L1 = [1.5 2.4; -0.3 3; 2 3; 2 2.5];     L2 = [2.5 1.4; -1.3 2; 3 3.1; 1.2 2.5];      weight init, step:               step = 0.1     w = [1 1];      Then L = [L1; L2]; Then we have this formula to calculate mean sum square sum_all - I mean sigma symbol which means sum all elements.               E = 1/N * sum_all(d(j) - o(j))^2      j - is index of vector inside L matrix. Then new weight equation:               w(new) = w + step/N * sum_all(d(j) - o(j) * o(j) * (1 - o(j)) * x(j))      Then we have this logical function (pseudo code as it may look clearer):               d(j) = 1 if x(j) belongs to L1     d(j) = 0 if x(j) belongs to L2      And I don't understand how I have to get o(j) or in other words output. I didn't find anywhere where it would say how you get it. It kind of not clear for me. Should I do something like this?               if w * x(j)' > 0 then o(j) = 1     elseif w * x(j)' <= then o(j) = 0      But then if I run my algorithm, it stops training after one epoch and of course does not train it at all (as actual output won't differ from desired output in a first epoch at all). If I try something like this: o(j) = w * x(j)', then numbers go so high that matlab stop recognizing it as number and it goes into infinite loop.. So how I should get actual output x(j) - o(j) using single layer perceptron and backpropagation algorithm?..


title: Binary classification to evaluate a object classifer. How to merge two samples/groups (different size)?

text: I'm evaluating a classifier to detect objects. For example:   I have a set1 with 50 pictures with the object to detect.   And set2 with 100 pictures without the object. I determine the true positves (0 in set2), true negatives (0 in set1), false negatives and false positives for those sets. If I want to merge set1 with set2 (to have a set with and without objects and real values for TP and TP) I have to decrease the size of set2. **Question:** Can I do this with a simple division? (divide the values of TN, TP etc by 2) or do I have to randomly pick 50 pictures of set2 and use these TN's etc.? Edit: **Example:**               Set1: TN = 43 FP = 6 FN = 1 TP = 0     Set2: TN = 0 FP = 13 FN = 0 TP = 87     Set2 divided by 2: TN = 0 FP = 7 FN = 0 TP = 44          Set12: TN = 43 FP = 13 FN = 1 TP = 44      I need an explanation to back this up. thx..


title: Comparing 2 different sets of features

text: So let's say I have two different sets of features A and B. I'm trying to determine which set of features is the best. I'm using leave-one-out cross validation as the final metric as my data set is small. I'm trying to figure out my experimental setup and I'm deciding between a few ways: 1) Give feature set A to my classifer (and optionally run feature selection) on it and give feature set B to a different classifier (also optionally run feature selection) and then compare the LOOCV error between these 2 classifiers? 2) Give feature set A AND B to the classifier and then definitely run feature selection it and then draw a higher level conclusion based on the features selected. (e.g. if more from A are selected state that feature set A appears to hold more prognostic value) 3) Some other way which I dunno


title: Total probability theorem with normal probability density functions

text: Let $X$ and $Y$ be two continuous random variables. Suppose that:   1. $X$ has normal pdf with mean mu_1 and variance $\sigma_1^2$   2. $Y|X$ has normal pdf with mean $mu_2=(a+b*X)$ and variance $\sigma_2^2$; where a and b are two constants from linear regression   3. $f(Y)=\int f(Y|X)f(X) dx$; where $f()$ is the pdf Question: How do I show that $f(Y)$ is a normal distribution? How do I determine its parameters based on the information I'm given? Note: I have tried to brute force multiply the densities and integrate, which yields the answer to my question. However, is there a faster way to see this because I'm interested in extending to the case where $X$ is a random vector?


title: How to take mobility into account: multiple regression

text: I have a panel including two years and I want to compare urban and rural wages: $\log \text{wage}_{it} =\beta_0 +\beta_1\text{urban}_{it} +\beta_2\text{educ}_{it}+\beta_3\text{exper}_{it}+ \beta_4\text{exper}_{it}^2+\delta_2+\nu_it$ urban is a dummy as well as $\delta_2$ for the second year How to I take into account that people might have changed location over time? Wealthy people might have moved to the city and thus raised the difference in urban and rural wages. Should I add a dummy variable or even an interaction variable..?   I'm working with R and I use fixed effects and first difference regression.


title: Asymmetrical Order Effects in Within Subject Design

text: I am currently running a study with two groups (controll and current psy. disorder). All participants have three conversations, are videotaped and watch the videos of their interactions afterwards, two times. While watching the video, I assess my dependent variable. The tree conversations are the three levels of the within-subject factor. It starts with a neutral conversation, followed by a medium emotional conversation and a strongly emotional conversation. Now we are considering, to counterbalance the order in which participants watch their videos to prevent learning and fatigue effects. However, I assume that the effect of for example the highly emotional video on any subsequent video is higher than all other combinations. It is not possible, to run a between subject design. Are there any other options I am not aware of? Perhaps just asking participants how tired they are and how comftable they feel with the tool? Best, A.


title: ANOVA with huge dataset - use only the mean for each condition?

text: I need advice about how to carry out an ANOVA. I studied some theory of ANOVA, but apparently it is not enough. Basically I collected around 300 reaction times for 12 subjects in my experiment. For each subject I had 2 conditions. The first has 3 levels, the second has 5 levels. Each subject went through all the possible level combinations, so I have 50 reaction times for each combination (15 different combinations), for each subject. I want to prove the statistical significance of the different reaction time with different condition. What is not clear to me is the following point: Should I use ALL the dataset, or only the MEAN RT for each subject / condition (that is, 15 data points for each subject) ? If I use the whole dataset, I have about 10,000 data points. If I do so, the degrees of freedom for the Error is 10,557, which confuses me. In the reaction time papers that I have around they report an F(x,y) where this x and y are usually really small, not more than 100, but their dataset is usually extremely big, more or less like mine. This make me think that they don't actually use all the reaction times, but only the mean reaction times (I suppose that, during the computation of the ANOVA, the mean of the means is then calculated). But does this make sense, statistically?


title: Standard Deviation and Variance in Sample and Population Formulas for all?

text: What are the differences between those Population and Standard Deviation and what are their fomulas? I was told that they are different formulas and which situation do you use population or sample? If a user gives 10 random numbers, would I be using population or sample?


title: difficulties evaluating if a ordinal multnomial model fitted with polr (in Rstudio) can be used for predicting over new data

text: I'm getting very different output when i perform a multinomial ordinal model depending on wether I scale or not my data first. Anyone could tell me why? Here are the results with raw data:                polr(formula = result ~ 1 + Speed + rsvl + rsvl:Speed + tailcontrst:Reltailsize +       strategy:rsvl, data = results)          Coefficients:                              Value Std. Error t value     Speed                   -0.1758606  0.0149032 -11.800     rsvl                     7.9648646  0.0345985 230.209     Speed:rsvl              -0.3342932  0.0259837 -12.865     tailcontrst:Reltailsize  0.0004315  0.0001678   2.572     rsvl:strategyburrowing  -0.8982385  0.0746702 -12.029          Intercepts:                 Value    Std. Error t value      fail|tail  -0.1186   0.0689    -1.7216     tail|body   0.4728   0.0669     7.0701          Residual Deviance: 20457.78      AIC: 20471.78       Here are the results with scaled data:               polr(formula = result ~ 1 + Speed + rsvl + rsvl:Speed + tailcontrst:Reltailsize +      strategy:rsvl, data = results)          Coefficients:                                 Value Std. Error  t value     Speed                   -0.97646    0.02083 -46.8841     rsvl                     0.30276    0.02726  11.1064     Speed:rsvl              -0.06810    0.02097  -3.2467     tailcontrst:Reltailsize  0.01092    0.01895   0.5760     rsvl:strategyburrowing  -0.02332    0.03841  -0.6069          Intercepts:                 Value    Std. Error t value      fail|tail  -0.5332   0.0208   -25.6952     tail|body   0.0516   0.0202     2.5480          Residual Deviance: 20608.95      AIC: 20622.95      The best AIC is given by the non transformed data so I would rather use them, if none here tells me the contrary. The model is suposed to be the best combination of factors among all the possible combinations. I obtained it using the function GLmulti. However, despite the t values are quite big and the graphs make sense with the existence of significant effects, the chi-squared test tells me the fit is very bad (value=0). Someone could tell me how do I evaluate if i can rely on this model to make predictions based on new data? Thanks in advance.


title: MAP estimate of posterior parameters

text: I have a setup where the joint posterior is written as: $$ P(w, \lambda, \phi \vert y) = P(\phi) \times P(w \vert \lambda) \times P(\lambda) \times \prod_{i=1}^{N}P(y_i \vert w_i, \phi, \lambda) $$ Now $\phi$ and $\lambda$ are modelled as Gamma distributions and the likelihood and prior on $w$ is normally distributed. So, taking advantage of the conjugate prior property, the posterior on $\lambda$ and $\phi$ will also be Gamma distributed. So, I can derive expressions for $P(w|\lambda, \phi, y)$, $P(\lambda|w, \phi, y)$, $P(\phi|\lambda, w, y)$ because of the conjugate property. What I would like to do is estimate the parameters $w$, $\phi$ and $\lambda$. So my idea was to start with some initial initialisation of $w$, $\phi$ and $\lambda$, and then while keeping the others fixed, update one of them in turn by maximising the log-posterior probability of these derived expressions for the marginal distribution using some simple optimization routine like conjugate conjugate gradient descent. The questions I have are the following:   * Is this the EM way of doing things or have I stated the EM algorithm? I am finding the parameters that correspond to a local maximum mode but in the discussion of EM I always read that this so called latent variables must be present in an EM formulation and I have never got my head around that? What is the difference between the method I have proposed and the EM algorithm?   * Will this approach find me the local optima corresponding to the maximum of the full joint distribution i.e. $P(w, \lambda, \phi \vert y)$. I am doing this iteratively and does the solution suffer because of this? I am guessing that this will get me the joint optimum but could not completely convince myself. These iterations, I guess, correspond to computing the partial derivatives w.r.t to each of the parameters to be estimated in turn, which is similar to optimisation of a multi variable function.   * Also, is this a reasonable way to solve the problem? I know I will only be getting a local solution but are there any other approaches to tackle this problem in this generative model setup?


title: Which is the formula from Silverman to calculate the bandwidth in a kernel density estimation?

text: I want to calculate a better bandwidh for my kernel density estimator, which is an Epanechnikov. I use Silverman's formula which involves the standard deviation of the sample, the sample size and a constant, but I'm getting a very smooth curve in most cases and I would prefer if it were more balanced. Thank you for any help you can give me.


title: Calculating the absenteeism improvement for cabin crew

text: I have the number of monthly flights, and the number of the active cabin crew in monthly basis, and the number of crew who were absence (sick leave or so) in monthly bases as well. How do I get an index or indicator about the absence rate? I want to know are things are improving or not by using all three parameters available. How to achieve that?


title: Finding most significant data set within data

text: Basically, I'm conducting research based on two types of data: Noise levels and the temperature of a room. I've recoded data for 2 days.. I am using the spearman's correlation methods to determine whether or not there is a correlation between temperature changes and the level of noise.. The problem is that there is a lot of data, some of which I don't need. (During the night, for example, there will be no noise since no-one will be in the lab) whereas, in the day time, and evening there will be noise levels. I'm wondering what method could I therefore apply to the data, in order to just extract the most significant parts to perform the analysis. I was thinking PCA and ICA to reduce dimensionality, but I'm interesting in knowing what others think.. Thanks


title: Which model for survey analysis across countries and time

text: I'm analyzing data from a political opinion survey. The survey is administered once a year to roughly 1000 citizens from 18 different European countries. In a first step, I am interested in gauging the relationship between various individual-level variables for all respondents over the past four years. My dependent variable is a dichotomous variable indicating whether the respondent approves or disapproves of the sitting government. Textbooks on nested data/multilevel analysis advise me to use a random intercept hierarchical model. However, it is never explicitly mentioned how I am supposed to account for the time dimension inherent in my data. As I understand it, there are three levels to my data: respondent, country, and year. How should I specify my model accordingly? Thanks a lot!


title: Quiz Competition Regression

text: If you were attempting to forecast the outcome of a general knowledge quiz competition with predictors such as age, education levels, IQ etc there seems to be several obstacles to using say conditional logistic regression, or other classifications into Y=1 (100% winner) or Y=0 100% loser, even if you came second in the quiz. The winner of a particular quiz depends on the questions asked. The winner depends not only on own knowledge but on the other competitors knowing fewer correct answers The quiz round may be very low scoring, the winner just gets a few more questions correct so has little winner ability. The winner may not win any subsequent rounds against the same competitors (is the winner truly a logistic Y=1 winner, in that case?) If the local winner goes on to compete against quizzers from other States then he/she may outclass them or be outclassed. Who knows until after the competition? I have been advised that conditional logistic regression can deal with this. I am puzzled as to how. Ranking the competitors in each round would not seem to work as they could be any position, winner to last, in subsequent rounds. You could add up the ranking each round and the lowest number would be the overall winner and so on and regress on the final outcome data but that might not tell you anything forecasting the overall winner of an inter state competition. Any suggestions?


title: Bernoulli mixture models for image classification, pathological cases

text: I'm trying to use a Bernoulli mixture model to classify MNIST images, and I'm running into pathological cases which screw up my calculations. The pdf of a multidimensional (let's say N dimensions) Bernoulli distribution is as such: $P(x|\mu) = \prod\limits_{i=0}^{N}\mu_{i}^{x_i}(1 - \mu_{i})^{(1 - x_i)}$ The pdf of a mixture of thoses being a weighted sum of them. Dealing with small probabilities, I naturally turned to log probabilities, and transformed the equation as such: $log(P(x|\mu) = \sum\limits_{i=0}^{N}log(\mu_{i}^{x_i}) + log((1 - \mu_{i})^{(1 - x_i)})$ My problem is, what happens when $\mu_i$ becomes either 0 or 1 during training? I end up with $log(0^1)$, which yields -Inf and ruins my classification. Are those case pathological, or is there a way to represent them in a sane way? I thought of considering those cases a pathological, because they typically represent corners of the image on which no pixel ever wanders, and running a PCA to reduce dimensionality, but I'm not sure how I can proceed from there: I can understand reducing my dataset to D' dimensions, and so having my centers of D' dimensions, but how can I use them with an image from my validation set, which is still D dimensions? I feel using the coefficients from my training set will introduce bias. Thanks for your help, and cheers!


title: Programming Multiple Variable PCA Ratios

text: I would like to generalize Paul Teetor's **A Better Hedge Ratio**, which uses prcomp() to determine a TLS ratio between two variables. I am hoping to extend this to multiple variables, but am having trouble finding accurate results. Any suggestions?


title: How to construct the feature weight vector (or Decision Boundry) from a linear SVM classifier from scikit?

text: I use the following code to train an svm classifier:               clf = svm.SVC(kernel='linear')     clf.fit(train_mat, train_labels)      that fit the data and save the info in the `clf` object. Now I know how theoretically the `w` vector is constructed. It is a sum over all support vectors multiplied by their labels and the corresponding alpha values. Problem is, I cant seem to find this info in `clf`. I have a `_coef` attribute but I not sure if it is what im looking for. Documentation in scikit is simply lacking. I couldn't find the info there. Does anyone knows how to do that? I need the `w` vector to check the wights given to each of the features and see which one are the most valuable.


title: strange behavior for LIBSVM during cross-validation

text: I am working on a 2-class classification problem and using LIBSVM. I am facing this following problem with cross validation (CV). Any help is highly appreciated. When I perform a 10 fold CV with -v 10 option, I get an accuracy of 90%. But when I do a training and testing on the same data with 90-10 split, the accuracy drops to 50%. How can this be explained? Thanks in advance. A Research student.


title: How to do cross validation for repeated measurement (proc mixed)

text: I had run the mixed model for repeated mesurement using SAS proc mixed, I got the model(s) already for my dependent variable. However, proc mixed doesn't provide R^2. Is there any cross validation methods that can apply to my data set. Or just simply tell people that I used proc mixed and no R^2 (which i think that people with not trust my model, esp. when need to publish paper). I often see people use LOOC for their data and compute mean error but my data is repeated (4 repeated). Anyone have any idea with this? Thank you in advance ^^ Good day:) KK


title: Misclassification Rate

text: The misclassfication rate in decision trees is defined as $$1-\max_{j} p_{j}$$ Suppose we want to classify people into republicans and democrats. In the training data, $\hat{p}_1 = 0.7$ (probability of being republican) and $\hat{p}_2 = 0.3$ (probability of being a democrat). If we split on $\text{age} \leq 50$ and have this has the terminal node, how would you compute the misclassification rate? If there are 100 people in the training data, then 70 are republicans and 30 are democrats. Then if 40 are less than or equal to age 50 and 60 are greater, would you count the number of republicans and democrats among the 40 people and 60 people to get a new $\hat{p}_1$ and $\hat{p}_2$?


title: Permutation test based in the Wilcoxon test

text: I am testing the statistical significance of 60 strong correlated quantitative variables in 2 groups, with small sample size, non normality and presence of outliers. I have implemented a permutation test based in the t-test statistic for deal correctly the correlation structure of the variables; but I think this permutation test could be also based in the statistic of the nonparametric Wilcoxon test, due to the presence of the outliers. Is it correct?


title: How many zeros in an independent variable are too many for regression?

text: I use simple regressions. My independent variable is a count data (# of a drug used per year), and it has too many zeroes. Depending upon a dependent variable, for some simple regressions, I have 13 observations, and there are only 3 non-zero values out of 13. For other simple regressions, I have 150 observations and only 12 are non-zeros out of 150. Is it appropriate to use simple regressions for these data? Is there a rule about how many zeros should be in independent variable to be able to model the relationship? My dependent variable is continuous (weight).


title: How interactions between explanatory variables can be found using binary recursive partitioning?

text: I would like to investigate interactions between my explanatory variables prior to building a statsitical model. Apparently it is possible to do it in R using a regression tree (library: tree). This method is based on binary recursive partitioning. The algorithm predicts a response variable by partitioning the data into subgroups based on the predictor variables. I do not understand, however, how this algorithm investigates interactions between predictors. I would be very grateful for an easy to understand explanation for non-expert.


title: Are there some risks connected with using k-NN with k=1 by F=0.77?

text: I am not very expirienced in the field of machine learning. For now I model some binary classifier using k-NN. I experimented a bit with different values of k. In some cases the best one was with `k=1`. I know that there is no majority voting, and my decision is made on single neighbour, what seems risky. F-measure for such case is `0.77` the second best option suggest for same data `k=25` with `F=0.63`               F=2*recall*precision / (recall + precision)


title: Approximating original formula from a dataset

text: _Disclaimer: I don't have sufficient statistics' background or vocabulary for this question, but I did my best to do a preliminary search both within the forum and the broader interwebs_ Is there a way for one to determine a formula or its simplified approximation from a set of results? My sole expectation of accuracy from such a "reverse-engineered" formula is the ability to reproduce the original set - or a set close enough to it - _not_ the ability to further expand the set or introduce intermediate values.


title: Who invented profile maximum likelihood estimation?

text: Could anyone give me some information on who invented profile maximum likelihood estimation or who first use profile maximum likelihood estimation and the short history of profile maximum likelihood estimation? I would like to know about any paper on this.


title: Laplace errors and ridge regression

text: Thinking about a normal linear regression, penalized by the L1 norm. That is the lasso, is there any literature on median ridge regression?. That is the residuals are calculated from an L1 norm instead of a squared norm. Then the penalization is by the square norm. My intuition tells me this will bring about variable selection too. Any literature on this?????


title: Statistical test for two distributions where only 5-number summary is known

text: I have two distributions where only the 5-number summary (minimum, 1st quartile, median, 3rd quartile, maximum) and sample size are known. Cntrary to the question here, not all data points are available. Is there any non-parametric statistical test which allows me to check whether the underlying distributions of the two are different? Thanks!


title: How can I regress household income on three factors

text: I have to regress family income (`faminc`; in dollars) onto **husband's educational attainment** (`he`; in years), **wife's educational attainment** (`we`; in years), and **number of children less than 6 years old in household** (`kl6`) using Stata. (the file only contains data of 4 above factors) I use OLS to estimate a model in the form: $$faminc = b_1 + b_2 * he + b_3 * we + b_4 * kl6 + \epsilon $$                     Source |       SS       df       MS              Number of obs =     430     -------------+------------------------------           F(  3,   426) =   28.77            Model |  1.4002e+11     3  4.6673e+10           Prob > F      =  0.0000         Residual |  6.9100e+11   426  1.6221e+09           R-squared     =  0.1685     -------------+------------------------------           Adj R-squared =  0.1626            Total |  8.3102e+11   429  1.9371e+09           Root MSE      =   40275          ------------------------------------------------------------------------------           faminc |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------               he |   3185.882   795.4493     4.01   0.000     1622.388    4749.376               we |   4637.415   1059.177     4.38   0.000     2555.551    6719.279              kl6 |  -8372.704   4343.059    -1.93   0.055     -16909.2    163.7893            _cons |  -5998.224   11161.51    -0.54   0.591    -27936.72    15940.27      I have some questions: 1) The regression yields $b4<0$. Is this true in fact? I mean that if the family has more children, the less income they gain? 2) Is this model good enough? Should I use natural logarithm or add dummy to make it better?


title: Using repeated measures anova: Am I right?

text: We have just finished an experiment, where we had samples distributed into 6 groups and each group receiving a different treatment. Each group had 5 samples, with 20 replicate measurements taken per sample (they do not represent timepont difference, but mainly replicates). The data looks something like this ![Data Sample](http://i.stack.imgur.com/6KmYi.png) I am stacking the data and making R chew on it, so that the final input looks something like this... ![DataInputToR](http://i.stack.imgur.com/hHA3Z.png) Now, We want to see if there is difference in TotalActivity (the values in the trial variable in the first table) Between the groups represented in FinalID (6 groups). I factor the relevant fields after importing the CSV and perform the repeated measures ANOVA in R like so               Data <- within(Data, {     FinalID <- factor(FinalID)     AnimalID <- factor(AnimalID)     Trial <- factor(Trial)     })          aov.out = aov(TotalActivity~FinalID*Trial + Error(AnimalID), data=Data)     PostHoc <- with(Data,pairwise.t.test(TotalActivity,FinalID,p.adjust.method="bonferroni"))     print(summary(aov.out))     print(PostHoc)      Okay, now the confession is that statistics is not my forte, and am a bit confused how to set up my repeated measures. I googled it and found a couple of formulas like so               aov.out = aov(TotalActivity~FinalID*Trial + Error(AnimalID), data=Data)          aov.out = aov(TotalActivity~FinalID + Error(AnimalID/FinalID), data=Data)      I have no idea what's the right formula to use, or even if any f the above are right statistical technique... suggestions please?


title: Variable order and accounted variability in linear mixed-effects modeling

text: Suppose that, in a study of 15 subjects, the response variable (res) is modeled with two explanatory variables, one (level) is categorical with 5 levels and the other (response time: RT) is continuous. With lmer in the lme4 package of R, I have:               fm1 <- lmer(res ~ level * RT + (level-1 | subject), data=mydata)     anova(fm1)                       Df  Sum Sq Mean Sq  F value     level        4  3974.9   993.7   9.2181     RT           1  1953.5  1953.5  18.1209     level:RT     4  5191.4  1297.9  12.0393      If I change the order of the two variables, I get slightly different results for the main effects:               fm2 <- lmer(res ~ RT * level + (level-1 | subject), data=mydata)     anova(fm2)                       Df  Sum Sq Mean Sq  F value     RT           1  1671.8  1671.8  15.5077     level        4  4256.7  1064.2   9.8715     RT:level     4  5191.4  1297.9  12.0393      Does such a difference come from the sequential (instead of marginal) approach in lme4 in accounting for data variability? In this case, the variable order change does not lead to a big difference, but previously I've seen dramatic differences. What does such a big difference mean? Does it mean that the model needs more tuning until big difference disappears? My second question is that, if I want to know which variable among the two (RT and level) accounts for more data variability, what would be a reasonable approach? Based on the relative magnitude of Sum Sq (or Mean Sq) of the two variables? Any statistical testing method to compare variability among explanatory variables?


title: ANOVA with time series as dependent variable

text: I have a randomized experimental dataset with six treatments with each approx. N=60. The outcome is a time-series, namely deforestation in a land-use simulation game over 40 rounds. I managed to show that the impact of the treatments on the state of the land (i.e. number of cummulative cells deforested) is highly significant in a single year, but I have a hard time finding the right method to show that the impact on the ENTIRE TIME SERIES is significant. I'm afraid testing in a single year is overestimating significance, as I can choose any year, adding "researchers degrees of freedom". I could (and successfully did) test every year seperatly, but that seems to be a very unelegant solution. In more general terms: I have a single independent variable from one 1 to 6, and my dependent variables is time series for every observation. What I want to do is basically an ANOVA, but feeding it with a whole time series as dependent variable instead of single values for each observation. If possible, it would be cool if the method also allows for controling for other independent fixed factors, such as player age, occupation etc., and ideally for more than one time series as dependent variable, as I also have data for intensification, savings, cows sold and some other values for every year in the game. Any expert insights? My data is a SQL database with a single entry for every round of the time series for every subject, linked to the subject-properties via a unique ID => I can bring it into any shape needed for the analysis. My problem is not to shape it but to find the right test. I'm using mySQL & R.


title: Sampling with non-uniform costs

text: Suppose that I have a population, each represented by a bit $b_i$ for $i \in \\{1,\ldots, n\\}$. I would like to compute an estimate $\hat{B}$ of the parameter $B = \sum_{i=1}^nb_i$ so that with high probability, the error $|\hat{B}-B| \leq k$ for some fixed $k$. However, I have to pay a cost $c_i$ to sample bit $b_i$, and this cost may be different for each $i$. I want to find the minimum-cost sample that satisfies my accuracy constraint. Clearly, uniform sampling is not necessarily optimal. Has this been studied? Is there a known optimal solution specifying the probability $p_i$ that I should sample each bit $b_i$ to compute $\hat{B}$?


title: Specifying conditional probabilities in hybrid Bayesian networks

text: I am trying to get a deeper understanding of the various types of Bayesian networks. Most of the literature/lectures I've come across use discrete random variables exclusively and only mention continuous random variables in passing. It seems if you want to mix discrete and continuous variables in a hybrid network, then you have a few different scenarios to handle:               Child Is       Parent(s) Are     ----------------------------------------     Discrete       All Discrete     Discrete       All Continuous     Discrete       Hybrid (Discrete and Continuous)     Continuous     All Discrete     Continuous     All Continuous     Continuous     Hybrid (Discrete and Continuous)      I understand how to specify conditional probabilities (i.e. $p(N=n_i|P)$) for the cases where the parents (P) are all discrete. The conditional probabilities are essentially a lookup table that maps a condition to either a probability mass function (when child is discrete) or probability density function (when child is continuous). The other cases are numbing my brain. How does one go about specifying conditional probabilities of a given node when any of its parents are continuous? Thanks **Edit** Terminology has been (hopefully) corrected based on feedback.


title: Soft and Hard EM (Expectation Maximization)

text: What is the difference between soft and hard expectation maximization? EDIT: ok, i've found out this paper: http://ttic.uchicago.edu/~dmcallester/ttic101-07/lectures/em/em.pdf that explain quite well the situations


title: How to predict max value of a heavy tailed variable

text: I have a response variable which is right skewed. There are lot of low values, less than 10 and then extremely large outliers. I want to create a model which can predict the max values that the variable can attain. Which techniques can I use to predict the max value


title: Generative vs discriminant models

text: What are the differences between generative and discriminant models (in the context of Bayesian learning and inference)? and what it is concerned with prediction, decision theory or unsupervised learning?


title: Most famous statisticians

text: What are the most important statisticians, and what is it that made them famous?   (Reply just one scientist per answer please.)


title: Seeking a distribution, perhaps uncommon, consistent with two data points and expert constraints?

text: I am trying to state a prior distribution for a Bayesian meta-analysis. I have the following information about a random variable:   1. Two observations: 3.0, 3.6   2. a scientist who studies the variable has told me that $P(X<2)=P(X>8)=0$, and that values as high as 6 have nonzero probability. I have used the following approach to optimization (the mode of log-N = $e^{\mu-\sigma^2)}$:               prior <- function(parms, x, alpha) {       a <- abs(plnorm(x[1], parms[1], parms[2]) - (alpha/2))       b <- abs(plnorm(x[2], parms[1], parms[2]) - (1-alpha/2))       mode <- exp(parms[1] - parms[2]^2)       c <- abs(mode-3.3)       return(a + b + c)     }     v = nlm(prior,c(log(3.3),0.14),alpha=0.05,x=c(2.5,7.5))     x <- seq(1,10,0.1)     plot(x, dlnorm(x, v$estimate[1], v$estimate[2]))     abline(v=c(2.5,7.5), lty=2) #95%CI      ![alt text](http://i.stack.imgur.com/hvKna.png) In the figure, you can see the distribution that this returns, but I would like to find something more like the red lines I have drawn in. This provides the same shaped distribution using the lognormal, gamma, or the normal, and it results in a distribution with $P(X=5)<0.05$ and $P(X=6)< 0.01$, i.e.:                plnorm(c(5,6), v$estimate[1],v$estimate[2])      Can anyone suggest alternatives? I would prefer to stick with a single distribution rather than a mixture. Thanks!