Data Management
- 載入資料
dta2 <- read.table("C:/Users/ASUS/Desktop/data/insomnia.txt", h=T)
str(dta2)
## 'data.frame': 478 obs. of 4 variables:
## $ case : int 1 1 2 2 3 3 4 4 5 5 ...
## $ treat : int 1 1 1 1 1 1 1 1 1 1 ...
## $ occasion: int 0 1 0 1 0 1 0 1 0 1 ...
## $ resp : int 1 1 1 1 1 1 1 1 1 1 ...
head(dta2)
## 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
- 重新命名
colnames(dta2) <- c("ID", "Trt", "Time", "Sleeptime")
str(dta2)
## 'data.frame': 478 obs. of 4 variables:
## $ ID : int 1 1 2 2 3 3 4 4 5 5 ...
## $ Trt : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Time : int 0 1 0 1 0 1 0 1 0 1 ...
## $ Sleeptime: int 1 1 1 1 1 1 1 1 1 1 ...
Data Analysis
library(repolr)
## Warning: package 'repolr' was built under R version 4.0.3
dta2_gee<-repolr(Sleeptime ~ Trt + Time + Trt:Time,
subjects="ID", data=as.data.frame(dta2), times=c(1,2), categories=4,
corr.mod ="independence", po.test=TRUE)
summary(dta2_gee)
##
## repolr: 2016-02-26 version 3.4
##
## Call:
## repolr(formula = Sleeptime ~ Trt + Time + Trt:Time, subjects = "ID",
## data = as.data.frame(dta2), 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
## Trt 0.0336 0.2369 0.1418 0.8872
## Time 1.0381 0.1564 6.6375 0.0000
## Trt: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)
Graph
# hand-made cumulative prop. against time to sleep
#
dta2$Trt <- factor(ifelse(dta2$Trt==1, "Active", "Placebo"))
#
dta2$Time<- factor(ifelse(dta2$Time==1, "Follow-up", "Initial"))
#
dta.tbl <- ftable(dta2$Trt, dta2$Time , dta2$Sleeptime)
#
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
# t: matrix transpose() 矩陣轉置?
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)
#
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()
# The end