[R-sig-ME] adding a kinship matrix to a GLMM

[mhorton at uchicago.edu]

Hello,

I'd like to test the mixed-model approach taken in EMMA (-X):
http://www.nature.com/ng/journal/v42/n4/abs/ng.548.html

in a GLMM framework to analyze non-normal data (e.g. count data; sometimes 
there are lots of zeros). I'm trying this with MASS' glmmPQL, which seems to 
allow the user to provide an "an optional correlation structure".

There are a few classes that extend corStruct, but if I try corSymm, things seem 
to work with:

K <- read.table("kinship.txt",...); # an n x n 'identity by state' matrix
cs.K <- corSymm(K[lower.tri(K)],fixed=T);
cs.K <- Initialize(cs.K,data=data.init);
test <- glmmPQL( response ~ snps_i + offset(log(offset)),
          random=~1|subject,correlation=cs.K, family="quasipoisson" );

# data.init is a 1-column matrix of unique subject ids.
# there are between 1-4 biological replicates per subject

I'm writing to ask is this the best model formulation?
Also, if I add a co-factor 'block', which has 6 levels (subjects do not occur in all 
blocks), I get a lot of NA errors:
test <- glmmPQL( response ~ snps_i + as.factor(block)+ offset(log(offset)), 
          random=~1|subject,correlation=cs.K, family="quasipoisson" );
iteration 1
iteration 2
iteration 3
Error in na.fail.default(list(subject = c(1L, 1L, 1L, 2L, 2L, 2L, 3L,  : 
  missing values in object

Again, is there a better model formula?  Is corSymm even doing what I think it is 
doing here? 

Thanks for your time and any comments!
Matt