R Stan examples: GLMM: Treatment for Migraines

York University

if(.Platform$OS.type == 'windows') windowsFonts(Arial=windowsFont("TT Arial")) 
library(spida2)
library(car)
library(lattice)
library(latticeExtra)

## Loading required package: RColorBrewer

library(effects)

## 
## Attaching package: 'effects'

## The following object is masked from 'package:car':
## 
##     Prestige

library(MASS)      # for glmmPQL

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:spida2':
## 
##     Null

library(rstan)

## Loading required package: ggplot2

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:latticeExtra':
## 
##     layer

## The following object is masked from 'package:spida2':
## 
##     labs

## Loading required package: StanHeaders

## rstan (Version 2.15.1, packaged: 2017-04-19 05:03:57 UTC, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

data(migraines)
dd <- migraines
?migraines

## starting httpd help server ...

##  done

Quick look

head(dd)

##   id ha time dos hatype female age airq medication
## 1  1  1  -11 753   Aura      1  30    9 Continuing
## 2  1  1  -10 754   Aura      1  30    7 Continuing
## 3  1  1   -9 755   Aura      1  30   10 Continuing
## 4  1  1   -8 756   Aura      1  30   13 Continuing
## 5  1  1   -7 757   Aura      1  30   18 Continuing
## 6  1  1   -6 758   Aura      1  30   19 Continuing

xqplot(dd)

xqplot(up(dd, ~id))

gicc(dd, ~ id)

##         id         ha       time        dos     hatype     female 
##  1.0000000  0.2998662  0.3155582  0.9958465  1.0000000  1.0000000 
##        age       airq medication 
##  1.0000000  0.3705792  1.0000000

A few summary stats:

dd$length <- capply(dd$id, dd$id, length)
tab(dd, ~length) # not quite what we want

## length
##     7    11    12    13    14    16    17    18    19    20    21    22 
##     7    11    48    13    28    64    68   180   171   140   168    66 
##    23    24    25    26    27    28    29    30    31    32    33    34 
##    23   120    50   130    81   112   116    90    31    96   198   170 
##    35    38    39    40    41    42    43    45    46    47    48    51 
##   105    76    39    80    41    84    43    90    92   141    48   102 
##    54    55    56    59    60    62    63    70    75    80    91   121 
##   108   110    56    59    60    62    63    70   150    80    91   121 
## Total 
##  4152

dd %>% up(~id) %>% tab(~length) # better: 1 sub with 7 records to 1 with 121

## length
##     7    11    12    13    14    16    17    18    19    20    21    22 
##     1     1     4     1     2     4     4    10     9     7     8     3 
##    23    24    25    26    27    28    29    30    31    32    33    34 
##     1     5     2     5     3     4     4     3     1     3     6     5 
##    35    38    39    40    41    42    43    45    46    47    48    51 
##     3     2     1     2     1     2     1     2     2     3     1     2 
##    54    55    56    59    60    62    63    70    75    80    91   121 
##     2     2     1     1     1     1     1     1     2     1     1     1 
## Total 
##   133

dd$ha_mean <- capply(dd$ha, dd$id, mean)
dd %>% up(~id) %>% xqplot 

dd %>% up(~id) %>% scatterplotMatrix

• time = 0 indicates the onset of treatment
• We need to create a Treatment variable by using time

This is slow but worth it if you are seriously studying this data set

dd$id.o <- dd$id %>% factor %>% reorder(dd$length)
xyplot( jitter(ha) ~ time | id.o, dd, strip = FALSE,
        layout = c( 8, 9 ),
        panel = function ( x , y , ...){
          panel.xyplot( x, y , ...)
          panel.loess (x , y , ..., family = 'gaussian', lwd = 2 , col = 'red')
          panel.abline(v=0)
        })

knitr::knit_exit()