Ordinal Regression with plssem

This vignette demonstrates how to fit an ordinal regression model with plssem using the titanic dataset.

In plssem, regression-style model syntax like y ~ x1 + x2 is supported. When the dependent variable (and/or predictors) are ordinal. Ordinal variables can be supplied via the ordered argument, or by making sure they are ordered in the dataset.

Data

head(titanic[, c("Survived", "Age", "Sex", "Female", "Pclass")])
#>   Survived Age    Sex Female Pclass
#> 1        0  22   male      0      3
#> 2        1  38 female      1      1
#> 3        1  26 female      1      3
#> 4        1  35 female      1      1
#> 5        0  35   male      0      3
#> 6        0  NA   male      0      3

Linear ordinal regression

This model predicts survival as a function of age and sex.

m_linear <- "Survived ~ Age + Female"

fit_linear <- pls(
  m_linear,
  data          = titanic,
  ordered       = "Survived",
  boot.R        = 50,
  bootstrap     = TRUE,
  boot.parallel = "multicore",
  boot.ncores   = 2
)

summary(fit_linear)
#> plssem (0.1.2) ended normally after 1 iterations
#>   Estimator                                    OrdPLSc
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                           714
#>   Number of iterations                               1
#>   Number of latent variables                         0
#>   Number of observed variables                       3
#> 
#> Fit Measures:
#>   Chi-Square                                     0.000
#>   Degrees of Freedom                                 0
#>   SRMR                                           0.000
#>   RMSEA                                            NaN
#> 
#> R-squared (latents):
#>   Survived                                       0.347
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Survived ~    
#>     Age            -0.043      0.041   -1.061    0.289
#>     Female          0.584      0.029   20.035    0.000
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Age ~~        
#>     Female         -0.093      0.037   -2.524    0.012
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .Survived        0.653      0.035   18.486    0.000
#>     Age             1.000                             
#>     Female          1.000

Optional: evaluate predictive performance.

pls_predict(fit_linear, benchmark = "acc")
#> PlsSemPredict object
#> Available fields: $Y, $X.cont, $X.cont.pred, $X.ord, $X.ord.pred, $benchmark
#> 
#> Y [714 x 3] (head)
#>   Survvd    Age Female
#> 1 -0.420 -0.530 -0.759
#> 2  0.743  0.571  1.317
#> 3  0.779 -0.255  1.317
#> 4  0.752  0.365  1.317
#> 5 -0.458  0.365 -0.759
#> 7 -0.515  1.673 -0.759
#> 
#> X.cont [714 x 3] (head)
#>   Survvd    Age Female
#> 1 -0.653 -0.530 -0.759
#> 2  0.955  0.571  1.317
#> 3  0.955 -0.255  1.317
#> 4  0.955  0.365  1.317
#> 5 -0.653  0.365 -0.759
#> 7 -0.653  1.673 -0.759
#> 
#> X.cont.pred [714 x 3] (head)
#>   Survvd    Age Female
#> 1 -0.420 -0.530 -0.759
#> 2  0.743  0.571  1.317
#> 3  0.779 -0.255  1.317
#> 4  0.752  0.365  1.317
#> 5 -0.458  0.365 -0.759
#> 7 -0.515  1.673 -0.759
#> 
#> X.ord [714 x 3] (head)
#>   Survvd    Age Female
#> 1  1.000 -0.530 -0.759
#> 2  2.000  0.571  1.317
#> 3  2.000 -0.255  1.317
#> 4  2.000  0.365  1.317
#> 5  1.000  0.365 -0.759
#> 7  1.000  1.673 -0.759
#> 
#> X.ord.pred [714 x 3] (head)
#>   Survvd    Age Female
#> 1  1.000 -0.530 -0.759
#> 2  2.000  0.571  1.317
#> 3  2.000 -0.255  1.317
#> 4  2.000  0.365  1.317
#> 5  1.000  0.365 -0.759
#> 7  1.000  1.673 -0.759
#> 
#> Benchmark summary
#> acc: n=1, mean=0.780, median=0.780, min=0.780, max=0.780
#>  variable value
#>  Survived 0.780

Non-linear ordinal regression (interaction)

To include a non-linear (interaction) effect, add an interaction term. With ordinal indicators and interactions, plssem automatically switches to the Monte-Carlo ordinal PLSc estimator.

m_int <- "Survived ~ Age + Female + Age:Female"

fit_int <- pls(
  m_int,
  data          = titanic,
  ordered       = "Survived",
  boot.R        = 50,
  bootstrap     = TRUE,
  boot.parallel = "multicore",
  boot.ncores   = 2
)

summary(fit_int)
#> plssem (0.1.2) ended normally after 77 iterations
#>   Estimator                                  MCOrdPLSc
#>   Link                                          PROBIT
#>                                                       
#>   Number of observations                           714
#>   Number of iterations                              77
#>   Number of latent variables                         0
#>   Number of observed variables                       3
#> 
#> Fit Measures:
#>   Chi-Square                                     0.066
#>   Degrees of Freedom                                 0
#>   SRMR                                           0.003
#>   RMSEA                                            Inf
#> 
#> R-squared (latents):
#>   Survived                                       0.564
#> 
#> Regressions:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Survived ~    
#>     Age            -0.071      0.034   -2.069    0.039
#>     Female          0.697      0.036   19.295    0.000
#>     Age:Female      0.249      0.054    4.637    0.000
#> 
#> Covariances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>   Age ~~        
#>     Female         -0.092      0.039   -2.356    0.018
#>     Age:Female      0.008                             
#>   Female ~~     
#>     Age:Female      0.001                             
#> 
#> Variances:
#>                  Estimate  Std.Error  z.value  P(>|z|)
#>    .Survived        0.436                             
#>     Age             1.000                             
#>     Female          1.000                             
#>     Age:Female      1.046

pls_predict(fit_int, benchmark = "acc")
#> PlsSemPredict object
#> Available fields: $Y, $X.cont, $X.cont.pred, $X.ord, $X.ord.pred, $benchmark
#> 
#> Y [714 x 4] (head)
#>   Survvd    Age Female Ag:Fml
#> 1 -0.367 -0.530 -0.759  0.495
#> 2  1.087  0.571  1.317  0.845
#> 3  0.875 -0.255  1.317 -0.242
#> 4  1.034  0.365  1.317  0.574
#> 5 -0.600  0.365 -0.759 -0.184
#> 7 -0.940  1.673 -0.759 -1.176
#> 
#> X.cont [714 x 3] (head)
#>   Survvd    Age Female
#> 1 -0.638 -0.530 -0.759
#> 2  0.933  0.571  1.317
#> 3  0.933 -0.255  1.317
#> 4  0.933  0.365  1.317
#> 5 -0.638  0.365 -0.759
#> 7 -0.638  1.673 -0.759
#> 
#> X.cont.pred [714 x 3] (head)
#>   Survvd    Age Female
#> 1 -0.367 -0.530 -0.759
#> 2  1.087  0.571  1.317
#> 3  0.875 -0.255  1.317
#> 4  1.034  0.365  1.317
#> 5 -0.600  0.365 -0.759
#> 7 -0.940  1.673 -0.759
#> 
#> X.ord [714 x 3] (head)
#>   Survvd    Age Female
#> 1  1.000 -0.530 -0.759
#> 2  2.000  0.571  1.317
#> 3  2.000 -0.255  1.317
#> 4  2.000  0.365  1.317
#> 5  1.000  0.365 -0.759
#> 7  1.000  1.673 -0.759
#> 
#> X.ord.pred [714 x 3] (head)
#>   Survvd    Age Female
#> 1  1.000 -0.530 -0.759
#> 2  2.000  0.571  1.317
#> 3  2.000 -0.255  1.317
#> 4  2.000  0.365  1.317
#> 5  1.000  0.365 -0.759
#> 7  1.000  1.673 -0.759
#> 
#> Benchmark summary
#> acc: n=1, mean=0.780, median=0.780, min=0.780, max=0.780
#>  variable value
#>  Survived 0.780