gmane.comp.lang.r.lme4.devel

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. Error t value
(Intercept) 2.877840   0.019648  146.47
x           0.001337   0.001264    1.06

Correlation of Fixed Effects:
  (Intr)
x 0.002

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 weights. Is this a valid approach?Alternatively I was wondering if I can do this directly in a mixed model. From an example in Gelman and Hill (2007) book I thought I could model the dependence of conspecific effects on species-level attributes by including an interaction between my conspecific effect and the species level attribute. This is, for the fixed part of my model:

ALIVE~CONSPp+TOTABN+ LANDSPABN+ CONSPp:LANDSPABN

A non-trivial question about the model implementation is whether the random term (1+CONSPp| fSPECIES) has to be included or not in the model. The inclusion of this term in the model is what makes the difference between the interaction being significant or not. In the traditional ANOVA framework (which I am still more familiar with) I would say that once I include the random slope in the model my single species-level moderator ( LANDSPABN) does not add anything to the model but I think this is not the way it works in the mixed modeling (likelihood) framework.
Is the interaction term (CONSPp:LANDSPABN) really answering the question on the dependence of the conspecific effect on the commoness of the species?Is including the random term (1+CONSPp| fSPECIES) needed in the model?Any light on these issues is highly appreciated.Best,

Edwin

[[alternative HTML version deleted]]

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 this?

Greetings,
Bertolt

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: num  225 225 225 225 225 225 225 225 225 225 ...
 $ [log transformed ] outcome : num  NA 2.887 0.963 4.006 2.06 ...

Using nlme package, my current model looks like

summary(mR01A<- lme(lnSICIrest~factor(agegroup)*time_point*time_point2*factor(stim.cond),data=tDCSrest, random=~1|session/ID,na.action=na.exclude,method="REML"))

and 
summary(mR01_corAR <- update(mR01, correlation = corAR1()))
respectively.


My questions are: 
1. Does the rather complex random effects structure is in any way reasonable or am I totally on the wrong path here?

2. For the model with corAR1 correlation structure, 95% CI cannot be produced - the following error message is returned: "Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance." 
Is it necessary and appropriate in addition to specifying the random effects as above to account for the autocorrelation structure of the longitudinal data? Using anova to compare the models, I get a BIC change of -8.621 with a df increase of 2 when corAR1 correlation structure is added to the model. And if I do not need to include corAR1, would it be better to use the lme4 package or is it just a matter of taste?


I would be very grateful for any thoughts and comments!

kirstin







On 18.05.2012, at 15:13, Kirstin-Friederike Heise wrote:



--
Pflichtangaben gemÃ¤Ã Gesetz Ã¼ber elektronische Handelsregister und Genossenschaftsregister sowie das Unternehmensregister (EHUG):

UniversitÃ¤tsklinikum Hamburg-Eppendorf; KÃ¶rperschaft des Ã¶ffentlichen Rechts; Gerichtsstand: Hamburg

Vorstandsmitglieder: Prof. Dr. Guido Sauter (Vertreter des Vorsitzenden), Dr. Alexander Kirstein, Joachim PrÃ¶lÃ, Prof. Dr. Dr. Uwe Koch-Gromus 

[[alternative HTML version deleted]]

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 / 99 97 97 00




2012/5/23 Edwin Lebrija Trejos <elebrija-PkbjNfxxIARBDgjK7y7TUQ< at >public.gmane.org>


[[alternative HTML version deleted]]

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 varying intercept and slope for the effect of conspecifics, i.e. (1+CONSPp|SPECIES) and to extract from the results the conditional modes (using terminology in Prof. Bates draft book) of the random effects (and their variances) to model them separately as a function other continuous predictors (using simple weighted linear regressions). The full model of this approach is: M1<-lmer(ALIVE~CONSPp+TOTABN+(1+CONSPp|SPECIES)+(1|PLOT)+(1|YEAR), data=oneyrseedl,family=binomial).  After further reading I found an example from Gelman and Hill 2007 in which he includes "group level predictors". Specifically in their example, a predictor X1 appears both as fixed (y~ X1) and random effect (1+ X1 | grouping factor) with the grouping factor predictor (X2) modelled as an interaction with X1. So their full model is y ~ X1 + X2 + X1:X2+ (1 + X1 | grouping factor).  I thought I may be able to specify a similar model to directly incorporate my continuous species-level variable (species abundance across the landscape) to explain the variation in the effect of conspecifics on seedling survival. One crucial difference between my case and Gelman and Hill´s case is that my continuous species-level variable has only one value per species (my grouping factor) while their X2 has a series of values per grouping factor. I thus think that if include the continuous species-level variable in the model I do not need to include a random slope for the effects of conspecifics (with SPECIES as a grouping factor) as I feel it would control for all the variation that could be more meaningfully explained by my species-level variable (or any other species-level variable I may come up later).  Thus, the model I am thinking about would  be: M2<-lmer(ALIVE~CONSPp+TOTABN+ LANDSPABN+ CONSPp:LANDSPABN + (1|PLOT)+(1|YEAR), data=oneyrseedl,family=binomial).  I am doubting if a random intercept by SPECIES is still needed (i.e. (1|fSPECIES).I hope I have exposed my problem clearly. I would highly appreciate any insight from the mixed model experts.Kind regards,Edwin        
[[alternative HTML version deleted]]

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, and you can do so.  If
you want a SAS independent solution, you can get the DF from lme() and
use those with gls().

Doug Bates has removed DF calculations from his more recent package,
lme4 which is sort of the new nlme, so the precedent is actually for
not providing any, leaving to the user to choose sensibly.  Finally
because new development is happening in lme4 not nlme, I rather doubt
any of the authors are going to want to put effort into altering code
in nlme (particularly because it could potentially break other code).

If you are doing this a lot and want a niceish interface, my
suggestion would be to hack the anova.gls() code to include
denominator DF, call it anovadf.gls() or something, and just put that
code in your .Rprofile so that function is always available to you.

Cheers,

Josh

On Tue, May 22, 2012 at 8:01 AM, Charles Determan Jr <deter088< at >umn.edu> wrote:

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

Charles Determan Jr — 2012-05-22T15:31:57

Jim,

Thank you for contributing.  I don't think it is ever a bad thing to
explain some SAS information here and I for one appreciate it.  I am well
aware of the different ways in which degrees of freedom (as seen in SAS)
can be calculated and that there is no true 'correct' way for all
situations.  However, I am only interested in this particular one for
consistency. I very much prefer using R for the flexibility in the analysis
and would like to be able to accomplish this particular component.

Regards,
Charles

On Tue, May 22, 2012 at 10:20 AM, Baldwin, Jim -FS <jbaldwin-SjfvJOW+wKv1P9xLtpHBDw< at >public.gmane.org>wrote:


[[alternative HTML version deleted]]

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

Baldwin, Jim -FS — 2012-05-22T15:20:52

Have you considered the other direction?  SAS offers 5 different options for constructing the denominator degrees of freedom:  BETWITHIN, CONTAIN, KENWARDROGER, RESIDUAL, SATTERTHWAITE.  Also reading the associated references for those methods will help with understanding why there isn't a single method that is "correct" for all situations.  (My apologies for showing part of a SAS manual on this list.)

Jim

Jim Baldwin
Station Statistician
Pacific Southwest Research Station
USDA Forest Service
Albany, California


-----Original Message-----
From: r-sig-mixed-models-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org [mailto:r-sig-mixed-models-bounces-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org] On Behalf Of Charles Determan Jr
Sent: Tuesday, May 22, 2012 8:01 AM
To: Douglas Bates
Cc: r-sig-mixed-models-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
Subject: Re: [R-sig-ME] 2nd attempt - conflict of dfs or f value in lme

Dr. Bates,

After another thought, would it be possible to modify the gls source code to incorporate the same denDF calculation used in the lme function from this package?  I have been looking through the source code of the nlme package and by adding the denDF to the gls function I would solve my initial problem of having consistent F values and dfs.  However, I would prefer to work with one of the authors regarding any modification.  Perhaps I should contact R core mailing list?

Regards,
Charles

On Mon, May 21, 2012 at 12:22 PM, Douglas Bates <bates-GX8I/T4BApV4piUD7e9S/g< at >public.gmane.org> wrote:


        [[alternative HTML version deleted]]





This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.

Re: Calculating repeatability with ordinal data

Samantha Patrick — 2012-05-22T15:12:57

Hi Pierre

Thank you for hints on priors - I will definitely try them out!

Unfortunately adding a global intercept does not fix the problems with 
the contrasts.

Thanks again

Sam

Le 22/05/2012 16:43, Pierre de Villemereuil a Ã©crit :

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

Charles Determan Jr — 2012-05-22T15:01:15

Dr. Bates,

After another thought, would it be possible to modify the gls source code
to incorporate the same denDF calculation used in the lme function from
this package?  I have been looking through the source code of the nlme
package and by adding the denDF to the gls function I would solve my
initial problem of having consistent F values and dfs.  However, I would
prefer to work with one of the authors regarding any modification.  Perhaps
I should contact R core mailing list?

Regards,
Charles

On Mon, May 21, 2012 at 12:22 PM, Douglas Bates <bates-GX8I/T4BApV4piUD7e9S/g< at >public.gmane.org> wrote:


[[alternative HTML version deleted]]

Re: Calculating repeatability with ordinal data

Pierre de Villemereuil — 2012-05-22T14:43:55

Hi !

You might like to look at the following reference:
Nakagawa, Shinichi, and Holger Schielzeth. âRepeatability for Gaussian 
and Non Gaussian Data: a Practical Guide for Biologists.â /Biological 
Reviews/ 85, no. 4 (November 1, 2010): 935â956.

Concerning the error you get, did you try to remove the "-1" in your 
formula ? Maybe this is the contrast MCMCglmm is complaining about:
model1c<-MCMCglmm(scoremax~trait, random =~Bird + year, rcov = 
~us(trait):units, data = Data, prior = priortest2, family = "ordinal")

Besides, I think you should try out a 'chi2' prior distribution for your 
variances and compare with the results obtained using 'priortest2':
priortest3 <- list(R=list(V = 1 , fix=1), G = list(G1=list(V=1, nu=1000, 
alpha.mu=0, alpha.V=1),G2=list(V=1, nu=1000, alpha.mu=0, alpha.V=1)))

I found that, at least in some cases, this prior behaves better than the 
one in 'priortest2' for binary data and a check  comparing the two 
priors is worth a try ! It is always better to check the sensitivity to 
prior in a Bayesian analysis, and especially for binary data, because 
they 'contain' less information and you have to fix the residual 
variance. As such, your model can be very sensitive to the choice of 
your prior !

Cheers,
Pierre.

Le 22/05/2012 16:22, Samantha Patrick a Ã©crit :


[[alternative HTML version deleted]]