tests/mcseTest.R

In glmm: Generalized Linear Mixed Models via Monte Carlo Likelihood Approximation

#some functions from glmm that we need
library(glmm)
getFamily <- glmm:::getFamily
getEk<-glmm:::getEk
addVecs<-glmm:::addVecs
tconstant <- glmm:::tconstant
addMats<-glmm:::addMats

#other functions we'll need

#here is el, likelihood
elR <- function(Y,X,eta,family.mcml,wts){
	family.mcml<-getFamily(family.mcml)
	neta<-length(eta)
	ntrials <- rep(1, neta)

	if(family.mcml$family.glmm=="bernoulli.glmm"){
		foo<-.C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), cumout=double(1))$cumout
		mu<-.C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(1), ntrials=as.integer(ntrials), cpout=double(neta))$cpout
		cdub<-.C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(1), ntrials=as.integer(ntrials), cppout=double(neta))$cppout
	}
	if(family.mcml$family.glmm=="poisson.glmm"){
		foo<-.C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(2), ntrials=as.integer(ntrials), wts=as.double(wts),cumout=double(1))$cumout
		mu<-.C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(2), ntrials=as.integer(ntrials), cpout=double(neta))$cpout
		cdub<-.C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(2), ntrials=as.integer(ntrials), cppout=double(neta))$cppout
	}

	value<-as.numeric(Y%*%eta-foo)
	gradient<-t(X)%*%(Y-mu)	
	cdubmat<-diag(cdub)
	hessian<-t(X)%*%(-cdubmat)%*%X
	
	list(value=value,gradient=gradient,hessian=hessian)
}

#for calc t
tdist2<-function(tconst,u, Dstarinv,zeta,myq){
	inside<-1+t(u)%*%Dstarinv%*%u/zeta
	logft<-tconst - ((zeta+myq)/2)*log(inside)
	as.vector(logft)
}


distRandGeneral<-function(uvec,mu,Sigma.inv){
	logDetSigmaInv<-sum(log(eigen(Sigma.inv,symmetric=TRUE)$values))
	umu<-uvec-mu
	piece2<-t(umu)%*%Sigma.inv%*%umu
	out<-as.vector(.5*(logDetSigmaInv-piece2))
	const<-length(uvec)*.5*log(2*pi)
	out<-out-const
	out
}


distRand <-
function(nu,U,z.list,mu){
	# T=number variance components
	T<-length(z.list)
	
	#nrandom is q_t
	nrand<-lapply(z.list,ncol)
	nrandom<-unlist(nrand)
	totnrandom<-sum(nrandom)
	
	mu.list<-U.list<-NULL
	if(T==1) {
		U.list[[1]]<-U
		mu.list[[1]]<-mu
		}

	if(T>1){
		U.list[[1]]<-U[1:nrandom[1]] 
		mu.list[[1]]<-mu[1:nrandom[1]]
		for(t in 2:T){
			thing1<-sum(nrandom[1:t-1])+1
			thing2<-sum(nrandom[1:t])
			U.list[[t]]<-U[thing1:thing2]
			mu.list[[t]]<-mu[thing1:thing2]
		}
	}
	
	val<-gradient<-Hessian<-rep(0,T)
	
	#for each variance component
	for(t in 1:T){
		you<-as.vector(U.list[[t]])
		mew<-as.vector(mu.list[[t]])
		Umu<-(you-mew)%*%(you-mew)
		val[t]<- -length(U)*log(2*pi)/2+as.numeric(-.5*nrandom[t]*log(nu[t])-Umu/(2*nu[t]))
		
		gradient[t]<- -nrandom[t]/(2*nu[t])+Umu/(2*(nu[t])^2)
		
		Hessian[t]<- nrandom[t]/(2*(nu[t])^2)- Umu/((nu[t])^3)
		
	}
		
	value<-sum(val)
	if(T>1) hessian<-diag(Hessian)
	if(T==1) hessian<-matrix(Hessian,nrow=1,ncol=1)
	
	list(value=value,gradient=gradient,hessian=hessian)		
}

mcseTEST <- function(mod){



	beta <- mod$beta
	nu <-mod$nu
	umat<-mod$umat
	m<-nrow(umat)

	x <- mod$x
	y <- mod$y
	Z=do.call(cbind,mod$z)
	T<-length(mod$z)
	nrand<-lapply(mod$z,ncol)
	nrandom<-unlist(nrand)
	
	if(is.null(mod$weights)){
	  wts <- rep(1, length(y))
	} else{
	  wts <- mod$weights
	}


	family.glmm<-getFamily(mod$family.glmm)
	if(family.glmm$family.glmm=="bernoulli.glmm"){family_glmm=1}	
	if(family.glmm$family.glmm=="poisson.glmm"){family_glmm=2}	
	if(family.glmm$family.glmm=="binomial.glmm"){family_glmm=3}	

	pvec <-mod$pvec
	p1 <- pvec[1]
	p2 <- pvec[2]
	p3 <- pvec[3]
	zeta <- mod$zeta
	beta.pql <- mod$beta.pql
	nu.pql <- mod$nu.pql



	#need to get D*
	eek<-getEk(mod$z)
	Aks<-Map("*",eek,nu.pql)
	D.star<-addVecs(Aks) 
	D.star<-diag(D.star)
	D.star.inv<-solve(D.star)

	#need to recreate the variance matrix of  imp sampling distribution
	
	eta.star<-as.vector(x%*%beta.pql+Z%*%mod$u.pql)
	cdouble<-bernoulli.glmm()$cpp(eta.star) #still a vector
	cdouble<-diag(cdouble)
	Sigmuh.inv<- t(Z)%*%cdouble%*%Z+D.star.inv
	Sigmuh<-solve(Sigmuh.inv)

	Dstinvdiag<-diag(D.star.inv)

	tconst<-tconstant(zeta,nrow(D.star.inv),Dstinvdiag)

	eta<-b<-rep(0,m)
	lfu<-lfyu<-list(rep(c(0,0,0),m))
	lfu.twid<-matrix(data=NA,nrow=m,ncol=4)

	#stuff for the numerator of v-hat
	npar <- length(beta)+length(nu)
	Gpiece<-rep(0,npar)
	squaretop <- b
	numpiece <- matrix(data=NA, nrow=m, ncol=npar*npar)
	
	#for each simulated random effect vector
	for(k in 1:m){
		Uk<-umat[k,]  #use the simulated vector as our random effect vec
		eta<-mod$x%*%beta+Z%*%Uk # calculate eta using it
		zeros<-rep(0,length(Uk))

		#log f_theta(u_k)
		lfu[[k]]<-distRand(nu,Uk,mod$z,zeros) 

		#log f_theta(y|u_k)
		lfyu[[k]]<-elR(mod$y,mod$x,eta,family.glmm, wts) 

		#log f~_theta(u_k)
		lfu.twid[k,1]<-tdist2(tconst,Uk,D.star.inv,zeta=zeta,myq=nrow(D.star.inv))
		lfu.twid[k,2]<-distRandGeneral(Uk,mod$u.pql,D.star.inv)
		lfu.twid[k,3]<-distRandGeneral(Uk,mod$u.pql,Sigmuh.inv)
		
		tempmax<-max(lfu.twid[k,1:3])
		blah<-exp(lfu.twid[k,1:3]-tempmax)
		pea<-c(p1,p2,p3)
		qux<-pea%*%blah
		lfu.twid[k,4]<-tempmax+log(qux)
		
		b[k]<-as.numeric(lfu[[k]]$value)+as.numeric(lfyu[[k]]$value)-lfu.twid[k,4]

		#things for the numerator now
		Gpiece <- c(lfyu[[k]]$gradient,lfu[[k]]$gradient)
		GGT <- Gpiece %*% t(Gpiece)	

		squaretop <- exp(	2*as.numeric(lfu[[k]]$value) + 2*as.numeric(lfyu[[k]]$value) - 2*lfu.twid[k,4] )
	
		numpiece[k,] <- as.numeric(GGT*squaretop)
	}

	#return to denominator of Vhat
	#bk are unnormalized log weights
	a<-max(b) #biggest logftheta(uk,y)-logftwiddle
	thing<-exp(b-a) 
	mygamma<- exp(a)*sum(thing)/m
	Vdenom <- mygamma^2 	#denom of V hat
	#now the denominator is finished

	
	#return to numerator of Vhat

	num<-apply( numpiece   ,2,sum)/m
	temp<-num/Vdenom
	Vhat<- matrix(temp,nrow=npar,ncol=npar)
	Vhat 

	Uhat <- mod$loglike.hessian
	Uhatinv <- qr.solve(Uhat)

	UinvVUinv <- Uhatinv %*% Vhat %*% Uhatinv

	MCSE <- sqrt(diag(UinvVUinv)/m)
	MCSE
}


data(BoothHobert)
clust <- makeCluster(2)
set.seed(123)
mod.mcml1<-glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=1000, doPQL=TRUE, cluster=clust)
all.equal(mcseTEST(mod.mcml1), as.numeric(mcse(mod.mcml1)))

stopCluster(clust)