How do I specify random effects using gamm models in R?

0
votes

I am using gamm models in mgcv package to analyse how specific diversity measures eg. Shannon vary over time and with environmental variables eg. temperature.

I have the initial model so far to analyse time series: 

modf<-gamm(y~ as.factor(year) + s(doy,bs='cc',k=kdy),method=mth,correlation=tcor,data=d,
           control=ctrl,random=NULL,gamma=1)

I want to include temperature as random effect and thought of doing something like:

modf<-gamm(y~ as.factor(year) + s(doy,bs='cc',k=kdy), + s(temp,bs="re"),method=mth,
           correlation=tcor,data=d,control=ctrl,gamma=1)

However, so far I have only seen this for gam not gamm. Does it still work this way?

An example of type of data structure:

  • $ total_abundance: num 6364161 1929775 7057036 1266342 3981198 ...
  • $ shannon : num 1.87 2.08 1.95 1.84 2.06 ...
  • $ turnover : num 0.613 0.475 0.525 0.556 0.429 ...
  • $ year : int 1985 1986 1987 1987 1987 1988 1989 1989 1991
  • $ month : int 8 12 3 7 8 5 1 8 1 9 ...
  • $ day : int 20 8 15 6 17 9 16 29 14 4 ...
  • $ temp : num 25.5 9.87 4.8 19.72 26.03 ...
  • $ doy : num 232 342 74 187 229 130 16 241 14 247 ...

where doy is 'day of year' and accounts for seasonality

Thanks

1
rmixed-modelsgammgcv
Can you elaborate on what you mean by "add in temperature as random effect"? Is temp a continuous variable or something else? What is the grouping structure you have in mind? Might I suggest including some sample data in your question -- you are more likely to get a helpful answer that way. – Weihuang Wong

1 Answers

0
votes

What you seem to be wanting to do doesn't make any sense unless you want the linear effect of temp to vary within the levels of a grouping factor.

Typically, you specify this random slope using random as

list(group = ~ x)

where group is a factor grouping variable and x is your temp.

What you are asking for with temp = ~ 1 is a random intercept for the unique values of temp, which is probably asking too much of the data.

The equivalent of s(time, bs = "re") requires you to remove the intercept from the random formula:

 list(group = ~ x - 1)

but you still need a group variable.

If you just want to control for temp, add it as a linear parametric effect (+ temp) or a smooth effect (+ s(temp)) in the gamm() model formula.