covaruber / sommer

36 stars 23 forks source link

Query regarding MEMMA output #34

Closed jrhadams closed 2 years ago

jrhadams commented 2 years ago

I had a question regarding MEMMA's output. Reading the documentation I got the impression that I could extract variance components and for each random effect specified, but in my attempt to create a genetic model with two random effects, I only see one random component.

Here is a minimal example where I attempt to perform a single random marker test with background genetic effects.

###R and package version
>packageVersion('sommer')
[1] ‘4.1.3’
> R.version$version.string
[1] "R version 4.0.4 (2021-02-15)"

####Define your design and covariance matrices####
Z1<- apply(as.matrix(dat[,c(11:13)]),2,
             function(x) as.integer(x*2)) #your qtl effects based on 3 founders of population
K1<- diag(ncol(Z1)) #Single variance with no covariance among founder effects
dimnames(K1)<-list(colnames(Z1),colnames(Z1)) 
Z2<-model.matrix(~-1+geno,data=dat) #Your background genetic effects
K2<- Gmat # additive covariance matrix
X<- model.matrix(~1+block,data=dat) #Fixed trial effects
Z<- list(qtl=list(Z=Z1,K=K1),
           poly=list(Z=Z2,K=K2))

testMod<-tryCatch({
   MEMMA(Y=dat$yield
       X=X,ZETA=Z,silent=F)
       },
   error=function(c){
      cat(c$message,file = "log.txt",append = T)
      return(NULL)
})

Examining the variances, I get the following:

> testMod$var.comp$Vu
         T.1
T.1 26.11816

And extracting the realised effects, I get an array the same size as the original dataset.

> nrow(testMod$u.hat)==nrow(dat)
[1] TRUE

I was expecting mod$u.hat to only contain the effects from Z1 and Z2 with the length of the output being the column ranks of Z1 and Z2. I suspect that the output of u.hat is actually the sum of all random effects for each listing in dat.

My question is whether I am mistaken in my understanding of MEMMA or have I made a mistake in how I specified these models?

Thanks in advance.

covaruber commented 2 years ago

Hi, no actually the biggest weakness of the EMMA algorithm (the one used in the rrBLUP package) is that can only estimate one variance component other than the error, that's why is not part of the mmer algorithms available. I just left it there for curious people but not very useful. On top of that is an R based algorithm whereas the ones that mmer calls are coded in Rcpp.

Cheers, Eduardo