adding a kinship matrix to a GLMM

Thanks David. I'll see if I can fit a kinship matrix using the other  
packages you show.
It would also be great to setup the same test case to see how  
quickly any of this will
actually run. I can't do that immediately, but I will try to post  
the results next week.

---- Original message ----

Date: Thu, 26 Jan 2012 10:33:21 +1000 (EST)
From: "David Duffy" <David.Duffy at qimr.edu.au>
Subject: Re: [R-sig-ME] adding a kinship matrix to a GLMM
To: <mhorton at uchicago.edu>
Cc: <r-sig-mixed-models at r-project.org>

On Tue, 24 Jan 2012, mhorton at uchicago.edu wrote:

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:

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?

I think so.  That error might be arising from incomplete data, as it says
(do debug(lme) to see where, as this is what glmmPQL calls).  PQL might
not be giving brilliant estimates of the random effects, but I presume you
are most interested in snps_1.  Here is one result from a simulation of
mine for a binomial trait, using R GLMM packages, augmented by your
PQL

head(x)
  ped id fa mo sex trait locus
1   1  1 NA NA   m  <NA>   1/2
2   1  2 NA NA   f  <NA>   1/2
3   1  3  1  2   f     n   1/2
4   1  4  1  2   f     y   1/1
5   1  5  1  2   f     n   1/1
6   2  6 NA NA   f  <NA>   1/2


library(pedigreemm)
ped <- pedigree(x$fa, x$mo, x$id)
pedigreemm(trait=="y" ~ (1|id), pedigree=list(id=ped), family=binomial(),
data=x)
pedigreemm(trait=="y" ~ (1|id), pedigree=list(id=ped),
family=binomial(link=probit), data=x)

library(AnimalINLA)
ped2 <- x[,2:4]
ped2[is.na(ped2)] <- 0
Ainv <- compute.Ainverse(ped2)
pheno <- data.frame(id=x$id, trait=(x$trait=="y"), Individual=x$id)
pheno <- pheno[complete.cases(pheno$trait),]
m1 <- animal.inla("trait", genetic="id", Ainverse=Ainv, data=pheno,
type.data = "binomial")
summary(m1)

library(MASS)
library(kinship)
K <- kinship(x$id, x$fa, x$mo)
observed <- !is.na(x$trait)
K <- K[observed, observed]
cs.K <- corSymm(2*K[lower.tri(K)],fixed=T)
id <- as.matrix(as.factor(x$id[observed]))
colnames(id) <- "id"
cs.K <- Initialize(cs.K,data=id)
pql1 <- glmmPQL(trait ~ 1, random=~1|id,
                correlation=cs.K,
                data=x[observed,], family="binomial
summary(pql1)

library(MCMCglmm)
pheno <- data.frame(animal=x$id, sire=x$fa, dam=x$mo, y=(x$trait=="y"))
g1 <- MCMCglmm(y~1, random=~animal, pedigree=pheno[,1:3],
               family="categorical", data=pheno, verbose=FALSE)
summary(g1)
g2 <- MCMCglmm(y~1, random=~animal, pedigree=pheno[,1:3],
               family="ordinal", data=pheno, verbose=FALSE)
summary(g2)

-----------------------------
logistic-normal
                 RE SD
pedigreemm       0.847
AnimalINLA       0.671
glmmPQL          1.185
MCMCglmm         1.389


--
| David Duffy (MBBS PhD)                                         ,-_|\
| email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v