1 ######################## RECAP ###############################
   2 
   3 load("../week4/data/births.rdata")
   4 
   5 ## gaussian model
   6 m.lm <- lm(bweight ~ hyp, data=births)
   7 m.glm <- glm(bweight ~ hyp, family=gaussian, data=births)
   8 
   9 ## recap anova
  10 m.bin1 <- glm(lowbw ~ hyp, family=binomial, data=births)
  11 summary(m.bin1)
  12 
  13 ## back transforming coefficients
  14 invlogit <- function(x){ exp(x)/(exp(x) + 1)}
  15 invlogit(coef(m.bin1)[1]) 
  16 
  17 ## effects
  18 require(effects)
  19 Effect("hyp",m.bin1)
  20 
  21 ## logistic regression
  22 m.bin5 <- glm(lowbw ~ gestwks, family=binomial, data=births)
  23 summary(m.bin5)
  24 
  25 ## understanding the coefficients
  26 invlogit(coef(m.bin5)[1])
  27 exp(coef(m.bin5)[2])
  28 
  29 
  30 Effect("gestwks",m.bin5)
  31 
  32 
  33 #### O-Ring Example
  34 library(faraway)
  35 data(orings)
  36 head(orings)
  37 
  38 m.or <- glm(cbind(damage,6-damage) ~ temp,
  39                       family=binomial, orings)
  40 
  41 summary(m.or)
  42 
  43 invlogit(coef(m.or)[1])
  44 exp(coef(m.or)[2])
  45 
  46 require(effects)
  47 allEffects(m.or)
  48 
  49 tf <- 20:100
  50 pd <- predict(m.or,newdata=list(temp=tf), type="response")
  51 
  52 plot(tf,pd,type="l",
  53      xlab=expression(paste("temperature (",degree,"F)",sep=" ")),
  54      ylab="probability of damage")
  55 
  56 
  57 ## ggplot
  58 require(scales)
  59 require(ggplot2)
  60 
  61 orings$trials <- 6
  62 ggplot(orings,aes(x=temp,y=damage/trials)) +
  63     geom_point() +
  64     geom_smooth(method = "glm", family = "binomial", aes(weight = trials)) +
  65     xlab(expression(paste("temperature (",degree,"F)"))) +
  66     ylab("probability of damage") +
  67     scale_y_continuous(labels = percent)
  68 
  69 ggplot(orings,aes(x=temp,y=damage/trials)) +
  70     geom_point() +
  71     geom_smooth(method = "glm", family = "binomial", aes(weight = trials), fullrange = T) +
  72     xlab(expression(paste("temperature (",degree,"F)"))) +
  73     ylab("probability of damage") +
  74     scale_x_continuous(limits = c(0,100)) +
  75     scale_y_continuous(labels = percent)
  76 
  77 
  78 ggsave("img/challenger2.png")
  79 
  80 ## Binary Ancova
  81 infection <- read.table("../week5/data/infection.txt",header=T)
  82 summary(infection)
  83 
  84 infection$sex <- factor(infection$sex, labels = c("male","female"))
  85 infection$infected <- factor(infection$infected, labels = c("non infected","infected"))
  86 
  87 
  88 m.inf <- glm(infected~age*sex,family=binomial,
  89              data=infection)
  90 
  91 summary(m.inf)
  92 
  93 coef(m.inf)
  94 invlogit(coef(m.inf)[1])
  95 
  96 invlogit(coef(m.inf)[1]+coef(m.inf)[3])
  97 
  98 solve(0.015657,3.000513)
  99 solve(0.02670685,2.883849)
 100 
 101 
 102 allEffects(m)
 103 
 104 ###############################################################
 105 ######################  Exercise  #############################
 106 ###############################################################
 107 
 108 ## change sex and infected to factors with the labels male and
 109 ## female, resp. infected/non infected if not already done
 110 
 111 
 112 ## Try to reproduce the plot! Hints:
 113 
 114 ## set up a ggplot object, think about the aesthetics aes().
 115 ## Which quality of the graph you wanna set to which variable?
 116 
 117 
 118 
 119 ## begin with the lines (geom_smooth())
 120 
 121 
 122 
 123 
 124 ## add the points (geom\_jitter());
 125 ## do not think about the symbols in the first place;
 126 ## try to adjust the width and height appropriately
 127 
 128 
 129 
 130 
 131 ## change the colour of the lines and points
 132 ## (scale_colour_manual()); I used midnightblue
 133 ## for male and deeppink for female
 134 
 135 
 136 
 137 ## change the symbols (scale_shape_manual()); use
 138 ## values = c("male" = "\u2642","female" = "\u2640")) as values
 139 
 140 
 141 
 142 ## set the axes titles
 143 
 144 
 145 ## change to text of the y axis to percentage, etc
 146 
 147 
 148 require(ggplot2)
 149 require(scales)
 150 
 151 
 152 ggplot(infection,aes(y = as.numeric(infected=="infected"), x = age, color=sex, shape=sex))+
 153   geom_smooth(method = "glm", family = "binomial")+
 154   geom_jitter(position = position_jitter(width = 0, height= 0.08),size=5,face="bold")+
 155   scale_color_manual(values = c("male"="midnightblue", "female"= "deeppink"))+
 156   scale_shape_manual(values = c("male" = "\u2642","female" = "\u2640"))+
 157   labs(title="Infection", x = "age in days", y= "infected in percent") +
 158   scale_y_continuous(breaks=seq(0,1,by=0.25),labels = percent)
 159 
 160 
 161 
 162 
 163 ###############################################################
 164 ###################  Exercise  contra##########################
 165 ###############################################################
 166 
 167 ## load the data using
 168 cuse <- read.table("http://data.princeton.edu/wws509/datasets/cuse.dat", header=TRUE)
 169 
 170 ## description here: http://data.princeton.edu/wws509/datasets/#cuse
 171 
 172 
 173 ## model the use of contraceptiva dependent on age, education,
 174 ## and the wish for more children
 175 lrfit <- glm( cbind(using, notUsing) ~ age + education + wantsMore ,
 176              family = binomial, data = cuse)
 177 
 178 ## what are the probabilities according to age group?
 179 ## what for the education group?
 180 ## wants more?
 181 
 182 require(effects)
 183 allEffects(lrfit)
 184 
 185 ## add an interaction effect age:wants More
 186 lrfit2 <- glm( cbind(using, notUsing) ~ education + wantsMore * age,
 187               family = binomial, data = cuse)
 188 
 189 lrfit2 <- glm( cbind(using, notUsing) ~ education + wantsMore + age + wantsMore:age,
 190              family = binomial, data = cuse)
 191 
 192 ## how does the effects change
 193 allEffects(lrfit2)
 194 
 195 ## does the interaction term improves the model significantly?
 196 anova(lrfit,lrfit2, test = "Rao")
 197 anova(lrfit,lrfit2, test = "LRT")
 198 
 199 
 200 require(ggplot2)
 201 require(scales)
 202 
 203 
 204 ggplot(infection,aes(y = as.numeric(infected=="infected"), x = age, color=sex, shape=sex))+
 205   geom_smooth(method = "glm", family = "binomial")+
 206   geom_jitter(position = position_jitter(width = 0, height= 0.08),size=5,face="bold")+
 207   scale_color_manual(values = c("male"="midnightblue", "female"= "deeppink"))+
 208   scale_shape_manual(values = c("male" = "\u2642","female" = "\u2640"))+
 209   labs(title="Infection", x = "age in days", y= "infected in percent") +
 210   scale_y_continuous(breaks=seq(0,1,by=0.25),labels = percent)
 211 
 212 
 213 ggplot(infection,aes(y = as.numeric(infected)-1, x = age, color=sex, shape=sex))+
 214   geom_smooth(method = "glm", family = "binomial")+
 215   geom_jitter(position = position_jitter(width = 0,height = 0.07))
 216 
 217 ggplot(infection,aes(y = as.numeric(infected)-1, x = age, color=sex, shape=sex))+
 218   geom_smooth(method = "glm", family = "binomial")+  
 219   geom_jitter(position = position_jitter(width = 0, height= 0.08),size=5,face="bold")+
 220   scale_color_manual(values = c("male"="midnightblue", "female"= "deeppink"))+
 221   scale_shape_manual(values = c("male" = "\u2642","female" = "\u2640"))+
 222   labs(title="Infection", x = "age in days", y= "infected in percent") +
 223   scale_y_continuous(breaks=seq(0,1,by=0.25),labels = percent)
 224 
 225 
 226 
 227 
 228 ###############################################################
 229 ########################## mixed models #######################
 230 ###############################################################
 231 
 232 
 233 require(lme4)
 234 require(nlmeU)
 235 require(lattice)
 236 
 237 head(armd)
 238 
 239 ## formula (explicitly)
 240 f1 <- formula(visual ~ visual0 + time + treat.f + ## main effects
 241                 treat.f:time + ## interaction 
 242                 (1|subject))  ## random effect
 243 
 244 
 245 m1 <- lmer(f1, data = armd)
 246 
 247 (m1.sum <- summary(m1))
 248 show(m1)
 249 
 250 ## estimation method
 251 isREML(m1)
 252 
 253 ## beta estimates (fixed effects)
 254 fixef(m1)
 255 
 256 ## beta estimates incl. se and t-test
 257 coef(m1.sum)
 258 
 259 ## Var beta
 260 vcov(m1)
 261 
 262 ## b estimates (random effects)
 263 ranef(m1)
 264 
 265 ## beta + coupled b
 266 coef(m1)
 267 
 268 
 269 ## confidence intervals
 270 m1.prall <- profile(m1)
 271 confint(m1.prall)
 272 
 273 confint(logProf(m1.prall))
 274 
 275 
 276 ## model with random intercept
 277 f1 <- formula(visual ~ visual0 + time + treat.f + ## main effects
 278                 treat.f:time + ## interaction 
 279                 (1|subject))  ## random effect
 280 
 281 m1 <- lmer(f1, data = armd)
 282 
 283 
 284 ### fixed effects
 285 fixef(m1)
 286 
 287 coef(summary(m1))
 288 
 289 ### random effects
 290 head(ranef(m1)[[1]])
 291 
 292 require(broom)
 293 glance(m1)
 294 head(tidy(m1))
 295 
 296 ## different kinds of predictions
 297 ## no random effects
 298 armd$pred_overall_gen <- predict(m1, type='response', re.form=NA)
 299 ## all random effects
 300 armd$pred_overall <- predict(m1, type='response')
 301 ## specified random effects
 302 armd$pred_subj <- predict(m1, type='response', re.form=~ (1|subject))
 303 
 304 ## keep only predictions
 305 armd_pred <- merge(armd, as.data.frame(table(subject=armd$subject)))
 306 
 307 
 308 
 309 ggplot(armd[armd$subject %in% 1:16,], aes(colour = treat.f)) +
 310   geom_smooth(aes(x = time, y = pred_overall_gen), method = "lm") +
 311   geom_smooth(aes(x = time, y = pred_overall), linetype = 2, method = "lm", fullrange=T) +
 312   geom_point(aes(x = time, y = visual)) +
 313   facet_wrap(~ subject)
 314 
 315 ggsave("img/armd1.png")
 316 
 317 
 318 ## profile plots
 319 require(lattice)
 320 xyplot(logProf(m1.prall),as.table = T)
 321 xyplot(logProf(m1.prall),absVal = T,as.table = T)
 322 
 323 
 324 ## add random slope
 325 f2 <- formula(visual ~ visual0 + time + treat.f + ## main effects
 326                 treat.f:time + ## interaction 
 327                 (1 + time |subject))  ## random effect
 328 
 329 
 330 m2 <- lmer(f2, data = armd)
 331 
 332 ## different kinds of predictions
 333 ## no random effects
 334 armd$pred_overall_gen2 <- predict(m2, type='response', re.form=NA)
 335 ## all random effects
 336 armd$pred_overall2 <- predict(m2, type='response')
 337 ## specified random effects
 338 armd$pred_subj2 <- predict(m2, type='response', re.form=~ (1|subject))
 339 
 340 
 341 ggplot(armd[armd$subject %in% 1:16,], aes(colour = treat.f)) +
 342   geom_smooth(aes(x = time, y = pred_overall_gen2), method = "lm") +
 343   geom_smooth(aes(x = time, y = pred_overall2), linetype = 2, method = "lm", fullrange=T) +
 344   geom_point(aes(x = time, y = visual)) +
 345   facet_wrap(~ subject)
 346 
 347 ggsave("img/armd2.png")
 348 
 349 
 350 ## residuals
 351 
 352 ## default residual plot
 353 plot(m1)
 354 savePlot("img/defresidplot.png")
 355 
 356 
 357 plot(m1, type = c("p","smooth"))
 358 savePlot("img/defresidplot2.png")
 359 
 360 
 361 ## residuals per timepoints
 362 p1 <- plot(m1, factor(time) ~ resid(., type = "pearson"),
 363            abline=c(v=0), lty=2, xlim = c(-30,30))
 364 p2 <- plot(m2, factor(time) ~ resid(., type = "pearson"),
 365            abline=c(v=0), lty=2, xlim = c(-30,30))
 366 
 367 require(gridExtra)
 368 grid.arrange(p1,p2,nrow = 2)
 369 
 370 savePlot("img/residualsbp.png")
 371 
 372 ## multcomp
 373 require(multcomp)
 374 require(ggplot2)
 375 require(broom)
 376 
 377 ## default behavior
 378 m1.mc <- glht(m1)
 379 m1.mc
 380 
 381 summary(m1.mc)
 382 
 383 m1.mc.tidy <- tidy(summary(m1.mc))
 384 m1.mc.tidy
 385 
 386 ## use resulting object for ggplot
 387 ggplot(m1.mc.tidy, aes(x = lhs, y = estimate)) +
 388   geom_point()
 389 
 390 ggsave("img/multcomp1.png")
 391 
 392 ## extract confidence intervals
 393 m1.mc.ci <- tidy(confint(m1.mc))
 394 
 395 ggplot(m1.mc.ci, aes(x = lhs,
 396                        y = estimate,
 397                        ymin = conf.low,
 398                        ymax = conf.high)) +
 399                          geom_point() +
 400                          geom_pointrange()
 401 
 402 ggsave("img/multcomp2.png")
 403 
 404 
 405 plot(m1.mc)
 406 savePlot("img/multcomp3.png")
 407 
 408 
 409 ## refit the model using time coded as factor
 410 f1 <- formula(visual ~ visual0 + time.f + treat.f + ## main effects
 411                 treat.f:time.f + ## interaction 
 412                 (1|subject))  ## random effect
 413 
 414 m1 <- lmer(f1, data = armd)
 415 
 416 ### tukey
 417 summary(glht(m1,linfct = mcp(time.f = "Tukey")))
 418 
 419 ### dunnett
 420 glht(m1,linfct = mcp(time.f = "Dunnett"))
 421 
 422 ### taking interaction into account
 423 summary(glht(m1,linfct = mcp(time.f = "Tukey", interaction_average = T)))
 424 
 425 
 426 ### extract contrast matrix
 427 glht(m1,linfct = mcp(time.f = "Tukey"))$linfct
 428 
 429 class(armd$time.f)
 430 
 431 glht(m1,linfct = mcp(time.f = "Tukey", interaction_average = T))$linfct
 432 
 433 ### recode main/interaction 
 434 armd$ia <- interaction(armd$treat.f,armd$time.f)
 435 f3 <- formula(visual ~ visual0 + ia +
 436                 (1|subject))  ## random effect
 437 m3 <- lmer(f3, data = armd)
 438 
 439 summary(glht(m3,linfct = mcp(ia = "Tukey")))
 440 
 441 ### set adjustment method
 442 summary(glht(m3,linfct = mcp(ia = "Tukey"),test = "fdr"))
 443 
 444 
 445 ### specify contrasts of interest
 446 mt3.1 <- glht(m3,linfct = mcp(ia = c("Active.52wks - Active.24wks = 0",
 447                          "Placebo.52wks - Active.4wks = 0")))
 448 
 449 mt3.2 <- glht(m3,linfct = mcp(ia = c("Active.52wks - Active.24wks = 1",
 450                             "Placebo.52wks - Active.4wks = 1")))
 451 
 452 
 453 mt3.3 <- glht(m3,linfct = mcp(ia = c("Active.52wks - Active.24wks >= 0",
 454                             "Placebo.52wks - Active.4wks >= 0")))
 455 
 456 mt3.3 <- glht(m3,linfct = mcp(ia = c("Active.52wks - Active.24wks <= 1",
 457                             "Placebo.52wks - Active.4wks <= 1")))
 458 
 459 summary(mt3.1)
 460 summary(mt3.2)
 461 summary(mt3.3)
 462 
 463 
 464 ### simul. testing for both factors
 465 mat1 <- glht(m1,linfct = mcp(time.f = "Tukey", interaction_average = T))$linfct
 466 mat2 <- glht(m1,linfct = mcp(treat.f = "Tukey", interaction_average = T))$linfct
 467 
 468 mat <- rbind(mat1,mat2)
 469 
 470 summary(glht(m1,linfct = mat))
 471 
 472 tidy(summary(glht(m1,linfct = mcp(time.f = "Tukey", interaction_average = T))))
 473 tidy(summary(glht(m1,linfct = mcp(treat.f = "Tukey", interaction_average = T))))
 474 
 475 
 476 ### using model matrix
 477 tmp <- expand.grid(time.f = levels(armd$time.f),
 478                    treat.f = levels(armd$treat.f))
 479 tmp$time.f <- factor(tmp$time.f,ordered = T)
 480 
 481 cm <- model.matrix(~ time.f * treat.f, data = tmp)
 482 cm <- cbind(cm, visual0 = 0)
 483 
 484 glht(m1,linfct = cm)
 485