# The GLMM
#Dr.ELDIRDIRI FADOL IBRAHIM FADOL
n <- 250
id <- seq(n)
day <- 1:14
d <- expand.grid(id = id, day = day)
set.seed(1)
trt <- sample(c("control", "treat"), size = n, replace = TRUE)
sex <- sample(c("female", "male"), size = n, replace = TRUE)
d$trt <- trt[d$id]
d$sex <- sex[d$id]
d <- d[order(d$id, d$day),]
rownames(d) <- NULL
head(d, n = 15)
## id day trt sex
## 1 1 1 control female
## 2 1 2 control female
## 3 1 3 control female
## 4 1 4 control female
## 5 1 5 control female
## 6 1 6 control female
## 7 1 7 control female
## 8 1 8 control female
## 9 1 9 control female
## 10 1 10 control female
## 11 1 11 control female
## 12 1 12 control female
## 13 1 13 control female
## 14 1 14 control female
## 15 2 1 treat male
rbinom(n = 5, size = 1, prob = 0.5)
## [1] 1 1 1 1 0
d$trtsex <- interaction(d$trt, d$sex)
probs <- c(0.40, 0.85, 0.30, 0.50)
names(probs) <- levels(d$trtsex)
d$p <- probs[d$trtsex]
set.seed(3)
r_probs <- rnorm(n = n, mean = 0, sd = 0.03)
d$random_p <- r_probs[d$id]
d$p <- d$p + d$random_p
d$y <- rbinom(n = nrow(d), size = 1, prob = d$p)
head(d[c("id", "day", "trt", "sex", "p", "y")])
## id day trt sex p y
## 1 1 1 control female 0.371142 0
## 2 1 2 control female 0.371142 0
## 3 1 3 control female 0.371142 1
## 4 1 4 control female 0.371142 0
## 5 1 5 control female 0.371142 0
## 6 1 6 control female 0.371142 0
head(d[d$id == 5, c("id", "day", "trt", "sex", "p", "y")])
## id day trt sex p y
## 57 5 1 treat female 0.8558735 1
## 58 5 2 treat female 0.8558735 1
## 59 5 3 treat female 0.8558735 1
## 60 5 4 treat female 0.8558735 1
## 61 5 5 treat female 0.8558735 1
## 62 5 6 treat female 0.8558735 1
library(lme4)
## Loading required package: Matrix
m <- glmer(y ~ trt * sex + (1|id), data = d, family = binomial)
summary(m, corr = FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt * sex + (1 | id)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 4132.4 4163.3 -2061.2 4122.4 3495
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6008 -0.8434 0.3774 0.9704 1.6878
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.03944 0.1986
## Number of obs: 3500, groups: id, 250
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.30840 0.07127 -4.327 1.51e-05 ***
## trttreat 2.18835 0.12423 17.615 < 2e-16 ***
## sexmale -0.63440 0.10859 -5.842 5.16e-09 ***
## trttreat:sexmale -1.21669 0.16556 -7.349 2.00e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plogis(-0.30840)
## [1] 0.4235053
plogis(1.87995)
## [1] 0.8676054
plogis(-0.9428)
## [1] 0.2803351
plogis(0.02886)
## [1] 0.5072145
predict(m, type = "response",
re.form = NA,
newdata = data.frame(sex = c("female", "female",
"male", "male"),
trt = c("control", "treat",
"control", "treat"))
)
## 1 2 3 4
## 0.4235043 0.8676050 0.2803350 0.5072137
exp(2.18835)
## [1] 8.920482
exp(confint(m, parm = "trttreat"))
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## trttreat 7.023558 11.44831
odds1 <- 0.001/(1 - 0.001)
odds2 <- 0.009/(1 - 0.009)
odds2/odds1
## [1] 9.072654
VarCorr(m)
## Groups Name Std.Dev.
## id (Intercept) 0.1986
confint(m, parm = "theta_", oldNames = FALSE)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## sd_(Intercept)|id 0 0.3453926
d$p2 <- plogis(-0.4054651 +
2.140066*(d$trt == "treat") +
-0.4418328*(d$sex == "male") +
-1.292768*(d$trt == "treat")*(d$sex == "male") +
rnorm(n = n, mean = 0, sd = 0.03)[d$id])
d$y2 <- rbinom(n = nrow(d), size = 1, prob = d$p2)
sim <- simulate(m, nsim = 100)
obs <- aggregate(d$y, by = list(trt = d$trt, sex = d$sex), mean)
obs
## trt sex x
## 1 control female 0.4242424
## 2 treat female 0.8659341
## 3 control male 0.2820823
## 4 treat male 0.5071429
s_p <- lapply(sim, aggregate,
by = list(trt = d$trt,sex = d$sex),
mean)
s_df <- do.call(rbind, s_p)
library(ggplot2)
ggplot() +
geom_point(aes(y = x, x = trt), data = obs,
size = 4, color = "blue") +
geom_jitter(aes(y = x, x = trt), data = s_df,
width = 0.2, height = 0, shape = 1, alpha = 1/2) +
facet_wrap(~sex) +
labs(x = "Treatment", y = "Predicted Probability",
title = "Model with interaction")
m2 <- glmer(y ~ trt + sex + (1|id), data = d, family = binomial)
sim2 <- simulate(m2, nsim = 100)
s_p2 <- lapply(sim2, aggregate,
by = list(trt = d$trt, sex = d$sex),
mean)
s_df2 <- do.call(rbind, s_p2)
ggplot() +
geom_point(aes(y = x, x = trt), data = obs,
size = 4, color = "blue") +
geom_jitter(aes(y = x, x = trt), data = s_df2,
width = 0.2, height = 0, shape = 1, alpha = 1/2) +
facet_wrap(~sex) +
labs(x = "Treatment", y = "Predicted Probability",
title = "Wrong model: without interaction")
#Eldirdiri Fadol= excuted