1 data input

dta <- read.table("insomnia.txt", h=T)

head(dta)
##   case treat occasion resp
## 1    1     1        0    1
## 2    1     1        1    1
## 3    2     1        0    1
## 4    2     1        1    1
## 5    3     1        0    1
## 6    3     1        1    1

2 repolr function: GEE model fitted in ordinal data

library(repolr)
dta_repolr <- repolr(resp ~ treat + occasion + treat:occasion, 
                     subjects="case", data=dta, times=c(1,2), categories=4, 
                     corr.mod="independence")

# Show repolr modle results (1)
summary.repolr(dta_repolr)
## 
## repolr: 2016-02-26 version 3.4 
## 
## Call:
## repolr(formula = resp ~ treat + occasion + treat:occasion, subjects = "case", 
##     data = dta, times = c(1, 2), categories = 4, corr.mod = "independence")
## 
## Coefficients: 
##                 coeff     se.robust  z.robust  p.value 
## cuts1|2          -2.2671    0.2091   -10.8422    0.0000
## cuts2|3          -0.9515    0.1769    -5.3787    0.0000
## cuts3|4           0.3517    0.1751     2.0086    0.0446
## treat             0.0336    0.2369     0.1418    0.8872
## occasion          1.0381    0.1564     6.6375    0.0000
## treat:occasion    0.7077    0.2458     2.8792    0.0040
## 
## Correlation Structure:  independence 
## Fixed Correlation:  0
# Show repolr modle results (2)
library(MASS)
sjPlot::tab_model(m0 <- MASS::polr(ordered(resp) ~ treat + occasion + treat:occasion, data=dta), show.obs=F, show.r2=F, show.se=T, transform=NULL)
## 
## Re-fitting to get Hessian
  ordered(resp)
Predictors Log-Odds std. Error CI p
1|2 -2.27 0.20 -2.73 – -1.80 <0.001
2|3 -0.95 0.18 -1.43 – -0.48 <0.001
3|4 0.35 0.17 -0.30 – 1.01 0.045
treat -0.03 0.24 -0.50 – 0.43 0.888
occasion -1.04 0.24 -1.51 – -0.57 <0.001
treat * occasion -0.71 0.33 -1.36 – -0.05 0.035

Occasion effect = 1.04 for placebo, 1.04 + 0.71 = 1.75 for treatment.

Odds ratios: exp(1.04) = 2.8, exp(1.75) = 5.7

Treatment effect exp(0.03) = 1.03 initially, exp(0.03+0.71) = 2.1 at follow-up.

3 data management

dta$treat <- factor(ifelse(dta$treat==1, "Active", "Placebo"))

dta$occasion <- factor(ifelse(dta$occasion==1, "Follow-up", "Initial"))

dta.tbl <- ftable(dta$treat, dta$occasion , dta$resp)

dta.tbl
##                     1  2  3  4
##                               
## Active  Follow-up  40 49 19 11
##         Initial    12 20 40 47
## Placebo Follow-up  31 29 35 25
##         Initial    14 20 35 51

4 hand-made cumulative prop. against time to sleep

#
obsp <- as.numeric(t(prop.table(dta.tbl, m=1)))

#
clp <- c(cumsum(obsp[1:4]), cumsum(obsp[5:8]), cumsum(obsp[9:12]), 
         cumsum(obsp[13:16]))
#
yt <-rep(c(20,25,45,60),4)

#
y <- data.frame(time2sleep=yt, clprop=clp)

5 plot

#
plot(y[,1], y[,2], type="n",
     xlab="Time to Sleep (min.)", ylab="Cumulative proportion")

#
lines(y[1:4,  1],  y[1:4,2], lty=1, col="steelblue")
lines(y[5:8,  1],  y[5:8,2], lty=5, col="steelblue")
lines(y[9:12, 1],  y[9:12,2], lty=2, col="indianred")
lines(y[13:16,1],  y[13:16,2], lty=4, col="indianred")

#
text(22, .12, "T,I")
text(22, .47, "T,F")
text(22, .20, "P,I")
text(22, .33, "P,F")
#
grid()

累積折線圖顯示，採用催眠藥物組和安慰劑組間起初皆以相似的時間入睡。但在兩週後，接受催眠藥物治療的人比起安慰劑組，能更迅速入眠。(TF>PF)