1 ## Exercise
   2 ## The excel file data.xlsx contains the worksheets mother and child
   3 ## containing respective parts of the births data set. Use the
   4 ## read_excel() command to read both data sets and
   5 
   6 
   7 
   8 
   9 
  10 
  11 ## use merge() to join them
  12 names(mdata)
  13 names(cdata)
  14 
  15 intersect(names(mdata),names(cdata))
  16 
  17 
  18 
  19 str(births)
  20 
  21 summary(births$hyp)
  22 births$hyp <- factor(births$hyp,levels = c("normal","hyper"))
  23 summary(births$hyp)
  24 
  25 summary(births$lowbw)
  26 births$lowbw <- factor(births$lowbw,levels = c("normal","low"))
  27 summary(births$lowbw)
  28 
  29 summary(births$sex)
  30 births$sex <- factor(births$sex,levels = c("M","F"))
  31 summary(births$sex)
  32 
  33 ## from lm() to glm()
  34 m1 <- lm(bweight ~ hyp, data=births)
  35 m2 <- glm(bweight ~ hyp, family=gaussian, data=births)
  36 
  37 ## binary response
  38 y <- sample(0:1,100, replace = T)
  39 
  40 plot(y)
  41 abline(h=mean(y))
  42 
  43 p <- seq(0,1,by=0.05)
  44 
  45 log(p/(1-p))
  46 
  47 ## inverse logit function
  48 
  49 invlogit <- function(x){
  50     exp(x)/(1+exp(x))
  51 }
  52 
  53 ## install and load the arm package, there is also a invlogit function
  54 ## low birth weight model
  55 m <- glm(lowbw ~ hyp, family=binomial, data=births)
  56 
  57 invlogit(coef(m)[1])
  58 
  59 table(births$lowbw,births$hyp)
  60 40/(388+40)
  61 
  62 ## proportion test
  63 prop.test(c(20,40),c(72,428))
  64 
  65 chisq.test(table(births$lowbw,births$hyp))
  66 
  67 ## effects plot
  68 require(effects)
  69 plot(Effect("hyp",m))
  70 
  71 ## effects
  72 Effect("hyp",m)
  73 
  74 
  75 ## controlling for sex
  76 m2 <- glm(lowbw ~ hyp+sex, family=binomial, data=births)
  77 summary(m2)
  78 
  79 
  80 m3 <- glm(lowbw ~ hyp*sex, family=binomial, data=births)
  81 summary(m3)
  82 
  83 
  84 
  85 Effect("hyp",m2)
  86 Effect("sex",m2)
  87 
  88 plot(Effect("hyp",m2))
  89 plot(Effect("sex",m2))
  90 
  91 Effect(c("hyp","sex"),m2)
  92 Effect(c("sex","hyp"),m2)
  93 
  94 plot(Effect(c("hyp","sex"),m2))
  95 plot(Effect(c("sex","hyp"),m2))
  96 
  97 
  98 
  99 ## Exercise
 100 ## the estimated probability for moms with hypertension to get
 101 ## a baby with low birth weight for all three models
 102 
 103 
 104 
 105 
 106 
 107 ## is their a difference in effects between boys and girls?
 108 ## Which model can answer this question?
 109 
 110 
 111 
 112 
 113 ## stratified model
 114 
 115 m4 <- glm(lowbw ~ sex + sex:hyp, family=binomial, data=births)
 116 summary(m4)
 117 
 118 m4 <- glm(lowbw ~ sex/hyp, family=binomial, data=births)
 119 m4
 120 
 121 ## Exercise
 122 
 123 
 124 
 125 ## understanding the coefficients
 126 ftable(births$hyp,
 127        births$sex,
 128        births$lowbw)
 129 
 130 ## male/normal bp
 131 15/(206+15)
 132 ## female/normal bp
 133 25/(25+182)
 134 ## male/high bp
 135 12/(12+31)
 136 ## female/high bp
 137 8/(8+21)
 138 
 139 prop.table(table(births$lowbw,
 140                  births$hyp,
 141                  births$sex),c(2,3))
 142 
 143 
 144 ## logistic regression
 145 m5 <- glm(lowbw ~ gestwks, family = binomial, data= births)
 146 summary(m5)
 147 
 148 ## understanding the effects
 149 ### intercept
 150 invlogit(coef(m5)[1])
 151 
 152 ### slope
 153 exp(coef(m5)[2])
 154 
 155 ## exercise
 156 
 157 Effect("gestwks",m5)
 158 Effect("gestwks",m5,xlevels = list(gestwks = c(20,30,40)))
 159 
 160 ## use the command to gain the estimated probability of
 161 ## low birth weight for a gestational age of 27 and 36 weeks
 162 
 163 
 164 
 165 
 166 ## example for calculating the odds
 167 Effect("gestwks",m5,xlevels = list(gestwks = c(27,28)))
 168 Effect("gestwks",m5,xlevels = list(gestwks = c(39,40)))
 169 
 170 p.from.odds <- function(odds) odds/(1+odds)
 171 odds.from.p <- function(p) p/(1-p)
 172 
 173 odds.from.p(0.01779725)/odds.from.p(0.04252149)
 174 
 175 
 176 require(ggplot2)
 177 ggplot(births,aes(x = gestwks, y = as.numeric(lowbw)-1)) +
 178     geom_smooth(method = "glm", family = "binomial",se = T,size = 2) +
 179     geom_point(shape="|")
 180 
 181 ggsave("img/glmggplot.png")
 182 
 183 
 184 ## try to change the axis titles (xlab() and ylab())
 185 ## add a title (ggtitle())
 186 
 187 ggplot(births,aes(x = gestwks, y = as.numeric(lowbw)-1)) +
 188     geom_smooth(method = "glm", family = "binomial",se = T,size = 2) +
 189     geom_point(shape="|")
 190 
 191 
 192 
 193 
 194 
 195 ## change the colour of the function to black
 196 ## change the colour of the points to red for the low birth weight
 197 ## and to green for the one with normal birth weight
 198 
 199 ggplot(births,aes(x = gestwks, y = as.numeric(lowbw)-1)) +
 200     geom_smooth(method = "glm", family = "binomial",se = T,size = 2) +
 201     geom_point(shape="|")
 202 
 203 
 204 
 205 
 206 
 207 ## change the position of the legend;
 208 ## place it somewhere near the upper right corner
 209 ## inside the plotting area
 210 
 211 
 212 ggplot(births,aes(x = gestwks, y = as.numeric(lowbw)-1)) +
 213     geom_smooth(method = "glm", family = "binomial",se = T,size = 2) +
 214     geom_point(shape="|")
 215