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