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
#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)
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
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
# 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()