[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
bbolker at gmail.com
Sat Apr 16 17:08:20 CEST 2011
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.
More information about the R-sig-mixed-models
mailing list