Dear all, is there existing code to fit a ordinal regression model with probit link where the variance can be modeled? For example, the polr function from the MASS package fits: probitP(Y <= k | x) = zeta_k - eta Instead of, probitP(Y <= k | x) = (zeta_k - eta)/scale, where scale = exp(b*z) which is what I need. I've seen these models go by the name of heteroskedastic ordered probit, location-scale glms, generalized nonlinear models, ... . I know the package gnlm (Jim Lindsey) has a function to fit this model using the logit link, so I'm wondering if this is already done in R using probit instead. Thanks. Regina