getLoglik: get log likelihood from an input variance - lian0090/SKAT2 GitHub Wiki

install_github("lian0090/SKAT2")
library(SKAT2)
data(mouse)
pheno=mouse.pheno
Z1=mouse.X[,1:40]
Z2=mouse.X[,41:60]
Z3=mouse.X[,61:80]
X=cbind(1,model.matrix(~pheno$GENDER)[,-1],pheno$CageDensity)
y=pheno$Obesity.BMI

Usage

getLoglik(Var,y,X,Z=list())

• Var: current values of variance. The first element of Var is for error variance, the rest are for the elements of Z.
• y: phenotype
• X: incidence matrix for fixed effects
• Z: a list of incidence matrix for random effects. If there are two random components, Z will be a list of two.

getLoglik(Var=c(0.1,0.1),y=y,X=X,Z=list(Z1))
[1] 2022.239
getLoglik(Var=c(0.1,0.1,0.1),y=y,X=X,Z=list(Z1,Z2))
[1] 2001.678
getLoglik(Var=c(0.1,0.1,0.1,0.1),y=y,X=X,Z=list(Z1,Z2,Z3))
[1] 1982.07

If there is one random component has a lot of columns, and it is repeated in a lot of computations, it is much better to directly provide it as an object from eigen decomposition to save computation. For example, if we want to use all the genotypes to control genetic background, we can do it like this:

eigenG=getEigenG(Zg=mouse.X)
getLoglik(Var=c(0.1,1e-6,0.1),y=y,X=X,Z=list(eigenG,Z1))
[1] 1995.731