Interpretación de los resultados de la regresión circular múltiple

Question

Soy nuevo en las estadísticas circulares y he estado jugando con circglmbayes y bpnreg . No estoy seguro de cómo interpretar el resultado. He consultado la base de datos CrossValidated y los documentos/tutoriales asociados a los paquetes mencionados. Sin embargo, mi limitada formación matemática me impide generalizar a mi caso.

Tengo una variable dependiente circular y dos variables independientes circulares. Mi objetivo es comparar la fuerza con la que las dos variables independientes predicen la variable dependiente:

$$outangle = b1*angle1 + b2*angle2$$

Entiendo que para los paquetes mencionados debo incluir dos predictores lineales por predictor circular:

$$outangle = b1*sin(angle1) + b2*cos(angle1) + b3*sin(angle2) + b4*cos(angle2)$$

Los modelos se ejecutan, pero no sé cómo interpretar sus resultados (ya que hay dos predictores por variable). Me gustaría poder decir algo así como "el ángulo 1 es el que más contribuye al ángulo de resultado, mientras que el ángulo 2 también contribuye, pero en menor medida". En otras palabras: ¿cómo puedo convertir la salida de los dos predictores lineales para representar un único coeficiente para la variable circular? Los resultados de los dos paquetes también son bastante diferentes. Entiendo que utilizan métodos diferentes, por ejemplo, bpnreg utiliza la distribución normal proyectada, pero en última instancia esperaría que produjeran la misma respuesta cualitativa.

A continuación se muestra un ejemplo simplificado de código y resultados:

library(bpnreg)
library(circglmbayes)

# test data
angles1  <- c(1,   1.1, 0.9, 2,   1.5, 2.5, 3.0, 0.5)
angles2  <- c(1.2, 0.7, 1.0, 2.3, 1.4, 2.8, 0.1, 0.2)
outangle <- c(1.1, 1,   0.8, 2.1, 1.4, 2.7, 3.1, 0.4)

# find linear components
sori1 <- sin(angles1)
cori1 <- cos(angles1)
sori2 <- sin(angles2)
cori2 <- cos(angles2)

# make dataframe
df  <- data.frame(outangle, sori1, cori1, sori2, cori2)

# run bpnr regression
circfit <- bpnr(outangle ~ 1 + sori1 + cori1 + sori2 + cori2, df)
circfit

# run circGLM regression
circfit2 <- circGLM(outangle ~ 1 + sori1 + cori1 + sori2 + cori2, df)
circfit2

salida (truncada) para bpnr:

Linear Coefficients 

Component I: 
                 mean      mode        sd     LB HPD    UB HPD
(Intercept) -31.38534 -34.43952 16.939204 -63.276156 -1.988348
sori1        53.25621  77.42373 29.152543  -2.583563 96.693508
cori1        12.73618  12.34990  5.416473   4.027589 24.099637
sori2       -21.23710 -27.35492 16.259473 -48.074359  9.390253
cori2        21.34205  28.25890 11.224285  -2.308050 36.577529

Component II: 
                 mean      mode        sd     LB HPD      UB HPD
(Intercept) -33.16929 -49.09227 19.125987 -58.474295   1.0513257
sori1       108.48108 156.23939 56.547787   5.689475 182.3388648
cori1       -12.65639 -17.02416  7.716536 -24.791198   0.5974747
sori2       -43.46557 -68.49125 29.332958 -83.259452   8.4257935
cori2        10.43962  16.47495  6.794084  -1.659991  20.2801103

Circular Coefficients 

Continuous variables: 
         mean ax    mode ax     sd ax      LB ax     UB ax
sori1  0.3582156  0.3453063 0.2579013  0.1491464 0.6097201
cori1 -0.3463878  0.4257769 0.9575427 -2.0946958 1.1603828
sori2 -0.6953270 -0.8023903 1.6901613 -2.0029269 1.5081406
cori2  1.6737354  1.7497827 1.2241280 -1.0348414 3.3354418

salida (truncada) para circGLM:

Coefficients:
          Estimate     SD     LB      UB
Intercept    1.604  0.061  1.493   1.712
Kappa       56.128 47.798  2.052 168.050
sori1       -0.150  0.099 -0.359   0.039
cori1       -0.522  0.092 -0.680  -0.344
sori2        0.108  0.077 -0.049   0.250
cori2        0.016  0.072 -0.129   0.189

Answer 1

La respuesta breve es que no hay forma de calcular un coeficiente singular a partir de un predictor circular, de modo que si queremos cuantificar el efecto de un único predictor circular, tenemos que utilizar algo como el ajuste del modelo o los criterios de información.

Por qué tenemos este problema

Interpretar los coeficientes cos y sin por separado no funcionará, ya que eso significaría que si rotamos el predictor, lo que no debería cambiar la fuerza de su efecto, obtendríamos intervalos de credibilidad diferentes para los componentes del predictor cos y sin. Además, en casi todos los casos, el modelo con sólo uno de los dos elementos del predictor (cos y sin) no tiene sentido y no debe interpretarse. Por lo tanto, deberíamos añadir sólo ambos o no añadir ninguno al modelo.

Por último, debo tratar de evitar confusiones: el bpnreg no sólo tiene un componente cos y otro sin del predictor; también divide el resultado circular en Componente I (cos) y Componente II (sin).

Cómo comparar modelos con criterios de ajuste o de información

El enfoque (que simplificaré para mayor claridad para tener sólo un predictor circular en lugar de los dos que tiene) podría ser algo como esto:

• Ejecute el modelo sin el predictor circular, llámelo M1.
• Ejecute el modelo con el predictor circular, llámelo M2.
• Compara los dos modelos M1 y M2. Podemos hacerlo de varias maneras, en complejidad creciente:
  • Compare el DIC de un modelo con el del siguiente y prefiera el modelo con el BIC más bajo. Esto se da por defecto en ambos paquetes.
  • Igual que arriba, pero con WAIC. También lo ofrecen ambos paquetes y, por lo general, es un poco mejor que el DIC.
  • Calcule la verosimilitud marginal de cada modelo y el factor de Bayes asociado, que es la relación entre las verosimilitudes marginales. Se trata de una opción del paquete circbayes pero ese paquete es todavía algo experimental. También puede hacerlo usted mismo con el paquete bridgesampling pero es más complicado. Es la mejor opción, pero la más compleja de aplicar.
  • Si el objetivo final es la predicción activa, la validación cruzada también es una opción.

Al tener dos predictores circulares, el proceso es el mismo, pero hay cuatro modelos para comparar.

(Revelación completa: soy el autor de circglmbayes .)