title: L1 norm linear function estimation

text: I have a question related to linear function estimation. Suppose the potential line function model can be expressed as `ax+by+c=0`, where `a, b, c` are the unknown parameters that needed to be estimated. Now a set of data is given `(x1,y1), (x2, y2), ... (xn, yn)`, and we want to make use of these data to estimate the parameters. As outliers (points that are not suited for model estimation) may exist in this data set, one way to avoid outliers is to minimize `sum(abs(a*xi+b*yi+c)) (i=1,2, ..., n)`. However, I could not figure out how to minimize this function step by step, could anyone give me some clues? Thanks!


title: Multiple comparisons with the glht function in R gives strange results

text: As a follow-up to this question, I have the following data:                  Site Treatment Survival     1   BED        DN      1.0     2   BED        DN      1.0     3   BED        DN      1.0     4   BED        MB      1.0     5   BED        MB      1.0     6   BED        MB      0.9     7   BED    Forest      0.4     8   BED    Forest      0.5     9   BED    Forest      0.4     10  BRO        DN      0.9     11  BRO        DN      1.0     12  BRO        DN      1.0     13  BRO        MB      1.0     14  BRO        MB      1.0     15  BRO        MB      1.0     16  BRO    Forest      1.0     17  BRO    Forest      1.0     18  BRO    Forest      1.0     19  LAP        DN      0.8     20  LAP        DN      0.4     21  LAP        DN      0.6     22  LAP        MB      0.5     23  LAP        MB      1.0     24  LAP        MB      0.7     25  LAP    Forest      0.2     26  LAP    Forest      0.2     27  LAP    Forest      0.4      on which I ran a binomial glm :               > glm.out <- glm(Survival~Site*Treatment, data=surv,         family="binomial", weights=rep(10, nrow(surv)))     > anova(glm.out, test="Chisq")          Analysis of Deviance Table     Model: binomial, link: logit     Response: Survival     Terms added sequentially (first to last)                         Df Deviance Resid. Df Resid. Dev P(>|Chi|)         NULL                              26    138.254                   Site            2   63.098        24     75.155 1.988e-14 ***     Treatment       2   42.991        22     32.164 4.620e-10 ***     Site:Treatment  4   13.874        18     18.290  0.007707 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       All effects are significant, so now I want to do post-hoc comparisons. I have discovered the `glht` function in the multcomp package, which seems to do what I want. I read somewhere that with 2 independent variables, you should set `interaction_average=TRUE` to get a result equivalent to a TukeyHSD for a linear model.               > Treat.comp <- glht(glm.out, mcp(Treatment="Tukey", interaction_average=TRUE))     > summary(Treat.comp)      which gives me these strange results:                    Simultaneous Tests for General Linear Hypotheses     Multiple Comparisons of Means: Tukey Contrasts     Fit: glm(formula = Survival ~ Site * Treatment, family = "binomial",          data = surv.san, weights = rep(10, nrow(surv.san)))          Linear Hypotheses:                      Estimate Std. Error z value Pr(>|z|)     DN - MB == 0       -0.202   3316.127   0.000        1     Forest - MB == 0   -1.886   3316.127  -0.001        1     Forest - DN == 0   -1.684   3316.127  -0.001        1     (Adjusted p values reported -- single-step method)      Why are my P-values = 1 for all comparisons even though the Treatment effect was highly significant? If I try using `glht` with `interaction_average=FALSE` I get more reasonable results but with a warning message:               Warning message:     In mcp2matrix(model, linfct = linfct) :       covariate interactions found -- default contrast might be inappropriate      Can someone help me to understand what I am doing wrong? Thank-you!!!


title: Tutorials for feature engineering

text: As is known to all, feature engineering is extremely important to machine learning, however I found few materials associated with this area. I participated to several competitions in Kaggle and believe that good features may even be more important than a good classifier in some cases. Does anyone know any tutorials about feature engineering, or is this pure experience?


title: Extending logistic regression for outcomes in the range between 0 and 1

text: I have a regression problem where the outcomes are not strictly 0, 1 but rather in the range of all real numbers from 0 to 1 included $Y = [ 0, 0.12, 0.31, ..., 1 ]$. This problem has already been discussed in this thread, although my question is slightly different. I can't use linear regression for the same reasons that logistic regression is normally used. In linear regression A) very large IVs values will skew the predicted outcome to 1 and B) the outcome of linear regression is not bounded to the 0,1 limits. Looking at this logistic cost function from my textbook $$\text{Cost} = -y \log(h(x)) - (1 - y) \log(1-h(x))$$ I gather that the equation is designed to calculate a cost greater than 0 only when $y$ and $x$ do not have the same value 0 or 1. Would it be possible to use logistic regression by modifying the cost function to measure all hypothesis errors?


title: Appropriateness of K-S test and Kruskal-Wallis for assessing medical data set

text: I am analyzing some stats for a paper I am writing. I have statistics help available through my faculty but am not there for a couple of months so hoping you can help out. n = 463, 2 patient groups( <70 & >= 70 years). Using SPSS. I have assessed a whole bunch of stats using Pearson's $\chi^2$ test - They are categorical such as Smoking History (T/F), Smoking Current (T/F), BMI>30, Diabetes,etc. So far so good (Unless anyone can tell me if this is an inappropriate test?) I have more variables such as Pre-operative Creatine, ICU hours stay, etc., which are not normally distributed (I used to know how to test for normality 'properly' but I have just drawn Q-Q plots and seen if they are on the line and they aren't. My intention was to use student's t but my rusty stats knowledge tells me I can't now because they aren't normally distributed. In many aspects I am using the statistical methods of a friend that wrote a similar paper using a similar dataset, and for his continuous variables he used **Kolmogorov-Smirnov** but I'm not sure if he just used this for normality. He then says that he used **Kruskal-Wallis** to assess these variables but I can't work out why this is appropriate. Would someone mind explaining to me if these tests are appropriate for comparing these samples?


title: What do I do when values of AIC are low and approximately equal?

text: Chris Chatfield, whose many quality books and papers I enjoyed reading, in (1) gives the following advice: > For example, the choice between ARIMA time-series models with low and > approximately equal values of the AIC should probably be made, not on which > happens to give the minimum AIC, but on which gives the best forecasts of > the most recent year's data. What is the rationale for such advice? If it is sound, why does forecast::auto.arima and other forecasting routines do not follow it? Yet to be implemented? It has already been discussed here that to look for models that just happened to give the minimum AIC is probably not a good idea. Why is the option to have $n\ge1$ ARIMA models with low but approximately equal (e.g. within 1 or 2 values of the minimal AIC) is not a default in much of the time series forecasting software? (1) Chatfield, C. (1991). Avoiding statistical pitfalls. _Statistical Science, 6_ (3), 240–252. Available online, URL: http://valjhun.fmf.uni- lj.si/~mihael/ul/vs/pdfseminar/chatfield.pdf.


title: Calculating standard error of a coefficient that is calculated from other estimated coefficient

text: I'm working on the Gompertz growth model to fit weight at age: $$ m(a)=m_{\infty}e^{-\gamma exp(-g{1}a)}$$ Where $m_{\infty}$ and $g_{1}$ are coefficients to be estimated. To deal with lack of normality I transformed my data using Box&Cox, and to deal with heterogeneity I used weights. After applying _Box &Cox_ I fit a model with the **R** function `nls()` and I used this model's residuals as weights $(w_{i}=1/StudentizedResiduals_{i})$ in a new fitting. I need to estimate the coefficient $\mu=m_{\infty}e^{-1}$ and its standard error. **Question 1:** As I've got $m_{\infty}$ and its standard error estimated with `nls()`, so Can I calculate the standard error of $\mu$ like this:$$StdErr(\mu)=StdErr(m_{\infty})e^{-1}$$ Or should I estimate this estandard error in a different way? **Question 2:** As a result of the _Box &Cox_ application my estimated coefficients are in a different scale and it's difficult to interpret them. From this question I understand that I can't transform-back the coefficients, but any suggestion how make easier the interpretation of 'transformed' coefficients? Thanks in advance


title: Understanding confidence bounds in logistic and multinomial regression

text: I used both logistic (binomial distribution) and multinomial regression on a range of different datasets. For simplicity, I only have one explanatory variable in the linear predictor. I also calculated the confidence bounds for the predicted probabilities, following this procedure. What I observed, is that for the large datasets my confidence intervals were much narrower than for smaller datasets, which I thought made sense. However, when I look at a diagnostic plot such as expected vs. observed counts, even for the large datasets there are some visible systematic deviations of the data from the fit which apparently are not reflected by the confidence intervals, which seem to underestimate such variations. What I mean is, I would expect the CIs to be wider to reflect these deviations. Any idea on why this would be the case? Or am I missing something in my interpretation / calculation of confidence intervals?


title: Annotating graphs in R

text: This is more of a "how to use R" question than an actual hardcore statistics question, but I think the concentration of R masters here makes this a good forum for it. I'm refreshing a time series graphing package that currently uses gnuplot. The first step is getting somewhere close to the current graphs, and then I hope to be able to add more statistical analysis from R later. Currently, the gnuplot-generated graphs look like this: ![alt text](http://i.stack.imgur.com/EkuiN.png) The corresponding graph from R looks like this: ![alt text](http://i.stack.imgur.com/mqWEr.png) I'm reasonably happy with this. What's missing is the legend stuff from below the gnuplot graph. I can add a legend saying "In" and "Out", but how to best present the average, maximal and minimal values? I could add horizontal lines to the graph and print the values inside the graph area, but the placement might interfere with the actual graph.


title: Different time scales in time-dependent cox proportional hazards

text: I'm doing survival analysis with time-dependent variables and I have two types of data: birds followed in a biweekly basis and others followed randomly. I think that analysing both types of data together is a major problem. Any suggestions in how to deal with it?


title: First passage time distribution in a irreducible transient discrete-time Markov chain (DTMC)

text: In a Markov chain, a state $j$ is transient if $f_{jj}<1$ ($f_{jj}$ is probability of ever visiting state $j$ starting from state $j$ ). Suppose, I have an irreducible transient DTMC (means all states are transient). Now, I want to prove that for any $i,j$ in $S$ ($S$ is DTMC state space), $f_{ij}$ (i.e probability of ever reaching state $j$ starting from state $i$) is less than 1. It is clear that $f_{ii}<1$ and $f_{jj}<1$. But, how to prove that $f_{ij}<1$ for any $i,j$. Thanks Prasenjit


title: Logistic Regression: Interpreting Continuous Variables

text: I had a couple of questions about interpreting odds ratios for continuous variables in logistic regression. I feel like these are basic questions about logistic regression (and probably about regression in general), and although I'm slightly ashamed that I don't know the answers, I'm gonna swallow my pride and ask them so I know them in the future! Here's my situation...I'm looking at a sample of adjudicated youth who, as part of their probation, were enrolled in a job/life skills training program. I wanted to see the extent to which the age at which they were released from the program predicted employment six months post-release from the program. (Also, keep in mind that there are other predictors in the model, but I've excluded them because they are not statistically significant and I want to keep this as clear as possible.) Predictor: Age of release from training program (Mean age = 17.4, SD=1.2, Range 14.3-20.5) Outcome: Employed or not (Employed=1, Not Employed=0) Result: Odds ratio 3.01 (p<.005) (I've excluded goodness of fit stats, etc. because I'm seeking answers about the interpretation of the odds ratio only; I feel comfortable w/ the evaluation of model fit, CI's, etc.) Putting it into words: As age increases by one year, the odds of being employed six months post-discharge increase by three units. Questions: 1) When I say, "As age increases by one year..." what is the starting point for age? Does age start at zero? For example, "As age increases from 0 [i.e., the lowest age if you were to place this model on a graph]...” Does age start at the the lowest age among the range of ages in the sample? For example, “As age increases from 14.3...” OR Does age start at the mean age of the sample? For example, “As age increases from 17.4...”, 2) Would centering help me interpret this result OR is that only effective in interpreting the y-int? If it would help, I was thinking of doing either mean centering or subtracting the lowest age in the range from all the other ages in the sample. Any suggestions? 3) Finally, is it appropriate to say that compared to a 14-year-old youth, a 17-year-old youth is nine times more likely to be employed? I ask because I know that logistic regression assumes a sigmoidal relationship, and I’m curious as to whether this 3 unit increase in odds remains consistent at any point along the regression line. Thanks so much! Aaron


title: PCA and degrees of freedom

text: I had thought that a principal components analysis gives us information about the number of degrees of freedom (i.e. independent variables) underlying a system which we may only be able to observe through derived (and possibly correlated) variables. But a simple experiment has made me doubt this: I generated 1,000 sets of 10 drawings from independent standard normal distributions ($\mathcal{N}(0,1)$) and calculated the PCA for that data matrix. The correlation matrix was very close to diagonal, as I'd hope. The eigenvalues, and corresponding proportion of variances were:                1.14   12%      1.11   11%      1.07   11%      1.04   10%      1.02   10%      0.98   10%      0.96   10%      0.93    9%      0.85    9%      0.81    8%      This chimed with my expectations: this system has ten independent variables, and the eigenvalues make it clear that all ten principal components carry a significant contribution to the total variance of the system. Next, for each set of observations, I derived a new set of variables, $\{y_1,...,y_{10}\}$, by taking linear combinations in the form of cumulative averages of the original variables $\{x_1,...,x_{10}\}$: $$ y_n = \frac{1}{n}\sum_{i=1}^{n}x_i$$ I've lost no information by doing this: all the $\{x\}$ can be exactly recovered from $\{y\}$. However, while the $x$ were independent, the $y$ are clearly correlated (and the correlation matrix confirms this). I calculated the PCA for $\{y\}$ and obtained this set of eigenvalues:                1.98   68%      0.57   19%      0.18    6%      0.09    3%      0.04    1%      0.03    1%      0.01    0%      0.01    0%      0.01    0%      0.00    0%      Which suggests that only five or six principal components are necessary to capture virtually all the variance in the system. Now, to some extent, I see that this makes sense: the $y$ are quite strongly correlated, so e.g. knowing $\{y_1,...,y_9\}$ places tight constraints on $y_{10}$. But the system is still driven by ten underlying independent variables, and all I've done is take linear combinations of them to form my new variables. Shouldn't the PCA be able to see through this? Clearly my intuition is a little awry. Can anyone help me to see this more clearly?


title: Interpretation of point-biserial correlation

text: How can I interpret point-biserial correlation? If the results give me a positive and significant correlation, how should I interpret it? Should I say that the variable category that I coded 1 is positively correlated with the outcome variable?


title: Does similarity of coefficients and p-values regardless of whether dependent variable is transformed  suggest  untransformed model is reliable?

text: Related to my earlier question, I need to perform regression on a skewed dependent variable (n = 500). Since the residuals weren't normally distributed, I was able to transform the DV non-linearly in a way that it now approaches normality. Residuals are normal when using this transformed variable as a dependent variable. For the two models, The p-values for the various predictors are very much alike, and the relative sizes of the coefficients are very similar as well.   * To what extent are those two facts are indicators (or not) of the reliability of the coefficients and p-values obtained in the first model (using raw data)?


title: How do I set up an unbalanced repeated measures analysis in R?

text: ### Context: I am analysing some impact assessment data (measuring invertebrate richness in response to pollution), but they are unbalanced - there are not data for every site at sampling occasion, and there were more datapoints recorded _after_ the impact than _before_ the impact. I am a new user to R, and have gathered through reading on this site and others that the standard anova packages `aov()` and `ezAnova()` can't deal with unbalanced designs. I assume I should instead be using a package like `lme4`. However, I am not sure how to structure my data, or program the analysis. One of the problems is that I'm not sure how to incorporate sampling dates as the repeated measures aspect of my design. My data has 5 columns Site code, Date, BeforeAfter, ControlImpact, Richness. ### Questions:   * How should I set up my data for conducting repeated measures analysis with unbalanced data in R?   * Should I use `lme4` or some other package?


title: Whether to use structural equation modelling to analyse observational studies in psychology

text: I've noticed this issue coming up a lot in statistical consulting settings and i was keen to get your thoughts. ### Context I often speak to research students that have conducted a study approximately as follows:   * Observational study   * Sample size might be 100, 200, 300, etc.   * Multiple psychological scales have been measured (e.g., perhaps anxiety, depression, personality, attitudes, other clinical scales, perhaps intelligence, etc.) The researchers have read the relevant literature and have some thoughts about possible causal processes. Often there will be some general conceptualisation of variables into antecedents, process variables, and outcome variables. They have also often heard that structural equation modelling is more appropriate for testing overall models of the relationships between the set of variables that they are studying. ### Question   * Under what conditions do you think structural equation modelling is an appropriate technique for analysing such studies?   * If you would not recommend structural equation modelling, what alternative techniques would you recommend?   * What advice would you give to researchers considering using structural equation modelling in such cases?


title: Weight variables for predictive model

text: I received a question today that I wasn't exactly sure how to answer. I have built a predictive model using a fairly basic logistic regression that works pretty well and fits our business needs. Recently, we purchased a CRM tool that allows us to build "probability" scores, but only allows the end users to give integer weights to various factors. Said differently, one can arbitrarily assign a weight of 10 points to one factor and -5 points to another with the sum of all weights representing the "probability" for a given entity in our database. What I am looking to do is translate my model to this new format such that the resulting score equals the calculated probability from my logistic model. This is not out of desire, but business needs. Admittedly I am not sure how to use the calculated coefficients and "adjust" them to these requirements. What is the best approach, if any? General thoughts on how to assign statistically valid integer weights to business criteria given these constraints? Any thoughts or insight will be very much appreciated.


title: Factor dependent correlation

text: If I take a set of measurements and test correlation of variable $A$ vs variable $B$ and get a significant correlation, that makes sense to me. But what if further analysis reveals that of those factors, there is only a significant positive correlation within one group, and that group is over- represented. Is the global correlation still valid, or is it, upon more detailed inspection a sample-bias effect? Here is some graphs to explain: The global correlation ![global correlation](http://i.stack.imgur.com/9tFJx.png) The group separated correlations ![grouping separated correlation](http://i.stack.imgur.com/cTPwA.png)


title: Superiority of LASSO over forward selection/backward elimination in terms of the cross validation prediction error of the model

text: I obtained three reduced models from a original full model using   * forward selection   * backward elimination   * L1 penalization technique (LASSO) For the models obtained using forward selection/backward elimination, I obtained the cross validated estimate of prediction error using `CVlm` in package `DAAG` available in `R`. For the model selected via LASSO, I used `cv.glm`. The prediction error for LASSO was less than than the ones obtained for the others. So the model obtained via LASSO seems to be better in terms of its predictive capacity and variability. Is this a general phenomenon that always occurs or is it problem specific? What is the theoretical reasoning for this if this is a general phenomenon?


title: Statistical testing: Multiple classifiers, 1 domain. Would rANOVA be appropriate?

text: When comparing the performance of two classifiers over a single domain, in the context of a classification problem in machine learning, it is common to use a paired t-test, using the 10 average results from 10x10-fold cross-validation as measurements, where the folds at each iteration are the same for the two classifiers. An obvious yet incorrect generalization of this t-test to the case where there are multiple classifiers would be to take all pairs, or just the interesting pairs if we only want to compare 1 classifier to all the others, and to work from there. The problem, however, is that the results of these multiple t-tests are not independent from one another. I've heard of the ANOVA as a generalization of the t-test to situations with more than 2 groups to compare. I thought it was the solution to this problem, but then I read more and stumbled upon the repeated measures ANOVA, which seems to be even closer to what I'm looking for. Could anyone confirm (or invalidate) that a repeated measures ANOVA is indeed what should be used in this situation? What is the difference between a regular ANOVA and the rANOVA in this case?


title: Objections to randomization

text: In Clinical trials - a methodologic perspective, Steven Piantadosi writes (ch.13, p. 334): > In Chapter 2, I noted the objections to randomization by Abel and Koch > (1997) and Urbach (1993), and indicated the worth of studying their concerns > and likely errors. They reject randomization as a > >   1. means to validate certain statistical tests, >   2. basis for causal inference, >   3. facilitation of masking, and >   4. method to balance comparisons groups. > According to me, (1)-(4) are benefits of randomization. So, why do Abel, Koch, and Urbach reject randomization on the basis of those arguments?


title: Inspecting assumption of homoscedasticity

text: Using a Fligner test to infer about the respect of the assumption of homoscedasticity is not very smart given that the Fligner test tests to the null that there is no difference of variance between the groups. This will wrongly favors small sample size. As it has been said by @Michael Mayer here. How can we further investigate if the assumption of homoscedasticity is respected? Is it worth plotting the model's residuals versus the fitted values? Below lines are R coded:               m = aov(myFormula, myData)     plot(residuals(m), m$fit)      I don't have much experience in statistics and it seems rather hard for me to decide from this plot whether the assumption of homoscedasticity is respected. What else can I do?


Does regression analysis measure cause and effect? If yes, then how? If no, then what is done? Please describe with an example.


title: Should I always order the data vector in ascending order before using "summary" function in R?

text: Should I always **order** the data vector in ascending order before using `summary` function in R? Or should I just apply `summary` function directly on raw data vector?


title: Is it necessary to check ergodicity in estimation of autocorrelation function?

text: Given a sample path from a process supposed to be stationary, I saw the sample autocorrelation function of the sample path is used to estimate the autocorrelation function of the process. But this requires that the stationary process is ergodic. So   * is ergodicity checked on the sample path before estimation of the autocorrelation function? If yes, how is ergodicity checked?   * is the estimation done assuming ergodicity without checking it, and after the estimation, testing/checking if the estimation fits the data? If yes, how is the testing done? A side question: how is stationarity is tested for estimation of autocorrelation function? Thanks and regards!


title: Prior of multivariate Polya distribution?

text: Anyone knows a prior (preferably conjugate) to the multivariate Polya distribution? I need it for Gibbs sampling. So if anyone has another idea, I am interested.


title: Log-likelihood ratio in document summarization

text: I initially asked this on stack overflow and was referred to this site, so here goes: I am implementing some unsupervised methods of content-selection/extraction based document summarization and I'm confused about what my textbook calls the "log-likelihood ratio". The book briefly describes it as such: "The LLR for a word, generally called lambda(w), is the ratio between the probability of observing w in both the input and in the background corpus assuming equal probabilities in both corpora, and the probability of observing w in both assuming different probabilities for w in the input and the background corpus." Breaking that down, we have the numerator: "the probability of observing w in both the input and in the background corpus assuming equal probabilities in both corpora" - How do I calculate what probability to use here? and the denominator: "the probability of observing w in both assuming different probabilities for w in the input and the background corpus". - is this as simple as the probability of the word occurring in the input times the probability of the word occurring in the corpus? ex: (count(word,input) / total words in input) * (count(word,corpus) / total words in corpus) I've been looking over a paper my book references, Accurate Methods for the Statistics of Surprise and Coincidence (Dunning 1993), but I'm finding it difficult to relate to the problem of calculating LLR values for individual words in extraction based summarization. Any clarification here would be really appreciated.


title: Feature / attribute selection for k-means or other clustering

text: It seems to me that in literature it is assumed that one knows which features / attributes to choose to characterize an item in clustering. If I have a database with items which have many attributes, how do I know which attributes to choose for a good clustering? Is there any guideline or literature which deals with this problem?


title: Proportional hazards assumption

text: The proportional hazards assumption basically says that the hazard rate does not vary with time. That is, $\text{HR}(t) \equiv \text{HR}$. When can we assume this? What if the hazard ratios at various times are: $2.4, 2.36, 2.27$ and $2.03$? Can we make the proportional hazards assumption? Also we have $$ \log[h(t|\textbf{x})] = \log[h_{0}(t)] + \beta_{1}x_1 + \dots + \beta_{p}x_{p}$$ Why do we need to estimate $h_{0}(t)$? If we have $h(t|\textbf{x})$, why can't we just put all the values of the predictors to zero to get $h_{0}(t)$? **Edit.** I want a means to assess whether the PH assumption is true.


title: Finding median survival time from survival function

text: Is the best way to find the median survival time from a survival plot just to draw a horizontal line from $p = 0.5$ to the curve and project down to the x-axis?


title: 5-way interaction

text: I am running a mixed factors ANOVA for a brain imaging study on language processing. The design includes four within-subject factors:   1. complexity: simple/complex;   2. agreement: correct/number violation/gender violation;    3. hemisphere: left/right;   4. region: anterior/posterior;  and one between-subject factor:   1. proficiency: low/high proficiency.  Factors 3 and 4 (hemisphere, region) only matter if they interact with the other factors. This is because the effects of complexity and agreement are predicted to emerge in specific areas of the scalp. The results reveal a significant 5-way interaction (complexity by agreement by hemisphere by region by group), and I find it almost impossible to figure out what is driving the interaction by just looking at the means. Am I justified in looking at the two proficiency groups separately? Most similar studies do this, but I'm unclear as to how to do this. For example, once I decide to examine proficiency learners separately, I find a 3-way interaction between agreement, hemisphere, and region in the high-proficiency learners (meaning that there is a brain response for agreement violations, which is captured in the left anterior portion of the scalp), but there is no such interaction for low-proficiency learners (p = .13) Is it licit to report this difference even if the original omnibus ANOVA (with group as a factor) showed that there was no agreement by hemisphere by region by proficiency interaction? In other words, once I decide to look at the two proficiency groups, can I run the analysis as if I had never compared the two groups or am I constrained by the original ANOVA as to what follow-ups I can do?


title: Choosing the right size of an out of sample data

text: I am beginner in forecasting, especially forecasting with R and I am really willing to improve my knowledge. Recently, I started practicing electricity consumption time series forecasting. The first barrier I faced is the choice of out of sample data for assessing the forecast accuracy of the forecast model i will be using (regression with ARIMA errors). I have data for 147 months and I want to forecast the next 24 months, for the period June 2013 to January 2015. Furthermore, I have read in @RobHyndman's online text book Hyndman, R.J. and Athanasopoulos, G. (2013),   _Forecasting: principles and practice_ ,   (accessed 28 May 2013), section 2.5   under ' **Training and test sets** ', that: > size of the test set is typically about 20% of the total sample, although > this value depends on how long the sample is and how far ahead you want to > forecast If I divide the dataset with 20% of it being the out-of-sample data, the forecast model applied to the in-sample data is not quite accurate, since I guess it fails to capture the recent trend, (which began in the middle of the last year), of decreased electricity consumption due to a significantly raised electricity tariff. What do I do about this? Can you possibly give me instruction on what would be considered appropriate size of the out of sample data. I also tried with 7 months of out-of-sample data, but I am afraid that there might be an overfitting issue. Is that right?


title: How to check for a relationship between water bacterial levels

text: I have a set of data: http://pastebin.com/KHLKD8XB which is based on TVC (total viable counts - i.e., numbers of bacteria) from water going into a machine (BM), water taken from the machine (IM) and water taken from a device which is washed in that machine (Scope). Whilst the counts are performed on a 200ml sample, results are reported as cfu/100ml, hence you appear to find 1/2 a bacteria in some samples.               library(ggplot2)     testscope<-read.csv2("http://pastebin.com/raw.php?i=KHLKD8XB")     testscope$date<-as.Date(testscope$date)      ggplot(data=testscope,aes(x=date, y=tvc, colour=class)) +          geom_step()  + facet_grid(class ~ ., scales="free")  +          labs(title="Machine TVC versus Scope TVC") + ylab("cfu/100ml")      Plotted using ggplot2: ![enter image description here](http://i.stack.imgur.com/pksjL.png) **My question is: How can I check whether TVC readings from the Scope are related to TVC readings from the machine ?** By eye, I can't discern a pattern between the machine level (IM) and the Scope levels. Note that although I have the BM data, that is not something we are particularly interested in at present - it is IM versus Scope that we need to understand. EDIT: Note that the dataset uses natural log, not base 10. My mistake.


title: R Code for Yeo-Johnson transformation

text: I have writen code for a Box-Cox transformation (see below). But now I want to do a Yeo-Johnson transformation because `datc$plot` contains zeros. I tried, but I didn't find a solution.               lambda.fm1 <- boxcox(datc$plot ~ datc$cond.evlot*datc$cond.dl*datc$version),                           family="yjPower")     lambda.max <- lambda.fm1$x[which.max(lambda.fm1$y)]     require(car)     datc$plott <- bcPower(datc$plot, lambda = lambda.max, jacobian.adjusted = FALSE)


title: Methods for comparing proportions in partially overlapping samples

text: We conducted 2 similar surveys in parallel to a group of health care facilities. One was sent to managers, one to staff. Some questions were identical, others asked about different areas of the facility. We received responses from groups of facilities that do not completely overlap. So our response rates were the following (approximately).   * survey a - to staff member at facilities: 70%   * survey b - to manager at facilities: 81%   * both surveys: 62% We are interested in reporting proportions, difference in proportions with (95% CI), and statistical test of difference between proportions for identical question on both surveys (answered by 2 non-overlapping populations) Martin Bland has a good paper on different approaches...does anyone have SAS code for similar analysis? http://www-users.york.ac.uk/~mb55/overlap.pdf


title: Find k in k-means, but only between two options

text: There are many algorithms for finding $k$ in the $k$-means algorithm that depend on finding the "kink" in a graph of the objective function. What if I know that the data have either $2$ or $3$ clusters? How do I find which $k$ is "better"? Of course, a higher $k$ will lower the objective function, but this in not necessarily the "correct" k. Any ideas?


title: Any public available evaluation frameworks for recommender systems?

text: I would like to know if there are any evaluation frameworks of recommender systems which are capable of evaluating rating prediction and topN recommendation (Precision and recall etc.). Maybe I need to find them in recommender frameworks? If so are they easy to plug in? As a matter of fact I am working on dynamic recommender systems which may need the items for training are earlier than item for testing/recommendations, which is also the real application scenario I think. In my point of view, as long as the data records have timestamps, the evaluation work should use training and testing items with time order. But many of the two sets of items are randomly selected. Am I right?


title: Survey, does this make sense?

text: I'm trying to learn a bit of statistics, and I thought I'd ask you guys about conducting surveys. Let's say I conduct a survey like the following https://docs.google.com/forms/d/1YrqcWQNl1CmbINaVZuL18ZseqBSWemdHjy_P4ES637A/viewform I am surveying a lot of cyclists and I want to know how they are motivated. Different questions describes different types of motivation (seven categories of motivation with 4 questions in each category). Does it make sense to do any statistics on this? For example I would think testing whether one type of motivation is significantly more prevalent than another. I guess I am comparing the different groups (of motivation). I'm not sure if this is a Likert scale. If you guys have any suggestions on which types of tests I should do, it would be really helpful. Thanks in advance.


title: Generalization error for classification with a nonconvex loss function

text: I've been working my way through Vapnik's 1998 _Statistical Learning Theory_ book and one thing that I'm still unsure of is if his risk bounds hold for nonconvex loss functions -- i.e., when we can't be sure we've minimized the empirical risk. Does anyone know of any references that discuss this?


title: independent or dependent

text: I have situation in which I query different sites on the genome (biologically independent) and for each site I get a measurement. I have a number of measurements that are close together on the genome and I call that a region. I then want to compare the measurements of a region between two biological specimens using a Mann Whitney U test (small number of measurements per region). I want to stress that I am not repeatedly measuring the same thing, but that all measurements concern separate loci on the genome. We do know however that loci that are in close proximity to each other show strongly correlated measurements. A reviewer pointed out to me that these measurements are not independent because they are correlated and that the subsequent test is therefore not valid. In my opinion however there is no reason to assume that the value of one locus influences the value of the adjacent loci. Therefore the measurements are independent, regardless of their correlation. On the internet I foundvarious opinions. If measurements influence each other, they are dependent, that much is clear. But the term correlation also pops up a lot. Indeed there will be correlation in the case of measurements influencing each other, but if that is not the case and there is still correlation present, would you agree that the measurements are independent regardless of their correlation? The answer from Peter Flom already provides some clues: Dependent vs. independent samples Thanks for any insights. Kind regards,


title: Difference in AIC and BIC values between sem and lavaan packages in R

text: I ran the same SEM model in _sem_ and _lavaan_. I got the same parameters and - generally - very close test values, with the exception of AIC and BIC which were immensely different between the two packages. The following is the resulting AIC and BIC from _sem_ :               AIC =  2913.849     BIC =  -1777.617      The following is the resulting AIC and BIC from _lavaan_ :               Akaike (AIC)                               37780.878     Bayesian (BIC)                             38178.999      Why is there such a huge difference in these values? Are they calculated differently in each package? edit: Here are pieces of code on how I got such values. In the case of package _sem_ , I used the following line:               options(fit.indices = c("GFI", "AGFI", "RMSEA", "NFI", "NNFI", "CFI", "RNI", "IFI", "SRMR", "AIC", "AICc", "BIC", "CAIC"))     fit = sem(mymodel, cov(mydata), nrow(mydata), data = mydata)     summary(fit)      As for the package _lavaan_ , the following line was used instead:               fit = sem(mymodel, data = mydata, estimator = "ML")     summary(fit,fit.measures=TRUE)      I am using _R_ version 3.0.3 (2014-03-06) "Warm Puppy" (x64 version), _sem_ package version 3.1-3 and _lavaan_ package version 0.5-16 . Thank you, Deuterium


title: Should the control variables in the first and second stage of 2SLS be different?

text: I first noticed this when I was running some models is Stata and the variables that I had used as controls in the second stage were showing up in the first stage (I think this is just Stata's default approach). A lot of these variables are in no way related to the dependent variable in the first stage. Similarly, there are a lot of variables that I would like to include in the first stage but don't think they should be included in the second stage. Would it be better to just do OLS on the first stage with a set of controls (including the instrument) relevant to that dependent variable (the endogenous variable I want to instrument), save the fitted values and use them in another OLS regression with another set of controls relevant to the final dependent variable? Rather than use the Stata command approach?


title: Imputation with mice: recode variables before or after imputation?

text: I am using `mice` in `R`, a chained equations (sequential regression) algorithm, to impute a series of polytomous variables (e.g. scales from 1 to 5 or 1 to 3). Ultimately, my interest is in estimating statistics on simple proportions of categories on these variables (e.g. 1 and 2 versus other). There are two options for doing this recoding, i.e. before and after imputation. Recoding after imputation would mean treating variables as polytomous in `mice` using a multinomial or ordinal regression model. On some categories (e.g. top or bottom, 1 or 5) the sample data is sparse, however. So I am not sure about the variance of imputing these categories by polytomous models. However, the model would be most adequate to the original data. Recoding before imputation might solve this problem, because sparsely populated categories are summarized to indicator variables. In addition sequential regression would be limitted to less categories, in most cases thus 2 (e.g., 1 and 2 versus other). A logistic model with far less parameters would be used which may increase accuracy of imputations. But I am not sure whether recoding before imputation might change the data to an extend that imputations are not correct anymore. So my questions are should I recode before or after imputation? Does this make any difference, is one or the other preferable?


title: direction of effect size

text: Should my effect size be positive or negative when the control group gains more weight than the experimental group (when it is hypothesised that the experimental group would lose more weight i.e., the difference between groups is due to more weight gain in the control rather than weight loss int the experimental group)?


title: Validity of an iteration for Goodness of fit to a specific application

text: I am not a statistics expert and would like to check the validity of a test I hav used on a survey results. **Survey Results** I asked people (n=264) to characterize a panoramic road they chose in 24 different categories. Each of them could choose as many categories he/she desired. Than I counted the number of times each category was mentioned for two subsets of my sample: experts (n=37) and public (n=227). I received a distribution of the magnitude (or frequencies) for each category in each group. **The Question** I would like to check whether the distribution of choices in m two subsets is the same. If it is different, I would like to know which categories are the cause fr this variance. **The test** I chose to use Goodness of Fit (Chi-Square) Using R, while x=(public counts) and p=(experts counts) with rescale.p = TRUE. I set my significance level to 5%. Than I wrote a while loop, that works untill test's P-Value > 5%. The loop omits variables in each iteration. Omitting variables is based upon a delta vector (delta<-abs(public_counts - expert_counts). For each iteration the variables that has the highest delta are omitted. That is untill H0 is accepted. I would like to know if this extension of the chi.square function using the while loop is a valid methd to answer my question. In particular, is the delta vector that is based on absolute counts is a correct specification. I would be happy to provide raw data and code for my multi.chi function if necessary.


title: How to compare the variance from published summary statistics with own data?

text: I have compiled a very small set of summary data from the literature, and I wish to compare the variances between aspects of the literature-based data, and to some of my own data. The summary data includes the mean, standard deviation and sample size. In earlier tests, I compared the variances of one continuous dependent variable among 2 age classes and 2 years. I used the Fligner-Killeen test since I have extreme values in the data and I'm not sure it was normal (can't remember now!). I followed up this broad test with pairwise multiple comparisons using F-tests in R `var.test(x~age)` What would be a good way to compare the variances of the literature-sourced data? I've searched through the help files in R, and came up with this method, which I believe generates a random set of numbers with the specified sample size, mean and standard deviation, and then I used those datasets to conduct the F-test:               herring_year1<-rnorm(10,mean=10.5,sd=0.51)     herring_year2<-rnorm(15,mean=10.9,sd=0.43)     var.test(herring_year1,herring_year2)      Would this be a good approach? If not, can you suggest what might be? If so, how can I then compare these variances to my own data set? Should I essentially use the summary data from my own set in the same manner? Or generate the random data for the summary data from the literature and stick it in a file to compare to my raw data? Also, do I _need_ a broad test initially, or can I go straight to the pairwise comparisons?


title: What can we say about the likelihood function, besides using it in maximum likelihood estimation?

text: I found that often in literature that likelihood values are often used to compare different estimation method for the same model. And I got the impression that is the only way likelihood values are used. However, I wonder what else we can say about the likelihood function. For example, can we compare two totally different models' likelihood functions? What's the minimum of the likelihood (obviously not zero)? What's a good guess about the maximum without going through into iterative method of estimation? (Can we build up inequality equations and claim boundaries for the likelihood function?)


title: Restricted mean residual life: Survival analysis: Discrete

text: I am hoping someone can help with this - perhaps trivial question of interpreting the formula below. It is the mean residual lifetime (restricted to the window of data) for the case of discrete time. ### Restricted Mean Residual Life > Remaining time until event or upper bound--whichever comes first $$\mu_{[\gamma]}(t|\mathbf{x})=E\left(\min(T, t+\gamma+1)-t|T\ge t, \mathbf{x}\right)=\sum_{j=t}^{t+\gamma}\prod_{k=t}^{j}\left(1-h(k|\mathbf{x})\right)$$ My question is how do you interpret the sigma and product? $t$ is the current time (e.g. "5") and $\gamma$ is the upper bound from my understanding. How do you compute this?


title: How to model probability of a web click using bayesian modelling with constraints?

text: I'm trying to model clicks in a Google AdWords campaigns. This is a function of impressions, the probability of a click (CTR) on that impression, the cost per click (CPC), and the budget of the advertiser. It's essentially $$\text{clicks} = \text{impressions} \times \text{CTR}$$ subject to $$\text{impressions} \times \text{CPC} \times \text{CTR} \leq \text{Budget}$$ **How should this constraint be incorporated into my model?** Happy to provide more info if it'd be helpful.


title: How to understand the label-bias problem in HMM?

text: How can I understand the label-bias problem in Hidden Markov Models? And why is CRF able to solve this problem?


title: Could prior elicitation by actually drawing the density of the prior be sensible? Has it been done/discussed?

text: Somebody mentioned, I don't remember who, that "there are many ways to specify the prior, you could even draw it!". It is clear to me that it is possible to actually draw the density of the prior using pen and paper and transform it into a histogram prior for use in an analysis. What I'm wondering is: **Has this approach been used/discussed in the literature? If not, are there reasons why this is generally not a good approach?**


title: Confusion related to convexity of a function

text: I was reading this paper http://www.umiacs.umd.edu/~jags/pdfs/KSforNLP.pdf where they say that $trace(K \pi^{T} L \pi)$ is a convex function of $\pi$ where $\pi$ is a permutation matrix and K and L are kernel matrices. I am not sure how this function is convex. Any clarifications?


title: Fixed effects in a linear probability model (LPM)

text: I have a dataset with about 6000 obs. in a cross sectional times series structure (countries and years). My dependent/response variable is binary and my independent/explanatory variables are a mix of binary and continuous. Now, I worry that the explanatory variable I'm interested in might not be independent of unobserved, country-specific, time-invariant factors, so I'd like to remove these with fixed effects (and probably do the same with year- specific, but country-invariant factors). Since my response variable is binary, I have been using conditional logit models grouped on country, but I find it undesirable for this project that they discard all information on countries without variation on the response variable (there are a lot of those in my sample). So, I'd like to supplement these with linear probability models (LPM's) with standard fixed effects. However - and here comes my problem - I figure that, besides the usual problems with LPM's, they must have an additional issue in this case since the binary response variable would imply an unnatural restriction on the unobserved effects. I'm just not sure how big of a problem this restriction is? And if it depends on whether or not the rest of the model violates the 0-1 bound of the response variable? Hope you guys can help me out! Best wishes Bertel


title: Which test to use when comparing variances for data that is not necessarily normally distributed

text: I am a biologist looking for guidance regarding which statistical test to use when comparing variances in data. I am comparing two theoretical distributions to the experimental data. The data in the distributions may not be normally distributed and there may not be an equal number data points. An example of my data is shown below: ![enter image description here](http://i.stack.imgur.com/zDW2d.png) So what I have done so for is to normalise the data on similar y-axis scales. I did this by taking the mean of each of the three distributions and plotting their difference from the mean. The distribution of the data is shown below: ![enter image description here](http://i.stack.imgur.com/bTuoy.png) The y-axis is the difference from the mean and the x-axis the type of data. Using the normalised distribution I have been using F-Tests in R (with command `var.test()`) to statistically analyse the variance in the data. But the F-Tests assumes a normal distribution. Is the F-Test the correct test to use in this case? If not can anyone suggest a more appropriate test?


title: How many observations do you need within each level of a random factor to fit a random effect?

text: I'm trying to analyse some data from a set of bird surveys. My response variable is "bird abundance", which is the number of birds counted over a five-minute period. These five-minute counts were conducted at ~200 sites. Counts were repeated three times at each site, although there are ~20 sites where only two counts were completed. I want to model bird abundance as a function of site-level attributes (habitat quality etc) as well as Count-level attributes (weather conditions at the time of the count etc). So, I have ~200 sites and ~600 individual counts, but only 2-3 counts per site. My question is: given that I only have 2-3 counts per site, can I include site as a random factor to account for the nonindependence of counts within sites? (note I can drop the sites that have only two counts if necessary). I've read conflicting information about the number of observations that you need within each level of the random factor. Ben Bolker’s paper on GLMMs in Trends in Ecology and Evolution says "5-6 random effect levels per random effect and 10-20 samples per treatment level or experimental unit", but then I've also read stuff that suggests using mixed models for repeated measures designs that take a pre-treatment, post-treatement, and follow up sample - ie only three observations within each level of the random effect. Thanks for the help! Jay


title: Marginalization from randomly-spaced grids of the joint distribution

text: I have the vectors of co-ordinates of the grid points $(x, y)$ and the joint probability values, $z$. Could you recommend a way to do marginalization for the below joint distribution? (In my problem, I have random sampled grid points from the joint distribution. Of course, if the points were lined up on straight lines then this is not so difficult). ![irregular grids](http://i.stack.imgur.com/GeO8m.png)


title: Ordinal data and zero inflated Poisson regression

text: Suppose we have ordinal data that can assume the following values: $0$, $1$, $2$, $3$ and $4$. Also suppose that there are a lot of data points that fall into category $1$. Can we use zero inflated Poisson regression to model this?


title: How to test whether the difference in difference between means is significantly different? (with example)

text: Basically I have   * Four groups (2-by-2) with sample means $\overline{x_{1,A}}$, $\overline{x_{2,A}}$, $\overline{x_{1,B}}$ and $\overline{x_{2,B}}$, standard deviations $s_{1,A}$, $s_{2,A}$, $s_{1,B}$ and $s_{2,B}$ with sample sizes $n_A$ and $n_B$ (i.e. $n_{1,A}=n_{2,A}$ and $n_{1,B}=n_{2,B}$). I calculate the differences between the two means as $d_{1-2, A} = \overline{x_{1,A}}-\overline{x_{2,A}}$ and $d_{1-2, B} = \overline{x_{1,B}}-\overline{x_{2,B}}$. My objective is to test whether $d_{1-2, A}$ is statistically different from $d_{1-2, B}$, corresponding to a two sided test. My suggestion would be: $\frac{d_{1-2, A}-d_{1-2, B}}{ [s_{1-2, A}^2/n_A] + [s_{1-2, B}^2/n_B] }$. The problem is that I don't know how to calculate the standard deviation of this difference, i.e. $s_{1-2, A}$ and $s_{1-2, B}$. So my questions are basically:   1. Do you agree with my procedure?   2. How do I calculate $s_{1-2, A}$ and $s_{1-2, B}$? p.s. I tried $s_{1-2, A} = \sqrt{s_{1,A}^2 + s_{2,A}^2}$ but it didn't work. * * * **Here is an example with numbers** $\overline{x_{1,A}} = 2.310$; $s_{1,A}=0.865$; $n_A = 16753$ $\overline{x_{1,B}} = 2.403$; $s_{1,B}=0.897$; $n_B = 5378$ $\overline{x_{2,A}} = 2.993$; $s_{2,A}=1.002$; $n_A = 16753$ $\overline{x_{2,B}} = 3.108$; $s_{2,B}=1.042$; $n_B = 5378$ The objective is to replicate a p-value for significane between $d_{1-2, A}$ and $d_{1-2, B}$ of **$0.032$** (two sided would be my guess).


title: regression model selection, back transformation of a logged response

text: EDITED: I want to select the a regression model (geographically weighted regression vs. OLS) based on a pseudo-R2 comparing the predicted values and the actual ones. I am aiming to regress the predicted values onto the actual values. My models have a logged response variable. My questions are as follows: 1) Is it necessary to back transform the response before calculating the pseudo R2? 2) If it is the case, is it necessary to consider the variance only in the regression output or also in the logged input variable? My R code looks like this:               lm( exp( predict(regressionmodel)+var(predict(regressionmodel))*0.5 ) ~ exp(data$loggedinput), data=data )      I appreciate your response very much. Thanks in advance!


title: R and EViews differences in AR(1) estimates

text: The main **problem** is: I cannot obtain similar parameter estimates with EViews and R. For reasons I do not know myself, I need to estimate parameters for certain data using EViews. This is done by picking the NLS (nonlinear least squares) option and using the following formula: `indep_var c dep_var ar(1)` EViews claims that they estimate linear AR(1) processes such as: $$ Y_t = \alpha + \beta X_t + u_t $$ where $u_t$ errors are defined as: $$ u_t = \rho \cdot u_{t-1} + \varepsilon $$ by using an equivalent equation (with some algebraic substitutions): $$ Y_t = (1 - \rho) \alpha + \rho Y_{t - 1} + \beta X_t - \rho \beta X_{t - 1} + \varepsilon_t $$ Furthermore, this thread over at the EViews forums suggests that their NLS estimations are generated by the Marquardt algorithm. Now, the go-to R function to estimate AR(1) processes is `arima`. However, there are two problems:   1. the estimates are maximum likelihood estimates;    2. the intercept estimate is not actually the intercept estimate (according to R.H. Shumway & D.S. Stoffer). Therefore, I turned to the `nlsLM` function from the minpack.lm package. This function uses the Marquardt algorithm to achieve nonlinear least squares estimates, which should yield the same results as the EViews implementation (or very similar ones, at least). Now the code. I have a data frame (`data`) with an independent variable and a dependent variable such as the one generated by the following code:               data <- data.frame(independent = abs(rnorm(48)), dependent = abs(rnorm(48)))      To estimate parameters in the equation EViews claims to estimate (3rd one on this post), I use the following commands:               library(minpack.lm)     result <-     nlsLM(dependentB ~ ((1 - theta1) * theta2) + (theta1 * dependentA) +                         (theta3 * independentB) - (theta1 * theta3 * independentA),     data = list(dependentB = data$dependent[2:48], dependentA = data$dependent[1:47],        independentB = data$independent[2:48], independentA = data$independent[1:47]),     start = list(theta1 = -10, theta2 = -10, theta3 = -10)     )      Unfortunately, the estimates output by `nlsLM` are not close to the ones output by EViews. Do you have any idea of what might be causing this? Or maybe my code is wrong? Finally, I would like to say that I personally am an R user - that is exactly why I'm trying to do this in R instead of EViews. I would also love to provide you the data I'm working with but it's impossible since it's confidential data.


title: Patterns in binary data sorted with rank indices

text: I have a sequence of binary data, each row of which has seemingly random bits and associated with an integer value (" _rank_ "). When we have a binary datum, we always can calculate the _rank_ related to it but not the data _from_ it's rank value. For an example, with 32-bit length and sorted by its rank:               rank     data     69538    10010110 01010100 00110101 11100000      57643    10010001 01110001 00111010 11001101     3478     01010001 11001001 10100110 11100011     278      11101000 11100000 00100100 01001000     ...      My question is how do I find a pattern in the data to reduce the rank? The final goal is to get data which have rank indices smaller than certain threshold value. I thought about a neural network but would there be any better approach? Thank you.


title: Explaining generalized method of moments to an non-statistician

text: How do I explain Generalized Methods of moments and how it is used to a non statistician? So far I am going with it is something we use to estimate conditions such as averages and variation based on samples we have collected. How do I explain the part where you estimate the parameter vector by minimize variance


title: Why only the mean value is used in (K-means) clustering method?

text: In clustering methods such as K-means, the euclidean distance is the metric to use. As a result, we only calculate the mean values within each cluster. And then adjustments are made on the elements based on their distance to each mean value. I was wondering why the Gaussian function is not used as the metric? Instead of using `xi -mean(X)`, we can use `exp(- (xi - mean(X)).^2/std(X).^2)`. Thus not only the similarity among the clusters are measured (mean), but the similarity within the cluster is also considered (std). Is this also equivalent to the Gaussian mixture model? It is beyond my question here but I think mean-shift may arise the same question above.


title: How do I regularize Python scipy's curve_fit when performing linear regression?

text: I have recently become proficient at using `Python/scipy` `curve_fit` to perform linear regression. However, with higher order polynomials my data is sometimes overfits. **How can I add regularization to reduce over-fitting?**


title: Confusion over interpreting regression coefficients

text: My lecturer's regression:               Variables:     *lprice*, the log of house prices     *y81*, which is a dummy for the year 1981 against 1978      *ldist*, log of distance from incinerator to houses      **These are comments made by my lecturer to help give context:** >   * From the lecture we know that in 1978 people did not know about the > incinerator site. >   * By 1981 people did know about the site. >   * If building the incinerator reduces the value of houses closer to the > site, the coefficient on y81*ln(dist), delta_1, will be positive. >   * If beta_1 is positive it indicates that the incinerator is being located > in an area where the houses are of lower value. >   * (If politicians/planners assume that poorer people have less lobbying > power, this would be evidence of strategic decision-making on their part.) > ![lprice](http://i.stack.imgur.com/phgkm.png) **After regression was done my lecturers comment was this:** >   * The coefficient on y81ldist is positive, but small and not significantly > different from zero at the 10% level (or the 5% level or the 1% level). >   * In contrast the coefficient on ldist is positive and significant at the > 0.1% level according to a two-tailed test. >   * Based on this preliminary analysis, I conclude that the incinerator has > no effect on house values. > **Now i changed _lprice_ to _price_ to see if it made a difference. Here is the regression I got:** ![enter image description here](http://i.stack.imgur.com/SZEuv.png) Here the coefficient of _y81ldist_ is quite large, and significant at 6.7% level, compared to a 55% in the first regression. Im curious as to why my lecturer concluded that building the incinerator had no effect on prices? Am I overlooking something? From my regression, it shows that building the incinerator reduced prices of houses that are nearer to it by an average of ~15,400 dollars. (Also, does the 19302 coefficient on ldist mean that before the incinerator was built, prices of houses that were further away from the proposed site of building were already 19302 dollars more expensive that the ones closer?) Are my interpretations right? Any insight would be appreciated.


title: Contrast to test significant interaction - why not include main effect?

text: Given a data set with a continuous outcome variable $Y$ and two nominal IV $X_1$ and $X_2$ which have 4 (A,B,C,D) and 3 (E,F,G) levels respectively, we can assign the first level of each as a reference level and create dummy or indicator variables: $X_1: X_1B$ , $X_1C$ , $X_1D $ $X_2: X_2F$ , $X_2G$ An ordinary linear regression of $Y$ on these variables can be written as (including the interaction terms): $Y = \beta_0 + \beta_1X_1B + \beta_2X_1C + \beta_3X_1D + \beta_4X_2F + \beta_5X_2G + \beta_6X_1B*X_2F + \beta_7X_1B*X_2G + \beta_8X_1C*X_2F + \beta_9X_1C*X_2G + \beta_{10}X_1D*X_2F + \beta_{11}X_1D*X_2G + e $ where $e $~$N(0,\sigma{^2})$, the error is normally distributed with constant variance. If we wished to test the significance of the interactions $X_1$ x $X_2$ to determine if, say, the effect of $X_1$ on $Y$ differs at different levels of $X_2$, one method would be to use a contrast by constructing the appropriate matrix $L$ to multiply against the regression beta vector with the null hypothesis $L\beta=0$. See for example here. My understanding is that we would be (with the example above) using null hypothesis involving just the interaction terms: $Ho:$ $\beta_{6}=0$ and $\beta_{7}=0$ and $\beta_{8}=0$ and $\beta_{9}=0$ and $\beta_{10}=0$ and $\beta_{11}=0$ but, why wouldn't testing $\beta_{1}=0$ and $\beta_{2}=0$ be included because these are jointly interpreted as the effect of $X_1$ on $Y$ when $X_2= $level E (the reference level for $X_2$)? **ADD After ahfoss response:** If the interpretation of $\beta_6$ through $\beta_{11}$ = 0 is the following: The effect of $X_1$ = B on Y is not different when $X_2$ = F The effect of $X_1$ = B on Y is not different when $X_2$ = G The effect of $X_1$ = C on Y is not different when $X_2$ = F The effect of $X_1$ = C on Y is not different when $X_2$ = G The effect of $X_1$ = D on Y is not different when $X_2$ = F The effect of $X_1$ = D on Y is not different when $X_2$ = G Isnt it possible that "The effect of $X_1$ = A on Y IS different when $X_2$ = E" and that this would suggest an interaction? I am missing the obvious I'm sure but it has not quite clicked for me. I am missing how the 6 statements above for $\beta_6$ through $\beta_{11}$ = 0 is equivalent to your equation 1.


title: Difference between loadings and correlations between observed variables and factor saved scores in factor analysis

text: I thought that the loadings in factor analysis were the correlations between the observed variables and the latent factors. However, when I do factor analysis in R using the psych package, this does not seem to be the case:                   library(psych)         set.seed(1)         X <- matrix(rnorm(200), ncol=10)         fa1 <- fa(X, nfactors=3, rotate="none", scores=TRUE)              cor(X, fa1$scores)  #correlations between original variables and factor scores                        MR2         MR1         MR3          [1,]  0.465509161  0.87299813  0.03241641          [2,] -0.010609644 -0.32714571  0.64968725          [3,] -0.219685860  0.47331827 -0.39132195          [4,] -0.815516983  0.22669390  0.42273446          [5,] -0.075178935 -0.40431701 -0.69661843          [6,] -0.204917832  0.07472006  0.05508017          [7,]  0.240675941  0.13027263  0.23238220          [8,]  0.756677687 -0.05621205  0.23746738          [9,]  0.004384459  0.12095273  0.55100943         [10,]  0.640507568 -0.67810600  0.18597947              fa1$loadings[1:10, 1:3]                        MR2         MR1         MR3          [1,]  0.433925641  0.82218385  0.02717957          [2,] -0.009889808 -0.30810366  0.54473104          [3,] -0.204780777  0.44576800 -0.32810435          [4,] -0.760186392  0.21349881  0.35444221          [5,] -0.070078250 -0.38078308 -0.58408054          [6,] -0.191014719  0.07037085  0.04618204          [7,]  0.224346738  0.12268990  0.19484113          [8,]  0.705339180 -0.05294013  0.19910480          [9,]  0.004086985  0.11391248  0.46199451         [10,]  0.597050885 -0.63863574  0.15593470              cor(fa1$scores)  # Check that factor scores are uncorrelated                   MR2          MR1           MR3         MR2  1.000000e+00 4.266996e-16 -1.299606e-16         MR1  4.266996e-16 1.000000e+00  1.961151e-16         MR3 -1.299606e-16 1.961151e-16  1.000000e+00      The loadings and correlations are similar, but I expected them to be the same. I tried looking at the source code for `fa` but had trouble understanding it. Could someone please tell me how the loadings differ from the correlations? Update: For each factor, the correlations with the observed variables are constant multiples of the loadings:               cor(X, fa1$scores)/fa1$loadings[1:10, 1:3]                MR2      MR1      MR3      [1,] 1.072786 1.061804 1.192675      [2,] 1.072786 1.061804 1.192675      [3,] 1.072786 1.061804 1.192675      [4,] 1.072786 1.061804 1.192675      [5,] 1.072786 1.061804 1.192675      [6,] 1.072786 1.061804 1.192675      [7,] 1.072786 1.061804 1.192675      [8,] 1.072786 1.061804 1.192675      [9,] 1.072786 1.061804 1.192675     [10,] 1.072786 1.061804 1.192675


title: Dumb Bayes theorem question

text: What's the difference between "probability of H given E" as calculated by Bayes theorem and the "positive predictive value" of E?


title: Pattern recognition with discrete data (preferably in R)

text: I am trying to get a grip on pattern recognition with discrete data. Currently I started to work through the classics of Duda and Hart. Is there something you could recommend particularly for applied use, preferably with R. Is there even a useR Springer book I haven't found yet or any other good tutorial? Something that helps with algorithm selection?


title: How to condition a marked point processes?

text: **Conditioning a point process is the concept being investigated in this post.** The marked point processes are of interest. Figure **A** demonstrates a phenomenon being investigated for modelling. We recognize a feature such can be modeled by a point as the center and an angle (Figures **B** and **C** ). Therefore we construct a point process based on the central points (Figure **C** ). Then we develop our model by adding a mark as representative of angle (Figures **C** and **D** ). Figure **(i)** demonstrates a field being studied in which only some samples in some regions are available. Some regions therefore are unsampled and unknown. As described above we generate a marked point process to get a model describing the phenomenon in entire space. Figure **(ii)** is the result. As it can be seen the left hand side of (ii) does not match the existing samples. There are indeed generated stochastically. figure **(iii)** demonstrates the final goal in which a full map of the phenomenon has been modeled in entire space conditioned to the existing samples at the left hand side.   **Question: How can we do it?**   **Note:** The conditional modelling referred to conditional simulations those are common in geostatistics cannot be applied for this case at least directly. For example, a grade estimation made by Kriging or IDW is already conditioned to the existing samples. This is because an interpolation (weighted averaging) is being used for estimation which makes sense for grade. However for the case being studied here, there is no such understanding for variation of angle (as one of variables of interest). Furthermore, the output of the estimations/simulations in geostatistics etc is often a full grid which means you have a point very close to the sampling point. Indeed the location itself is not considered as variable while for PP stochastic generated locations are the base. The marks also are generated by means of stochastic models where again cannot be adjusted regarding to nearby samples. **Updates:** of course an answer is of interest if you don't mind! ![enter image description here](http://i.stack.imgur.com/rd2BY.jpg)


title: Estimate beta binomial distribution

text: I have a dataset as follows: ni-xi 76 48 38 19 20 14 26 29 39 45 38 36 32 34 26 21 23 12 24 18 14 61 35 21 32 20 89 30 34 35 xi 0 0 16 26 1 8 23 19 24 16 16 19 0 0 0 6 8 3 5 9 5 0 34 17 13 0 1 11 13 22 I want to model this dataset as a beta-binomial distribution: xi ~ Binomial(p, ni) where p ~ Beta(α, β). I use VGAM package to de the fitting. fit=vglm(cbind(ni-xi, xi) ~ 1, betabinomial.ab, trace = TRUE) coef(fit) VGLM linear loop 1 : loglikelihood = -95.9867 VGLM linear loop 2 : loglikelihood = -95.6633 VGLM linear loop 3 : loglikelihood = -95.6483 VGLM linear loop 4 : loglikelihood = -95.6477 VGLM linear loop 5 : loglikelihood = -95.6477 VGLM linear loop 6 : loglikelihood = -95.6477 However, I get the results as follows: (Intercept):1 (Intercept):2 0.9006006 -0.3794863 Can anyone tell me why I get the negative parameters? How can i solve this problem.


title: Training Hidden Markov Models for multiple input observations

text: I'm working with Hidden Markov Models and I have a dataset composed by independent phrases, where each word is an observation. Hence, the best way to adjust my parameters (via Baum-Welch algorithm) is considering each phrase per time and not all phrases concatenated. I would like to know if there is an algorithm that do the training in this way. If not, what are the strategy to avoid transitions created by the concatenation (last word of to first word).


title: Which test(s) for analyzing yes-no answers for control vs treated conditions

text: To keep it simple, I have sets of data in which I have a control and a treated sample. The data I obtain is "yes" or "no" (0 or 1) for a particular behaviour. I have 150 data points for each condition (150 for control, 150 treated), obtained across 3 experimental runs, which I intend to pool. They are not paired in any way. Which test would you recommend? I am currently thinking that binary logistic regression is the best choice, but I'm not sure. I am a biologist and not a mathematician, so please keep that in mind :)


title: Tukey's Hinges: Grouping Data

text: I am having trouble figuring out how to group data using the inter-quartile ranges calculated with a box and whiskers plot as well as looking at Tukey's Hinges. I understand the IQR and Tukey's hinges are not the same thing and that there are different interpretations of Tukey's Hinges. Basically, I did a calculation using SPSS and the output provided me the weighted quartiles and below that, Tukey's hinges. My question is about whether or not you can use the hinges or the weighted quartile values within the data to group the data. For example. Let's say, this is your set of data (this is not the data set I'm working, but just a simple example to illustrate my question): 1 2 2.5 2.5 2.5 3 3 4 5 5 6 7 7.5 7.5 8 9 9 10 Let's then say your IQR and Tukey Hinges are 25% = 2.5, 50% = 5, 75% = 7.5 (This is just random, this may not be the actual case for the data but I am using these values just to explain my question). Now let's say you want to divide the data into 4 groups using the IQR and/or Tukey's Hinges. Which group do the values which are the hinge or "division point" between the groups go? When I was reading up on Tukey's Hinges, Tukey stated that the hinge is a point of division in the data but that is still vague when trying to group data. So, would a groups be like this: Group A 1 2 2.5 2.5 2.5 Group B 3 3 4 5 5 Group C 6 7 7.5 7.5 Group D 8 9 9 10 Or can you exclude the "hinge" values and group the data like this? Group A 1 2 Group B 3 3 4 Group C 6 7 Group D 8 9 9 10 I am having a tough time finding good explanations and research papers I could use to back up why I choose one methodology over another. This project is for an internship, none of the data I've written out in this e-mail is being used in my analysis. I am simply confused about dividing groups based on Tukey's Hinges. I would appreciate your thoughts. Thanks!


title: Significance test in linear regression with only one predictor

text: As far as I know, I have two options for tests in linear regression. The $F$-test for the model (if it explains more variance than it has error variance) and the $t$-test (to see if the slope is not zero). With more than one predictor, I can see why both tests are there. But in my case, I only have one predictor. In my opinion, the $F$-test and the $t$-test do the same in this case because if the slope is zero, it is exactly the model it is compared to in the $F$-test. Where is the flaw in my logic?


title: How do hidden units in a Boltzmann Machine enforce this parity constraint?

text: In Ackley and Hinton's paper "A Learning Algorithm for Boltzmann Machines", they write that > A hidden unit would be needed, for example, if the environment demanded that > the states of three visible units should have even parity-a regularity that > cannot be enforced by pairwise interactions alone. Could someone explain how a hidden unit enforces this parity constraint? I'm having a hard time seeing what the structure and weights of the network would be. (In general, I see intuitively why hidden units add power, but I don't have a rigorous understanding.)


title: The prestige magician paradox

text: You probably know the trick in the movie The Prestige : **[MOVIE SPOILER]** A magician has found an impressive magic trick: he goes into a machine, close the door, and then disappears and reappears in the other side of the room. But the machine is not perfect : instead of just teleporting him, it duplicates him. The magician stays where he is, and a copy is created at the other side of the room. Then, the magician in the machine falls discreetly in a water tank under the floor and is drowned. **Edit:** The probability of the new copy of the magician being drowned is 1/2 (in other words, the new copy has 1/2 chances of being drowned, and 1/2 chances of popping into the room). Also, the water tank never fails and the chances are 1 that the magician dropping in the tank dies. So the magician doesn't really like doing this trick, because "you never know where you are going to be, on the other side of the room or drowned". Now, the paradox is the following : **Imagine the magician does the trick 100 times. What are his chances of surviving ?** Edit, additional question: **What are the chances of the magician of keeping his physical brain and not having a new one ?** * * * Quick analysis: One one hand, there is one magician alive, and 100 drowned magicians, so his chances are 1 out of 100. On the other hand, each time he does the trick, he has 1/2 chances of staying alive, so his chances are $(1/2)^{100}=1/(2^{100})$ of staying alive. What is the right response and why ?


title: Matching training and testing datasets

text: This is more a general question on machine learning in general, but I have two datasets ~75k rows of patients and 100 columns (call it `training` and `testing`)--each column is either a numeric predictor corresponding to how many times they've visited a specific physician, or a factor (two examples are Gender and 10-year age bins). There are also two outcome arrays `training.obs` and `testing.obs` that correspond to how many days each patient has spent in the hospital. Let's say I wanted to train a linear model `model <\- lm(training.obs ~ ., data=training)`, and cross-validate it with the testing set by               testing.pred <- predict(model, testing)      and take the RMSE of `testing.pred` and `testing.obs`. The two datasets are greatly different, in terms of factor ratios (i.e. a larger proportion of males in the training set), and distribution of numeric predictors. Would it make sense to match these two datasets by sampling `training`, or will I lose model accuracy because I'm excluding rows? How would I go about doing this?


title: How do I determine how well a dataset approximates a distribution?

text: Quite simple, I have some probability distribution p(x), how can I measure whether one empirical density (set of delta masses) is a better approximation than another. I know that KL-divergence is a well accepted measure between two continuous densities, but it's not clear how to apply that to a set of samples.


title: Finding groups of similar people

text: After asking 16 questions (yes or no) to 75 people I have a table of their answers coded like `00110011110101010` ('0'=no and '1'=yes). Now I would like to find groups of people who answered in the most similar way (ideally to find all possible groups with similitudes from 16/16 to 12/16). How to do this in SPSS?


title: Changes in R-squared

text: I was reading online, and I found that If the Variance of X increases then the value of R-squared increases If the Variance of the residuals increases then the value of R-squared decreases Can someone explain the intuition behind these two statements. I am generally confused about the impact of variance of X or the residuals on the explanatory power of the regression.


title: Weird looking binned residual plot

text: ![Binned Residual Plot](http://i.stack.imgur.com/jgduc.png) My binned residual plot is quite strange looking, the 95% confidence lines are so very jagged, with points between. I have colored this "inside" of the 95% confidence interval because it is really hard otherwise to tell which points are inside and which are outside. This plot was made in R. This is from a generalized estimating equation-logistic regression, with many correlated clusters and a AR-1 correlation structure. I am not sure if this is a reasonable way for this plot to look or if I have done something wrong. Has anyone else had a plot that looks like this? Thanks in advance, Jesse


title: Bayesian: Sampling from Truncated Distributions

text: When would the rejection sampling method be preferred to the inverse CDF method for sampling truncated random variables? And when would the inverse CDF method be preferred to the rejection sampling method for sampling truncated random variables?


title: computing KL divergence: M projections for arbitrary distributions

text: **Background** I have a generative model for a process that can be described as follows: $$ y = t(x, w) + e $$ where $x$ and $y$ observations of a set of random variables which are related by a non-linear transformation function $t$, parameterised by the unknown parameters to be estimated given by $w$. $e$ is the normally distributed error term with 0 mean and a diagonal covariance matrix given by $\sigma^{-1} I$ where $I$ is the identity matrix and $\sigma$ is the global noise precision. So, assuming independence along each pixel of the image, I have the following likelihood term as a product over the individual pixels $i$: $$ P(y|x, w, \sigma) = \prod_{i} (\frac{\sigma}{2\pi})^{\frac{1}{2}} \exp^{-\frac{1}{2}e_i \sigma e_i} $$ Each $w_i$ parameter is a 3-dimensional parameter (for each spatial dimension) to be estimated. **KL-Divergence and M-Projections:** I am using Expectation propagation (EP) to estimate the posterior distribution and EP has a M-projection step which projects a distribution onto a simpler approximating distribution which in my case is a multivariate normal distribution over my parameters $w$. The way EP works is by exploiting the factored form of the likelihood term i.e. it starts with an approximating distribution $Q$ (multivariate normal over $w$) and then replaces each ith factor by the exact term and projecting this distribution on the current estimate of the Q. For example, in this problem I replace one ith term from my likelihood expression to form the distribution to project as: $P(y|x_i, w_i, \sigma) q_{j \neq i}(w_j)$ where the first term is the exact factor and $q_{j \neq i}(w_j)$ is the approximated distribution with the influence of the ith term removed. Let us call this distribution $U$ According to the literature, I need to find the moments $E_U[w]$. So, I need to match the first and second order moments i.e. compute these expectations for my parameters using the distribution $U$. However, I am completely lost as to how to do this. Can someone give me some suggestion on how I should proceed?


title: What to do when heteroskedasticity is expected?

text: I've got a dataset where I'm attempting to predict when an individual will develop a particular disease based on a set of biomarkers. I'm able to find a pretty good fitting model, but it has a high degree of heteroskedasticity. However, this heteroskedasticity is expected--it makes sense that the model will have smaller residuals as the individual nears diagnosis. I began thinking about various "fixes," but wasn't sure if I should fix it. Any thoughts on this?


title: Difference between t-test and ANOVA in linear regression

text: I wonder what differences are between t-test and ANOVA in linear regression?   1. Is a t-test to test whether any one of the slopes and intercept has mean zero, while ANOVA to test whether all slopes have mean zero? Is this the only difference between them?   2. In simple linear regression i.e. where there is only one predictor variable, there is only one slope to estimate. So are t-test and ANOVA equivalent, and if yes, how, given that they are using different statistics (t-test is using t-statistic and ANOVA is using F-statistic)? Thanks and regards!


title: How many data points to fit an ex-Gaussian distribution?

text: I want to fit an ex-Gaussian distribution (which has three parameters: $\mu$, $\sigma$, and $\nu$) to reaction time data. **How many data points are needed as a minimum to reasonably fit this distribution?** I am aware that an answer to this question depends on many factors (e.g., how much data actually deviate from the theoretical distribution) and that "more is better" in general. Anyway, are there any rules-of-thumb, or publications which could point me to the needed minimum of data points?


title: Why do people use the term "weight of evidence" and how does it differ from "pointwise mutual information"?

text: Here, "weight of evidence" (WOE) is a common term in the published scientific and policy-making literature, most often seen in the context of risk assessment, defined by: $$w(e : h) = \log\frac{p(e|h)}{p(e|\overline{h})}$$ where $e$ is evidence, $h$ is hypothesis. Now, I want to know what is the main difference with PMI (pointwise mutual information) $$pmi(e,h)=\log\frac{p(e,h)}{p(e)*p(h)}$$


title: Unbiased estimator for the smaller of two random variables

text: Suppose $X \sim \mathcal{N}(\mu_x, \sigma^2_x)$ and $Y \sim \mathcal{N}(\mu_y, \sigma^2_y)$ I am interested in $z = \min(\mu_x, \mu_y)$. Is there an unbiased estimator for $z$? The simple estimator of $\min(\bar{x}, \bar{y})$ where $\bar{x}$ and $\bar{y}$ are sample means of $X$ and $Y$, for example, is biased (though consistent). It tends to undershoot $z$. I can't think of an unbiased estimator for $z$. Does one exist? Thanks for any help.


title: What to make of explanatories in time series?

text: Having worked mostly with cross sectional data so far and very very recently browsing, scanning stumbling through a bunch of introductory time series literature I wonder what which role explanatory variables are playing in time series analysis. I would like to _explain a trend_ instead of de-trending. Most of what I read as an introduction assumes that the series is stemming from some stochastic process. I read about AR(p) and MA processes as well as ARIMA modelling. Wanting to deal with more information than only autoregressive processes I found VAR / VECM and ran some examples, but still I wonder if there is some case that is related closer to what explanatories do in cross sections. The motivation behind this is that decomposition of my series shows that the trend is the major contributor while remainder and seasonal effect hardly play a role. I would like to explain this trend. Can / should I regress my series on multiple different series? Intuitively I would use gls because of serial correlation (I am not so sure about the cor structure). I heard about spurious regression and understand that this is a pitfall, nevertheless I am looking for a way to explain a trend. Is this completely wrong or uncommon? Or have I just missed the right chapter so far?


title: How to generate uniformly distributed points on the surface of the 3-d unit sphere?

text: I am wondering how to generate uniformly distributed points on the surface of the 3-d unit sphere? Also after generating those points, what is the best way to visualize and check whether they are truly uniform on the surface $x^2+y^2+z^2=1$?


title: Can we compute bivariate from marginal distributions?

text: Can we do that? If yes, then what are the conditions which should be met? Thanks.


title: How to create a dependent variable in regression when it does not exist in data?

text: ![enter image description here](http://i.stack.imgur.com/upzR9.jpg) In ques 5, I have 10 factors which could affect purchasing decision. Now i wanted to find the co-relation between 10 factors which i will treat as Independent factors and Purchase decision which i intend to use as a dependent variable but couldn't quantify it with anything. Now in order to apply regression, i need an actual dependent variable which can be used to co-relate with other factors. How can i get around this problem or is there any other test which i can use.


title: sample variance of vector variable

text: I've been stuck on this problem for several days: In a sample of size $I$, each sample member $i=1,\ldots,I$, is described by some feature $m(i)$. Value $B$ is a certain function of $\mbox{avg}(m(i))$, that can be approximated by: $B=1/\mbox{avg}(m(i))$. Now I have to analyze the following expression: ![](http://i.stack.imgur.com/ERwh3.png) If $m(i)$ is a scalar, then this expression equals $B*\mbox{Var}(m(i))$, where $\mbox{Var}(\cdot)$ is a sample variance. But what if $m(i)$ is a $1 \times K$ row vector (which is a second case that I have to consider), and $B$ equals inverse of of diagonal matrix with $\mbox{avg}(m(j))$, for $j=1,\ldots,K$ on it's diagonal and zeros elswhere, i.e. $B$ is positive definite matrix, so there exists unique square root of $B$. But I don't see how to express it using variance of $m(i)$ or something of that sort. What I hoped for was to express it as a quadratic form with some variance- covariance matrix associated to it, simmilar to standard expression for e.g. portfolio variance, just with other arguments. Or to manage to represent it as some quadratic form. It would be more than appreciated if someone could actually help me with this!


title: Incremental one-way ANOVA

text: Imagine you have a set of four elements (A-D) with some numeric values of a measured property (several observations for each element):               A: 26 25 29 21     B: 24 17 16     C: 32 34 29 19 25 27 28     D: 23 29 26 20 14      I have to detect if there are significant differences on the average levels. So I run a one way ANOVA to determine if differences are found. It works fine but, when I get a new sample, I need to execute the ANOVA again for just one new sample. Is there any way to do an "incremental one-way ANOVA" that works **without redoing the whole computations** ? I have a lot of data and the process is expensive in terms of time and memory consumption.


title: Principal components of spatial variables

text: I have a number of rasters of environmental data (~10) which may be important predictors for modelling species presence and abundance at ~10 different locations. I would like to know which of the rasters are important in explaining the variance in observed results. Is it appropriate to look for principal components of the rasters? Is this done across the entire extent, or only with respect to my observed sites? Any reference on this would be appreciated. EDIT: It turns out that PCs are not useful in this case. Even though two rasters may have very large correlation across the landscape, if they are different at the narrow range of sites occupied by the modelled species, this is what really matters.


title: nMDS in vegan for soil data

text: I am working with abiotic soil data such as bulk density, moisture levels and soil chemistry as response data (some quantitative some as percentages) and a mix of abiotic and biotic data as environmental data. To compare different sites I want to use multivariate statistics. For this PCA is an accepted method in this field, but with PCA I get an arch/horse shoe effect in my ordination plots.   I have been pointed to nMDS Non-metric multidimensional scaling as an alternative indirect ordination method used in ecology. I intent to use the metaMDS function from the vegan package in R, which gives me reasonable results with the standard settings. But there are some questions that I couldn’t find clear answers to: The help says that metaMDS is designed to be used with community data. Is in this context only biotic data meant? What would other data in this context mean? And what would my soil data fall into? In case my soil data is not community data a range of settings has to be changed. This includes the dissimilarity index. As I am dealing with abiotic data, I generally would assume a linear gradient (as opposed to unimodal), should I in this case chose euclidean distance as dissimilarity index?               rankindex(env, soil, c("euc","man","bray", "jac", "kul"))      gives me the lowest score for euclidean (man and bray the highest). **EDIT** I found several publications that point towards using the Gower metric dissimilarity index for soil chemical data. See the following book page 320 for more info: Legendre, Pierre, and Loic FJ Legendre. Numerical ecology. Vol. 20. Elsevier, 2012.. **END EDIT** Do I have to change the name from nMDS to MDS if I use euclidean distance, as euclidean would be a metric method?). Should the data be transformed before the calculations? Would autotransform = TRUE be sufficient?


title: Recommendations for MRI classification in R of large dataset (n=100, p=20000)

text: I am working on a magnetic-resonance imaging dataset which includes about 100 observations (= subjects) and 20000 predictors (=voxels). I would like to conduct classification in R using methods like support-vector machines or random-forests. I believe that analysis of large datasets will require special setup & would like to know about recommendations, if thats possible in R and what hardware is recommended.


title: Fourier transform of distributions

text: Which distributions are their own Fourier transform besides the _normal distribution_ and the _generalized arcsine distribution_?


title: What is a good way to compare similarity between datasets with little variance?

text: Let's say I have a list of 100 MLB pitchers and 5 statistics for each of them. The difference between, for example, an ERA of 3.5 and 3.1 might not look like a lot to a naive similarity algorithm, but is a lot in baseball. Given that a lot of the player statistics that I'm looking at have this little variance, a lot of a little variance like this, what is the best way to calculate similarity between two players? An example of the data might look like this:               Player | ERA | Wins | Strikeouts     --------------------------------     A      | 3.5 | 15   | 180     B      | 3.1 | 12   | 210     C      | 3.4 | 13   | 150      I've used cosine similarity, and the results are all too similar, down to a thousandth of a decimal place.


title: Expectation of the variance of the sampling set without replacement

text: Select $n$ numbers without replacement from the set $\{1,2,...,m\}$, and generate the set $S=\{a_1,a_2,...,a_n\}$. I want to calculate the expectation of the variance for the sampling set $\mathbb{E}[Var(S)]$ and the maximum variance among all samples : $\max{Var(S)}$. Besides, what's the distribution of the sample variance?


title: Calculation confidence interval for estimated logit in logistic regression in R

text: I have read from here and understand how to calculate the estimated logit from a fitted logistic regression model, but how to work on the confidence interval? As it involved a variance-covariance matrix and I think it is better to have a program to do the calculation, rather then doing it by myself. Thanks. ### Edit 01 I have added a script here:               chdage.dummy <- data.frame(chd=c(rep(1,50),rep(0,50)),                                race=c(rep("white",5),rep("black",20),rep("hispanic",15),rep("other",10),                                       rep("white",20),rep("black",10),rep("hispanic",10),rep("other",10)),                                stringsAsFactors=FALSE)     chdage.dummy[,"race"] <- factor(chdage.dummy[,"race"],levels=c("white","black","hispanic","other"))     chdage.lr.02 <- glm(chd~race,data=chdage.dummy,family="binomial")     predict(chdage.lr.02,newdata=data.frame(race="white"))      `predict` function can give me an estimate, but I can't use `confint` outside `predict`, so what can I do?


title: Average of percentages / prove causal relationship between sale dates & margin sales

text: I just wanted to make sure I'm right here. I have a situation where I need to prove that on days where we sell more than 10 items, our margin (sale price vs. suggested price) goes up. I'm not really sure it's possible to prove this, but this is what I'm working with. The table at the top is an itemized list of each sale we did every day. I notice that if I take a sum of column R on days with more than 10 sales (so 5/1/2011 & 5/4/2011) and divide it by the corresponding totals in column Q ($36,018/$73,920) I get 48.7% (I then take the total of all everything in each column and subtract the corresponding totals, divide it by the appropriate suggested price total- this is shown at the very bottom in column R). Is this the right approach to this? I just need to be sure. The table at the bottom is a separate table I have that aggregates all the info from the above table. Previously I was taking an average of the corresponding percentages (so average 44.2% & 56.3%) but I realized this was likely very incorrect. I'm not completely clear on why, but I avoid taking averages of averages and this situation seemed very similar to me, so that's why I'm posting here to double check. This calculation gives me very different results (shown in column W). Why is this the wrong way? I hope this all makes sense. I know this is one of those fundamental rules of math, yet my brain is fried right now so I need to understand fully why the second method is wrong. I'm not really seeing any analogous examples online. (*Please also note that this is just a subset of my data. The real dataset consists of 3+ years of records and a few other data points, but this was the main one I was asked to look into, though I'm not convinced a causal relationship can be proven from this data alone.) Thank you. ![Example](http://i.stack.imgur.com/SO0ob.png)