[R-sig-ME] Fwd: from general r-help to mixed help group: please help: extract p-value from mixed model in kinship package

[Ben Bolker]

Ram H. Sharma <sharma.ram.h at ...> writes:


> > Dear R experts
> >
> > I was using kinship package to fit mixed model with kinship matrix.
> > The package looks like lme4, but I could find a way to extract p-value
> > out of it. I need to extract is as I need to analyse large number of
> > variables (> 10000).
> >
> > Please help me:
> >
  
require(kinship)

##Generating random example  data

#********************pedigree data*****************************

 
id <- 1:100
 
dadid <- rep(c(0,seq(1,21,by=2)),rep(c(5,10),c(4,8)))
momid <- rep(c(0,seq(2,22,by=2)),rep(c(5,10),c(4,8)))

ped <- data.frame(id, dadid, momid)
 
# *****************************kmatrix**************************************

cfam <- makefamid(ped$id,ped$momid, ped$dadid)

kmat <- makekinship(cfam, ped$id, ped$momid, ped$dadid)

#*****************************************x and y variables

set.seed(3456)

dat <- sample(c(-1,0,1), 10000, replace = TRUE)

snpmat<- data.frame(matrix(dat, ncol = 100))

names(snpmat) <- c(paste ("VR",1:100, sep='' ))

yvar <- rnorm(100, 30, 5)
covtrait <-  rnorm(100, 10, 5)

mydata <- data.frame(id, yvar, covtrait, snpmat)

#******************************mixed model in lmekin


fmod <- lmekin(yvar ~ covtrait , data= mydata, random = ~1|id,
               varlist=list(kmat))

## note warning message ???

## $coefficients[2,4] # does not work

Your best bet in this situation is to look at (1) print.lmekin
and (2) str(fmod) to see what is present in the object and how
the package is extracting the information when it prints it.
This will show you that the information you want (the table
of parameter estimates and Wald p-values) is contained in
a list element called "ctable" so fmod$ctable[2,4]
is what you want.

  I do think it is most helpful when package authors make
this information accessible: the (or at least a) standard
convention is to have a summary()  method that returns an
object of class summary.xxx (in this case summary.lmekin)
which in turn has a coef() method to extract the coefficient
matrix with standard errors and p-values, so that
coef(summary(fmod)) would get what you wanted (that is
not implemented here).


#*******************************ultimate target: to put in

>

nvars <- ncol(mydata)-2
P <- numeric(nvars)
for (i in seq(nvars)) {
  P[i] <- lmekin(yvar~ mydata[,i+2] , data= mydata, random = ~1|id,
     varlist=list(kmat))$ctable[2,4]
}

You might want to watch out for the warning messages: not
being familar with the package, I don't know what they signify.
But ignoring warning messages is always dangerous.