x<- sample(c("M","F"),100,replace=TRUE)
y<- rnorm(100) + (x=="M")*1.0
v<- cov(matrix(rnorm(10^6),ncol=100))
# no sex effect
o<- estVC(y, v = list(AA=v,DD=NULL,HH=NULL,AD=NULL,
MH=NULL,EE=diag(100)),method="BFGS")
o
# sex as fixed effect
fo<- estVC(y, x, v = list(AA=v,DD=NULL,HH=NULL,AD=NULL,
MH=NULL,EE=diag(100)),method="BFGS")
fo
2*(fo$value-o$value) # log-likelihood test statistic
# sex as random effect
SM<- rem(~x,data.frame(x=x))
ro<- estVC(y, v = list(AA=v,DD=NULL,HH=NULL,AD=NULL,
MH=NULL,SE=SM$x,EE=diag(100)),method="BFGS")
ro
2*(ro$value-o$value) # log-likelihood test statistic
Run the code above in your browser using DataLab