(To the list moderator: I just subscribed to the list. Apologies for not
having done so longer before trying to post.)
Hi all,
I am currently using lmer to analyze data from an experiment with a
single fixed factor (treatment, 6 levels) and a single random factor
(block). I've been trying to follow the online guidance for estimating
p-values for parameter estimates on these and other help threads:
https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
http://pidgin.ucsd.edu/pipermail/r-lang/2007-August/000057.html
but have not yet been successful. This is my code:
lnmass <- lmer(log.mass ~ treatment + (1|block), data=exp1)
summary(lnmass)
samp <- rnorm(n=10000)
mcmcpvalue <- function(samp)
{std <- backsolve(chol(var(samp)),
cbind(0,t(samp)) - colMeans(samp),
transpose = TRUE)
sqdist <- colSums(std*std)
sum(sqdist[-1] > sqdist[1]/nrow(samp) }
markov1 <- mcmcsamp(lnmass, 10000)
HPDinterval(markov1)
mcmcpvalue(as.matrix(markov1[,1]))
mcmcpvalue(as.matrix(markov1[,2]))
mcmcpvalue(as.matrix(markov1[,3]))
mcmcpvalue(as.matrix(markov1[,4]))
mcmcpvalue(as.matrix(markov1[,5]))
mcmcpvalue(as.matrix(markov1[,6]))
The first time I tried it HPDinterval generated CIs for each treatment
level, etc. but after each of the last 6 lines of code I received the
following error messages:
Error in markov[,1] : object is not subsettable
Error in as.matrix(markov1[,1]) : error in evaluating the argument 'x'
in
selecting a method for function 'as.matrix'
Error in chol(val(samp)) : error in evaluating the argument 'x' in
selecting
a method for function 'chol'
When I came back to it today HPDinterval(markov1) also no longer worked and
gave me the following error message:
Error in UseMethod("HPDinterval"): no applicable method for
"HPDinterval"
I imagine at least some of this has to do with the fact that I don't fully
understand what exactly the mcmcpvalue function is doing. Any help on
either problem would be very much appreciated.
Thanks, Maureen
*****
Maureen Ryan
Department of Evolution and Ecology
University of California Davis
One Shields Avenue
Davis, CA 95616
(530)304-2266
[[alternative HTML version deleted]]
Hi Maureen, There is a specialist mailing list for mixed-models, where you are likely to get a well-informed answer to this, quite possibly from Prof. Bates himself (the author of lmer). https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models Regards, Mark. PS: R Hompage -> Mailing Lists for other specialist list. Maureen Ryan wrote:> > (To the list moderator: I just subscribed to the list. Apologies for not > having done so longer before trying to post.) > > Hi all, > I am currently using lmer to analyze data from an experiment with a > single fixed factor (treatment, 6 levels) and a single random factor > (block). I've been trying to follow the online guidance for estimating > p-values for parameter estimates on these and other help threads: > > https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html > http://pidgin.ucsd.edu/pipermail/r-lang/2007-August/000057.html > > but have not yet been successful. This is my code: > > lnmass <- lmer(log.mass ~ treatment + (1|block), data=exp1) > summary(lnmass) > > samp <- rnorm(n=10000) > mcmcpvalue <- function(samp) > {std <- backsolve(chol(var(samp)), > cbind(0,t(samp)) - colMeans(samp), > transpose = TRUE) > sqdist <- colSums(std*std) > sum(sqdist[-1] > sqdist[1]/nrow(samp) } > > markov1 <- mcmcsamp(lnmass, 10000) > HPDinterval(markov1) > > mcmcpvalue(as.matrix(markov1[,1])) > mcmcpvalue(as.matrix(markov1[,2])) > mcmcpvalue(as.matrix(markov1[,3])) > mcmcpvalue(as.matrix(markov1[,4])) > mcmcpvalue(as.matrix(markov1[,5])) > mcmcpvalue(as.matrix(markov1[,6])) > > The first time I tried it HPDinterval generated CIs for each treatment > level, etc. but after each of the last 6 lines of code I received the > following error messages: > > Error in markov[,1] : object is not subsettable > Error in as.matrix(markov1[,1]) : error in evaluating the argument 'x' in > selecting a method for function 'as.matrix' > Error in chol(val(samp)) : error in evaluating the argument 'x' in > selecting > a method for function 'chol' > > When I came back to it today HPDinterval(markov1) also no longer worked > and > gave me the following error message: > > Error in UseMethod("HPDinterval"): no applicable method for "HPDinterval" > > I imagine at least some of this has to do with the fact that I don't fully > understand what exactly the mcmcpvalue function is doing. Any help on > either problem would be very much appreciated. > > Thanks, Maureen > > ***** > Maureen Ryan > Department of Evolution and Ecology > University of California Davis > One Shields Avenue > Davis, CA 95616 > (530)304-2266 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- View this message in context: http://www.nabble.com/lmer%2C-estimation-of-p-values-and-mcmcsamp-tp22177307p22178938.html Sent from the R help mailing list archive at Nabble.com.
For getting the p-values from a lmer object you need to use the pvals.fnc() function. do the following: library(lme4) library(languageR) lnmass <- lmer(log.mass ~ treatment + (1|block), data=exp1) pvals.fnc(lnmass, nsim=10000) See the documentation for futher details. Best regards, Ofir. Maureen Ryan wrote:> > (To the list moderator: I just subscribed to the list. Apologies for not > having done so longer before trying to post.) > > Hi all, > I am currently using lmer to analyze data from an experiment with a > single fixed factor (treatment, 6 levels) and a single random factor > (block). I've been trying to follow the online guidance for estimating > p-values for parameter estimates on these and other help threads: > > https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html > http://pidgin.ucsd.edu/pipermail/r-lang/2007-August/000057.html > > but have not yet been successful. This is my code: > > lnmass <- lmer(log.mass ~ treatment + (1|block), data=exp1) > summary(lnmass) > > samp <- rnorm(n=10000) > mcmcpvalue <- function(samp) > {std <- backsolve(chol(var(samp)), > cbind(0,t(samp)) - colMeans(samp), > transpose = TRUE) > sqdist <- colSums(std*std) > sum(sqdist[-1] > sqdist[1]/nrow(samp) } > > markov1 <- mcmcsamp(lnmass, 10000) > HPDinterval(markov1) > > mcmcpvalue(as.matrix(markov1[,1])) > mcmcpvalue(as.matrix(markov1[,2])) > mcmcpvalue(as.matrix(markov1[,3])) > mcmcpvalue(as.matrix(markov1[,4])) > mcmcpvalue(as.matrix(markov1[,5])) > mcmcpvalue(as.matrix(markov1[,6])) > > The first time I tried it HPDinterval generated CIs for each treatment > level, etc. but after each of the last 6 lines of code I received the > following error messages: > > Error in markov[,1] : object is not subsettable > Error in as.matrix(markov1[,1]) : error in evaluating the argument 'x' in > selecting a method for function 'as.matrix' > Error in chol(val(samp)) : error in evaluating the argument 'x' in > selecting > a method for function 'chol' > > When I came back to it today HPDinterval(markov1) also no longer worked > and > gave me the following error message: > > Error in UseMethod("HPDinterval"): no applicable method for "HPDinterval" > > I imagine at least some of this has to do with the fact that I don't fully > understand what exactly the mcmcpvalue function is doing. Any help on > either problem would be very much appreciated. > > Thanks, Maureen > > ***** > Maureen Ryan > Department of Evolution and Ecology > University of California Davis > One Shields Avenue > Davis, CA 95616 > (530)304-2266 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- View this message in context: http://www.nabble.com/lmer%2C-estimation-of-p-values-and-mcmcsamp-tp22177307p22180956.html Sent from the R help mailing list archive at Nabble.com.