Re: [R] no solution yet, please help: extract p-value from mixed model in kinship package

From: Juliet Hannah <juliet.hannah_at_gmail.com>
Date: Mon, 18 Apr 2011 14:44:28 -0400

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.
>



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. Received on Mon 18 Apr 2011 - 18:47:52 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 18 Apr 2011 - 19:40:31 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive