Generalized ordered regression with clustered standard error using the vglm function

Question

I am working with an ordinal outcome variable (state_prev) measuring sexual violence prevalence, using a generalized ordered regression. I fitted the model: sup_state <- vglm(state_prev ~ ext_sup + ext_f + ANIpC + GDPpC + unemployment + population, family = cumulative(parallel = FALSE, reverse = FALSE), data = SVAC_ext_supp)

As I am doing an analysis in a conflict context I want to add clustered standard errors (variable: conflictid). But i was unable to do so. Can anyone help implementing clustered standard errors in a generalized ordered regression?

Thanks so much!

• method: generalized ordered regression
• model: sup_state <- vglm(state_prev ~ ext_sup + ext_f + ANIpC + GDPpC + unemployment + population, family = cumulative(parallel = FALSE, reverse = FALSE), data = SVAC_ext_supp)

coeftest() (package sandwich) is not comptible with a vglm
se <- coeftest(base_state, vcov = vcovCL, cluster = ~conflictid)

ordLORgee() (package multgee) gets an error message: vector memory exhausted (limit reached?)

se <- ordLORgee(formula = state_prev ~ factor(ext_sup) + factor(ext_f), link = "logit",
                    data = SVAC_ext_supp, id = conflictid, LORstr = "uniform").

robcov does not work either (package rms). i get:

Error in h(simpleError(msg, call)) : 
 error in evaluating the argument 'x' in selecting a method for function 'as.matrix': $ operator not defined for this S4 class

robcov(base_state, cluster=SVAC_data$conflictid) 

I don't want to use polr() (package mass) because i cannot satisfy the proportional odds assumption. If i do, i get the error:

Re-fitting to get Hessian
Error in svd(X) : infinite or missing values in 'x'

base_state <- polr(as.factor(state_prev) ~ ext_sup + ext_f + ANIpC + GDPpC + poverty + unemployment + population, data = SVAC_ext_supp)

se <- coeftest(base_state, vcov = vcovCL, cluster = ~conflictid) 

lm.cluster() is not adequate for ordinal outcomes

the only thing that works is a binary sexual violence outcome variable and using glm()

mod <- glm(SV_bin ~ ext_sup + ext_f, family = binomial(link = "logit"), data = SVAC_ext_supp)

vcov_cluster <- vcovCL(mod, cluster = ~conflictid)

coeftest(m, vcov. = vcov_cluster)

but has anyone an idea how i can implement clustered standard errors into a generalized ordered regression?