Dear all,
I am trying to simulate data sets from a model fitted with glmmPQL, in
order to compute the distribution of a summary statistics. My data are
binomial and I have a correlation term in my model. My model is
structured in the following way
m <- glmmPQL( fixed = cbind(sucess,failure) ~ x1 + x2 + ... ,
random = ~ 1 | bidon,
correlation = corGaus(form=~ longitude + latitude),
family = binomial)
My simulation model would be something like
#the success probability
y <- rmvnorm(1,predict(m),var.covar.matrix)
#the binomial sampling
k <- rbinom(length(y),size=n,prob=y)
In order to run these simulations, therefore, I need the predict of my
model, which I extract easily, and the residual variance-covariance
matrix (the var.covar.matrix object in the previous code). I can extract
with no problem the correlation matrix, from the corStruct object that
glmmPQL returns. I guess I then need to multiply this matrix by a
residual variance to obtain what I need. And here, I do have a problem:
I cannot find a way to estimate the variance that is solely due to the
random effect. The sigma that glmmPQL returns seems to be the sum of the
variance I need and the variance that is produced by the binomial
sampling.
There is probably a simple way to do this, but I just can't find it. Any
help on this would be welcome !
A separate issue: I have no idea how I should run my simulation in case
I need to set family to quasibinomial, instead of binomial... And this
is something I might consider in a next future.
Thanks in advance for your help,
--
Jean-Baptiste Ferdy
Institut des Sciences de l'?volution de Montpellier - UMR 5554
Universit? Montpellier 2
t?l. (0)4 67 14 42 27
fax ?(0)4 67 14 36 22
--
Jean-Baptiste Ferdy
Institut des Sciences de l'?volution de Montpellier - UMR 5554
Universit? Montpellier 2
t?l. (0)4 67 14 42 27
fax ?(0)4 67 14 36 22
