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.