title: Outliers change the sign of linear relation

text: I have a scatter plot with Hgb A1C as the explanatory variable ('x' axis) and FPG ('y' axis)as the result variable. There is a positive linear relationship but it is weak. There is 1 outlier in the y direction. There is 1 outlier in the x direction. When each of these points is removed, r changes in opposite directions. I am trying to explain why this is happening, in simple language. Would it make sense to say: When the outlier in the y direction is removed, r increases because an outlier that normally falls a distance away from the regression line would decrease the size of the correlation coefficient. When the outlier in the x direction is removed, r decreases because an outlier that normally falls near the regression line would increase the size of the correlation coefficient. Does this make sense?


title: Order statistics of absolute value of bivariate normal distribution

text: Suppose $X_1$ and $X_2$ are bivariate normal and let $|X|_{(1)}$ and $|X|_{(2)}$ be the ordered version of their absolute value. I am interesting in finding the following probabilities or some bounds on it. \begin{align*} Pr(|X|_{(1)} < c_1, |X|_{(2)} < c_2) \end{align*} and \begin{align*} Pr(|X_1| < c_1, |X|_{(1)} < c_2) \end{align*} Can anyone familier with this give some hint or some references? Thank you very much. Hanna


title: Which statistical tool do I apply to measure the association of Brand Identity and Brand Image?

text: I am pursuing my Ph.D. I am stuck with the last research objective. It talks about exploring the association of Brand Identity and Brand Image. Sample Size for measuring Brand Image (Customer Respondents): 1200 (3 sectors with 4 brands each; Total Brands 12; Sample size of each brand 100) Sample Size for measuring Brand Identity (Brand Managers are the respondents): 48 (3 sectors with 4 brands each; Total Brands 12; Sample size of each brand 4) Which statistical tool do I apply to measure the association of Brand Identity and Brand Image (The sample size for both the entities are unequal; Ex Brand Image responses (customer responses) for Brand Lux= 100 and Brand Identity responses (Brand Managers responses) for Brand Lux = 4)?


title: Probability that 2 OH NFL teams go 31 weeks w/o wins on the same day

text: I did this the obvious way, and my friend came back with a better idea. Can you guys adjudicate or improve on both? My way: > The Cincinnati Bengals and the Cleveland Browns both won on Sunday for the > first time in 46 weeks (says ESPN). That seemed way too improbable... > > 46 is too high though. If we account for bye weeks, MNF nights, head-to-head > games, etc., we get 31 weeks where they each had a chance to win. > > Now we can take their respective records since 2009 (CLE: 11-31, CIN: 18-24) > to calculate win probabilities for any given week. This gives a probability > of 11% that both win in the same week (assume independence). > > So...the probability of this 31-week drought? 2.5%...statistically > significant but not earth-shattering. For reference, if these teams had even > odds of winning any given week, the probability would plummet to 0.01%! My friend's response: > man this is the most thought provoking post i've ever read…now i've spent 30 > minutes thinking about it. anyway…i'm probably making myself sound like an > idiot right now, but i'm not sure that independence assumption is correct. i > think a more accurate way to think about this problem is the classical jar > and ball problem. so if we disregard head to head games and mnf games and > all that, the bengals and browns have played 42 games. now let's first put > the 11 wins by the browns into separate jars. so if we now put the 18 wins > by the bengals one by one into the jars, there is a 31/42 chance that the > first win won't end up in a jar with a browns win…the second one has a > probability of 30/41 (since it doesn't have the option to end up in the same > jar as a previous bengals win)…third one has a 29/40 chance of not ending up > in a jar with a bengals win…so on and so forth. if we think about it that > way, the chances that a bengals win and browns win does not end up in the > same jar after the 18 bengals wins and 11 browns wins have all been put into > jars is ~.058%. > > anyway…just figured it might be a slightly more accurate way to think about > it since the probability of two team with winning percentages above .500 not > having wins in the same weekend over the span which they were both above > .500 is 0%…i think. That all makes sense except for the intuitive feeling that the games were independent events (assuming neither team thought about the other). Who's right? Thanks!


title: Comparing fits in simulation for power analysis

text: I have the following data (just an small example here). I want to know effect of time and ecotype over height. `subdata <- id time ecotype height 1 17 a 73.692 1 22 a 213.010 1 25 a 343.700 1 28 a 663.030 109 32 a 1267.300 109 17 b 47.445 109 22 b 148.050 109 25 b 280.570 109 28 b 509.000 109 32 b 954.050 121 17 c 33.972 121 22 c 71.235 121 25 c 130.650 121 28 c 220.190 121 32 c 452.740 133 17 d 38.365 133 22 d 84.068 133 25 d 137.560 133 28 d 255.390 133 32 d 426.970` I have fitted the following model `m1 <- lme(height ~ ecotype*time + (time | id))` and now I want to know the Power of this analysis in order to do that I use simulation `time <- (subdata$time) ecotype <- as.factor(subdata$ecotype) id <- subdata$id height=subdata$height` `m2 <- lmer(height ~ ecotype * time + (time| id), subdata) s2 <- simulate(m2) beta.hat=fixef(m2) se=sqrt(diag(vcov(m2)))` `k <- c() B=1000 tstar=rep(0,B) set.seed(781) for(b in 1:B) { ystar=drop(simulate(m2)) ostar=lmer(ystar$sim_1~ecotype +(time|id)) for(i in 2:19) { k <- append(k, ((fixef(ostar)[i]-beta.hat[i])/sqrt(vcov(ostar)`[i,i]))) }}` The problem I have is I don't know how to calculate an overall p-value to say this model's power is low or high. Normally, the examples on the web always have variables (e.g ecotype) with two classes so the results of the simulation will always produce one value. However, in my case I have n number of ecotypes therefore more than one output: `fixef(ostar) (Intercept) ecotypea ecotypeb ecotypec ecotyped 78.98846 -22.39385 -37.27676 -38.66168` All I can think is to extract all the fixed effects estimates and the cv and compare them against, get a mean and compare it against a threshold (e.g. 0.05). However, I don’t know whether that’s right. Please can you give me a hand on this Please can you suggest how can I do that? Thanks


title: What is the Australian Calculator?

text: i need to know how to use the Australian Calculator to get my sample size. i was instructed to carry out my research using the Australian Calculator, but i don't know what it means.


title: Is skewness always bad?

text: In my experiment, I hypothesised that individuals in one treatment condition would give higher values on a likert scale than individuals in the other treatment condition. It was a one tailed hypothesis. Histograms and Q plots show that my data is skewed in the direction I would imagine for the first treatment condition - i.e. all clustering around the higher values in the Likert scale. This isn't the case, again as expected, in the second treatment condition. Even though I predicted this, is this skew still bad? Should I always use non- parametric tests in this situation (the results of parametric and nonparametric tests are the same - I just want to make sure I am using the right tests!). Any help would be really gratefully received! Thanks!


title: Significance test based on precision/recall/F1

text: Is it possible to do a significance test based solely on precision/recall/F1 scores? For example, if you come across 2 systems in a paper for which only P/R/F1 are reported (on the same dataset, etc.), can you then perform a statistical significance test? If yes, how is that done?


title: Multivariate data analyis of compositional data

text: Suppose I have a multivariate, compositional dataset that depicts the concentration of different elements. However, the data are not available on a single scale; i.e., some are of form 0.00x while others are integers. Should I apply any kind of normalization / standardization technique before applying, or do I need to do anything with the data at all before apply transformations (isometric log ratio [ilr], centered log ratio [clr], etc.) and start the data analysis, imputation on missing value using robust methods, robust PCA and data clustering? Some pointers for understanding compositional data analysis are also welcomed. [Update] For example: Consider two vectors:               [ 0.016, 71.2, 0.123, 1.74, 14.0, 0.002, 2310, 0.064, 0.29, 0.32,5.63, 96.5, 0.044,        4360, 1110, 585, 0.052, 62.9, 4.45, 1110, 1.50, 15.10, 783, 0.015,78.9, 5.61, 0.007,        0.022, 0.007, 0.53, 29.3 ]       [ 0.073, 245.0, 0.299, 2.77, 17.4, 0.039, 2460, 0.145, 0.85, 0.99, 20.40, 359.0 0.062,        4040, 1530, 148, 0.113, 217.0, 18.10, 1310, 4.61, 4.56, 880, 0.069, 230.0, 12.20,        0.028, 0.025, 0.013, 9.92, 34.1]      These two represent concentrations of different elements in soil samples collected from two different positions. If I were to analyse them using robust methods, should I preprocess them standardize / normalize in any form, or should I simply transform them into Aitchison geometry and start my analysis?


title: Computing probability from the Poisson distribution

text: I have a question regarding the `ppois` function in R:   * `ppois(x,lambda,lower.tail=TRUE)` gives $P[X>x]$   * `ppois(x,lambda,lower.tail=FALSE)` gives $P[X\le x]$ How can I calculate $P[X=n]$? I mean is there any function in R for this?


title: Measuring accuracy of a logistic regression-based model

text: I have a trained logistic regression model that I am applying to a testing data set. The dependent variable is binary (boolean). For each sample in the testing data set, I apply the logistic regression model to generates a % probability that the dependent variable will be true. Then I record whether the acutal value was true or false. I'm trying to calculate an $R^2$ or Adjusted $R^2$ figure as in a linear regression model. This gives me a record for each sample in the testing set like:               prob_value_is_true         acutal_value        .34                          0        .45                          1        .11                          0        .84                          0         ....                        ....                I am wondering how to test the accuracy of the model. My first attempt was to use a contingency table and say "if `prob_value_is_true` > 0.80, guess that the actual value is true" and then measure the ratio of correct to incorrect classifications. But I don't like that, because it feels more like I'm just evaluating the 0.80 as a boundary, not the accuracy of the model as a whole and at all `prob_value_is_true` values. Then I tried to just look at each prob_value_is_true discrete value, as an example, looking at all samples where `prob_value_is_true`=0.34 and measuring the % of those samples where the acutal value is true (in this case, perfect accuracy would be if the % of samples that was true = 34%). I might create a model accuracy score by summing the difference at each discrete value of `prob_value_is_true`. But sample sizes are a huge concern here, especially for the extremes (nearing 0% or 100%), such that the averages of the acutal values are not accurate, so using them to measure the model accuracy doesn't seem right. I even tried creating huge ranges to ensure sufficient sample sizes (0-.25, .25-.50, .50-.75, .75-1.0), but how to measure "goodness" of that % of actual value stumps me. Say all samples where `prob_value_is_true` is between 0.25 and 0.50 have an average `acutal_value` of 0.45. Is that good since its in the range? Bad since its not near 37.5% (the center of the range)? So I'm stuck at what seems like should be an easy question, and hoping someone can point me to a resource or method to calculate an accuracy stastic for a logistic regression model.


title: got stuck in Cluster analysis, way forward?

text: I have a situation, where I need to classify items into groups (lets say 6). When I ran k-means 90% of my data fall in 1 group remaining 10% fall in other groups. What's next step? In order to further group the data, I have taken the 90% data group and once again I ran k-means.This time I have 15 new groups within this new dataset. But now again 76% fell in one group remaining in 14 groups? How to deal in such situation?


title: repeated measures ANOVA - non-experimental data

text: For every year, for eight years, I have a set of projects (n=332). The number of projects that exist in each year varies (for example, some years there are 10 projects others there are 60). In each year, projects are either approved or not approved based on a score. I want to compare the means for approved and not approved projects for each year as well as compare the scores for approved and not approved projects over time. A repeated measures ANOVA with contrasts lets me do this but because the projects are not the same every year I believe this is not the right test. Thoughts? Also, the scores used to approve or not approve a project is a sum of several sub-scores. I want to compare the standard deviations about the project score for both approved and not approved project and over time. PS a mixed model, projects within years in STATA won't converge


title: Marginally Uncorrelated Responses in Multivariate Regression

text: Consider the multivariate regression with random predictors as described in Izenman's _Modern Multivariate Statistical Techniques_ , that is, $X : r \times 1$ and $Y: s \times 1$ are jointly normally distributed with respective mean vectors $\mu_X$ and $\mu_Y$ and they have joint covariance matrix $$\left( \begin{array}{cc} \Sigma_{XX} & \Sigma{XY} \\\ \Sigma_{YX} & \Sigma_{YY} \end{array} \right)$$ In the general regression problem, we are concerned with the covariance matrix of $(Y|X=x)$, which is $\Sigma_{YY} - \Sigma_{XY}^T \Sigma_{XX}^{-1}\Sigma_{XY}$. However, my question is: in a real data setting, would we expect our response $Y$ to be marginally uncorrelated? That is, would we expect $\Sigma_{YY}$ to be the identity matrix, or would there be some marginal correlation between the responses? I have my own opinion on this but am interested in a fresh perspective. Simple examples or anecdotes about your own experience with this data would be much appreciated.


title: What math subjects would you suggest to prepare for data mining and machine learning?

text: I'm trying to put together a self-directed math curriculum to prepare for learning data mining and machine learning. This is motivated by starting Andrew Ng's machine learning class on Coursera and feeling that before proceeding I needed to improve my math skills. I graduated from college a while ago so my algebra and statistics (specifically from political science/psychology classes) are rusty. The answers in the thread Is a strong background in maths a total requisite for ML? only suggest books or classes directly related to machine learning; I have already looked into some of those classes and books and do not know exactly what math subject to study (for instance: what field[s] of math address deriving an equation to "minimize a cost function"?). The other thread suggested (Skills & coursework needed to be a data analyst) only mentions broad categories of skills needed for analyzing data. The thread Introduction to statistics for mathematicians does not apply because I do not already have a degree in math; a similar thread Mathematician wants the equivalent knowledge to a quality stats degree has an incredible list of stats books, but again, I'm looking at starting math from a rusty recollection of algebra and moving up from there. So, for those that work in machine learning and data mining, _what fields of math do you find essential to do your job?_ What math subjects would you suggest to prepare for data mining and machine learning, and in what order? Here is the list and order I have so far:   * Algebra   * Pre-calculus   * Calculus   * Linear algebra   * Probability   * Statistics (many different sub-fields here, but don't know how to break them out) As for the data mining and machine learning, through my current job I have access to records on website/app activity, customer/subscription transactions, and real estate data (both static and time-series). I'm hoping to apply the data mining and machine learning to these datasets. Thank you! EDIT: For posterity's sake, I wanted to share a helpful math self-assessment for Geoffrey Gordon's/Alex Smola's Intro to Machine Learning class at CMU.


title: Predicting one time series from another if they are related - algorithm in R

text: I am new to time series analysis; I am not even sure if this is even a TS problem. I have looked at other TS posts, but I have a hard time to translate the responses to my needs. For now I am hoping for help to define the procedure on how to solve my problem. For the future I hope to be able to "fish" by myself. Thus it would be great to get help with some example R code, hints on textbooks to look at, and some idea on how to interpret the results from what analysis will be suggested. The problem is that I don't have data yet. So this is what I expect. I hope that this makes sense. I am trying to keep this more abstract, because I tend to be too applied, and I like to understand the theory better. I expect to have two time series A and B. A and B are related; A causes B, but in an unknown association and with some delay. Visually, I am expecting to see a profile for A and B that appear to be similar, but are substantially off from each other. I hope that the number of peaks are similar. However, aside from a lag I expect a lot of differences between peak height, peak width, area under the peak (assuming some normalization) etc. I expect the differences to be mix of systematic and random effects. I need to predict A from B profiles only. For several locations I will have a test setup with at least 10 pairs of A profiles. I will have approximately 5 groups of 5 B profiles for one A profile. I am expecting to use information that I gain from this analysis to apply to a set of B profiles and retrospectively predict A profiles. What would be a good approach to predict A profiles from a B profile? Any applied literature that does deal with natural processes ie geology to help with this problem? Any other thoughts? Thank you, Boris


title: Can McNemar's test be improved upon by adjustments for zeros like those in a sign test?

text: McNemar's test is a special example of the binomial sign test, but the ('vanilla') sign test suffers from bias due to ignoring differences equal to zero. McNemar's test statistic is given by: $\chi^{2} = \frac{\left(|r-s|-1\right)^{2}}{r+s}$, where $r$ and $s$ are the counts of discordant pairs (0,1) versus (1,0), distributed $\chi^{2}$ with 1 degree of freedom under the null hypothesis. I am having a hard time parsing Sribney on the sign test: > The test statistic for the sign test is the number $n_{+}$ of observations > greater than zero. Assuming that the probability of an observation being > equal to zero is exactly zero, then, under the null hypothesis, $n_{+} \sim > \text{Binomial}(n, p=\frac{1}{2})$, where $n$ is the total number of > observations. But what do we do if we have some observations that are zero? > > **Fisher’s Principle of Randomization** > > We have a ready answer to this question if we view the test from the > perspective of Fisher’s Principle of Randomization (Fisher 1935). Fisher’s > idea (stated in a modern way) was to look at a family of transformations of > the observed data such that the a priori likelihood (under the null > hypothesis) of the transformed data is the same as the likelihood of the > observed data. The distribution of the test statistic is then produced by > calculating its value for each of the transformed “randomization” data sets, > considering each data set equally likely. > > For the sign test, the “data” are simply the set of signs of the > observations. Under the null hypothesis of the sign test, $P(X_{i}>0)= > P(X_{i}<0)$, so we can transform the observed signs by flipping any number > of them and the set of signs will have the same likelihood. The $2^{n}$ > possible sign changes form the family of randomization data sets. If we have > no zeros, this procedure again leads to $n_{+} \sim \text{Binomial}(n, > p=\frac{1}{2})$. > > If we do have zeros, changing their signs leaves them as zeros. So if we > observe $n_{0}$ zeros, each of the $2^{n}$ sign-change data sets will also > have $n_{0}$ zeros. Hence, the values of $n_{+}$ calculated over the sign- > change data sets range from 0 to $n-n_{0}$, and the “randomization” > distribution of $n_{+}$ is $n_{+} \sim \text{Binomial}(n-n_{0}, > p=\frac{1}{2})$. Because this seems to be saying go ahead and ignore zeros. But then later in the paper, Sribney provides an adjustment for the sign-rank test that accounts for zeros just along the lines I am wondering about: > The adjustment for zeros is the change in the variance when the ranks for > the zeros are signed to make $r_{j}=0$; i.e., the variance is reduced by > $\frac{1}{4}\sum_{i=1}^{n-{0}}{i^{2}}=n_{0}\frac{\left(n_{0}+1\right)\left(2n_{0}+1\right)}{24}$. Should I instead be asking whether or not to apply the signed-rank test to individually-matched case control data? A simple made up example will illustrate why ignoring zeros presents a problem. Imagine you've paired data with no differences equal to zero (this would correspond to data for a McNemar's test with only discordant pairs present). With a sample size of, say, 20, you find 15 positive signs of differences and 5 negative signs of differences, and conclude significant difference. Now imagine that you have 1000 observed differences equal to zero in addition to those 15 positive and 5 negative signs of differences: now you conclude difference is not significant. If McNemar's test is conducted on 1020 pairs, 1000 of which are zeros, and with discordant pairs of 15 and 5, we should not reject the null hypothesis (e.g. at $\alpha = 0.05$). There is an adjustment to the sign test to correct for observed zero differences based upon Fisher’s "Principle of Randomization" (Sribney, 1995). Is there a way of improving on McNemar's test that addresses the effect of observed zero differences (i.e. by accounting for number of concordant pairs relative to number of discordant pairs)? How? What about for the asymptotic _z_ approximation for the sign test?    **References** Sribney WM. (1995) Correcting for ties and zeros in sign and rank tests. _Stata Technical Bulletin_. 26:2–4.


title: Measure topic centrality

text: Does anyone know any metric that measure **_topic centrality_** or topic stickiness (I've seen someone used this term too). For example _Cluster 1_ $d_1$ = { a, b, c }, $d_2$ = { b, c, d } and _Cluster 2_ $d_3$ = { a, b, c }, $d_4$ = { d, e, f } How can we compute a numerical value that tells us cluster 1 is better than cluster 2. In this case I guess we can measure **cluster purity** by calculating pairwise similarity of documents in the cluster. However, if the clusters are big (contain many documents) then it becomes computationally expensive. Is there any other more commonly used metric for this?


title: Logistic Regression not quite working

text: I'm playing around with logistic regression (using scikit-learn, which uses liblinear) I created some example data sets, and it often works well. Since I created the data sets, I already know the optimal bounds, and a 100% correct linear classifier is possible. Fairly often (~1 in 10) the predicted bound is off (see figure). It seems obvious to me that moving the bound a fraction higher would be better. What might have gone wrong? Are there parameters I should try changing, or better algorithms to use? ![example of failing](http://i.stack.imgur.com/Yfw0p.png) Thanks!


title: Multilevel regression problem: some level 2 cases are part of more than one level 1 category

text: I am new to multilevel models and I am basically wondering whether I can use a multilevel model at all, given the structure of my data. I have collected data on membership of lobby groups in committees used by the Directorates-General (DGs) of the European Commission (the DGs are a bit like ministries in a national context). I have data on the number of committee seats occupied by each group for all DGs. Level 2 of the analysis would be lobby groups, level 1 would be the DGs. Now, my problem is that many lobby groups are present in more than one DG. It is my understanding that multilevel regression should normally only be used if the cases on the lower level are nested in those of the higher level (i.e. students are parts of a class, classes are part of a school, etc.). In effect, the "students" (lobby groups) in my data are members of several "classes" (DGs) each. Is it still possible to do a multilevel model in this case? I am thinking of "forcing" the lobby groups to be nested in the DGs by adding data rows, so that one lobby group present in 5 DGs effectively becomes 5 identical lobby groups present in 1 DG each. Is this advisable?


title: Computing the pre-study probability of a finding being true

text: I read a fascinating free article about why most published research findings might be false. I am able to follow along almost completely, except for one tiny part. In it, the author, Ioannidis, writes               Let R be the ratio of the number of “true relationships” to “no relationships”     among those tested in the ﬁeld.     ...     The pre-study probability of a relationship being true is R⁄(R + 1).          Why is `P(T)=R/(1+R)`? He writes in the same paragraph:               Let us also consider, for computational simplicity, circumscribed ﬁelds where     either there is only one true relationship (among many that can be     hypothesized) or the power is similar to ﬁnd any of the several existing true     relationships.      Perhaps this can help explain it, but I do not see how.


title: Goodness of fit of an approximation of a PDF

text: I'd like to evaluate the performance of an algorithm which learns to approximate a discrete PDF for a latent variable given some noisy input. Now, is there a standard test to evaluate the goodness of that approximation? I would like to train the algorithm with synthetic data and then see how close the approximation is to the PDF I can compute knowing how the data was generated. I thought of computing the Kullback-Leibler divergence, but   * I don't know what's a good value. KL divergence is the information gain, so could I argue that eg. losing only <= 5% of the total information when using the approximation is good?   * How do I combine KL divergences from multiple trials with different input? Can I average? I suppose I might, considering that all the computed PDFs can be seen as one PDF over a multi-dimensional latent variable.


title: Is there a standard technique to debug MCMC programs?

text: Debugging MCMC programs is notoriously difficult. The difficulty arises because of several issues some of which are: (a) Cyclic nature of the algorithm We iteratively draw parameters conditional on all other parameters. Thus, if a implementation is not working properly it is difficult to isolate the bug as the issue can be anywhere in the iterative sampler. (b) The correct answer is not necessarily known. We have no way to tell if we have achieved convergence. To some extent this can be mitigated by testing the code on simulated data. In light of the above issues, I was wondering if there is a standard technique that can be used to debug MCMC programs. **Edit** I wanted to share the approach I use to debug my own programs. I, of course, do all of the things that PeterR mentioned. Apart from those, I perform the following tests using simulated data:   1. Start all parameters from true values and see if the sampler diverges too far from the true values.   2. I have flags for each parameter in my iterative sampler that determines whether I am drawing that parameter in the iterative sampler. For example, if a flag 'gen_param1' is set to true then I draw 'param1' from its full conditional in the iterative sampler. If this is set to false then 'param1' is set to its true value. Once I finish writing up the sampler, I test the program using the following recipe:   * Set the generate flag for one parameter to true and everything else to false and assess convergence with respect to true value.   * Set the generate flag for another parameter in conjunction with the first one and again assess convergence. The above steps have been incredibly helpful to me.


title: Does the cross validation implementation influence its results?

text: As you know, there are two popular types of cross-validation, K-fold and random subsampling (as described in Wikipedia). Nevertheless, I know that some researchers are making and publishing papers where something that is described as a K-fold CV is indeed a random subsampling one, so in practice you never know what is really in the article you're reading.   Usually of course the difference is unnoticeable, and so goes my question -- can you think of an example when the result of one type is significantly different from another?


title: Why do AIC and BIC show inversed outputs?

text: I am comparing three relatively simple GLMs having a Gamma distribution with AIC and BIC. The aim is to identify the effects of fertilizers (fdung), year and site on biomass of a specific grass species. Hence, the aim is not to predict new values but merely to identify the effects of the three factors. Here are the models used:               res1 <-  glm((Biomass..g.m².) ~  fdung * fyear * fblock, family=Gamma(link="identity"))     res2 <-  glm((Biomass..g.m².) ~  fdung * fyear + fblock, family=Gamma(link="identity"))     res3 <-  glm((Biomass..g.m².) ~  fdung + fyear + fblock, family=Gamma(link="identity"))      I expect the third model to be the most simplistic one and want to confirm this by an information criterion. However, when looking at AIC and BIC I get this output.               AIC(res1,res2,res3)     BIC(result1,res2,res3)             df      AIC            df      BIC     res1    49 5271.617     res1   49 5465.198     res2    16 5334.234     res2   16 5397.44     res3    10 5331.253     res3   10 5370.760      For AIC, the most complex model is "best" and for BIC the one with fewest df is best. I am thinking that with regard to my aim (identify effects on biomass) I should trust BIC. Am I wrong here with my conclusion? I already tried mixed effect models with the fblock as random factor but then the model with the Gamma distribution did not work any more and also I could not use fblock as fixed effect any more (leading to NAs for fblock), but this is not part of my question.


title: $E(x^k)$ under truncated $\mathcal{N}(\mu,1)$

text: There is a similar question in $E(x^k)$ under a Gaussian. However, it doesn't seem to be trivial when $\mu\ne0$. As mentioned in the previous question **$k$ is not an integer**. The integral that I need to evaluate is as follows: $$\int_0^\infty x^k\exp\left(-\frac{(x-\mu)^2}{2}\right)dx$$ If it helps for the case that $\mu=0$ the answer is $\frac{2^{(k-2)/2}}{\sqrt{\pi}}\Gamma(\frac{k+1}{2})$


title: PCA and PLS: testing variables for significance

text: I'm trying to understand the process for statistical testing for principal component analysis or partial least squares. **Step 1. PCA:** I feel that I have a not-terrible understanding of PCA: You find the ellipsoid described by the covariance matrix of the data, and then successively take the largest axis of variation (principal component 1), then the second largest (principal component 2), and so on. If the ellipsoid is long and stretched, then the variation is mostly along the first principal component (the eigenvector corresponding to the largest eigenvalue of the ellipsoid). If the ellipsoid is a planar "disc", then the variation in the data is explained well by two principal components, etc. I also understand that after choosing to use (for example) only the first two principal components, then all of the data points can be plotted on a "Scores" plot that shows, for each data point $D^{(i)}$, the projection of $D^{(i)}$ into the plane spanned by the first two principal components. Likewise, for the "Loadings" plot (I think) you write the first and second principal components as linear combinations of the input variables and then for each variable, plot the coefficients that it contributes to the first and second principal components. **Step 2. PLS or PLS-DA:** If there are labels on the data (let's say binary classes), then build a linear regression model to use the first and second principal components to discriminate class 0 (for data point $i$, that means $Y^{(i)}=0$) from class 1 (for data point $i$, that means $Y^{(i)}=1$) by first projecting all data to only lie in the plane spanned by the first and second principal components, and then regressing the projected input data $X_1', X_2'$ to $Y$. This regression could be written as (first step) the affine transformation (i.e. linear transformation + bias) that projects along $PC_1, PC_2$ (the first and second principal components), and then (second step) a second affine transformation that predicts $Y$ from $PC_1, PC_2$. Together these transformations $Y \approx Affine(Affine(X))$ can be written as a single affine transformation $Y \approx C (A X + B) + D = E X + F$. **Step 3. Testing variables from $X$ for significance in predicting the class $Y$:** This is where I could use some help (unless I'm way off already, in which case tell me!). How do you test whether an input variable (i.e. a feature that has not yet been projected onto the principal components (hyper)plane), and decide if it has a statistically significant coefficient in the regression $Y \approx E X + F$? Qualitatively, a coefficient in $E$ that is further from zero (i.e. positive and negative values with large magnitude) indicates a larger contribution from that variable. I remember seeing linear regression t-tests for normally distributed data (to test whether the coefficients were zero). Is this the standard approach? In that case, I would guess that ever variable from $X$ has been transformed to have a roughly normal distribution in Step 0 (i.e. before any of these other steps are performed). Otherwise, I could see performing a permutation test (by running this entire procedure thousands of times and each time permuting $Y$ to shuffle the labels, and then comparing each single coefficient in $E$ from the un-shuffled analysis to the distribution of coefficients from shuffled analyses). Can you help me see anywhere my intuition is failing? I've been trying to look through papers using similar procedures to see what they did, and as is often the case, they're clear as mud. I'm preparing a tutorial for some other researchers, and I want to do a good job.


title: What's the correct way to test the significance of classification results

text: There are many situations where you may train several different classifiers, or use several different feature extraction methods. In the literature authors often give the mean classification error over a set of random splits of the data (i.e. after a doubly nested cross-validation), and sometimes give variances on the error over the splits as well. However this on its own is not enough to say that one classifier is significantly better than another. I've seen many different approaches to this - using Chi-squared tests, t-test, ANOVA with post-hoc testing etc. What method should be used to determine statistical significance? Underlying that question is: What assumptions should we make about the distribution of classification scores?


title: Weak law of large numbers - redundant?

text: I might be missing something basic - but it appears that the strong law of large numbers covers the weak law. If that case, why is the weak law needed? Thanks!   Liran


title: Statistical properties of parameters estimated by method of Lagrange multipliers in R

text: I have time series data that can be fitted by a constrained non-linear function (non-linear in the parameters to be estimated). I'd like to use the method of Lagrange Multipliers (link), specifically R's `Rsolnp` package (link), because the Lagrange Multiplier itself has economic meaning. I would like the statistical properties (p values or confidence intervals) of the parameters estimated by `solnp` in `Rsolnp`. If my problem didn't involve constraints, I could use R's `nls` class (link) and the `summary.nls` method (link) to retrieve the statistical properties of the fitted parameters. Questions:   * Is there an equivalent to `summary.nls` for `solnp`?   * If not, how can I obtain the statistical properties of the parameters estimated by the method of Lagrange Multipliers? A reference to a paper or textbook here will be sufficient. Here are some details of my function (a CES production function from economics) and constraints: Function:               iGDP ~ ((1-zeta) * (exp(lambda_L*iYear) * iCapStk^(alpha) * iLabor^beta)^phi + zeta*(exp(lambda_E*iYear) * iQ)^phi)^(1/phi)      Constraints:               alpha + beta = 1        (Note that I don't want to substitute 1-beta for alpha, because the numerical value of the Lagrange Multiplier has a valuable economic interpretation. I want to use the both the constraint equation and the method of Lagrange Multipliers.) Known data (on an annual basis):               iGDP, iYear, iLabor, iCapStk, iQ      Parameters to be estimated (these are the parameters for which I would like statistical properties: p values, confidence intervals, etc.):               zeta, lambda_E, lambda_L, alpha, beta, phi


title: Variable entered in logistic regression model is part of another variable entered in the same model

text: I´m trying to find variables predicting a disease by using first logistic regression for each variable on the disease and then entering the significant variables into a multiple logistic regression model. However, one of the variables in the multivariate model is a clinical score, which contains some of the variables adjusted for (among others bmi<18). My question is: does it make sense to have both in the multivariate model? And if yes, what is the interpretation of an insignificant clinical score in the multivariate model? Is it insignificant because I controlled for some of its variables?


title: How to do Regression with many replicates on the dependent variable at each level of the independent variable?

text: How do I do regression with many replicates on the dependent variable at each level of the independent variable? Pardon my ignorance on this as its my first time posting here.   I'll try to be as precise as can be.   I have data that looks as follows:                 X,    y1,       y2,         y3,         y4        0,   54.241,   127.728,  127.73802,  127.73802         31,   65.132,   127.729,  127.73787,  127.73792         59,   65.364,   127.729,  127.73782,  127.73789      where X is the independent variable and y1-y41 are replicates that were measured at the particular level of X. This is clearly different from the usual case of regression where we have one column of X and another column for y! I had the idea that I should just average the values y1-y41 and do my regression on the mean of those values but I would like my analysis to be optimal. Thus, I have tried to fit random effects model with the replicate as random in `lmer` package in R. At first I transformed the data to the long for so it looked as follows:                 X,  replicate,  value       0,    y1,      54.241       31,   y1,      65.132        59,   y1,      65.364        90,   y1,      66.465      120    y1,      68.439      I fitted the following model: `lmer(value~X+(1|replicate),data=long_form)` I obtained the following random effects variance components estimates:               Random effects:     Groups Name            Variance   Std.Dev.     replicate(Intercept)    16.292     4.0364     Residual               100.557    10.0278      Now, I have the following questions: 1) Is my model correct? or is there a better model for such data? 2) Is the variance for replicate indicative of random variation between replicates at the same level for the entire experiment i.e what does a variance of 16.292 tell me here? 3) By looking at the distribution of the averages of the y values, it seems like the data is heavily skewed to the left and box-cox function suggested a power transformation of 2 i.e y^2. I have transformed my data and tested for normality using shapiro.test function but the data still isn't normal after transformation. Can I still use a random effects model (which I believe assumes normality) or is there another way of going around this without having to transform the data? I tried to be as complete as I can so please excuse my ignorance and let me know if there is something I need to further clarify. I really would appreciate your help and suggestions with this.


title: Data visualization from hierarchical clustering

text: I have data from questionnaire from school. 35 questions are various questions (influence of friends etc.) Possible answers for 35 questions are "definitely yes", "mostly yes", "mostly no" and "definitely no". I did hierarchical clustering using `hclust` in R. Then I used `cutree` for cut the dendrogram. How to visualize data about clusters from `cutree`? I wrote function for export information about clusters to CSV, but I want to display graphical information. Thanks


title: Is there a statistical test to show that sugar liking flies tend to belong to bins where there are flies that also likes sugar and vice-versa?

text: I have about 900 flies that either like sugar or are neutral (about 50% in each group). However, on closer inspection I do notice a trend that sugar liking flies are from bins 1-10 (each bin has about 45 flies), while the neutral ones tend to be from 11-20 bins. I'm not sure why, but is there a statistical test I can do to show that sugar liking flies tend to belong to bins where there are flies that also likes sugar and vice-versa?


title: How to reference R's automatically generated dummy variables in a linear model

text: R is creating dummy variables for the color variable (color is either white or red) in my linear regression model. For example, color:pH returns the interaction term colorred:pH. Some questions: Why is R only returning `colorwhite` for some interaction terms and `colorred` with other terms? Why doesn't it do `colorwhite` and `colorred` for each term interacting with color? Next, if I want to rebuild a model using only some of the dummy variables R generated, how can I reference them in a linear model? If I retype the model manually and include `colorwhite:pH` as a term, R returns an error that the object `colorwhite` isn't found. ![enter image description here](http://i.stack.imgur.com/dVR9H.png)


title: counfounded batch effects in microarray dataset - can I do partial experiment redesign?

text: I'm working with a microarray data set where the batches are completely confounded treatment time, i.e. time t1 is all in batches b1 and b2, and time t2 is all in batches b3, b4, and b5. I know this means that I can't do much abound looking for differential expression over time. However, we do have the original tissue samples and the option to re- extract/reprocess some samples in a new batch. Due to the study size, rerunning all samples with a proper randomized design is out of the question. Are there any studies describing how to rerun a small subset of samples to enable a statistically valid comparison across time points? Or does anyone have any advice? For example, can I reprocess 10 samples from each time point in a new batch - could that be sufficient? Here's the sample breakdown with the relevant treatments/covariates - you can see that timepoints t1 and t2 are completely confounded with batch, while timepoints t2 and t3 are not. I also included another covariate that we want to account for, Ethnicity, and the number of samples in each group:               Batch  Timepoint   Ethnicity  Num_Samples     b1         t1         A         21     b1         t1         B         54     b2         t1         A         20     b2         t1         B         56     b3         t2         A         10     b3         t2         B         35     b4         t2         B         49     b5         t2         A         28     b5         t2         B         4      If we do reprocess samples, do I need to ensure even representation of the two ethnicities in the new batch?


title: post-hoc test after logistic regression with interaction. Risk higher for type 1 or type 2 error?

text: I have a large data set (>1000 obs) and i'm performing regressions tests, both linear and logistic, on a series of clinical outcomes. In this test I verify the effect of interactions between two cat variables by creating a dummy interaction variable with the two vars joined in one (with R is done using the interaction() fuction). In the equation are present also other vars for adjustment. After every test I perform a post-hoc to verify if there is difference between every level of this interaction variable. The test is performed via glht() of the multcomp package in R. I believe, (i'm not 100% sure) that the post-hoc is performed using Tukey methodology. Of course applying the post-hoc procedure I loose significance. Often appears that results that were significant in the regression test become not significant. It's ok. I always think that p<0.05 it's a very shallow threshold, especially when you have large numbers like I have. Some times it even happens that some relevant effects between levels that were not show in the regression become visibile in the post hoc and this is a good thing. My question is therefore, when you perform a post hoc of a regression analysis, using Tukey, and with these large numbers, it's still easier to make a type 1 error and accept no real results, or a type 2 error, being too hard on data and loosing interesting results? Or it's impossible to know? (probably this is the correct answer) Thanks!


title: Transition matrix in left-right hidden semi-Markov model

text: i'm developing a hidden semi-Markov model left-right . In a left-right model a sequence of $M$ states starts in state 1 and ends in state M, with no repetition of states. Since the model is semi-Markovian a self transition $a_{ii}=0$ is not allowed. I think the last row of the matrix should be 0, but the rule is $\sum_{j}a_{ij}=1$. How can i solve this problem?


title: Correcting for restriction of range in multiple regression

text: I've conducted a multiple linear regression with 3 predictors, A, B, and C, all three on the same scale. One of the three predictors displays a smaller range than the other two (range A=0.53; B=0.98; C=0.86), as well as a smaller standard deviation (SD A=0.18; B=0.30; C=0.25).This is precisely the predictor that does not show a significant effect, while the other two (B and C) do. I would like to correct the coefficient (or standardized coefficient) of A for its restricted range. I have only found corrections for the correlation coefficient, such as Thorndike's case 2 (see, e.g., here: http://pareonline.net/pdf/v14n5.pdf ).   1. How can coefficients in multiple regression be corrected? Could I perhaps apply Thorndike's correction to the partial correlation between A and Y, while controlling for B and C? (ideally, however, I would like to have a correction of the regression coefficient)   2. The correction requires knowing the unrestricted SD. This is not possible, but I can use the SD of one of the other predictors. The argument would go: "The lack of an effect for A, as opposed to B and C might be due to its restricted range; with this correction we show that assuming for A as large an SD as for B, then A's coefficient would still be small (or is now as large as that for B)". Is this at all sound? Thank you!


title: Selecting features and estimating their out-of-sample performance with cross-validation

text: I have only a small dataset. I want to 1. select the most predictive features out of a large candidate pool and 2. get an estimate of their expected predictive performance. In _the elements of statistical learning (page 245ff)_, the authors stress the importance of including variable selection within the cross-validation loop for obtaining unbiased estimates of expected out-of-sample performance. However, the estimates of model performance obtained in this manner will not be for one defined set of features, as each cross-validation fold or repetition may lead to the selection of different features. I am, however, not interested in model performance averaged over different sets of features. I want to obtain one set of "best" features and the performance I can expect conditional on them in independent datasets. Do I have any options in simultaneously getting both 1. and 2.? Related question


title: Non-Orthogonal Planned Comparison: How to correct for alpha?

text: I have a data set containing three groups (A, B, C). Group A and Group B have undergone two different types of training designed to improve the same thing. Group C is my control group with no training. I have specific hypotheses stating that: Group A will be better than group C Group B will be better than group C and group A This is because the training used in group A has already been shown to improve over control (so I am trying to replicate that effect), while group B is a new training method which I expect to exceed both groups A and C. How would I go about conducting this analysis? I have conducted a one-way ANOVA with planned comparisons, but I realized that my comparisons would have to be non-orthogonal since comparing B to A and C in one contrast would be pooling variance from an experimental and control condition together. Optimally, I want these comparisons: A B C 1 0 -1 0 1 -1 -1 1 0 What is the proper way to conduct such a comparison? I would greatly appreciate any input!


title: repeated measures survival analysis

text: I have various samples on which three survival time variates have been collected. How can I compare the difference of such variates taking into account within subject correlation (possibly in R)?


title: impossible distribution statistics?

text: I'm currently reviewing an article where authors presented distribution statistics that look erroneous to me. But I'm not able to find a way to ascertain it. The article presented results with a mean of 95% and standard deviation of 25%. Maximum value can't be more than 100% and minimum value can't be less than 0%. I don’t have the sample number. I tried to generate lognormal random numbers with these statistics without success:               # R Code     require("Runuran")     d1 <- urlnorm(n = 1000, meanlog = log(95), sdlog = log(25), lb = 0, ub = 100)      How can I generalize this conclusion whatever the distribution?


title: Holtz-eakin, Whitney, Rosen Panel VAR

text: I'm trying to write a bit if R code to run the estimator in "Estimating Vector Autoregressions with Panel Data" Econometrica 56-6 (1988). I'm stuck on equation (3.12) I feel like $Z_{ir}$ and $Z_{is}$ should be of different lengths, and therefore the dot product $Z'_{ir}Z_{is}$ isn't a valid operation. If someone could take a look at the paper and explain what I'm missing, that would be really awesome.


title: Using decision trees with very infrequent outcomes

text: I am working on decision tree model and a value in dependent variable (churn) is very less. we have 1.5 lac records and only 1700 records have churn = 1. While using decision tree model, tree is not created, quality of model is okay but it is not able to predict where churn should be 1. Please let us know if any other model can be used to tackle the issue. I have used logistic regression but getting similar results. Model is not able to predict churn =1 while validating.


title: R package for marketing

text: I'm taking an intro to marketing class this semester and the course requires us to buy a software which looks nothing more than a bunch of excel macros. From the software's page > Marketing Engineering for Excel provides tools to help address the most > common marketing problems Segmentation, Targeting and Positioning (STP); New > Product Decision; Sales Forecasting; Advertising and Communication > Decisions; Salesforce and Channel Decisions; Pricing; and Sales Promotion > Decisions. Users will be able to harness the power of world-class analytics > to run quick and simple analyses... all from within Excel. I'll buy the license for the software (40 bucks for 6 months)... but can anyone tell me if there are R packages that do similar stuff?


title: Sense of correlation coefficient matrix of particle filter's parameters

text: I am using a particle filter to estimate the parameters($\Phi_{n\times1}$) of a non-linear model. Say my input (observations) is $t=1:k$, I will have a vector of length $k$ for each of the parameter estimates. Thus, at the end of the particle filter, I will end up with a parameter estimate matrix ($\Theta$) of $k\times n$ size, where the $k$ rows are the estimates at each time a new observation is made with the $k^{th}$ observation as my final estimate. _My question_ : Does it make sense to find the correlation coefficient matrix ($C_{n \times n}$) of this parameter estimate matrix $\Theta$? If so, what will this represent? \--AM _edit_ :   1. The non-linear model is a physiological signal model with 4 states and 7 unknown parameters. These 7 parameters are treated as state variables, thus yielding a 11 state vector which needs to be estimated. This is done in a particle filter framework.   2. Yes, the final estimates are the only ones of interest but I came across a report where the correlation of the estimates of the particles in the filtering timeline (as I mentioned above) was measured. Seeing how the estimates evolve tells a great deal about the filtering process. So I am interesting in knowing what (if any) does the correlation matrix of this evolution tell.


title: LASSO assumptions

text: In a LASSO regression scenario where $y= X \beta + \epsilon$, and the LASSO estimates are given by the following optimization problem $ \min_\beta ||y - X \beta|| + \tau||\beta||_1$ Are there any distributional assumptions regarding the $\epsilon$? In an OLS scenario, one would expect that the $\epsilon$ are independent and normally distributed. Does it make any sense to analyze the residuals in a LASSO regression? I know that the LASSO estimate can be obtained as the posterior mode under independent double-exponential priors for the $\beta_j$. But I haven't found any standard "assumption checking phase". Thanks in advance (:


title: How to accurately quantify forecast uncertainty in a special case of robust linear regression?

text: If I'm using OLS linear regression, and I want to know the uncertainty of my forecasts I can quantify it using residuals (MSE, median absolute deviation, etc). But if I'm using robust linear regression that down-weighs some observations as outliers and I don't have access to those weights, how can I accurately quantity the forecast uncertainty?


title: How to carry out repeated measures ANOVA?

text: I have the following data: IV = 2 levels (Experimental Group + Control Group). DV = I'm using signal detection variables i.e. (Hits, CRs, FAs & Misses) Both groups were tested twice so therefore I have 2 sets of DVs for each group (i.e. Hits 1 & 2, CRs 1 & 2) and so on. I want to compare between-group differences as well as within subject differences on both tests. How do I carry a repeated measures analysis? Thanks!


title: Bagging with oversampling for rare event predictive models

text: Does anyone know if the following has been described and (either way) if it sounds like a plausible method for learning a predictive model with a very unbalanced target variable? Often in CRM applications of data mining, we will seek a model where the positive event (success) is very rare relative to the majority (negative class). For example, I may have 500,000 instances where only 0.1% are of the positive class of interest (e.g. the customer bought). So, in order to create a predictive model, one method is to sample the data whereby you keep all positive class instances and only a sample of the negative class instances so that the ratio of positive to negative class is closer to 1 (maybe 25% to 75% positive to negative). Over sampling, undersampling, SMOTE etc are all methods in the literature. What I am curious about is combining the basic sampling strategy above but with bagging of the negative class.Something simply like:   * Keep all positive class instances (e.g. 1,000)    * Sample the negative classe instances in order to create a balanced sample (e.g. 1,000).   * Fit the model    * Repeat Anyone hear of doing this before? The trouble it seems without bagging is that sampling only 1,000 instances of the negative class when there are 500,000 is that the predictor space will be sparse and you may well not have a representation of possible predictor values/patterns. Bagging seems to help this. I looked at rpart and nothing "breaks" when one of the samples does not have all the values for a predictor (does not break when then predicting instances with those predictor values:               library(rpart)     tree<-rpart(skips ~ PadType,data=solder[solder$PadType !='D6',], method="anova")     predict(tree,newdata=subset(solder,PadType =='D6'))      Any thoughts? **UPDATE:** I took a real world data set (marketing direct mail response data) and randomly partitioned it into training and validation. There are 618 predictors and 1 binary target (very rare).               Training:     Total Cases: 167,923     Cases with Y=1: 521          Validation:     Total Cases: 141,755     Cases with Y=1: 410      I took all positive examples (521) from the training set and a random sample of negative examples of the same size for a balanced sample. I fit an rpart tree:               models[[length(models)+1]]<-rpart(Y~.,data=trainSample,method="class")      I repeated this process 100 times. Then predicted the probability of Y=1 on the cases of the validation sample for each of these 100 models. I simply averaged the 100 probabilities for a final estimate. I deciled the probabilities on the validation set and in each decile calculated the percentage of cases where Y=1 (the traditional method for estimating the ranking ability of the model).               Result$decile<-as.numeric(cut(Result[,"Score"],breaks=10,labels=1:10))       Here is the performance: ![enter image description here](http://i.stack.imgur.com/zMI5M.jpg) To see how this compared to no bagging, I predicted the validation sample with the first sample only (all positive cases and a random sample of the same size). Clearly, the data sampled was too sparse or overfit to be effective on the hold out validation sample. Suggesting the efficacy of the bagging routine when there is a rare event and large n and p. ![enter image description here](http://i.stack.imgur.com/X9v97.jpg)


title: Likelihood ratio test in R

text: I'm not sure if I'm asking something stupid or off topic here, but I can't think where can I ask this question. suppose I am going to do a univariate logistic regression on several independent variables, like this:               mod.a <- glm(x ~ a, data=z, family=binominal("logistic"))     mod.b <- glm(x ~ b, data=z, family=binominal("logistic"))      I did a model comparison (likelihood ratio test) to see if the model is better than the null model by this command               1-pchisq(mod.a$null.deviance-mod.a$deviance, mod.a$df.null-mod.a$df.residual)      Then I built another model with all variables in it               mod.c <- glm(x ~ a+b, data=z, family=binomial("logistic"))      In order to see if the variable is statistically significant in the multivariate model, I used the `lrtest` command from `epicalc`               lrtest(mod.c,mod.a) ### see if variable b is statistically significant after adjustment of a     lrtest(mod.c,mod.b) ### see if variable a is statistically significant after adjustment of b      I wonder if the `pchisq` method and the `lrtest` method are equivalent for doing loglikelihood test? As I dunno how to use `lrtest` for univate logistic model. (Please kindly let me know if I asked the wrong question here, thanks!)


title: How to model the probability that a tweet is spam using normal distribution?

text: I am trying to process a set of tweets, by using its length to determine the probability of spam. The data looks like this, given 2000 tweets, n is the tweet length (number of words)               when n = 1, relevant tweets = 0     when n = 2, relevant tweets = 10     when n = 3, relevant tweets = 20     when n = 4, relevant tweets = 100     when n = 5, relevant tweets = 200,     ...     ...     when n = 9, relevant tweets = 10     when n = 10, relevant tweets = 0.      My question is, how do I build a normal distribution based on this set of data, and calculate the probability of a tweet being spam, given n? Thanks very much. This is of my own research work and it is not homework. Sorry if the questions are not asked in proper statistical way as I have limited knowledge on probability and normal distribution.


title: How to distinguish these two linear regressions?

text: I have two linear regressions. The linear coefficients of both of them are very close to $1$. The first plot seems reasonable linear regression, while the second one is tricky since if the bottom left part or the top right part of the data (they are quite random) are fitted, we can not get a close-to-1 coefficient. The question is how to distinguish the two regressions after I get their very similar regression coefficient values? Instead of by looking at the plots, a statistic is preferred to describe the difference of the two regressions. Any ideas? ![enter image description here](http://i.stack.imgur.com/nqFjm.png) ![enter image description here](http://i.stack.imgur.com/TfdGK.png)


title: Analysing a rating question

text: In a survey done to find the characteristics consumers prefer while choosing a retail store. Respondents were given 9 options and to scale them from 1-10 one for the most preferred value. Now we have responses in which respondents have given same ratings to multiple characters. Example given:               price             1     quality           1     variety,          3     customer service  1     store design      9      Like wise mentioned above. How do we analyse this?


title: Plotting summary statistics with mean, sd, min and max?

text: I am from an economics background and usually in the discipline the summary statistics of the variables are reported in a table. However, I wish to plot them. I could modify a box plot to allow it to display the mean,sd and min and maximum but I don't wish to do so as box plots are traditionally used to display medians and Q1 and Q3. All my variables have different scales. It would be great if somebody could suggest a meaningful way by which I could plot these summary statistics? I can work with R or Stata Many Thanks


title: When wouldn't I use LASSO for model selection?

text: Assume that you need to build a linear model to make predictions for new observations, and that there is uncertainty about which subset of variables should be included in the model. You are only interested in making predictions, there is no theory available to guide variable selection, and you aren't interested in drawing causal inferences. In this scenario, when _wouldn't_ you use LASSO? Everything I read makes it sound like the perfect tool in this situation. However, when I have a nice new hammer in my hand, all problems look like nails, so I would like some examples of reasonable scenarios where variable selection and avoiding overfit are important, but when LASSO is either entirely inappropriate, or unlikely to be the best option. There are other threads on why you should use LASSO; I want to know why I shouldn't.


title: Co-occurrence Poisson binomial distribution

text: I'm using a Poisson binomial distribution for co-occurrence testing. I've 2 variables: $X \sim Bin(p_x,n)$ and $Y \sim Bin(p_y,n)$, where $p_x$ and $p_y$ are the probabilities, and both variables have the same number of data $n$. Then, let $Z_{\rm th} = X+Y \sim PoisBin(p_x,p_y,n)$. From this theoretical variable I have 3 probabilities for values 0, 1 and 2. I also have the probability of the values according to the observation variable $Z_{\rm obs} = \sum x_i + y_i$. How can I compare probability value of 2 in $Z_{\rm obs}$ and in $Z_{\rm th}$?


title: Expected Value of cumulative distribution function

text: Let $\varepsilon$ be a Gaussian distributed random variable with mean $\mu_0$ and standard deviation $\sigma_0$. Is it possible to compute/approximate the expected value $$ \begin{eqnarray} & &\mathbb{E}\left[\Phi\left(\frac{\varepsilon-\mu_1}{\sigma_1}\right)\mid\varepsilon<c\right]=...\\\ & &... =\int_{-\infty}^{c}d\varepsilon \frac{1}{\sqrt{2\,\pi\,\sigma_0^2}}\,\exp\left(-\frac{\left(\varepsilon-\mu_0\right)^2}{2\,\sigma_0^2}\right)\,\int_{-\infty}^{\frac{\varepsilon-\mu_1}{\sigma_1}}dt\,\frac{1}{\sqrt{2\,\pi}}\,\exp\left(-\frac{t^2}{2}\right) ? \end{eqnarray} $$ Note that while the integral $$\int \exp\left(-x^2\right)\,\textrm{erf}\left(x\right)\,dx$$ exists (check on http://integrals.wolfram.com/index.jsp), any integral of the kind $$ \int \exp\left(-(x-a)^2/b\right)\,\textrm{erf}\left((x-c)/d\right)\,dx $$ does not exist (at least according to Wolfram integrals).


title: Why do parametric tests need normal distribution of sample means?

text: Every parametric test has the assumption that the sample means are following a normal distribution. This is the case if the sample itself is normal distributed or if approximately if the sample size is big enough. but why do they have that assumption? I thought the test statistic is used to make a decision and they have their own distribution (F, t, etc).


title: How to identify the distribution of data?

text: I am trying to pick up and learn some statistics. I am not good at it, so please forgive me if I'm asking something naive. Suppose you have a set of data values, say 10k of them, and the format looks like: 1.256, 9.264, 4.483, 460.357.... How to identify the distribution of it? Thank you very much.


title: Count event using Poisson random variable

text: An event occurs with a mean of 15 times per hour and a standard deviation of 10 times/h. From the population a sample of 100 individuals is drawn. What is the probability that the sample mean will be between 14 and 16? I want to solve this problem assuming that the event follows a Poisson distribution (I can get the right answer when assuming a normal distribution). My R input is:               > 1-ppois(16, lambda=15, lower.tail=FALSE)-ppois(14, lambda=15)     [1] 0.1984695      The first `ppois` call gives the right tail area and the second call the lower tail area, which I both subtract from 1. What is wrong with my calculation?


title: Confidence intervals on 'cross plots'

text: I have a large set of data comprising triples of points (t, y1, y2). At each t, I have a distribution of y1 and y2 points, each of which is approximately normal. I would like to create a plot of y1 versus y2 (mean values) and add a 90% or 95% confidence interval envelope. I would appreciate any advice on how to calculate the confidence intervals in this case such that I can subsequently construct a plot. Thank you. **Update and clarification:** In my field, if y1 is a function of x and y2 is a function of x, a plot of y1 vs y2 over the domain of x is called a 'cross plot'. Additional information:   * In my application, t is the independent variable and y1 and y2 are dependent variables.   * There is a complex (differential equation based) predictive model that produces many values of y1 and y2 for every value of t.    * There is no explicit model relating y1 and y2.    * I can already create separate confidence interval plots (t vs y1, t vs y2).    * For each value of t, I can estimate the mean and standard deviation for y1 and y2.  My aim: For each value of t, say ti, I would like to plot mean(y1i) vs mean(y2i) as a point on my y1 vs y2 plot and then include the confidence interval envelope.


title: Research methodology suggestions on Electric utility services

text: I have a data set of 30 million consumers of an electric utility service. \- consumption \- billing \- smart meter readings sampled every 5 seconds \- market share data area wise Any suggestions on how to compute weights for underlying factors that can predict customer satisfaction? This is an initial plan I have, but need help understanding its merits and limitations.   1. Identify a pool of "satisfied customers" based on business rules. Someone who has been using the service for a period of 2 years and has an increasing consumption trend.   2. Identify a pool of dissatisfied customers, ones who ended the relationship prematurely and may have exhibited late payment of fees.   3. Add cross sectional variables to account for difference in customers Use logistic regression to identify significant factors. How do I better this or are there serious flaws in this approach?


title: Determining the best correlated time series

text: Before asking, I read similar questions, but none of them lead to satisfying answers for my specific interest. I want to homogenize a climate time series of precipitation of the Dominican Republic over 64 years (1940-2003). For that, it is really important to select a reference series among a group of candidates. Let's say `sjo` is the base series, for which I want to find a good reference series; `bani`, `plc` and `ra` are reference candidates, because they are close to `sjo`. In the following map, the red point is the base station, and the green ones are the reference candidates: ![](http://i.stack.imgur.com/AxrWL.jpg) I performed three correlation analyses (done in R, function `cor()`), considering these monthly variables: raw precipitation value, normalized difference, and transformed values with Box-Cox. Those variables correspond, respectively, to fields that begin with `p`, `dian` and `pnorm`. Normalized difference comes from the first difference series method (FDM), which was proposed by Peterson, consisting of: $[Pm_t - Pm_{t-1}] / [Pm_t + Pm_{t-1}]$, where $Pm_t$ is the precipitation value for the month $m$, and $Pm_{t-1}$ is the precipitation for the same month 1 year before. I followed Peterson and colleagues' (1998) remark, which says that FDM applied to precipitation might work better using normalized difference. As can be seen in page 1 of this PDF file, correlation was calculated for the whole time series (1940-2003). For raw precipitation and Box-Cox transformed values, `bani` is the best correlated with `sjo` (yellow background cells shows the maximum correlation index). Notice that for raw precipitation, `bani` is significantly more correlated than others. For normalized difference, `ra` is only a bit more correlated than the rest. However, each candidate station has a statistically significant correlation index with `sjo` at the $\alpha=.05$ significance level, suggesting ANY of them could be used as a reference series. This is a bit confusing, so I was unsatisfied and decided to make a more detailed analysis, splitting the series in 5-year period intervals, and evaluating correlation between series for the same 3 variables: raw precipitation, normalized difference and Box-Cox transformed. Tables from page 2 to 8 in the PDF show the results of these partial correlations; the last page summarizes the times each station has had the maximum correlation value for each variable. As can be seen, `bani` is the most frequently correlated value for the 3 variables analyzed (in all cases, more than 7 times of the twelve 5-year periods analyzed). With these results, I think that `bani` is the best candidate as a reference series of `sjo`, but I'm not sure about it. Is the five-year period analysis OK? Should I perform some other analysis?


title: equivalent ways to write hypotheses

text: What are some additional ways to write the null and alternative hypotheses below (i.e., mathematically equivalent)? Null hypothesis: H0: µ1 = µ2 the two population means are equal Alternative hypothesis: H1: µ1 ≠ µ2 the two population means are not equal OR Null hypothesis: H0: µ1 - µ2 = 0 the difference between the two population means is 0 Alterative hypothesis: H1: µ1 - µ2 ≠ 0 the difference between the two population means is not 0 OR ...?


title: Checking Independence for functions in R

text: Given $$f_{(U_1,U_2)}(u_1,u_2)=\begin{cases} 1/2& -u_1<u_2<u_1 \text{ and } u_1 - 2 < u_2 < 2 - u_1 \text{ and } 0 < u_1 <2\\\ 0& \text{otherwise}\end{cases}$$ I found that $$f_{(U_1)}(u_1) = \begin{cases} u_1 & 0<u_1<1\\\2-u_1 &1<u_1<2 \\\0 & \text{otherwise}\end{cases}$$ And $$f_{(U_2)}(u_2) = \begin{cases} 1+u_2 & -1<u_2<0\\\1-u_2 &0<u_2<1 \\\0 & \text{otherwise}\end{cases}$$ I want to justify whether $U_1$ and $U_2$ are independent or not by using R.. I am wondering what is the command is like...Could anyone please give me some directions as for how to check independence for functions in R?


title: How to generate a probability of x given a value of y and a distribution?

text: I have a bivariate data set where x and y values, both continuous variables, correlate well (visually speaking) when x and y are small. As the values of x and y increase, the correlation decreases to an eventual random state when x and y are large. My question is, given a y value, I would like to know the probability of the value of a given x. Should I generate a probability distribution for this data? I would very much appreciate any practical advice for Python in specific.


title: Nonparametric chi-square test

text: I have two independent populations (up and down), each population contains three categories 1,0,-1. For example:                    up = 1,1,1,1,0,0,1,1,-1,-1,1          down = 1,-1,1,-1,0,1,-1      These populations are not the same size as each other and are much larger than the example above (around 3,000 datapoints) I want to know how similar these two distributions are to one another (distribution meaning density of -1,0,1) When I graph these points the data is not normal so I am unable to run a chi- square test.


title: Nested or Split plot

text: I have a design with four factors Frame (positive negative neutral) Graphic ( graph no graph) Color (red grey) nested in graphic Position of color (nested in color, which is nested in graphic) I am having difficulty understanding if this is an incomplete nested design or a split plot design. Any thoughts would be great, thanks.


title: How should I present linear regression output in a table format?

text: Is there any convention regarding the minimum information that should be included in a table showing a linear regression output? The reason for this question is that I want to make the output as easy to understand as possible. I have graphed the main findings but I want to include the other outputs in a detailed table. This is in a general social science context (i.e. a presentation; not journal publication) and my audience is not trained in statistics, so I want to keep everything as simple as possible (while still retaining the critical information).


title: Percentile Ranks in Matlab

text: I have 2 vectors of raw values. I need to do a percentile rank correlation and an ordinal rank correlation between them in Matlab.               Vec1 = [25    28    29    29    31    32    33    33    33    35    37] ;      I'm using `corr` function in Matlab with `type=Spearman` for the ordinal part. I'm confused about the percentile rank part. Until now I was using `tiedranks*100/length` and would feed the output to `corr` function with `type=Pearson` (courtesy of a forum answer). However, the correct way is at http://www.psychstat.missouristate.edu/introbook/sbk14m.htm.               Website Ranks    = [4.6  13.6  27.3  27.3  40.9  50    68.2  68.2  68.2  86.4  95.4] ;     TiedRanks*100/len= [9.09 18.18 31.81 31.81 45.45 54.54 72.72 72.72 72.72 90.90 100];      Please let me know which is correct? Can I use Matlab's `prctile` function here? Thanks.


title: Rao-Cramér inequality

text: Suppose $p({\bf x},\theta), \theta \in {\bf \Theta} \subset {\bf \mathbb R^k}$ is a (regular) family of distributions. The score function is defined as $s({\bf x},\theta)=\frac{\partial}{\partial \theta}log(p({\bf x},\theta))=grad(p({\bf x},\theta))$, the Fisher information matrix as ${\bf I}(\theta)=var({\bf s})$ and $E_\theta=({\bf t}({\bf X}))=\theta$. Deduce the Rao-Cramér inequality from the following expression: $$corr^2({\bf a}^T{\bf t}, {\bf c}^T{\bf s})=\frac{({\bf a}^T{\bf c})^2}{{\bf a}^Tvar(t){\bf a}\cdot{\bf c}^Tvar(s){\bf c}}$$ where ${\bf a}$ and ${\bf c}$ are fixed ${\bf k}$-dimensional vectors and ${\bf s}$ and ${\bf t}$ are random variables. The above stated expression comes from the fact that $E({\bf s})={\bf 0}$ and $E({\bf s}{\bf t}^T)={\bf I}$. I proved this equation, but I don't know how to maximize it over vector ${\bf c}$ and deduce the Rao-Cramér inequality to prove that the correlation coefficient is always less or equal 1.


title: Size of cluster as a predictor in a glmm

text: I have pupils clustered within schools, with a binary outcome (pupil experiences bullying or does not). So far all predictors are pupil-level data, and I run a the same model (same predictors) on each year of data. My teacher was talking about adding a school-level predictor: the size of the school, which will simply be the count of pupils in each school, that is the number of observations in each school. Is there any problem with doing this ? I tried to look for some references on this but I didn't find anything useful so far, and looking in some textbooks it doesn't seem like a normal thing to do ? Edit: I am interested in how the residual variance at the school level changes over time.


title: How can I use standard deviation/SEM to assess the appropriateness of replacing missing values with the mean?

text: I have a data set consisting of reaction time (RT) measurements. These will be used to calculate the duration of experimental trials in an MRI study. In each block (experimental condition) there are 10 trials. Now, due to various issues, there are some missing RTs. Due to the nature of MRI analysis, I need values for all 10 trials per block. Where more than 5 values are missing I will probably discard it from the analysis completely, but where only 1 or 2 are missing I plan on using the mean RT value for that condition in place of the missing value. However, I want to make sure that this is a principled decision to make, as the RT values in some conditions, by the same participant, can be quite variable. How can I use the standard deviation or standard error of the mean to ensure that it is 'fair' to use the mean in place of the missing value? For example, see the data below. Block 1 - Missing values: 2; Mean: 740; SD: 519; SEM: 196. Block 2 - Missing values: 1; Mean: 2245; SD: 292; SEM: 97. I'm trying to figure out an honest, consistent way of deciding whether the decision to replace the missing value with the mean is sound. Where it is not, I would rather leave the block out of the analysis than skew the data. Any advice? I hope this makes sense.


title: Predictions using glmnet in R

text: I am trying to model some data using glmnet package in R. Let say I have the following data               training_x <- data.frame(variable1=c(1,2,3,2,3),variable2=c(1,2,3,4,5))     y <- c(1,2,3,4,5)      (this is a simplification; my data is much more complicated). Then I used the following code to create the glmnet model.               x <- as.matrix(training_x)          library(glmnet)     GLMnet_model_1 <- glmnet(x ,y             ,family = "gaussian"             ,alpha = 0.755             ,nlambda = 1000             ,standardize = FALSE             ,maxit=100000)      I am using standardize=FALSE because my real life data is already standardized. Then I want to make prediction over a new set of data. Let's say my new data is               newdata <- data.frame(variable1=c(2,2,1,3),variable2=c(6,2,1,3))          newx <- as.matrix(newdata)     results <-predict(object = GLMnet_model_1, newx, type="response")      I would expect results to contain 4 elements (predictions of the newdata), but instead it gives me a 4x398 matrix... What am I doing wrong?


title: Weka: how do I score an arff file using a NaiveBayes model?

text: I trained a model with               java -Xmx6g -cp /usr/share/java/weka.jar weka.classifiers.bayes.NaiveBayes -t train.arff -d foo.nb      and I want to score a test file with it. I tried               java -Xmx6g -cp /usr/share/java/weka.jar weka.classifiers.bayes.NaiveBayes -T test.arff -l foo.nb -p N      which produces no files but writes some text pasts of which can be construed to be scores to stdout (and is unbelievable slow):                inst#     actual  predicted error prediction ()          1        1:0        1:0       0.549           2        1:0        1:0       0.55           3        1:0        1:0       0.531           4        1:0        1:0       0.515           5        2:1        1:0   +   0.552           6        1:0        1:0       0.532           7        1:0        2:1   +   0.519       If I read this correctly, the last column is the score and 3rd is the prediction based on it. Why does score 0.55 correspond to 0 (instance#2) but a smaller score 0.519 correspond to 1 (instance#7)? Where is the output documented? Is there a way to produce the csv score file? Thanks!


I have data on 26 participants (13 from computing and remaining 13 from non computing) who have participated in my research. Each participant is treated with a lab module (Hands on Robotics Session). Now each participant will be evaluated using a rubric on scale of 1 to 4. This experiment has both pre and post test. **Is data on a 1 to 4 scale interval or ordinal?**


title: Infinite fourth moment and maximum entropy

text: Alright, I expect this is a silly question, but I don't actually know, so. Suppose there is some random variable that's distributed on the reals, and all I know about the distribution is its mean $\mu$ and its variance $\sigma^2$. Maximum Entropy would have me say that I should choose the Normal distribution as the one that best fit my knowledge about it. However, suppose I have an utility function about this variable that goes $U(x)=x^4$. Then, my expected utility for drawing from that variable's distribution is $EU(x) = \int\limits_{-\infty}^{+\infty}\\!\\!x^4\ p(x)\ \mathrm{d}x$. If I choose to use MaxEnt, then that's just $3\sigma^4$. However, if the "true" distribution actually followed by that random variable is, say, the Student's t with $\nu\leq 4$, then my Expected Utility would diverge to infinity. If I treat the distribution as if I don't know what it is, then I'd have that $p(x)=\int\limits_{D\in\mathcal D}\\!\\!p(x|D)\ p(D)\ \mathrm{d}D$ where $\mathcal D$ is the space of all possible distributions. In that case, my expected utility is $EU(x) = \int\limits_{-\infty}^{+\infty}\\!\int\limits_{D\in\mathcal D}\\!\\!x^4\ p(x|D)\ p(D)\ \mathrm{d}D\ \mathrm{d}x$ which I expect is _not_ the same as the one I get by using the Normal except for very few priors over $\mathcal D$. My questions then are: how do I deal with _that_? What prior distribution over $\mathcal D$ gives me the same result as MaxEnt would? In what sense is the Normal my best guess for my state of ignorance? How do I guarantee that this doesn't diverge? How does having the Normal as the "representative" of my subjective uncertainty take other possible distributions such as the Student's t into account? Or, rephrasing this: I don't feel like the Gaussian actually represents my state-of-knowledge here, because nowhere in it is it represented that it's possible that the distribution has infinite fourth moment. In what sense is the Gaussian the best description of my state-of-knowledge? Is it, really? If not, what is?


title: Perl advantages in data management/statistical analysis?

text: I've frequently read that knowledge of Perl or Python complements SAS or R, and seen that Perl is a desirable language in job advertisements. Can anyone clarify this a bit more? Is Perl limited to text mining & database updates or does it have other applications in the statistics world? Thanks!


title: Multiplicity in sample size calculation for stratified estimation problem

text: We want to determine the public opinion about a recently administered intervention in a health care network. There's a brief questionnaire of 10 Y/N questions in which we'd like to estimate proportions within a specific margin of error (+/- 3%). Based on historical data, we feel we can reliably estimate response probabilities within specific groups. We're interested in stratified sampling in which we will use Horvitz-Thompson weighting to obtain the estimate of the population proportions of positive responses to the questions. My question is: how does one account for the multiplicity in calculation of the margin of error for individual questions? Intuiviely, estimation leads to the same issues in multiplicity that inference does. We might expect that this specific survey obtained, for least one question, an estimate of the population proportion and its estimated standard error which is inconsistent with that of the actual long-run average one would obtain from an infinite number of independent replications of the experiment. Hence, you would like a more conservative estimate of the required sample size to achieve a "95% confidence level" uniformly across all questions, or that there is a "grand" margin of error of 5% or so. Is this consistent with any existing survey methodology? Is multiple testing strictly considered only in the world of formal inference and p-values? How exactly would one go about calculating sample sizes or power with any existing multiplicity adjustment for estimation?


title: Regression analysis when the covariables is a sample from a population of potential variables

text: This question comes from trying to analyze my recent exam (exam I have given and corrected) statistically. I have a list of questions (20 in total) and each question is given a score from 0 to five, so total possible score is 100. Then I define pass/not pass via some cutoff (I give notes on the A-B-C-D-E-F scale, but here we only discuss F/notF). A logistic regression of pass on the 20 score variables obviously give a perfect fit, and the usual assumptions behind logistic regression are obviously not fulfilled (nothing special here about logistic regression , the question really are about regression modelling, not any specific kind of regression model). But in reality what we have sampled here is the variables, not the students! The interesting analysis is about the influence on notes of specific question types, for the actual students we have, which are not in any sense a sample. But the variables is a sample from a huge population of potential questions! So what could we do in direction of a formal analysis? Some form of cross- validation seems natural, but we should really subsample **variables** , not objects ( **students** ). Any good ideas, or good references? I tried google scholar, but couldn't find anything. EDIT as answer to comments. I am mostly after seeing which question have most influence on the pass/fail decision. My exam is a calculus exam, with questions as (simple example) "Find the derivative of $f(x) = e^x \cos(x^2+3) $." (there are also longer questions, but divided inparts) counting this parts I have 20 sub-questions. This questions might be seen as a sample (at least conceptually) from a much larger set of potential questions. My students are really a fixed set, and I want to evaluate the exam as a test for this specific set of students. (so resampling by bootstrapping students do not seem natural). I sum the scores for each question and then decide the note by some cutoff. Now, if I use the scores on the 20 questions as variables, and fit a logistic regression model as $$ P(\text{fail}) = \frac1{1+\exp(-\beta_0 -\beta_1 Q_1 - \dots -\beta_20 Q_{20})} $$ I do obviously get a perfect fit! since the pass/fail decision was taken on the basis of those 20 variables, in a linear fashion. But still, the decision is not perfect, obviosly, there are error sources. One source of error is things like the day form of the student. A very different kind of error source is the match between the questions and the knowledge of the student, in that with a different sample of 20 questions, the student might have got a different result. So the question is about how to model this situation (and then analyze according to that model).


title: Updating "support" or "range" in transformation variable

text: I submitted the final in a theoretical stats class and had trouble with several problems, but part of one is bugging me because it seems like it should be simple. Lets say we have two random independent variables $X$ and $Y$, each with their respective pdf $f(x)$ and $f(y)$. The support for each is: $-\infty < X < \infty$, and $-\infty < Y < \infty$. How do I determine what the support should be for a transformation to a new random variable $Z$, that utilizes both $X$ and $Y$? Several examples might be:   1. $Z = X-Y$   2. $Z = X\times Y$   3. $Z = X/Y$ What is the range of $Z$?   * $-\infty < Z < \infty$ ?   * $0 < Z < \infty$ ? I think the lower bound would be 0 for case #2, but how can you determine this for 1, 2, and/or 3 given above? I have Casella and Berger's (2002) _Statistical Inference_ text for the class, but I don't understand how to "update the support" from this book. I've had virtually 0 luck Googling this. Thanks for any suggestions either directly explaining how, or for a source.


title: Using subset to plot certain observations

text: I am trying to use the subset argument to restrict which observations are plotted. The doc for xyplot says the subset argument should be "An expression that evaluates to a logical or integer indexing vector." Here are two examples trying both, but I can not get either to work.               # Generate some data with obvious outliers, want to plot without these observations     library('xts')     a = as.xts(ts(rnorm(20), start=c(1990,1), freq=4))     b = as.xts(ts(rnorm(20), start=c(1990,1), freq=4))     a["1991"] = 20 + rnorm(4) # 1991Q2 - 1992Q1, observations 6-9     b["1991"] = 20 + rnorm(4) # 1991Q2 - 1992Q1, observations 6-9          library('lattice')     xyplot(a~b) # can see outliers     xyplot(a~b, subset=(time(a) %in% c("1991 Q2", "1991 Q3")) ) # logical vector is all false?     xyplot(a~b, subset=(1:20)[-c(6:9)] ) # produces a plot but changes the actual data plotted?


title: are "non-parametric" methods well-defined?

text: I have trouble understanding the delineation of "non-parametric" methods. It seems to me most non-parametric methods are parametric, just on a different space than the "standard ones". Graph-based inferences are often described as non-parametric, but graphs are fully described by the set of adjacency matrix parameter values. Likewise rankings can be considered a transformation of values onto a different parameter space. Non-parametric seems to get tossed around whenever something is parametric in a high dimensional way on a space we don't typically use. If the definition is a question of the use of a distribution, coming from a Bayesian background, I guess my intuition is often that usually when we don't put a distribution on something we're often neglecting to capture its uncertainty, rather than something intrinsic to the method. Hence there seem to be bayesian formulations of "non-parametric" methods which involve parameter distributions, suggesting that a lack of a distribution isn't inherent to the representation or method, but more of a modeling decision.


title: Recurrent neural networks in R

text: I've heard a bit about using neural networks to forecast time series, specifically recurrent neural networks. I was wondering, is there a recurrent neural network package for R? I can't seem to find one on CRAN. The closest I've come is the nnetTs function in the tsDyn package, but that just calls the nnet function from the nnet package. There's nothing special or "reccurant" about it.


title: Usual regression vs. regression when variables are differenced

text: I am just trying to understand what the relationship is between a normal multiple/simple regression vs. multiple/simple regression when the variables are differenced. For example, I am analyzing the relationship between deposit balance ($Y_T$) vs. market rates ($R_T$) If I run a simple linear regression, the correlation is negative and pretty significant (around -.74) However, if I take the log and difference of the dependent variable and the difference of the independent variable, so my equation is now $d\, \ln(Y_T)$ being regressed with $d\, R(T)$, my correlations and R^2's are not significant at all ($R^2 = .004$). I was just wondering whether this low $R^2$ even mean anything? Does it mean my model isn't a good fit, or do I ignore the $R^2$ when I am looking at differenced data? I know from the data there is a significant correlation between the original two variables, yet for my model I need to look at the variables differenced, so just wondering how to go about this.


title: Calculate frequency of 1D time series using autoregressive model parameters

text: I'm modeling some periodic data with a second-order autoregressive model, as follows: $$ x_3 = a_{1}x_1 + a_{2}x_2 $$ $$ x_4 = a_{1}x_2 + a_{2}x_3 $$ $$ ... $$ $$ x_n = a_{1}x_{n-2} + a_{2}x_{n-1} $$ I'm modeling the raw data here (imagine $x_i$ as a scalar pixel value, say, for $n$ time points) so I'm not worried about error terms, at least not yet. I can set this up as a system of equations and solve for the parameters $a_1$ and $a_2$. **Is there any way to determine the frequency of this 1D datasets using the calculated AR parameters, since the system intrinsically models harmonic oscillators?** I'm actually trying to avoid using fourier transforms, if possible. Thanks!


title: How to write a multinomial logit model as a generalized linear model?

text: Logistic, linear and probit regression can be written in terms of Generalized Linear Model (GLM). Can the multinomial logit regression be written in terms of GLM as well?


title: Combining variables

text: I'd like to combine several variables into one variable. Here is some context: Let's say I have two variables `Red.Beads` and `Blue.Beads` and each observation of each variable is a number representing the number of beads for that observation. Say each red bead has a mass of 5 grams, and each blue bead has a mass of 1 gram. I can a single variable representing overall bead mass for observation by doing the following: 5*Number.Red.Beadsi \+ 1*Number.Blue.Beadsi. But, if I do not have access to the mass of each type of bead, is there a meaningful way to combine the variables to get some proxy for mass of beads present. My initial thought was to standardize each variable (i.e. compute z-score), and then combine (i.e. z-score.Redi \+ z-score.Bluei). I recognize that this approach suffers from some problems, and I wonder what other options exists for combining the variables? In my situation, I want to combine different numbers of variables from 2 or 3 up to about 40 variables. Thanks in advance for any thoughts. Best Regards.


title: How to compare 37 dichotomous variables (lab tests) to 31 dichotomous variables (diagnosis)

text: I've constructed a database based on chart review of patients. Patients had a variety of tests performed (37 different types, which were either positive, negative, or not performed aka 1/0/blank - I presume this is consistent with a dichotomous/nominal variable). A total of 31 possible diagnoses could have resulted, and I have encoded these as dichotomous variables (each diagnosis as 1/0 aka yes/no). (There were different tests that assessed for the same condition, hence 37>31). Some collaborators have suggested creating a correlation matrix, which would doubly serve as input for data visualization. I am struggling with SPSS as to how best achieve this. In essence, I am unsure if Analyze>Correlate>Bivariate is the right path to take, or alternatively if a separate approach should be undertaken. For instance, certain tests were "clustered" together and I am not sure how to account for this occurrence. Further, I am not sure how to adjust for patient-level characteristics (which are not dichotomous variables, such as age, gender, and so on). Thanks for your help in advance


title: Determining the number of levels in a Repeated Measures ANOVA for multivariate or univariate within subjects effects

text: Maxwell and Delaney (1990) wrote that the multivariate approach should not be used for a repeated measures ANOVA for within subjects effects if n is less than "a + 10" (a is the number of levels for repeated measures). My data has 7 time levels and 2 levels for location and 2 levels for condition (control vs. non-control); therefore, I have 28 unique combinations. The total number of observations is 40, which is repeated 7 times. I am wondering what "a" is actually referring to…what would "a" equal in my case, 7, 11, or 28? Is my "n" 40 or 10? Also, for a repeated measures ANOVA, I read that I am supposed to test for normality beforehand. Do I just test all the combined data for normality?…meaning the number of "established seeds" column regardless of the levels? Then, if it fails normality, then try to log transform the data and test it again? Also, univariate tests are good for small sample sizes, but is there a sample size that is too small? I have another dataset that has only 6 observations. 2 types of traps were put in 3 different locations to collect ants every 2 weeks for a year. The stats run, but are they valid?