---
output: pdf_document
---
\large

#Multi-center clinical trial GLMM

##WARNING: glmer (version 1.1-10) MLEs are correct, but model assessment criteria are calculated erratically (glmmML used at bottom for comparison)

The true value of the log-likelihood of the random-clinic model is -37.0313 
(verified by ADMB using importance sampling).
Using the Laplace approximation the log-likelihood of this model is evaluated to be 
-37.1006 (ADMB using Laplace approx).

Under the model that also includes the subject-level random effect (i.e., overdispersion), 
the loglikelihood is -37.02 (from ADMB using importance sampling)

## Enter the data
```{r}
fav=c(11,10,16,22,14,7,2,1,6,0,1,0,1,1,4,6)
unfav=c(25,27,4,10,5,12,14,16,11,12,10,10,4,8,2,1)
y=fav; n=fav+unfav

clinic=factor(rep(1:8,rep(2,8)))
subj=factor(1:16)
trmt=factor(rep(c("Drug","Control"),8))

#Log-likelihood of saturated model
logsatdmodel=sum(log(dbinom(y,n,y/n))) 
```

Log-likelihood of saturated model is `r round(logsatdmodel,3)`

## Clinic random effect: Laplace approx via glmer
```{r}
library(lme4)
f=glmer(y/n~trmt+(1|clinic),family="binomial",weight=n)
summary(f)
cat("logLik(f)",logLik(f),
    "\ndeviance(f)",deviance(f),
    "\nAIC(f)",AIC(f),
    "\nPearson chisq",sum(resid(f,type="pearson")^2))
```
NOTES: 

- logLik and AIC are correct

- deviance in summary output is simply -2loglikelihood. 

- deviance(f) returns a modified deviance value (possibly suitable for model checking?)

## Clinic random effect: Gaussian quadrature approx via glmer
```{r}
#Using quadrature, deviance is wrto saturated model (and logLik=Deviance/2 is wrong)
fQ=glmer(y/n~trmt+(1|clinic),family="binomial",weight=n,nAGQ=9)
summary(fQ)
cat("logLik(fQ)",logLik(fQ),
    "\ndeviance(fQ)",deviance(fQ),
    "\nAIC(fQ)",AIC(fQ),
    "\nPearson chisq",sum(resid(fQ,type="pearson")^2))
```
NOTES: 

- deviance in summary output is now correct

- logLik in summary output, and from logLik(fQ), is -half of deviance, 
and is incorrect, as is AIC

- deviance(f) returns a modified deviance value (possibly suitable for model checking?)

## Clinic and subject (i.e., overdispersion) random effects: Laplace approx via glmer
```{r}
fOD=glmer(y/n~trmt+(1|clinic)+(1|subj),family="binomial",weight=n)
summary(fOD) #Model log-likelihood= -37.1007
cat("logLik(fOD)",logLik(fOD),
    "\ndeviance(fOD)",deviance(fOD),
    "\nAIC(fOD)",AIC(fOD),
    "\nPearson chisq",sum(resid(fOD,type="pearson")^2))
```
As per earlier use of Laplace for clinic only random effect.

## Clinic and trmt-within-clinic random effects: Laplace approx via glmer
## Model includes correlation between the two random effects
```{r}
fCT=glmer(y/n~trmt+(trmt|clinic),family="binomial",weight=n) #With re correlation
summary(fCT)
cat("logLik(fCT)",logLik(fCT),
    "\ndeviance(fCT)",deviance(fCT),
    "\nAIC(fCT)",AIC(fCT),
    "\nPearson chisq",sum(resid(fCT,type="pearson")^2))
```

\pagebreak

# Using glmmML
## Clinic random effect: Laplace approx via glmmML
```{r}
library(glmmML)
glmmML(y/n~trmt,cluster=clinic,family="binomial",weight=n,method="Laplace")
```
NOTE: deviance is correct, but AIC is not

## Clinic random effect: Gaussian quadrature approx via glmmML
```{r}
library(glmmML)
glmmML(y/n~trmt,cluster=clinic,family="binomial",weight=n,method="ghq",n.points=9)
```