1 load("birthsweights.rdata")
   2 
   3 ## Exercises
   4 ## Remember: We used hypertension of the mom to explain variation
   5 ## in the birth weight of the kid. Without looking in the material
   6 ## of the last session, try to redo the model.
   7 
   8 
   9 
  10 
  11 
  12 require(effects)
  13 
  14 
  15 
  16 
  17 res <- allEffects(m.hyp, se = T)
  18 summary(res)
  19 
  20 invlogit <- function(x){
  21     exp(x)/(1+exp(x))
  22 }
  23 
  24 
  25 ## What is the relationship between the coefficients of the model
  26 ## (from the model summary) and the effects?
  27 
  28 
  29 
  30 
  31 ## binomial/logistic regression
  32 m.wks <- glm(lowbw ~ gestwks, family=binomial, data=births)
  33 summary(m.wks)
  34 
  35 
  36 ## interpreting the coefficients
  37 invlogit(coef(m.wks)[1])
  38 
  39 exp(coef(m.wks)[2])
  40 
  41 
  42 ## effects
  43 Effect("gestwks",m.wks)
  44 
  45 
  46 ## Exercise
  47 Effect("gestwks",m.wks,xlevels = list(gestwks = c(20,30,40)))
  48 
  49 
  50 ## visualiziation
  51 require(ggplot2)
  52 ggplot(births,aes(x = gestwks, y = as.numeric(lowbw)-1)) +
  53     geom_smooth(method = "glm", family = "binomial",se = F,size = 2) +
  54     geom_point(shape="|")
  55 
  56 
  57 ## Exercises
  58 
  59 ## try to change the axis titles (xlab() and ylab())
  60 ## add a title (ggtitle())
  61 
  62 
  63 
  64 
  65 
  66 
  67 ## change the colour of the function to black
  68 ## change the colour of the points to red for the low birth weight
  69 ## and to green for the one with normal birth weight
  70 
  71 
  72 
  73 
  74 
  75 
  76 
  77 ## change the position of the legend;
  78 ## place it somewhere near the upper right corner
  79 ## inside the plotting area
  80 
  81 
  82 
  83 
  84 
  85 
  86 
  87 ## the famous Challenger DAta
  88 
  89 require(faraway)
  90 data(orings)
  91 
  92 ## model
  93 m.oring <- glm(cbind(damage,6-damage) ~ temp,family=binomial, orings)
  94 summary(m.oring)
  95 
  96 ## coefficients
  97 invlogit(coef(m.oring)[1])
  98 invlogit(coef(m.oring)[2])
  99 
 100 exp(coef(m.oring)[2])
 101 
 102 Effect("temp",m.oring,
 103        xlevels = list(temp = c(50:60)),
 104        transformation = list(link = logit,inverse = exp))
 105 
 106 
 107 
 108 
 109 require(scales)
 110 demo(plotmath)
 111 
 112 orings$trials <- 6
 113 ggplot(orings,aes(x=temp,y=damage/trials)) +
 114     geom_point() +
 115     geom_smooth(method = "glm", family = "binomial", aes(weight = trials)) +
 116     xlab(expression(paste("temperature (",degree,"F)"))) +
 117     ylab("probability of damage") +
 118     scale_y_continuous(labels = percent)
 119 
 120 ggsave("img/challenger.png")
 121 
 122 ggplot(orings,aes(x=temp,y=damage/trials)) +
 123     geom_point() +
 124     xlim(20,90) +
 125     geom_smooth(method = "glm", family = "binomial",
 126                 aes(weight = trials), fullrange =T) +
 127     xlab(expression(paste("temperature (",degree,"F)"))) +
 128     ylab("probability of damage") +
 129     scale_y_continuous(labels = percent)
 130 
 131 
 132 ggsave("img/challenger2.png")
 133 
 134 
 135 ## Infection Example
 136 
 137 load("infection.rdata")
 138 summary(infection)
 139 
 140 m.inf <- glm(infected~age*sex,family=binomial,
 141                                data=infection)
 142 
 143 summary(m.inf)
 144 
 145 ## coefficients
 146 ### intercepts
 147 invlogit(coef(m.inf)[1])
 148 invlogit(coef(m.inf)[1] + coef(m.inf)[3])
 149 
 150 ### slopes
 151 exp(coef(m.inf)[2])
 152 exp(coef(m.inf)[2] + coef(m.inf)[4])
 153 
 154 exp(30 * coef(m.inf)[2])
 155 exp(30 * (coef(m.inf)[2] + coef(m.inf)[4]))
 156 
 157 solve(0.015657,3.000513)
 158 solve(0.02670685,2.883849)
 159 
 160 ## all effects
 161 
 162 allEffects(m.inf)
 163 
 164 allEffects(m.inf,
 165            xlevels = list(age = seq(0,200,by = 50)))
 166 
 167 
 168 
 169 
 170 
 171 
 172 
 173 
 174 
 175 
 176 
 177 
 178 
 179 
 180 
 181 
 182 
 183 ## comparing different link functions
 184 ggplot(births,aes(x = gestwks, y = as.numeric(lowbw)-1)) +
 185     geom_smooth(method = "glm", family = "binomial",se = T,size = 2) +
 186     geom_smooth(method = "glm", family = binomial(link = "probit"),se = T,size = 2, colour = "red") +
 187     geom_smooth(method = "glm", family = binomial(link = "cloglog"),se = T,size = 2, colour = "green") +
 188     geom_point(shape="|")
 189