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