Notes on Analysing Experiments in R - Revision history

Seth: /* Good conversation snippets from the threads */

2012-04-24T13:58:35Z

‎Good conversation snippets from the threads

Seth: /* Good conversation snippets from the threads */

2012-04-24T13:58:16Z

‎Good conversation snippets from the threads

Seth at 13:54, 24 April 2012

2012-04-24T13:54:50Z

Seth at 13:54, 24 April 2012

2012-04-24T13:54:13Z

Show changes

140.182.235.250 at 01:31, 29 April 2010

2010-04-29T01:31:12Z

Seth at 18:44, 22 April 2010

2010-04-22T18:44:51Z

140.182.227.189: /* documentation */

2010-04-22T18:43:41Z

‎documentation

Seth: /* documentation */

2010-04-19T01:46:14Z

‎documentation

Seth: /* documentation */

2010-04-19T01:45:54Z

‎documentation

Seth: /* documentation */

2010-04-19T01:45:44Z

‎documentation

@@ Line 30: / Line 30: @@
 ==Good conversation snippets from the threads==
 * There is now an anova() method for lmer() and lmer2() fits performed using method="ML".  You can compare different models and get p-values for p-value obsessed journals using this approach. [https://stat.ethz.ch/pipermail/r-help/2007-August/138747.html]
-*This is Bate's answer to people who already know the answer:
+*This is Bate's answer to people who already know the answer: """With lmer fits I recommend checking a Markov chain Monte Carlo sample from the posterior distribtuion of the parameters to determine which are signification (although this is not terribly well documented at the present time).""" [http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3565.html]
-"""With lmer fits I recommend checking a Markov chain Monte Carlo sample from the posterior distribtuion of the parameters to determine which are signification (although this is not terribly well documented at the present time).""" [http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3565.html]
+*Similarly """Try using mcmcsamp() to sample from the posterior distribution of the parameter estimates. You can calculate a p-value from that, if that is  your desire.""" [http://www.opensubscriber.com/message/r-help@stat.math.ethz.ch/5544126.html]
 *Similarly
 """Try using mcmcsamp() to sample from the posterior distribution of the parameter estimates. You can calculate a p-value from that, if that is  your desire.""" [http://www.opensubscriber.com/message/r-help@stat.math.ethz.ch/5544126.html]
 *More for people who already know the answer: [http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/736?set_blog_all=yes]
 *Here is what I've been doing (snippet below. limit: only really works with infinite data): """My general advice to those who are required to produce a p-value for a particular fixed-effects term in a mixed-effects model is to use a likelihood ratio test.  Fit the model including that term using maximum likelihood (i.e. REML = FALSE), fit it again without the term and compare the results using anova. The likelihood ratio statistic will be compared to a chi-squared distribution to get a p-value and this process is somewhat suspect when the degrees of freedom would be small. However, so many other things could be going wrong when you are fitting complex models to few observations that this may be the least of your worries."""[http://markmail.org/message/56c4ck4mmjyouqfo]

@@ Line 31: / Line 31: @@
 * There is now an anova() method for lmer() and lmer2() fits performed using method="ML".  You can compare different models and get p-values for p-value obsessed journals using this approach. [https://stat.ethz.ch/pipermail/r-help/2007-August/138747.html]
 *This is Bate's answer to people who already know the answer:
- With lmer fits I recommend checking a Markov chain Monte Carlo sample from the posterior distribtuion of the parameters to determine which are signification (although this is not terribly well documented at the present time). [http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3565.html]
+"""With lmer fits I recommend checking a Markov chain Monte Carlo sample from the posterior distribtuion of the parameters to determine which are signification (although this is not terribly well documented at the present time).""" [http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3565.html]
 *Similarly
- Try using mcmcsamp() to sample from the posterior distribution of the parameter estimates. You can calculate a p-value from that, if that is  your desire. [http://www.opensubscriber.com/message/r-help@stat.math.ethz.ch/5544126.html]
+"""Try using mcmcsamp() to sample from the posterior distribution of the parameter estimates. You can calculate a p-value from that, if that is  your desire.""" [http://www.opensubscriber.com/message/r-help@stat.math.ethz.ch/5544126.html]
 *More for people who already know the answer: [http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/736?set_blog_all=yes]
 *Here is what I've been doing (snippet below. limit: only really works with infinite data): """My general advice to those who are required to produce a p-value for a particular fixed-effects term in a mixed-effects model is to use a likelihood ratio test.  Fit the model including that term using maximum likelihood (i.e. REML = FALSE), fit it again without the term and compare the results using anova. The likelihood ratio statistic will be compared to a chi-squared distribution to get a p-value and this process is somewhat suspect when the degrees of freedom would be small. However, so many other things could be going wrong when you are fitting complex models to few observations that this may be the least of your worries."""[http://markmail.org/message/56c4ck4mmjyouqfo]

← Older revision		Revision as of 13:54, 24 April 2012
Line 1:		Line 1:
−	'''NOTE:''' ''I wrote this in 2010 and touched it up recently in 2012. I can use bootstrapping now, but I still prefer to use the very simple ~~<pre>~~anova(mymodel, mymodelminusoneparameter)~~</pre>~~. There is one especially exciting approach, completely different, still very technical, that I'm not at all facile with, called [[Weblog:Bayesian_data_analysis_in_English\|Bayesian analysis]]. ''	+	'''NOTE:''' ''I wrote this in 2010 and touched it up recently in 2012. I can use bootstrapping now, but I still prefer to use the very simple "anova(mymodel, mymodelminusoneparameter)." There is one especially exciting approach, completely different, still very technical, that I'm not at all facile with, called [[Weblog:Bayesian_data_analysis_in_English\|Bayesian analysis]]. ''
		+

	I've been inferring my knowledge of statistical analysis from the more and less patient attacks that Douglas Bates unleashes on non-statisticians who need p-values in the language R. Douglas Bates is the author of R's mixed effects package, which lets you do tricky things. It turns out that, for complex enough designs, experimentalists are at the cutting edge of statistics, and a lot isn't known. Specifically, this article may be relevant to you if you use random effects (e.g. continuous-ish covariates), or nested experimental designs (e.g. involving phrases like "within subject" and "repeated measure"), a mix of ordinal/categorical and continuous variables and other non-vanilla flavors of analysis of variance.		I've been inferring my knowledge of statistical analysis from the more and less patient attacks that Douglas Bates unleashes on non-statisticians who need p-values in the language R. Douglas Bates is the author of R's mixed effects package, which lets you do tricky things. It turns out that, for complex enough designs, experimentalists are at the cutting edge of statistics, and a lot isn't known. Specifically, this article may be relevant to you if you use random effects (e.g. continuous-ish covariates), or nested experimental designs (e.g. involving phrases like "within subject" and "repeated measure"), a mix of ordinal/categorical and continuous variables and other non-vanilla flavors of analysis of variance.

@@ Line 51: / Line 51: @@
 *All the way back:  R is a programming language for statistical analysis.  It is free in both the bird and pocketbook senses.  It is almost identical to a very expensive program called S.  It is also very good, and in some ways it may already have eclipsed  S.
 *Less far back:  A t-test is a very simple test of whether two lists of numbers are not ''really'' different.  Comparing samples gets way more complicated.  Analysis of variance (ANOVA) is a step on the way, and is the most popular manifestation of The Linear Model.  Besides coursework, the best way that I know to get from knowing nothing to understanding this post is at sportsci.org .  It has the most thorough free resource that I've found for explaining it all to dummies [http://www.sportsci.org/resource/stats/].  I read the whole thing.
-*in R, aov() works most of the time.  If not, lm() should work.  If not, things get tangled and the learning curve gets worse.  There is glm(), lme(), nlme(), and lmer() (of which glmer(), nlmer() and lmer2() are variants), all used with anova() or glht() or MCMC/bootstrapping techniques involving pvals.func(), merMCMC-class(), mcmcsamp() or your own code.  Relevant libraries are nlme, lme4, languageR, and multcomp.  This page refers mostly to the use of lmer() in the package lme4 linked in the documentation section.
+*in R, aov() works most of the time.  If not, lm() should work.  If not, things get tangled and the learning curve gets worse.  There is glm(), lme(), nlme(), and lmer() (of which glmer(), nlmer() and lmer2() are variants), all used with anova() or glht() or MCMC/bootstrapping techniques involving pvals.func(), merMCMC-class(), mcmcsamp() or your own code.  Relevant libraries are nlme, lme4, languageR, RLRsim, and multcomp.  This page refers mostly to the use of lmer() in the package lme4 linked in the documentation section.
 *Hypothesis testing gets hard after lm().  You can't just type <pre>summary(yourlm)</pre> like with lm().  You have to either sample the model and infer test statistics from the simulated results (bootstrapping) or calculate the log likelihood ratio of the model against a null model*.  You do this by writing your own test, using the one above, or just typing <pre>anova(yourlm, yourotherlm)</pre>. Any of these could not-work, depending on the complexity of the model.
 <nowiki>*</nowiki> ideally one that is the same-but-for-one-factor (opaquely, the distribution of the log of the square of the ratio of two models is approximately chi-square, and I think it errs conservative.  sameness-but-for-one makes df=1 which is the easiest to interpret in this context)

@@ Line 15: / Line 15: @@
 *to google anything about stuff in R, "r stuff" won't work. Try "r cran stuff" or "r help stuff"
 *[http://zoonek2.free.fr/UNIX/48_R/14.html useful walk through the usefulness of mixed models, examples in R]. also [http://zoonek2.free.fr/UNIX/48_R/]
-*[worldwideweb.unconventionallylonguniformresourcelocator.com/forever/wikirefs/R_lmer_Dummies_Supplement.pdf a few slides on lmer] a simple thing I made.  More to document my confusion, but maybe helpful with lmer syntax
+*[http://worldwideweb.unconventionallylonguniformresourcelocator.com/forever/wikirefs/R_lmer_Dummies_Supplement.pdf a few slides on lmer], a simple thing I made.  More to document my confusion, but maybe helpful with lmer syntax
 *http://cran.r-project.org/web/packages/lme4/vignettes/Implementation.pdf
 *http://cran.r-project.org/web/packages/lme4/vignettes/