[R-sig-ME] profile and/or bootstrapped CIs for GAMM random effects

Wed Sep 2 15:48:01 CEST 2015

Hi Ben
Thanks for looking into this for me. It looks like you've been able to develop a reproducible example, but below is a subset of my data for reference. The full data are radio telemetry-based distances (longdist) for 28 fish on a flood plain, located approximately every 2 weeks. We're interested in exploring the temporal movement patterns of fish, and then identify potential environmental drivers of this. Fish leave their home pool in the wet season, move onto the floodplain for a period, then return pretty much back to their home pool

In regards to the Gamma/log-link model: I decided to go down this path as my response variable is distance from 'home' (longdist) which is strictly positive and continuous. I realise I could log-transform the data, but thought it good to just use an appropriate distribution. I'm happy to be convinced to go down the transformation path as it would help with other response variables that are continuous but include zeros (here I've had to split the analysis into two steps: 1) a binomial GLMM for zero or >0 distance; 2) a gamma GLMM for all distances >0)

Cheers
John

mydata <- structure(list(CODE = structure(c(4L, 3L, 8L, 3L, 8L, 4L, 8L,3L, 4L, 8L, 4L, 3L, 8L, 3L, 4L, 8L, 4L, 8L,
          4L, 8L, 1L, 3L, 7L,2L, 5L, 6L, 4L, 8L, 1L, 2L, 3L, 7L, 8L, 3L, 5L, 6L, 8L, 3L, 1L,7L, 2L, 5L, 6L, 1L, 3L, 8L, 5L,
          2L, 7L, 6L, 8L, 3L, 1L, 2L, 5L,6L, 7L, 8L, 3L, 1L, 5L, 2L, 6L, 7L, 3L, 1L, 5L, 2L, 7L, 8L, 1L,8L, 7L), 
         .Label = c("113.08", "113.09", "113.18", "113.2", "142","142.02", "142.11", "142.16"), class = "factor"),
          longdist = c(915.3204902,293.3651772, 121.817523, 926.8937579, 1657.624408, 634.3312583,192.829031, 258.3748562,
          612.3309929, 1369.839623, 515.8210633,243.1597688, 3955.710259, 9535.900065, 6882.407512, 1643.480296,9098.188882,
          1151.279039, 6526.611439, 2590.840377, 2872.818688,4482.909503, 277.6635567, 4188.271695, 8267.951169,
          3319.834396,11015.65632, 26744.47918, 1540.673734, 3623.761144, 5309.040932,452.5968432, 3043.91542, 4473.26457,
          7251.069327, 10195.13599,1455.050153, 4270.682608, 2341.194885, 1413.999748, 4296.287262,1870.788559, 5058.781157,
          1796.509776, 3564.458943, 321.5355384,2613.835935, 3616.955753, 167.6945737, 3768.178122, 1323.04003,3626.427229,
          342.0831861, 4025.501297, 1759.654454, 3575.230119,222.7191619, 1936.035776, 2087.192677, 66.80421153,
          3317.121218,3718.282656, 6090.502215, 328.044963, 3603.478304, 339.1607537,3139.444897, 4376.683207, 136.4894718, 1421.734876, 
         477.0550425,174.7315492, 1347.494408),
          time = c(6L, 6L, 6L, 16L, 16L, 16L,27L, 27L, 27L, 41L, 41L, 41L, 55L, 55L, 55L, 70L, 70L, 92L, 92L,105L, 105L, 105L,
          105L, 105L, 105L, 105L, 105L, 119L, 119L, 119L,119L, 119L, 133L, 133L, 133L, 133L, 147L, 147L, 147L, 147L, 147L,147L,
          147L, 164L, 164L, 164L, 164L, 164L, 164L, 164L, 176L, 176L,176L, 176L, 176L, 176L, 176L, 189L, 189L, 189L, 189L,
          189L, 189L,189L, 204L, 204L, 204L, 204L, 204L, 204L, 216L, 216L, 216L)), .Names = c("CODE","longdist", "time"),
          class = "data.frame", row.names = c(43L,48L, 53L, 80L, 82L, 84L, 101L, 109L, 111L, 129L, 132L, 134L,153L, 158L,
          159L, 174L, 189L, 205L, 208L, 241L, 247L, 252L, 253L,258L, 260L, 262L, 265L, 272L, 281L, 283L, 287L, 291L, 296L,
          301L,306L, 311L, 323L, 329L, 331L, 337L, 338L, 346L, 349L, 354L, 359L,365L, 367L, 369L, 373L, 377L, 380L, 381L,
          388L, 393L, 398L, 399L,401L, 407L, 409L, 415L, 421L, 424L, 432L, 433L, 436L, 437L, 449L,451L, 455L, 458L, 464L, 471L, 486L))

f2<-gamm4((longdist)~s(time,k=4),random=~(1|CODE),data=mydata,family=Gamma(link=log))
confint(f2$mer,method="profile")
confint(f2$mer,method="boot")

----------------------------------------------------------------------

Message: 1
Date: Tue, 1 Sep 2015 20:44:39 -0400
From: Ben Bolker <bbolker at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] profile and/or bootstrapped CIs for GAMM
	random effects
Message-ID: <55E64677.5030905 at gmail.com>
Content-Type: text/plain; charset=windows-1252

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

  Is there any chance you could post a reproducible example ... ?

  (Out of curiosity, do you have strong reasons to use a Gamma/log-link rather than a log-Normal model?)

  Ben Bolker

On 15-09-01 06:51 PM, John Morrongiello wrote:
> Hi all I thought I'd re-post this message from last week as it may 
> have got lost in the weekend traffic:
> 
> I'm having trouble estimating CIs for random effects (i.e. CODE) from 
> a gamm4 model fit with a gamma distribution, My model is:
> 
> c1<-gamm4((longdist)~s(time,k=4),random=~(1|CODE),data=catfish1,family
> =Gamma(link=log))
>
> 
confint(c1$mer,method="profile")
> ##this returns Error in mkMerMod(rho = environment(devfun), opt = opt, 
> reTrms = b$reTrms,  : unused argument (devFunOnly = TRUE)
> 
> confint(c1$mer,method="boot") ##this returns Error in 
> mkNewReTrms(object, newdata, compReForm, na.action = na.action,  :
> random effects specified in re.form that were not present in original 
> model Are these errors arising because I'm using a gamma distribution? 
> Is there another way to get the random effect CIs?
> 
> I can get random effect CIs for the following model: 
> c2<-gamm4((longdist)~s(time,k=4),random=~(1|CODE),data=catfish1,family
> =gaussian)
>
>  confint(c2$mer,method="profile")##works
> confint(c2$mer,method="boot")##same error message as above

--
Dr. John R. Morrongiello
School of BioSciences
University of Melbourne
Victoria 3010, Australia
T: +61 3 8344 8929
M: +61 403 338 554
E: john.morrongiello at unimelb.edu.au
W: morrongiellolab.com