mvord with autorregresive structure

Question

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

Answer 1

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