inst/test_RSKAT2.R

In lian0090/SKAT2: SKAT test for two or more random components

library(SKAT2)
data(mouse)
n=500
p=20
pheno=mouse.pheno[1:n,]
Z1=mouse.X[1:n,1:p]
X=cbind(model.matrix(~pheno$GENDER)[,-1],pheno$CageDensity)
X1=cbind(1,model.matrix(~pheno$GENDER)[,-1],pheno$CageDensity)
y=pheno$Obesity.BMI
G=tcrossprod(scale(mouse.X[1:n,],T,F))
eigenG=getEigenZd(Kd=G)
W=colmult(X,Z1)$Z

X=X1
Zt=Z1
tauRel=NULL
eigenZd=eigenG
U1=eigenZd$U1
d1=eigenZd$d1
windowtest=c("SKAT","Score")
tU1X=crossprod(U1,X)
tU1y=crossprod(U1,y)
tXX=crossprod(X)
tXy=crossprod(X,y)
tyy=crossprod(y,y)
tU1W=crossprod(U1,W)
tXW=crossprod(X,W)
tWW=crossprod(W)
tWy=crossprod(W,y)
tZtZt=crossprod(Zt)
tU1Zt=crossprod(U1,Zt)
tXZt=crossprod(X,Zt)
tyZt=crossprod(y,Zt)
tWZt=crossprod(W,Zt)
logVar=T
optimizer="bobyqa"
  out=getDL(0.5,0.5,d1=d1,n=n,tU1y=tU1y,tU1X=tU1X,tXX=tXX,tXy=tXy,tyy=tyy,tauw=0.5,kw=ncol(W),tU1W=tU1W,tXW=tXW,tWW=tWW,tWy=tWy,getQ=T,getS=T,get.tU1ehat=T,tZtZt=tZtZt,tU1Zt=tU1Zt,tXZt=tXZt,tyZt=tyZt,tWZt=tWZt,get_tSNP=F)
  
  Var=rep(0.5,3)
  names(Var)=c("var_e","taud","tauw1")
  neg2Log(Var,tU1y,tU1X,tXX,tXy,tyy,d1,n=length(y),tU1W=tU1W,tXW=tXW,tWW=tWW,tWy=tWy,kw=ncol(W),logVar=F)