1 input data

data(knee, package="catdata")
dta3<-knee
head(dta3)
##   N Th Age Sex R1 R2 R3 R4
## 1 1  1  28   1  4  4  4  4
## 2 2  1  32   1  4  4  4  4
## 3 3  1  41   1  3  3  3  3
## 4 4  2  21   1  4  3  3  2
## 5 5  2  34   1  4  3  3  2
## 6 6  1  24   1  3  3  3  2
str(dta3)
## 'data.frame':    127 obs. of  8 variables:
##  $ N  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Th : int  1 1 1 2 2 1 2 2 2 1 ...
##  $ Age: int  28 32 41 21 34 24 28 40 24 39 ...
##  $ Sex: num  1 1 1 1 1 1 1 1 0 0 ...
##  $ R1 : int  4 4 3 4 4 3 4 3 4 4 ...
##  $ R2 : int  4 4 3 3 3 3 3 2 4 4 ...
##  $ R3 : int  4 4 3 3 3 3 3 2 4 4 ...
##  $ R4 : int  4 4 3 2 2 2 2 2 3 3 ...

2 data management

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dta_f <- dta3 %>% 
  mutate(ID = factor(N),
         therapies =factor(Th, 1:2, labels = c('A', 'B')),
         Sex =factor(Sex, 0:1, labels = c("M", "F")))
str(dta_f)
## 'data.frame':    127 obs. of  10 variables:
##  $ N        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Th       : int  1 1 1 2 2 1 2 2 2 1 ...
##  $ Age      : int  28 32 41 21 34 24 28 40 24 39 ...
##  $ Sex      : Factor w/ 2 levels "M","F": 2 2 2 2 2 2 2 2 1 1 ...
##  $ R1       : int  4 4 3 4 4 3 4 3 4 4 ...
##  $ R2       : int  4 4 3 3 3 3 3 2 4 4 ...
##  $ R3       : int  4 4 3 3 3 3 3 2 4 4 ...
##  $ R4       : int  4 4 3 2 2 2 2 2 3 3 ...
##  $ ID       : Factor w/ 127 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ therapies: Factor w/ 2 levels "A","B": 1 1 1 2 2 1 2 2 2 1 ...

3 Long form

library(tidyr)
dtaL <- gather(data=dta_f, key=time, value=Painlevel, R1:R4, factor_key=TRUE)  
dtaL2<- dtaL[,c(5,6,4,3,7,8)]
head(dtaL2)
##   ID therapies Sex Age time Painlevel
## 1  1         A   F  28   R1         4
## 2  2         A   F  32   R1         4
## 3  3         A   F  41   R1         3
## 4  4         B   F  21   R1         4
## 5  5         B   F  34   R1         4
## 6  6         A   F  24   R1         3
str(dtaL2)
## 'data.frame':    508 obs. of  6 variables:
##  $ ID       : Factor w/ 127 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ therapies: Factor w/ 2 levels "A","B": 1 1 1 2 2 1 2 2 2 1 ...
##  $ Sex      : Factor w/ 2 levels "M","F": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Age      : int  28 32 41 21 34 24 28 40 24 39 ...
##  $ time     : Factor w/ 4 levels "R1","R2","R3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Painlevel: int  4 4 3 4 4 3 4 3 4 4 ...
dtaL3 <- dtaL2 %>%
  mutate(
    Painlevel=dplyr::recode(Painlevel, 'Nopain', 'Very mild', 'Mild', 'Moderate','Severe'),
    Painlevel=factor(Painlevel, levels=c('Nopain', 'Very mild', 'Mild', 'Moderate','Severe')))

head(dtaL3)
##   ID therapies Sex Age time Painlevel
## 1  1         A   F  28   R1  Moderate
## 2  2         A   F  32   R1  Moderate
## 3  3         A   F  41   R1      Mild
## 4  4         B   F  21   R1  Moderate
## 5  5         B   F  34   R1  Moderate
## 6  6         A   F  24   R1      Mild
str(dtaL3)
## 'data.frame':    508 obs. of  6 variables:
##  $ ID       : Factor w/ 127 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ therapies: Factor w/ 2 levels "A","B": 1 1 1 2 2 1 2 2 2 1 ...
##  $ Sex      : Factor w/ 2 levels "M","F": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Age      : int  28 32 41 21 34 24 28 40 24 39 ...
##  $ time     : Factor w/ 4 levels "R1","R2","R3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Painlevel: Factor w/ 5 levels "Nopain","Very mild",..: 4 4 3 4 4 3 4 3 4 4 ...

3.1 Cumulative of therapies

#knitr::kable(apply(with(dta_f, table(mental, ses)), 2, cumsum))
t(with(dtaL3, table(Painlevel, therapies)))
##          Painlevel
## therapies Nopain Very mild Mild Moderate Severe
##         A     56        17   65       87     27
##         B     63        45   67       59     22

3.2 freq figure

HH::likert(t(with(dtaL3, table(Painlevel, therapies))), main="", ylab="therapies")

4 model and comparison

library(VGAM)
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:tidyr':
## 
##     fill
m0 <- VGAM::vglm(ordered(Painlevel) ~ therapies, data=dtaL3, family=cumulative(parallel=TRUE))
m1 <- VGAM::vglm(ordered(Painlevel) ~ therapies, data=dtaL3, family=cumulative)
summary(m0)
## 
## Call:
## VGAM::vglm(formula = ordered(Painlevel) ~ therapies, family = cumulative(parallel = TRUE), 
##     data = dtaL3)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  -1.4361     0.1365 -10.518  < 2e-16 ***
## (Intercept):2  -0.8367     0.1255  -6.669 2.57e-11 ***
## (Intercept):3   0.2491     0.1201   2.074  0.03806 *  
## (Intercept):4   2.0308     0.1662  12.216  < 2e-16 ***
## therapiesB      0.4554     0.1594   2.857  0.00428 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), 
## logitlink(P[Y<=3]), logitlink(P[Y<=4])
## 
## Residual deviance: 1547.132 on 2027 degrees of freedom
## 
## Log-likelihood: -773.5662 on 2027 degrees of freedom
## 
## Number of Fisher scoring iterations: 4 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
## therapiesB 
##   1.576808
VGAM::lrtest(m0, m1)
## Likelihood ratio test
## 
## Model 1: ordered(Painlevel) ~ therapies
## Model 2: ordered(Painlevel) ~ therapies
##    #Df  LogLik Df Chisq Pr(>Chisq)  
## 1 2027 -773.57                      
## 2 2024 -767.93 -3 11.27    0.01035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(rms)
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
## 
## Attaching package: 'rms'
## The following objects are masked from 'package:VGAM':
## 
##     calibrate, lrtest
m0_po <- rms::lrm(Painlevel ~ therapies, data=dtaL3)
par(mfrow=c(1,2))
plot.xmean.ordinaly(Painlevel ~ therapies, data=dtaL3, subn=FALSE)
#partial propotional odds model