# --------------------------------------------------------------
# Estimation of POD-curve by probit binomial GLMM
# 
# Assumptions:
# The response occurs when the quantity (q1) of some material becomes larger than a threshold quantity (t1).
# The quantity (q1) fluctuates by following a normal distribution with a scaled variance (sigma^2).
# The threshold fluctuates among laboratories by following a normal distribution.
# 
# --------------------------------------------------------------

### Data from Kodama et al. (2011; J AOAC Int 94:224-231) ###

sink("Kodama.txt")
cat("
Obs x Lab r n
1  0.0005 1 6 6
2  0.0005 2 6 6
3  0.0005 3 5 6
4  0.0005 4 6 6
5  0.0005 5 6 6
6  0.0005 6 6 6
7  0.0005 7 5 6
8  0.0005 8 6 6
9  0.0005 9 6 6
10 0.0005 10 6 6
11 0.0005 11 5 6
12 0.0005 12 5 6
13 0.0005 13 6 6
14 0.0005 14 5 6
15 0.001 1 6 6
16 0.001 2 6 6
17 0.001 3 6 6
18 0.001 4 6 6
19 0.001 5 4 6
20 0.001 6 6 6
21 0.001 7 6 6
22 0.001 8 6 6
23 0.001 9 6 6
24 0.001 10 6 6
25 0.001 11 6 6
26 0.001 12 6 6
27 0.001 13 6 6
28 0.001 14 6 6
29 0 1 0 6
30 0 2 0 6
31 0 3 0 6
32 0 4 0 6
33 0 5 0 6
34 0 6 0 6
35 0 7 0 6
36 0 8 0 6
37 0 9 0 6
38 0 10 0 6
39 0 11 0 6
40 0 12 0 6
41 0 13 0 6
42 0 14 0 6
",fill=TRUE)
sink()

Kodama <- read.table("Kodama.txt", header=TRUE)
Kodama$response <- cbind(Kodama$r,Kodama$n-Kodama$r)

### GLMM ESTIMATION (threshold fluctuates among laboratories) ###

library(lme4)
Kodama.glmer <- glmer(response ~ x + (1|Lab), family=binomial(link = "probit"), data = Kodama)

sum.Kodama.glmer <- summary(Kodama.glmer)
eta <- sum.Kodama.glmer$coefficients[,1]
s2u <- sum.Kodama.glmer$varcor$Lab[1]
sum.Kodama.glmer

### MARGINAL ESTIMATION (threshold fluctuates among laboratories) ###

eta.Marginal <- eta/sqrt(1+s2u)
Estimates <- data.frame(c(round(eta[1],6),round(eta[2],6)), c(round(eta.Marginal[1],6),round(eta.Marginal[2],6)))
colnames(Estimates) <- c("GLMM estimates","Marginal estimates")
Estimates