gmane.comp.lang.r.lme4.devel

Fwd: error while trying to get predictions from a lmeobject

Antonio P. Ramos — 2012-05-26T00:29:35

---------- Forwarded message ----------
From: Antonio P. Ramos <ramos.grad.student-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org>
Date: Fri, May 25, 2012 at 5:29 PM
Subject: error while trying to get predictions from a lme object
To: Antonio Pedro Ramos <ramos.grad.student-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org>


Another example of the error I am taking about:

+              my.new.time.one.transition.low.and.middle + ttd +
+              maternal_educ+ log(IHME_id_gdppc) + hiv_prev-1,
+              merged0,na.action=na.omit,method="ML",weights=varPower(form
= ~ time),
+              random= ~ time| country.x,
+              correlation=corAR1(form = ~ time),
+              control=lmeControl(msMaxIter = 200, msVerbose = TRUE))


 expand.grid(time=20:29,country.x=c("Poland","Brazil","Argentina"))
as.data.frame(cbind(my.new.time.one.transition.low.and.middle=c(0,0,0,0,0,0,1,2,3,4),
+       ttd=c(0,0,0,0,0,0,1,0,0,0),
+       maternal_educ=seq(from=10.0, to=12.0, length.out=10),
+       IHME_id_gdppc=log(seq(fro

error while trying to get predictions from a lme object

Antonio P. Ramos — 2012-05-25T23:09:52

Hi all,

I am trying to get predictions for individual level observations for an lme
object. Yet, since I am get different types of errors for different trials,
it seems to me I am missing something.  My model is the following:

model <- lme(log(child_mortality) ~ as.factor(cluster)*time +
             my.new.time.one.transition.low.and.middle + ttd +
             maternal_educ+ log(IHME_id_gdppc) + hiv_prev-1,
             merged0,na.action=na.omit,method="ML",weights=varPower(form =
~ time),
             random= ~ time| country.x,
             correlation=corAR1(form = ~ time),
             control=lmeControl(msMaxIter = 200, msVerbose = TRUE))


It runs fine and the results make sense. Now my attempts to get
predictions. For example:

 data.frame(time=c(10,10,10,10),country.x=c("Poland","Brazil","Argentina","France"),
+
my.new.time.one.transition.low.and.middle=c(1,1,1,0),ttd=c(0,0,0,0),
+
maternal_educ=c(10,10,10,10),IHME_id_gdppc=c(log(5000),log(8000),log(8000),log(15000)),
+                          hi

Re: Error in lme4 0.999902344-0: "Object 'multResp'notfound" ?

Ben Bolker — 2012-05-25T14:40:44


  In principle yes; profiling a merMod (new-lme4) object creates
a profile object of class "thpr", which has a confint method.
The confidence intervals produced in this way are more accurate than
the Wald intervals implied by summary() (analogous to the difference
between confint.default() and confint.glm() applied to a glm() fit),
but are still based on Z-score/likelihood ratio test cutoffs, rather
than trying to do any kind of finite-size correction.  If you want
more accurate confidence intervals that incorporate the effects of
the sample size, then you'll need bootMer.

   However ... for GLMMs, we are still wrestling a bit with the
robustness of the fitting, which tends to be more of a problem
when profiling -- so it may not be practical to get profile
confidence intervals on GLMMs quite yet ...

  Ben Bolker

Re: Error in lme4 0.999902344-0: "Object 'multResp'notfound" ?

Bertolt Meyer — 2012-05-25T07:09:00


That does indeed solve the problem, thanks for the quick fix! A quick question: A colleague told me that this newest version of lme4 can do CIs for the parameter estimates of the fixed effects. Does that require bootMer() or is there another way?

Thanks again,
Bertolt

Re: Error in lme4 0.999902344-0: "Object 'multResp' notfound" ?

Ben Bolker — 2012-05-24T20:36:10


  thanks for the notification.  I've just posted 0.999902345-0
to the repository (it may take up to 24 hours to appear), should
address this issue.
   I've also updated the source version on the repository
(i.e. http://lme4.r-forge.r-project.org/repos ) to this
version; the Windows version is a bit older, but current
source and Windows versions (unlike MacOS) should be available
from the standard package directory at http://r-forge.r-project.org
(i.e. install.packages("lme4",repos="http://r-forge.r-project.org"))

  Ben Bolker

Re: Error in lme4 0.999902344-0: "Object 'multResp' not found" ?

Douglas Bates — 2012-05-24T18:17:16

Not a problem on an Ubuntu system.  It may be something to do with the
version of the Mac OS X package on R-forge.  Ben is the person who
creates those and he is very busy right now with teaching obligations.
 Do you have the capability of compiling the source package?

Here are the results on an Ubuntu 12.04 system

Loading required package: lattice
Loading required package: Matrix
+                       group <- factor(group)))
'data.frame':386 obs. of  4 variables:
 $ X    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ group: Factor w/ 97 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ y    : num  2.91 2.75 3.17 2.98 2.93 ...
 $ x    : num  -8.17 -13.5 -1.5 -4.84 -8.17 ...
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | group)
   Data: my_data

REML criterion at convergence: 142.4097

Random effects:
 Groups   Name        Variance Std.Dev.
 group    (Intercept) 0.02075  0.1440
 Residual             0.06639  0.2577
Number of obs: 386, groups: group, 97

Fixed effects:
            Estimate Std

Re: Group level predictors in mixed models

Edwin Lebrija Trejos — 2012-05-24T17:16:00




Hi again to all,
 
Honoré, thanks for your reply. I feel that I did not explain my main doubt properly (plus I realized that I made a mistake about Gelman and Hill´s book example: their grouping factor predictor also has one value per level of the grouping factor and not a series of them as I thought before). 
My main doubt is whether I can answer my question using a mixed model approach.  I want to understand whether the effect on survival of the proportion of conspecifics in the seedling neighborhood is moderated by species-level attributes. One of such attributes is how common (abundant) is the species across the landscape (my LANDSPABN variable). 
I have been thinking that I can do this in two ways. One is conducting two separate analysis: first using lmer to obtain (random) species level coeficients for the effect of conspecifics and then use simple weighted linear regressions to explain those coeficients as a function of species commonnes. I would use the variance of the estimated coeficients as

Error in lme4 0.999902344-0: "Object 'multResp' notfound" ?

Bertolt Meyer — 2012-05-24T15:30:07

Dear list,

The simplest random intercept model that I try to fit with the newest version of lme4 (0.999902344-0) throws an "Object 'multResp' not found" error on my Mac (Mac OSX 10.7.3, German R 2.15.0). Here is a reproducible example:

install.packages(c("minqa", "Rcpp"))
install.packages("lme4", repos="http://lme4.r-forge.r-project.org/repos")
library(lme4)

my_data <- read.csv("http://dl.dropbox.com/u/5384027/my_data.csv")
my_data$group <- as.factor(my_data$group)

# the data set consists of 97 four-person teams
# with two variables x and y observed on the individual
# level in each team

head(my_data)
#   X group        y          x
# 1 1     1 2.914286  -8.170984
# 2 2     1 2.746269 -13.504318
# 3 3     1 3.171429  -1.504318
# 4 4     1 2.978723  -4.837651
# 5 5     2 2.928571  -8.170984
# 6 6     2 2.987013   8.495682

mlmodel1_ri <- lmer(y ~ x + (1 | group), data = my_data)

# Fehler in lmer(y ~ x + (1 | group), data = my_data) : 
#   Objekt 'multResp' nicht gefunden

Does anyone know how to fix thi

Re: fitting model for repeated measures cross-overdesign?

Malcolm Fairbrother — 2012-05-24T14:23:50

Dear Kirstin,

If I've understood correctly, I'm pretty sure you want "ID/session" not "session/ID", because person-sessions are nested within persons.

And I don't think it makes sense to interact your linear and quadratic time effects (what would such an interaction mean?)… I would think this would make more sense:

factor(agegroup)*time_point*factor(stim.cond) + factor(agegroup)*time_point2*factor(stim.cond)

Even this will generate results that will be hard to interpret, with so many interactions… you'll need lots of plots of the expected outcome for different combinations of agegroup and stim.cond, versus time.

Good luck,
Malcolm

Re: fitting model for repeated measures cross-overdesign?

Kirstin-Friederike Heise — 2012-05-23T08:48:45

Dear list members, 

since I am still struggling with building a linear mixed effects model, I will try to rephrase and specify my questions:

I am analyzing a transient treatment effect (treatment vs. placebo) in a cross-over experiment (2 sessions) in two groups (old, young participants) with repeated measurements (before, during, after treatment in minutes) for both treatment conditions on metric outcome variables. Since I expect a non-linear development of a possible treatment effect over time, I tried to include a quadratic polynomial for time. 

This is my data frame: 

data.frame':5760 obs. of  21 variables:
 $ ID  [=subject level]       : Factor w/ 20 levels "OMI_01","OMI_02",..: 11 11 11 11 11 11 11 11 11 11 ...
 $ agegroup   : Factor w/ 2 levels "OLD     ","YOUNG   ": 2 2 2 2 2 2 2 2 2 2 ...
 $ session    : num  1 1 1 1 1 1 1 1 1 1 ...
 $ stim.cond [=treatment]: Factor w/ 2 levels "sham","tDCS": 1 1 1 1 1 1 1 1 1 1 ...
 $ time_point : num  -15 -15 -15 -15 -15 -15 -15 -15 -15 -15 ...
 $ time_point2

Re: Group level predictors in mixed models

Samadori Honoré Biaou — 2012-05-23T08:22:32

Hi Edwin,

You can perhaps test between alternate models using the R built in anova()
function. It gives a p-value indicating if models are significantly
different (if not, since the addition of the new term does not
significantly improve the model fit, you would prefer the simpler model).
For example, M1 and M2 below differ only by the inclusion of (1|fSPECIES),
and you can probably test if this addition improve the model fit by using
anova():

M1<-lmer(ALIVE~CONSPp+TOTABN+ LANDSPABN+ CONSPp:LANDSPABN +
(1|PLOT)+(1|YEAR), data=oneyrseedl,family=binomial)

M2<-lmer(ALIVE~CONSPp+TOTABN+ LANDSPABN+ CONSPp:LANDSPABN +
(1|PLOT)+(1|YEAR) + (1|fSPECIES), data=oneyrseedl,family=binomial)

anova (M1, M2)                     # The test for the alternate models


Please use "help (anova)" for more details and documentation about this
function.

Hope it helps,

Best regards,

Dr. Ir. Samadori Honoré BIAOU
Web pages:
http://sites.google.com/site/hbiaou/
www.fa-up.bj/staff/Biaou.html
Skype: hbiaou
Tel +229 94 55 81 46 /

Group level predictors in mixed models

Edwin Lebrija Trejos — 2012-05-23T03:35:21




Dear R-sig-mixed-models subscribers, I am analizing survival (ALIVE) of one year-old seedlings of over a hundred species as a function of the total number of seedlings (TOTABN) and the proportion of conspecifics seedlings (CONSPp) in their neighborhood (sampling quadrat). The data comes from monitoring seedling dynamics trough time (YEAR) in a series of forest plots (PLOT). A crucial part of the analysis is to understand the variation in responses between species (SPECIES). For now, I have one continuous species-level variable, abundance across the landscape (LANDSPABN), that I would like to use to explain the variation between species. After reading a number of printed and online references on advanced statistical models it seems to me that my problem could be properly analyzed using hierarchical bayesian modelling but I am just beggining to understand lmer (and R) and I was wondering if using such tools I can also perform a meaningfull analysis.My first approach in lmer was to specify a model with vary

Re: How to obtain t-values from mer objects

Ben Bolker — 2012-05-22T20:27:39


  Generally I think

  coef(summary(fm1))[,3]

 should allow you to get the t-statistics in either case,
*without* delving into the internal structure of the fitted model.
If you find that you consistently need to dig into the internals
of the fitted objects, you should ask on the list and/or ask
the maintainers to provide an appropriate accessor method ...

  (I will admit that I haven't tested this)

  Ben Bolker

Re: How to obtain t-values from mer objects

Douglas Bates — 2012-05-22T20:25:32

The preferred idiom is

coef(summary(fm1))

 num [1, 1:3] 1527.5 19.4 78.8
 - attr(*, "dimnames")=List of 2
  ..$ : chr "(Intercept)"
  ..$ : chr [1:3] "Estimate" "Std. Error" "t value"

because it works in lme4.0 and the development lme4 (and for many
other types of fitted models too).

On Tue, May 22, 2012 at 1:48 PM, Joehanes, Roby (NIH/NHLBI) [F]
<roby.joehanes-2zaOuxCdfhg< at >public.gmane.org> wrote:

Re: How to obtain t-values from mer objects

Gang Chen — 2012-05-22T19:40:34

Yes, that works well. Thanks a lot, Roby!

Gang


On Tue, May 22, 2012 at 2:48 PM, Joehanes, Roby (NIH/NHLBI) [F]
<roby.joehanes-2zaOuxCdfhg< at >public.gmane.org> wrote:

Re: How to obtain t-values from mer objects

Joehanes, Roby (NIH/NHLBI) [F] — 2012-05-22T18:48:07

How about

summary(fm1)< at >coefs[,3]

Or if you are using the next-generation lme4:
summary(fm1)$coef[,3]

Roby


On 5/22/12 2:35 PM, "Chen, Gang (NIH/NIMH) [C]" <gangchen-2loH/HJHQuifRvmTrFJqzg< at >public.gmane.org>
wrote:

Re: How to obtain t-values from mer objects

Gang Chen — 2012-05-22T18:35:04

Sorry I should have stated a little clearer.

I do know how to see or print the t-statistic values at the R prompt.
What I meant to ask is how to extract the t-statistic values from the
mer slots.

Gang

On Tue, May 22, 2012 at 2:21 PM, Andrzej <agalecki-63aXycvo3TyHXe+LvDLADg< at >public.gmane.org> wrote:

Re: How to obtain t-values from mer objects

Andrzej — 2012-05-22T18:21:02

Try summary(fm1)

Andrzej

On 5/22/2012 1:56 PM, Gang Chen wrote:

I know how to obtain the fixed effects from a mer object:



fixef(fm1)

or,


But how can I obtain the t-statistic values for the fixed effects?

Thanks,
Gang

How to obtain t-values from mer objects

Gang Chen — 2012-05-22T17:56:16

I know how to obtain the fixed effects from a mer object:



or,


But how can I obtain the t-statistic values for the fixed effects?

Thanks,
Gang

Re: 2nd attempt - conflict of dfs or f value in lme

Charles Determan Jr — 2012-05-22T15:50:13

Understood, thank you to all who have helped me in this endeavor.

Regards,
Charles

On Tue, May 22, 2012 at 10:42 AM, Joshua Wiley <jwiley.psych-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org>wrote:


[[alternative HTML version deleted]]

Re: 2nd attempt - conflict of dfs or f value in lme

Joshua Wiley — 2012-05-22T15:42:46

Hi Charles,

Unless you have a maniacal employer who demands that you replicate the
fixed effects from SAS, I do not think that adding DFs to the gls
output will really solve your problem.  You are still talking about
rather different models from PROC MIXED to gls().  You may be able to
get equivalent tests of the fixed effects, but I can promise the
random effects are not identical.  So you want to modify stable R
code, in order to get output from R to match a portion of output from
a different model altogether in SAS?  I appreciate that many people
are only interested in the fixed effects and just consider the
nonindependence in their data a nuissance factor, so you may have no
strong preference for random effects versus gls, but we are not
talking about most people here.  We are talking about R which is used
by all sorts of people.

If you need that in your work, you already showed in your original
email you can do it.  You have the F-values, you just want to
calculate p-values based on the DF from SAS, a