GEE Cacah/Poisson

library(gee)
## Warning: package 'gee' was built under R version 4.3.3

Input Data

library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.3.3
data <- read.xlsx("C:\\Users\\user\\Downloads\\Data GEE.xlsx")
head(data)
##   id trt base y1 y2 y3 y4
## 1  1   0    5  4  4  3  3
## 2  2   0    6  5  5  4  4
## 3  3   0    4  4  3  3  3
## 4  4   0    7  6  6  5  5
## 5  5   0    5  5  4  4  4
## 6  6   1    6  4  3  2  2

#Goodness of Fit Model #QIC

QIC.gee = function(model.R)
{
  library(MASS)
  model.indep = update(model.R, corstr = "independence")
  
  mu.R=model.R$fitted.values
  y = model.R$y
  type = family(model.R)$family
  
  quasi.R = switch(type,
                   poisson = sum((y*log(mu.R)) - mu.R),
                   gaussian = sum(((y - mu.R)^2)/-2),
                   binomial = sum(y*log(mu.R/(1 - mu.R)) + log(1 - mu.R)),
                   Gamma = sum(-y/(mu.R - log(mu.R))),
                   stop("Error"))
  
  omegaI = ginv(model.indep$naive.variance) 
  Vr = model.R$robust.variance
  trace.R = sum(diag(omegaI %*% Vr))
  px = length(mu.R)
  
  QIC = 2*trace.R -2*quasi.R
  output = c(QIC, quasi.R, trace.R, px)
  names(output) = c('QIC', 'Quasi Lik', 'Trace', 'px')
  return(output)
}

#====== Ubah Format Data dari Wide ke Long =========

data_long = reshape(data,
                    varying=list(c("base","y1", "y2", "y3", "y4")),
                    v.names="count", times=0:4, direction="long")

data_long = data_long[order(data_long$id, data_long$time),]
data_long$t = ifelse(data_long$time == 0, 8, 2)
data_long$visit = ifelse(data_long$time == 0, 0, 1)

Model

m11 = gee(count ~ offset(log(t)) + visit*trt,
          id = id, data = data_long,
          family = poisson, corstr = "exchangeable")
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)       visit         trt   visit:trt 
##  -0.3930426   1.1349799   0.1698990  -0.3960232
summary(m11)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logarithm 
##  Variance to Mean Relation: Poisson 
##  Correlation Structure:     Exchangeable 
## 
## Call:
## gee(formula = count ~ offset(log(t)) + visit * trt, id = id, 
##     data = data_long, family = poisson, corstr = "exchangeable")
## 
## Summary of Residuals:
##    Min     1Q Median     3Q    Max 
## 0.3250 1.3250 2.1125 3.3250 7.2000 
## 
## 
## Coefficients:
##               Estimate Naive S.E.   Naive z Robust S.E.  Robust z
## (Intercept) -0.3930426 0.10308861 -3.812667  0.08445744 -4.653736
## visit        1.1349799 0.08059462 14.082577  0.02742916 41.378596
## trt          0.1698990 0.13997862  1.213750  0.11050423  1.537489
## visit:trt   -0.3960232 0.11506311 -3.441791  0.05069255 -7.812256
## 
## Estimated Scale Parameter:  0.2869361
## Number of Iterations:  1
## 
## Working Correlation
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 1.0000000 0.5448543 0.5448543 0.5448543 0.5448543
## [2,] 0.5448543 1.0000000 0.5448543 0.5448543 0.5448543
## [3,] 0.5448543 0.5448543 1.0000000 0.5448543 0.5448543
## [4,] 0.5448543 0.5448543 0.5448543 1.0000000 0.5448543
## [5,] 0.5448543 0.5448543 0.5448543 0.5448543 1.0000000
m21 = gee(count ~ offset(log(t)) + trt*visit,
          id = id, data = data_long,
          family = poisson, corstr = "AR-M", Mv = 1)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)         trt       visit   trt:visit 
##  -0.3930426   0.1698990   1.1349799  -0.3960232
summary(m21)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logarithm 
##  Variance to Mean Relation: Poisson 
##  Correlation Structure:     AR-M , M = 1 
## 
## Call:
## gee(formula = count ~ offset(log(t)) + trt * visit, id = id, 
##     data = data_long, family = poisson, corstr = "AR-M", Mv = 1)
## 
## Summary of Residuals:
##       Min        1Q    Median        3Q       Max 
## 0.2740834 1.2740834 2.0700695 3.3556630 7.3203028 
## 
## 
## Coefficients:
##                Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept) -0.48265322 0.10300992 -4.6855024  0.09463243 -5.1002944
## trt          0.09654528 0.14228682  0.6785258  0.11843759  0.8151575
## visit        1.24062532 0.07097824 17.4789518  0.04111950 30.1712184
## trt:visit   -0.30875911 0.10355607 -2.9815647  0.05823357 -5.3020814
## 
## Estimated Scale Parameter:  0.3109166
## Number of Iterations:  5
## 
## Working Correlation
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 1.0000000 0.8087352 0.6540526 0.5289554 0.4277849
## [2,] 0.8087352 1.0000000 0.8087352 0.6540526 0.5289554
## [3,] 0.6540526 0.8087352 1.0000000 0.8087352 0.6540526
## [4,] 0.5289554 0.6540526 0.8087352 1.0000000 0.8087352
## [5,] 0.4277849 0.5289554 0.6540526 0.8087352 1.0000000

uji parsial

2*(pnorm(abs(summary(m11)$coefficients[,5]), lower.tail = FALSE))
##  (Intercept)        visit          trt    visit:trt 
## 3.259748e-06 0.000000e+00 1.241736e-01 5.617328e-15
2*(pnorm(abs(summary(m21)$coefficients[,5]), lower.tail = FALSE))
##   (Intercept)           trt         visit     trt:visit 
##  3.391257e-07  4.149822e-01 5.652140e-200  1.144898e-07
QIC.gee(m11)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)       visit         trt   visit:trt 
##  -0.3930426   1.1349799   0.1698990  -0.3960232
##       QIC Quasi Lik     Trace        px 
## 19.715963 -3.745524  6.112457 50.000000
QIC.gee(m21)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
## (Intercept)         trt       visit   trt:visit 
##  -0.3930426   0.1698990   1.1349799  -0.3960232
##       QIC Quasi Lik     Trace        px 
## 30.122300 -8.833056  6.228094 50.000000