library(readr)
## Warning: package 'readr' was built under R version 3.4.1
resp <- read_csv("~/BiostatCourses/PublicHealthComputing/Lectures/LongitudinalData/resp.csv")
## Parsed with column specification:
## cols(
##   id = col_integer(),
##   center = col_integer(),
##   treat = col_integer(),
##   sex = col_integer(),
##   age = col_integer(),
##   bl = col_integer(),
##   time = col_integer(),
##   status = col_integer()
## )

To fit a model with GEE, we can use the gee function in the R package gee.

library(gee)
## Warning: package 'gee' was built under R version 3.4.2
resp.mod1 <- gee(status ~ center + treat + sex + age + time + bl, data= resp,
                 family=binomial(link="logit"),id = id, corstr="independence")
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)      center       treat         sex         age        time 
## -2.83272325  0.67230658  1.30055737  0.11936743 -0.01818538 -0.06425079 
##          bl 
##  1.88406115
summary(resp.mod1)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Independent 
## 
## Call:
## gee(formula = status ~ center + treat + sex + age + time + bl, 
##     id = id, data = resp, family = binomial(link = "logit"), 
##     corstr = "independence")
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.93354161 -0.29377341  0.08661524  0.32842660  0.85563094 
## 
## 
## Coefficients:
##                Estimate  Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept) -2.83272325 0.683780391 -4.1427383  0.89810798 -3.1541010
## center       0.67230658 0.241080689  2.7887202  0.35722785  1.8820106
## treat        1.30055737 0.238356076  5.4563634  0.35097961  3.7055069
## sex          0.11936743 0.296519345  0.4025620  0.44369076  0.2690329
## age         -0.01818538 0.008920227 -2.0386674  0.01302443 -1.3962509
## time        -0.06425079 0.100074417 -0.6420302  0.08159845 -0.7874021
## bl           1.88406115 0.242876076  7.7572941  0.35018472  5.3801924
## 
## Estimated Scale Parameter:  1.011563
## Number of Iterations:  1
## 
## Working Correlation
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1
exp(coef(resp.mod1)[3])
##    treat 
## 3.671342
resp.mod2 <- gee(status ~ center + treat + sex + age + time + bl, data= resp,
                 family=binomial(link="logit"),id = id, corstr="unstructured")
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)      center       treat         sex         age        time 
## -2.83272325  0.67230658  1.30055737  0.11936743 -0.01818538 -0.06425079 
##          bl 
##  1.88406115
summary(resp.mod2)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Unstructured 
## 
## Call:
## gee(formula = status ~ center + treat + sex + age + time + bl, 
##     id = id, data = resp, family = binomial(link = "logit"), 
##     corstr = "unstructured")
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.93491199 -0.29761236  0.08394397  0.32853864  0.85504000 
## 
## 
## Coefficients:
##                Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept) -2.81444328 0.93082833 -3.0235901  0.90495875 -3.1100238
## center       0.68139125 0.34295483  1.9868251  0.35468060  1.9211405
## treat        1.27368865 0.33917818  3.7552199  0.34897488  3.6498004
## sex          0.10385184 0.42220560  0.2459746  0.44432393  0.2337300
## age         -0.01708821 0.01270377 -1.3451293  0.01291311 -1.3233224
## time        -0.07109647 0.08588739 -0.8277870  0.08137098 -0.8737325
## bl           1.93201249 0.34665147  5.5733573  0.34801055  5.5515917
## 
## Estimated Scale Parameter:  1.020586
## Number of Iterations:  3
## 
## Working Correlation
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.3272443 0.2036405 0.3010760
## [2,] 0.3272443 1.0000000 0.4339166 0.3559958
## [3,] 0.2036405 0.4339166 1.0000000 0.4014093
## [4,] 0.3010760 0.3559958 0.4014093 1.0000000
exp(coef(resp.mod2)[3])
##    treat 
## 3.574012
resp.mod3 <- gee(status ~ center + treat + sex + age + time + bl, data= resp,
                 family=binomial(link="logit"),id = id, corstr="exchangeable")
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)      center       treat         sex         age        time 
## -2.83272325  0.67230658  1.30055737  0.11936743 -0.01818538 -0.06425079 
##          bl 
##  1.88406115
summary(resp.mod3)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Exchangeable 
## 
## Call:
## gee(formula = status ~ center + treat + sex + age + time + bl, 
##     id = id, data = resp, family = binomial(link = "logit"), 
##     corstr = "exchangeable")
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.93336064 -0.29731249  0.08709546  0.32904495  0.85591042 
## 
## 
## Coefficients:
##                Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept) -2.83714023 0.92581792 -3.0644689  0.89532598 -3.1688349
## center       0.68091583 0.34137967  1.9945998  0.35676012  1.9086097
## treat        1.29222142 0.33738314  3.8301303  0.35053559  3.6864200
## sex          0.13076523 0.42009688  0.3112740  0.44409802  0.2944513
## age         -0.01841147 0.01263567 -1.4571028  0.01301183 -1.4149795
## time        -0.06419992 0.08142748 -0.7884306  0.08150945 -0.7876377
## bl           1.87781449 0.34368304  5.4637973  0.35009816  5.3636799
## 
## Estimated Scale Parameter:  1.010477
## Number of Iterations:  2
## 
## Working Correlation
##          [,1]     [,2]     [,3]     [,4]
## [1,] 1.000000 0.336918 0.336918 0.336918
## [2,] 0.336918 1.000000 0.336918 0.336918
## [3,] 0.336918 0.336918 1.000000 0.336918
## [4,] 0.336918 0.336918 0.336918 1.000000
exp(coef(resp.mod3)[3])
##    treat 
## 3.640865
resp.mod4 <- gee(status ~ center + treat + sex + age + time + bl, data= resp,
                 family=binomial(link="logit"),id = id, corstr="AR-M",Mv=1)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)      center       treat         sex         age        time 
## -2.83272325  0.67230658  1.30055737  0.11936743 -0.01818538 -0.06425079 
##          bl 
##  1.88406115
summary(resp.mod4)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     AR-M , M = 1 
## 
## Call:
## gee(formula = status ~ center + treat + sex + age + time + bl, 
##     id = id, data = resp, family = binomial(link = "logit"), 
##     corstr = "AR-M", Mv = 1)
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.93470602 -0.31201913  0.08454742  0.32964616  0.86032474 
## 
## 
## Coefficients:
##                Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept) -2.86519283 0.88267297 -3.2460412  0.91516891 -3.1307804
## center       0.74721383 0.31853870  2.3457553  0.35689307  2.0936630
## treat        1.24611300 0.31509343  3.9547413  0.35214254  3.5386608
## sex          0.11596705 0.39319669  0.2949339  0.45030751  0.2575286
## age         -0.01697868 0.01180962 -1.4376992  0.01295510 -1.3105784
## time        -0.08300192 0.10255013 -0.8093790  0.08208583 -1.0111601
## bl           1.91250297 0.32077029  5.9622198  0.35040176  5.4580290
## 
## Estimated Scale Parameter:  1.019906
## Number of Iterations:  3
## 
## Working Correlation
##            [,1]      [,2]      [,3]       [,4]
## [1,] 1.00000000 0.3901160 0.1521905 0.05937193
## [2,] 0.39011597 1.0000000 0.3901160 0.15219047
## [3,] 0.15219047 0.3901160 1.0000000 0.39011597
## [4,] 0.05937193 0.1521905 0.3901160 1.00000000
exp(coef(resp.mod4)[3])
##    treat 
## 3.476802

GLMM models can be fit with the glmer in the R package lme4.

library(lme4)
## Loading required package: Matrix
resp.mod5 <- glmer(status ~ center + treat + sex + age + time + bl + (1|id), data= resp,
                 family=binomial(link="logit"))
summary(resp.mod5)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: status ~ center + treat + sex + age + time + bl + (1 | id)
##    Data: resp
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.3    477.1   -214.2    428.3      436 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8574 -0.3590  0.1427  0.3693  2.2393 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  id     (Intercept) 3.89     1.972   
## Number of obs: 444, groups:  id, 111
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.80831    1.51838  -3.167 0.001542 ** 
## center       1.04669    0.54784   1.911 0.056060 .  
## treat        2.16328    0.55644   3.888 0.000101 ***
## sex          0.20247    0.67270   0.301 0.763425    
## age         -0.02546    0.02014  -1.264 0.206147    
## time        -0.10133    0.12518  -0.810 0.418221    
## bl           3.07840    0.60272   5.108 3.26e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) center treat  sex    age    time  
## center -0.360                                   
## treat  -0.703  0.058                            
## sex    -0.475 -0.147  0.219                     
## age    -0.155 -0.223 -0.050 -0.263              
## time   -0.178 -0.015 -0.031 -0.003  0.009       
## bl     -0.294 -0.150  0.301  0.102  0.015 -0.041
exp(coef(resp.mod5)$id[1,3])
## [1] 8.699627