1 Introduction

The data come from a randomized, double-blind clinical trial comparing active hypnotic drug with a placebo in patients who have insomnia problems.

The response is the patient’s reported time (in minutes) to fall asleep after going to bed. Patients responded before and after a two-week treatment period. The patients receiving the two treatments were independent samples.

Source: Agresti, A. (2002). Categorical Data Analysis. 2nd Ed. p. 462

Column 1: Subject ID Column 2: Treatment, Active = 1, Placebo = 0 Column 3: Time, Initial = 0, Follow-up = 1 Column 4: Time to sleep, < 20 min. = 1, 20-30 = 2, 30-60 = 3, > 60 min. = 4

2 Data management

#input data
dta <- read.table("C:/Users/Ching-Fang Wu/Documents/data/insomnia.txt", h=T)
#show first 6 lines
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
# 重新命名
names(dta) <- c("ID", "Treatment", "Time", "Sleeplatency")
str(dta)
## 'data.frame':    478 obs. of  4 variables:
##  $ ID          : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ Treatment   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Time        : int  0 1 0 1 0 1 0 1 0 1 ...
##  $ Sleeplatency: int  1 1 1 1 1 1 1 1 1 1 ...
#load packages
pacman::p_load(repolr,tidyverse)

3 GEE

summary(dta_repolr <- repolr(Sleeplatency ~ Treatment + Time + Treatment:Time, subjects="ID", data = as.data.frame(dta),
                 times = c(1,2),
                 categories = 4,
                 corr.mod = "independence",
                 po.test=TRUE))
## 
## repolr: 2016-02-26 version 3.4 
## 
## Call:
## repolr(formula = Sleeplatency ~ Treatment + Time + Treatment:Time, 
##     subjects = "ID", data = as.data.frame(dta), times = c(1, 
##         2), categories = 4, corr.mod = "independence", po.test = TRUE)
## 
## 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
## Treatment         0.0336    0.2369     0.1418    0.8872
## Time              1.0381    0.1564     6.6375    0.0000
## Treatment:Time    0.7077    0.2458     2.8792    0.0040
## 
## Correlation Structure:  independence 
## Fixed Correlation:  0 
## PO Score Test: 8.6394 (d.f. = 6 and p.value = 0.1949)
# polr
sjPlot::tab_model(m0 <- MASS::polr(ordered(Sleeplatency) ~ Treatment + Time + Treatment:Time, data=dta, Hess = T),
show.obs=F, show.r2=F, show.se=T, transform=NULL)
  ordered(Sleeplatency)
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
Treatment -0.03 0.24 -0.50 – 0.43 0.888
Time -1.04 0.24 -1.51 – -0.57 <0.001
Treatment * Time -0.71 0.33 -1.36 – -0.05 0.035

Time 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.

summary(dta_repolr$gee)
## Length  Class   Mode 
##      0   NULL   NULL

4 hand-made cumulative prop. against time to sleep

dta <- dta %>% 
  mutate(Treatment = factor(Treatment, 
                            levels=0:1, 
                            labels=c("Placebo","Active")),
         Time = factor(Time, 
                       levels=0:1, 
                       labels=c("Initial","Follow-up")),
         Sleeplatency = factor(Sleeplatency))

dta.tbl <- ftable(dta$Treatment, dta$Time, dta$Sleeplatency)

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

5 Cumulative prop.plot

# prop.table: matrix 1 indicates rows
# t: matrix transpose()
obsp <- as.numeric(t(prop.table(dta.tbl, m=1)))

# Cumulative Sums
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)

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=1, 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=2, col="indianred")

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

grid()

6 The end