rm(list=ls()) #remove list library(statmod) #for gauss quadrature calculation library(lme4) nloop=1000 H=200 F=100 sigma = 0.5 beta=c(1,0.5,-1) mu=log(3/7) #grand mean total=F*H n_ij = rpois(total, lambda = 6) #sample size: a lot of 0 n_ij[n_ij>3]<-0 #sum=7500 a1=seq(-1.2,1.2,by=.1) sid = Sys.getenv('SLURM_ARRAYID') sid = as.numeric(sid) gamma_1 <- a1[sid] output.data = NULL output.data1 = NULL loop=0 repeat{ NQ = 20 # number of quadrature points alpha = rnorm(H, mean=0, sd=sigma) #alpha is hospital effect gamma = rnorm(F-1, mean=0, sd=0.4) #gamma is facility effect gamma <- c( gamma_1,gamma) temp=expand.grid(gamma,alpha) # Create a data frame from all combinations of the supplied vectors or factors. gamma=temp[,1] alpha=temp[,2] alpha= rep(alpha,n_ij) gamma= rep(gamma,n_ij) temp=expand.grid(1:F,1:H) facility = temp[,1] hospital = temp[,2] try1=cbind(temp, n_ij) m=as.vector(sapply(split(n_ij,facility),sum)) m2=as.vector(sapply(split(n_ij, hospital), sum)) facility = rep(facility, n_ij) hospital = rep(hospital, n_ij) Z = cbind( rnorm(sum(n_ij),mean=0,sd=1), rnorm(sum(n_ij),mean=0,sd=1), rnorm(sum(n_ij),mean=0,sd=1) ) prob=exp(mu+alpha+gamma+Z%*%beta) prob=prob/(1+prob) Y = rbinom(n=sum(n_ij), size=1, prob=prob ) count <- as.vector(sapply(split(Y,facility),sum)) if( (sum(count == 0) + sum(count == m)) == 0 ) { #if1 loop <- loop + 1 ################# end of generating data ############################ GAMMA.HAT = NULL BETA.HAT = NULL SIGMA.HAT=NULL ALPHA.HAT=NULL V.HAT=NULL gamma.hat = rep(0,F) GAMMA.HAT= cbind(GAMMA.HAT,gamma.hat) gamma.hat=expand.grid(gamma.hat,1:H)[,1] gamma.hat= rep(gamma.hat,n_ij) beta.hat = c(1,1,1) ##### BETA.HAT=cbind(BETA.HAT,beta.hat) sigma.hat= 1 ###### SIGMA.HAT=cbind(SIGMA.HAT,sigma.hat) ptm <- proc.time() ########### start loop repeat{ ghq=gauss.quad.prob(NQ,"normal", sigma=sigma.hat) Z.beta.hat = Z%*%beta.hat ########### ghq.M=NULL for(i in 1:NQ){ ghq.q = 1/(1+exp(ghq$nodes[i]+gamma.hat+Z.beta.hat)) #nodes: at which evaluate the function ghq.p = 1-ghq.q ghq.p.or.q = Y*ghq.p + (1-Y)*ghq.q ghq.M = cbind(ghq.M, sapply( split( log(ghq.p.or.q), hospital ), sum) ) ##end cbind } ghq.M=ghq.M-apply(ghq.M,1,max) ghq.M=exp(ghq.M) alpha.hat = ghq.M %*% (ghq$nodes*ghq$weights) / ghq.M %*% ghq$weights ALPHA.HAT=cbind(ALPHA.HAT,alpha.hat) alpha.hat = expand.grid(1:F,alpha.hat)[,2] alpha.hat = rep(alpha.hat,n_ij) v.hat = ghq.M %*% (ghq$nodes^2*ghq$weights) / ghq.M %*% ghq$weights - (ghq.M %*% (ghq$nodes*ghq$weights) / ghq.M %*% ghq$weights)^2 V.HAT=cbind(V.HAT,v.hat) v.hat = expand.grid(1:F,v.hat)[,2] v.hat = rep(v.hat,n_ij) sigma.hat=sqrt( mean( ghq.M %*% (ghq$nodes^2*ghq$weights) / ghq.M %*% ghq$weights ) ) SIGMA.HAT=cbind(SIGMA.HAT,sigma.hat) q.hat = 1/(1+exp(alpha.hat+gamma.hat+Z.beta.hat)) p.hat = 1-q.hat u.gamma = Y-p.hat - 0.5*v.hat*p.hat*q.hat*(q.hat-p.hat) i.gamma = p.hat*q.hat + .5*v.hat*p.hat*q.hat*(p.hat^2+q.hat^2-4*p.hat*q.hat) gamma.update = sapply( split(u.gamma,facility), sum) / sapply( split(i.gamma,facility), sum) GAMMA.HAT = cbind(GAMMA.HAT, GAMMA.HAT[,ncol(GAMMA.HAT)] + gamma.update) gamma.update = expand.grid(gamma.update,1:H)[,1] gamma.update = rep(gamma.update,n_ij) gamma.hat = gamma.hat+gamma.update q.hat = 1/(1+exp(alpha.hat+gamma.hat+Z.beta.hat)) p.hat = 1-q.hat u.beta = t( Y-p.hat- 0.5*v.hat*p.hat*q.hat*(q.hat-p.hat) ) %*% Z i.beta = t(Z)%*%( Z* c( p.hat*q.hat +0.5*v.hat*p.hat*q.hat*(p.hat^2+q.hat^2-4*p.hat*q.hat) ) ) beta.hat=beta.hat+ solve(i.beta)%*%t(u.beta) BETA.HAT=cbind(BETA.HAT,beta.hat) beta.distance = BETA.HAT[,ncol(BETA.HAT)-1]-BETA.HAT[,ncol(BETA.HAT)] gamma.dis = GAMMA.HAT[,ncol(GAMMA.HAT)-1]-GAMMA.HAT[,ncol(GAMMA.HAT)] if(max(max(abs(beta.distance)), max(abs(gamma.dis)))<1e-9) break } sd.v.hat=sqrt(v.hat) B <- 10000 get.sl <- function(gamma) { alpha.mean<-alpha.hat[facility==1] alpha.sd<-sd.v.hat[facility==1] Z_facility<-cbind(Z[,1][facility==1], Z[,2][facility==1], Z[,3][facility==1]) prob_B<-NULL kk<-0 repeat{ kk<-kk+1 pre_prob_B=plogis(gamma +rnorm(m[1], mean=alpha.mean, sd=alpha.sd)+ Z_facility%*%beta.hat) prob_B=rbind(prob_B,pre_prob_B) if(kk>=B) break } Y.star <- rbinom(n=m[1]*B, size = 1, prob =prob_B ) Y.star <- matrix(Y.star,ncol=B) Y.star.sum <- apply(Y.star,2,sum) Y.sum = sum(Y[facility==1]) return ( 2*min( mean(Y.star.sum >= Y.sum), mean(Y.star.sum <= Y.sum) ) ) } proposed.gamma <- mu+gamma_1 sl <- sapply(proposed.gamma, get.sl) gamma_1_true=mu+gamma_1 CP_EM <- (sl>=0.05) # t(ranef(lmer.fit)$h)[1] proc.time() - ptm my.data=data.frame(Y=Y,z1=Z[,1],z2=Z[,2],z3=Z[,3],f=facility,h=hospital) lmer.fit<-lmer(Y~factor(f)+z1+z2+z3+(1|h)-1, family=binomial(link=logit),data=my.data, REML=FALSE) proc.time() - ptm gamma1_PQL<-coef(summary(lmer.fit))[1,1] gamma1_PQL_sd<-coef(summary(lmer.fit))[1,2] c.i.r.PQL<-gamma1_PQL+1.9645*gamma1_PQL_sd c.i.l.PQL<-gamma1_PQL-1.9645*gamma1_PQL_sd gamma_1_true=mu+gamma_1 CP_PQL <- ( gamma_1_true <= c.i.r.PQL ) & ( gamma_1_true >= c.i.l.PQL ) #fixef(lmer.fit) #GAMMA.HAT[,ncol(GAMMA.HAT)] #BETA.HAT[,ncol(BETA.HAT)] #round(ranef(lmer.fit)$h - ALPHA.HAT[,ncol(ALPHA.HAT)], digits = 2 ) # 1 indicates the approximate EM method EM.est1=cbind(t(GAMMA.HAT[1,ncol(GAMMA.HAT)]), CP_EM) # 2 indicates the lmer method LMER.est1=cbind(t(fixef(lmer.fit)[1]), gamma1_PQL_sd, CP_PQL) output.data1_pre=cbind(EM.est1,LMER.est1) output.data1 = rbind(output.data1, output.data1_pre) } #end if1 if(loop>=nloop) break } output.data1=as.data.frame(output.data1) colnames(output.data1) = c("Gamma1", "CP_EM","Gamma2_1", "gamma1_PQL_sd", "CP_PQL") gamma1_bias_EM=mean(output.data1$Gamma1-gamma_1_true) gamma1_bias_PQL=mean(output.data1$Gamma2_1-gamma_1_true) gamma1_CP_EM=mean(output.data1$CP_EM) gamma1_ASE_PQL=mean(output.data1$gamma1_PQL_sd) gamma1_CP_PQL=mean(output.data1$CP_PQL) gamma_1_true=mu+gamma_1 gamma1_sd_EM=sd(output.data1$Gamma1) gamma1_sd_PQL=sd(output.data1$Gamma2_1) mixed_summary=cbind(gamma_1_true, gamma1_bias_EM,gamma1_bias_PQL, gamma1_sd_EM,gamma1_sd_PQL, gamma1_ASE_PQL, gamma1_CP_EM, gamma1_CP_PQL) file.name = paste("/home/kevinhe/mix_100/mixed-", sid, ".csv",sep="") write.csv(mixed_summary, file=file.name) file.name2 = paste("/home/kevinhe/mix_100/mixed-data-", sid, ".csv",sep="") write.csv(output.data1, file=file.name2)