[R] extracting residual variance from glmmPQL

[Jean-Baptiste Ferdy]

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