library(compiler) library(nnet) #Sample code for Wang Y, and Chen H (2012). Generating type I error for scenario (a) in Table 1. # critical value of the F test seed.vec<-c(1:50) thresF=function(nsim=10000, n=50, nsub=10, p=2, K=18, seed){ set.seed(seed.vec[seed]) ###MODIFIED by YW ##computate the egivenvalues xi xi=matrix(rep(0,nsub*ngrid),nrow=ngrid) In=diag(rep(1,n)) one=rep(1,nsub) for(i in 1:ngrid){ lam1=as.double(lam.set[i]) inVlam1=(In-lam1*Z1%*%solve(diag(rep(1,K))+lam1*t(Z1)%*%Z1)%*%t(Z1))/sigma2 Plam1=(inVlam1-inVlam1%*%X%*%solve(t(X)%*%inVlam1%*%X)%*%t(X)%*%inVlam1) xi[i,]=as.double(svd(t(Z2)%*%Plam1%*%Z2)$d) } P0=diag(rep(1, n))-X%*%solve(t(X)%*%X)%*%t(X) mu=svd(t(Z1)%*%P0%*%Z1)$d vect=log((1+lam.set%*%t(mu)))%*%rep(1,K) Ft=rep(0, nsim) REML=matrix(rep(0,ngrid^2), nrow=ngrid) print(date()) for(k in 1: nsim){ u1=rchisq(nsub,1) u2=rchisq(1, df=n-nsub-p) #choose lam1 and lam2 via REML for(j in 1:ngrid){ lam2=as.double(lam.set[j]) REML[,j]=fREMLc(u1, u2, xi, lam2, one, n, p, vect) } mcol=max.col(-REML) minvect=cbind((1:ngrid),mcol) temp=REML[minvect] i=which.is.max(-temp) j=mcol[i] lam2=lam.set[j] Ft[k]=as.double((sum(lam2*xi[i,]*u1/(1+lam2*xi[i,])))*n/(sum(u1/(1+lam2*xi[i,]))+u2)) } thres=quantile(Ft, 0.95) print(date()) print(thres) return(thres) } fREML=function(u1, u2, xi, lam2, one, n, p, vect){ as.double((n-p)*log((1/(1+lam2*xi))%*%u1+u2)+(log(1+lam2*xi))%*%one+vect) } fREMLc=cmpfun(fREML) ## knots function default.knots <- function(x,num.knots){ if (missing(num.knots)) num.knots <- max(5,min(floor(length(unique(x))/4),35)) return(quantile(unique(x),seq(0,1,length=(num.knots+2))[-c(1,(num.knots+2))])) } ## test the random intercept with a nuisance smooth function ## n is the number of subjects; m is the number of measurements in each subject; K is the number of knots; ## p-1 is the order of the spline basis; ## ngrid is the number of grid points that used to search the smoothing parameter ## tsd is the standard deviation of the variance component corresponding to the smooth function (nuisance variance component) ## tsd2=0 is the standard deviation of the random intercept (variance component to be tested) ## nsim is the number of simulations library(nlme) library(RLRsim) ### settins to be changed seed=935493; tsd=1; n=10; m=5; set.seed(seed) N=n*m; K=min(N/10, 35); p=3; ngrid=100; tsd2=0; lam.set=c(0,exp(seq(-10,8,length=(ngrid-1)))) print(date()) nsim=5000 beta=c(1,-.2,.2) ID=sort(rep((1:n),m)) power=powerL=powerm=0 step1=0 while(step1<50){ T=runif(N) X1=runif(N) for(i in 1:n){T[ID==i]=sort(T[ID==i])} Z2=matrix(0, N, n) for(i in 1:n){ Z2[((i-1)*m+1):(i*m),i]=rep(1,m)} knots1 <- default.knots(T, K) X2=cbind(rep(1,N), T) Z1 <- outer(T, knots1,"-") Z1 <- (Z1*(Z1>0)) X=cbind(X1,X2) Y=X%*%beta+Z1%*%(rnorm(K,sd=tsd))+rnorm(N) # critical value of the F test group1=group=rep(1,N) fitd=try(lme(Y~-1+X, random=list(group=pdIdent(~-1+Z1), group=pdIdent(~-1+Z2)))) if(class(fitd)!="try-error"){ sigma2=as.double(fitd$sigma^2) lam1=as.double(as.numeric(nlme:::VarCorr(fitd)[2,1])/sigma2) thres=thresF(nsim=5000, n=N, nsub=n, p=p, K=K, seed=step1+1) ###MODIFIED by YW RLRT <- RLRTSim(X=cbind(X1,X2), Z=Z2, qrX=qr(X), sqrt.Sigma=diag(ncol(Z2)), nsim=1e5) (thresL=quantile(RLRT, prob=0.95)) print(thresL) # critical vaule of the mixture Chisq thresm=2.7 step2=0 while(step2<100){ Y=X%*%beta+Z1%*%(rnorm(K,sd=tsd))+rnorm(N) fit=try(lme(Y~-1+X, random=list(group=pdIdent(~-1+Z1), group1=pdIdent(~-1+Z2)))) fit0=try(lme(Y~-1+X, random=list(group=pdIdent(~-1+Z1)))) if(all(c(class(fit), class(fit0))!="try-error")){ sigma2=fit$sigma^2 lam1=as.numeric(VarCorr(fit)[2,1])/sigma2 var2=as.numeric(VarCorr(fit)[K+3,1]) #compute the observed test statistics inVlam1=(diag(rep(1,N))-lam1*Z1%*%solve(diag(rep(1,K))+lam1*t(Z1)%*%Z1)%*%t(Z1))/sigma2 inV=inVlam1-var2*inVlam1%*%Z2%*%solve(diag(rep(1,n))+var2*t(Z2)%*%inVlam1%*%Z2)%*%t(Z2)%*%inVlam1 H1=X%*%solve(t(X)%*%inV%*%X)%*%t(X)%*%inV H0=X%*%solve(t(X)%*%inVlam1%*%X)%*%t(X)%*%inVlam1 Yt1=Y-H1%*%Y RSS1=t(Yt1)%*%inV%*%Yt1 Yt0=Y-H0%*%Y RSS0=t(Yt0)%*%inVlam1%*%Yt0 obsFt=(RSS0-RSS1)*N/RSS1 # power of the F test power=power+1*(obsFt>thres) # compute the LRT test statistics LRTobs=anova(fit,fit0)$L.Ratio[2] # power of the LRT powerL=powerL+(LRTobs>thresL)*1 # power of the mixture Chisq powerm=powerm+(LRTobs>thresm)*1 step2=step2+1 } } step1=step1+1 } print(c(step1, power/100/step1,powerL/100/step1, powerm/100/step1)) } c(tsd, tsd2, n, power/nsim, powerL/nsim, powerm/nsim) seed binom.test(drop(power), nsim, p=power/100/step1) binom.test(drop(powerL), nsim, p=powerL/100/step1) binom.test(drop(powerm), nsim, p=powerm/100/step1)