R Under development (unstable) (2024-09-21 r87186 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #going to test el for the Booth and Hobert example
> 
> library(glmm)
Loading required package: trust
Loading required package: mvtnorm
Loading required package: Matrix
Loading required package: parallel
Loading required package: doParallel
Loading required package: foreach
Loading required package: iterators
> data(BoothHobert)
> clust <- makeCluster(2)
> set.seed(1234)
> mod.mcml1<-glmm(y~0+x1,list(y~0+z1),varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=21, doPQL=TRUE, debug=TRUE, cluster=clust)
> 
> mod.mcml<-mod.mcml1$mod.mcml
> z<-mod.mcml$z[[1]]
> x<-mod.mcml$x
> y<-mod.mcml$y
> 
> stuff<-mod.mcml1$debug
> beta.pql<-stuff$beta.pql
> nu.pql<-stuff$nu.pql
> u.pql<-u.star<-stuff$u.star
> umat<-stuff$umat
> 
> family.glmm<-bernoulli.glmm
> 
> objfun<-glmm:::objfun
> getEk<-glmm:::getEk
> addVecs<-glmm:::addVecs
> 
> ############################################
> #this should be the same as el
> getFamily<-glmm:::getFamily
> 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)
+ }
> 
> #compare elR and el.C for a value of eta
> neta <- 150
> eta<-rep(2,neta)
> ntrials <- rep(1,neta)
> that<-elR(mod.mcml$y,mod.mcml$x,eta,family.mcml=bernoulli.glmm,wts=rep(1,length(mod.mcml$y)))	
> this<-.C(glmm:::C_elc, as.double(mod.mcml$y), as.double(mod.mcml$x), as.integer(nrow(mod.mcml$x)), as.integer(ncol(mod.mcml$x)), as.double(eta), as.integer(1), as.integer(ntrials), wts=as.double(rep(1,length(mod.mcml$y))), value=double(1), gradient=double(ncol(mod.mcml$x)), hessian=double((ncol(mod.mcml$x)^2)))
> all.equal(as.numeric(that$value),this$value)
[1] TRUE
> all.equal(as.numeric(that$gradient),this$gradient)
[1] TRUE
> all.equal(as.numeric(that$hessian),this$hessian)
[1] TRUE
> 
> #compare to elval
> elvalout<-.C(glmm:::C_elval, as.double(mod.mcml$y), as.integer(nrow(mod.mcml$x)), as.integer(ncol(mod.mcml$x)), as.double(eta), as.integer(1), as.integer(ntrials), wts=as.double(rep(1,length(mod.mcml$y))), value=double(1))
> all.equal(as.numeric(that$value),elvalout$value)
[1] TRUE
> 
> #compare to elGH
> elGHout<-.C(glmm:::C_elGH,as.double(mod.mcml$y),as.double(mod.mcml$x),as.integer(nrow(mod.mcml$x)),as.integer(ncol(mod.mcml$x)),as.double(eta),as.integer(1), as.integer(ntrials), wts=as.double(rep(1,length(mod.mcml$y))), gradient=double(ncol(mod.mcml$x)),hessian=double((ncol(mod.mcml$x)^2)))
> all.equal(as.numeric(that$gradient),elGHout$gradient)
[1] TRUE
> all.equal(as.numeric(that$hessian),elGHout$hessian)
[1] TRUE
> 
> stopCluster(clust)
> 
> 
> proc.time()
   user  system elapsed 
   3.09    0.43    5.18