I have a simple script in R and I’m trying to run a bi-variate fixed effect ordered probit model using mvord:

# Load the data
safety_comfort_full <- read.csv("/Users/..../Data/roundabout_safety_comfort_responses.csv")

# Inspect the structure and summary of the data
str(safety_comfort_full)
summary(safety_comfort_full)

# Calculate z-scores and create new variables
safety_comfort_full <- safety_comfort_full %>%
  mutate(
    Violations = round(rowMeans(scale(select(., Q1_1, Q1_2, Q1_3)))),
    Errors = round(rowMeans(scale(select(., Q1_4, Q1_5, Q1_6, Q1_7, Q1_8, Q1_9)))),
    Positive_behavior = round(rowMeans(scale(select(., Q1_10, Q1_11, Q1_12))))
  )

# Convert `ResponseId` to a factor
safety_comfort_full$ResponseId <- as.factor(safety_comfort_full$ResponseId)

# Make sure 'Comfort' and 'Safety' are ordered factors
safety_comfort_full$Comfort <- as.ordered(safety_comfort_full$Comfort)
safety_comfort_full$Safety <- as.ordered(safety_comfort_full$Safety)


# Fit the multivariate ordinal regression model with fixed effects
fit <- mvord(
  formula = MMO2(Safety,Comfort) ~   Separation_from_Vehicle_Lane + Bike_Yielding_Priority +     Art_Ads_in_Center + Uni_Directional + Bike_Volumes + Facility_Style + X4_Legged 
  + Q23_2 + Q3 + Q4 + Q6 + Q32 + Q33 + Q5 + Q23_1 + Q23_3 + Q23_4 + Q23_5 + 
  Q23_6 + Q23_7 + Q23_8 + Violations + Errors + Positive_behavior + Roundabout_Number + ResponseId,  
  data = safety_comfort_full,
  error.structure = cor_general(~ 1),
  link = mvprobit(),
  threshold.constraints = NULL,
  control = mvord.control()
)

summary(fit)

However I receive an error when running the model:

Error in chol.default(derivs_for_se$H) : the leading minor of order 373 is not positive definite

I found this question from a couple of years ago but they were using a different library and I’m not sure how to solve it using the mvord library.

I’ve tried to implement regularizing means for the independent variables but it doesn’t seem to work.