mvord with autorregresive structure-CodePudding

This is my model:

library (mvord)
roedores2 <- mvord(datos(IndicAD) ~ Trat   Tiempo   
                    Trat:Tiempo (1|Tiempo)   (1|Granja), data =  datos, 
                    error.structure = cor_ar1(~Tiempo), link = mvprobit(), 
                    control = mvord.control(solver = "BFGS"))

"IndicAD" is a ordinal variable with 0,1,2,3,4 as posible responses

"Trat" is a factor with 3 levels

"Tiempo" is a factor (repeated measurement over time) with 11 levels (times)

"Granja" is a factor with two levels (two sites of sample)

But I keep getting this error:

Error in seq_len(rho$ndim) : 
  argument must be coercible to non-negative integer
In addition: Warning message:
In seq_len(rho$ndim) : first element used of 'length.out' argument

CodePudding user response：

The mvord package is working with multivariate ordinal regression with the meaning that response variable should contain more than one factor. In your case it is only one indcAD so mvord function fails.

In case you have two ordered factors indcAD1 and indcAD2 the algorith would work OK. In case you need to use univariate ordinal regression where response value is one factor you can use e.g. polr by MASS package.

library (mvord)
n <- 50
datos <- data.frame(
  IndiCAD1 = factor(sample(0:4, n, replace = TRUE), ordered = TRUE),
  IndiCAD2 = factor(sample(0:4, n, replace = TRUE), ordered = TRUE),
  Trat = factor(sample(3, n, replace = TRUE), ordered = TRUE),
  Tiempo = factor(sample(3, n, replace = TRUE)),
  Granja = factor(sample(2, n, replace = TRUE)))

res_mvord <-  mvord(MMO2(IndiCAD1, IndiCAD2) ~ Trat   Tiempo   
                 Trat:Tiempo, data =  datos, 
               error.structure = cor_ar1(~Tiempo), link = mvprobit(), 
               control = mvord.control(solver = "BFGS"))


res_mvord
# Call:
#   mvord(formula = MMO2(IndiCAD1, IndiCAD2) ~ Trat   Tiempo   Trat:Tiempo, 
#         data = datos, error.structure = cor_ar1(~Tiempo), link = mvprobit(), 
#         control = mvord.control(solver = "BFGS"))
# 
# Thresholds:
#   $IndiCAD1
# 0|1       1|2       2|3       3|4 
# 0.0000000 0.7113659 1.0971497 1.7800162 
# 
# $IndiCAD2
# 0|1      1|2      2|3      3|4 
# 0.000000 1.088273 1.486173 2.135159 
# 
# Coefficients:
#   (Intercept) 1    (Intercept) 2         Trat.L 1         Trat.L 2         Trat.Q 1 
# 0.7482449        1.2190449       -0.4892506        0.3545706        0.7432251 
# Trat.Q 2        Tiempo2 1        Tiempo2 2        Tiempo3 1        Tiempo3 2 
# -0.5710266        0.4844226       -0.6121860       -0.2728416       -0.3183369 
# Trat.L:Tiempo2 1 Trat.L:Tiempo2 2 Trat.Q:Tiempo2 1 Trat.Q:Tiempo2 2 Trat.L:Tiempo3 1 
# -0.5196205       -0.4058554       -0.4065644        0.6272309        0.3727869 
# Trat.L:Tiempo3 2 Trat.Q:Tiempo3 1 Trat.Q:Tiempo3 2 
# -1.6488805       -0.1636811        0.1608613 
# 
# Parameters of the error structure:
#   p25        p26        p27 
# 0.3447173 -0.3866346 -0.8193181 

res_polar <- polr(IndiCAD1 ~ Trat   Tiempo   
                    Trat:Tiempo, data = datos)

res_polar
# Call:
#   polr(formula = IndiCAD1 ~ Trat   Tiempo   Trat:Tiempo, data = datos)
# 
# Coefficients:
#   Trat.L         Trat.Q        Tiempo2        Tiempo3 Trat.L:Tiempo2 Trat.Q:Tiempo2 
# -0.84715296     1.26605745     0.72589772    -0.60879038    -0.74903326    -0.69633321 
# Trat.L:Tiempo3 Trat.Q:Tiempo3 
# 0.76901029    -0.09709166 
# 
# Intercepts:
#   0|1         1|2         2|3         3|4 
# -1.28052686 -0.05745217  0.58771283  1.71795395 
# 
# Residual Deviance: 150.1981 
# AIC: 174.1981

CodePudding user response：

Thanks for your answer. I thought that "mvord" might not work.

For ordinal models I use "clm" or "clmm" (if I have random factors), both from the "ordinal" package. The problem is that in both cases I don't see how to include the autoregressive correlation structure. In the example you give me, I don't see how you're including it either.

Does anyone know of any other alternatives?