Ram H. Sharma
2011-Apr-15 13:30 UTC
[R] no solution yet, please help: extract p-value from mixed model in kinship package
I am making the question clear. Please help.> 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 <- c(rep(0, 5), rep(1, 5), rep(3, 5), rep(5, 5), rep(7, 10), rep(9, 10), rep(11, 10), rep(13, 10), rep(15, 10), rep(17, 10), rep(19, 10), rep(21, 10)) momid <- c(rep(0, 5), rep(2, 5), rep(4, 5), rep(6, 5), rep(8, 10), rep(10, 10), rep(12, 10), rep(14, 10), rep(16, 10), rep(18, 10), rep(20, 10), rep(22, 10) ) 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 ~ mydata[,3] , data= mydata, random = ~1|id, > varlist=list(kmat)) $coefficients[2,4] # does not work > > # **************************************error > message********************************************> Error message:Error in lmekin(yvar ~ mydata[, 3], data = mydata, random = ~1 | id, varlist = list(kmat))$coefficients[2, : incorrect number of dimensions In addition: Warning message: In coxme.varcheck(ncluster, varlist, n, gvars, groups, sparse, rescale, : Diagonal of variance matrix is not constant #**************************************ultimate target: to put in loop*******************************> Ultimately I want to put into the loop: > > for(i in 3:length(mydata)) { > > P <- vector (mode="numeric", length = 1000) > > P[i] <- lmekin(yvar~ mydata[,i] , data= mydata, random = ~1|id, > varlist=list(kmat)) $coefficients[2,4] > > } >> Same errors: I tried lme4 conventioned but did not work !> I can extract fixed effects as well as I do in lme4 > b <- fixef(fit1) >> Error in UseMethod("fixef") : > no applicable method for 'fixef' applied to an object of class "lmekin" > > > -- > > Ram H >-- Ram H [[alternative HTML version deleted]]
Juliet Hannah
2011-Apr-18 18:44 UTC
[R] no solution yet, please help: extract p-value from mixed model in kinship package
Maybe the pedigree is not set up correctly. If this is the case, the kinship matrix will not be constructed correctly. I see that in this example, the diagonal terms differ. diag(kmat) lmekin runs fine for me, and I can extract p-values with: lmekinfit <- lmekin(...) pval <- lmekinfit$ctable; On Fri, Apr 15, 2011 at 9:30 AM, Ram H. Sharma <sharma.ram.h at gmail.com> wrote:> ?I am making the question clear. Please help. > > > >> 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 <- c(rep(0, 5), rep(1, 5), rep(3, 5), rep(5, 5), rep(7, 10), > rep(9, 10), rep(11, 10), rep(13, 10), rep(15, 10), rep(17, 10), > rep(19, 10), rep(21, 10)) > > momid <- c(rep(0, 5), rep(2, 5), rep(4, 5), rep(6, 5), rep(8, 10), > rep(10, 10), rep(12, 10), rep(14, 10), rep(16, 10), rep(18, 10), > rep(20, 10), rep(22, 10) ) > > 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 ~ mydata[,3] , data= mydata, random = ~1|id, >> varlist=list(kmat)) $coefficients[2,4] # does not work >> >> # **************************************error >> message******************************************** > > > >> Error message: > > Error in lmekin(yvar ~ mydata[, 3], data = mydata, random = ~1 | id, varlist > = list(kmat))$coefficients[2, ?: > ?incorrect number of dimensions > In addition: Warning message: > In coxme.varcheck(ncluster, varlist, n, gvars, groups, sparse, rescale, ?: > ?Diagonal of variance matrix is not constant > > #**************************************ultimate target: to put in > loop******************************* > >> Ultimately I want to put into the loop: >> >> for(i in 3:length(mydata)) { >> >> P <- vector (mode="numeric", length = 1000) >> >> P[i] <- lmekin(yvar~ mydata[,i] , data= mydata, random = ~1|id, >> varlist=list(kmat)) $coefficients[2,4] >> >> } >> > > >> Same errors: I tried lme4 conventioned but did not work ! > > > >> I can extract fixed effects as well as I do in lme4 >> ?b <- fixef(fit1) >> > > >> Error in UseMethod("fixef") : >> ? no applicable method for 'fixef' applied to an object of class "lmekin" >> >> >> -- >> >> Ram H >> > > > > -- > > Ram H > > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >