[R] Question about fitting a periodic model in glmm

Marc Girondot marc_grt at yahoo.fr
Tue May 7 20:24:42 CEST 2013


I would like to fit a period (annual) model in glmm. Here is the script 
I do:
# Generate "dummy" periodic counts with effect of a covariate co
# of course I plan to use this script on my own data !
d <- 1:500
co <- rnorm(500, 10, 2)
yco <- (1+sin(2*pi*(d+100)/365))*10*co/10+co
y <- floor(rnorm(500, yco, 10))
df1 <- data.frame(days=d, number=y, covariate=co, ID=1)
df1[df1$number<0, "number"] <- 0

# Just to look that all is ok:
plot(df1$days, df1$number, type="l", ylim=c(0,80), bty="n")
plot(df1$number, df1$covariate, bty="n")

# days is a fixed effect (I choose the days of observations)
# covariate is a random effect
# I fit periodic effect according to days as: sin( 2*3.1416 * (days) / 
365) + cos( 2*3.1416 * (days) / 365)
library(MASS)
fit <- glmmPQL ( number ~ covariate +
                    sin( 2*3.1416 * (days) / 365) +
                    cos( 2*3.1416 * (days) / 365),
                  family=quasipoisson(link = "log"), data=df1,
                  random = ~ 1+covariate | ID)

# test for the effects
library(spida)
wald( fit, list("covariate", "days"))

# predictions: all is good !
plot(df1$days, df1$number, type="l", ylim=c(0, 80), bty="n", main="ID=1")
par(new=TRUE)
newd1 <- data.frame(days=d, covariate=5, ID=1)
p1 <- predict(fit, newd1)
plot(d, exp(p1), type="l", col="red", ylim=c(0, 80),
      bty="n", axes=FALSE, xlab="", ylab="")
par(new=TRUE)
newd1 <- data.frame(days=d, covariate=10, ID=1)
p1 <- predict(fit, newd1)
plot(d, exp(p1), type="l", col="green", ylim=c(0, 80),
      bty="n", axes=FALSE, xlab="", ylab="")
par(new=TRUE)
newd1 <- data.frame(days=d, covariate=15, ID=1)
p1 <- predict(fit, newd1)
plot(d, exp(p1), type="l", col="blue", ylim=c(0, 80),
      bty="n", axes=FALSE, xlab="", ylab="")
legend("topleft", legend=c("covariate=5", "covariate=10", 
"covariate=15"), lty=1, col=c("red", "green", "blue"))

Now my questions:
- The periodic effect has two components, sin and cos dependent on 
"days". After the Wald test, I have then two p values for "days" effect 
(one for sin and one for cos). How can I combine these two p-values to 
get a global effect of periodic effect ?
- If I want to setup interaction between periodic effect of "days" and 
"covariate", how I can do as "days" appears in two effects (sin and cos) ?

A third question related, is it possible to use a glmm with negative 
binomial distribution ? I don't find still...

... or perhaps you have a better way to do all of that !

Thanks a lot.

Marc Girondot



More information about the R-help mailing list