glmmEP: Approximate frequentist inference for generalized linear...

Description Usage Arguments Author(s) References Examples

View source: R/glmmEP.r

Description

The machine learning paradigm known as expectation propagation, typically used for approximation inference in Bayesian graphical models contexts, is adapted to the problem of frequentist inference in generalized linear mixed model analysis. The essence is rapid approximation evaluation of the log-likelihood surface by using expectation propagation projections onto the Multivariate Normal family that circumvent the need for numerical integration. The approach applies to any reasonable dimension of the random effects vector, but currently is limited to probit mixed models with one level of nesting. The reference below provides full details.

Usage

1
glmmEP(y,Xfixed=NULL,Xrandom,idNum,control)

Arguments

y

Vector containing the (0/1) response data.

Xfixed

Matrix with columns consisting of predictors that enter the model as fixed effects. Each entry of the first column equal the number 1 corresponding to the fixed effects intercept.

Xrandom

Matrix with columns consisting of predictors that enter the model as random effects. Each entry of the first column must equal the number 1 corresponding to the random effects intercept.

idNum

Vector containing group membership identification numbers.

control

Function for controlling convergence and other specifications.

Author(s)

Matt Wand[email protected] and James Yu[email protected]

References

Hall, P., Johnstone, I.M., Ormerod, J.T., Wand, M.P. and Yu, J. (2018). Fast and accurate binary response mixed model analysis via expectation propagation. <arXiv:1805.08423v1>.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
library(glmmEP)

# Simulate data corresponding to a
# a simple random intercept model:

set.seed(1)  ; m <- 100 ; n <- 10
beta0True <- 0.37 ; beta1True <- 0.93 ; sigmaTrue <- 0.73
uTrue <- rnorm(m)*sigmaTrue
x <- runif(m*n)
y <- rbinom(m*n,1,pnorm(beta0True + beta1True*x 
            + crossprod(t(kronecker(diag(m),rep(1,n))),uTrue)))
idNum <- rep(1:m,each=n)
Xfixed <- cbind(1,x)
Xrandom <- matrix(1,length(y),1)
colnames(Xfixed) <- c("intercept","x")

# Obtain and summarise glmmEP() fit:

fit <- glmmEP(y,Xfixed,Xrandom,idNum)
summary(fit)

# Obtain simulated data corresponding to the
# simulation study in Section 4.1.2. of 
# Hall et al. (2018):

dataObj <- glmmSimData(seed=54321)
y <- dataObj$y  
Xfixed <- dataObj$Xfixed
Xrandom <- dataObj$Xrandom  
idNum <- dataObj$idNum

# Obtain the probit mixed model fit and summaries:


fitSimData <- glmmEP(y,Xfixed,Xrandom,idNum)
summary(fitSimData)
plot(fitSimData$randomEffects) ; abline(h=0) ; abline(v=0)


if (require("mlmRev"))
{
   # Example involving the data frame `Contraception'
   # in the package `mlmRev', starting with
   # setting up the input vectors and matrices
   # for a random intercepts model:

   data(Contraception)
   y <- as.numeric(as.character(Contraception$use)=="Y")
   Xfixed <- cbind(1,as.numeric(as.character(Contraception$urban)=="Y"),
                   Contraception$age,
                   as.numeric(as.character(Contraception$livch)=="1"),
                   as.numeric(as.character(Contraception$livch)=="2"),
                   as.numeric(as.character(Contraception$livch)=="3+"))
   colnames(Xfixed) <- c("intercept","isUrban","age","livChEq1",
                         "livChEq2","livChGe3")
   Xrandom <- as.matrix(rep(1,length(y)))
   colnames(Xrandom) <- "intercept"
   idNum <- as.numeric(as.character(Contraception$district))

   # Obtain random intercepts probit mixed model fit and summaries:

   
   fitContracRanInt <- glmmEP(y,Xfixed,Xrandom,idNum)
   summary(fitContracRanInt)
   hist(as.numeric(fitContracRanInt$randomEffects),
        xlab="random intercepts predicted values",probability=TRUE,
        col="dodgerblue",breaks=15,main="")
   abline(v=0,col="slateblue",lwd=2)
   

   # Change `Xrandom' to correspond to a random intercepts and 
   # slopes model:

   Xrandom <- cbind(1,as.numeric(as.character(Contraception$urban)=="Y"))
   colnames(Xrandom) <- c("intercept","isUrban")

   # Obtain random intercepts and slopes probit mixed model fit and summaries:

   
   fitContracRanIntAndSlp <- glmmEP(y,Xfixed,Xrandom,idNum)
   summary(fitContracRanIntAndSlp)
   plot(fitContracRanIntAndSlp$randomEffects,
        col="dodgerblue",xlab="random intercepts predicted values",
        ylab="random slopes predicted values",bty="l")
   abline(v=0,col="slateblue",lwd=2); abline(h=0,col="slateblue",lwd=2)
   
}

glmmEP documentation built on May 29, 2018, 9:04 a.m.