¿Qué prueba puedo utilizar para comparar las pendientes de dos o más modelos de regresión?

Question

Me gustaría probar la diferencia en la respuesta de dos variables a un predictor. He aquí un ejemplo mínimo reproducible.

library(nlme) 
## gls is used in the application; lm would suffice for this example
m.set <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "setosa")
m.vir <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "virginica")
m.ver <- gls(Sepal.Length ~ Petal.Width, data = iris, 
               subset = Species == "versicolor")

Puedo ver que los coeficientes de las pendientes son diferentes:

m.set$coefficients
(Intercept) Petal.Width 
  4.7771775   0.9301727
m.vir$coefficients
(Intercept) Petal.Width 
  5.2694172   0.6508306 
m.ver$coefficients
(Intercept) Petal.Width 
   4.044640    1.426365 

Tengo tres preguntas:

1. ¿Cómo puedo comprobar la diferencia entre las pendientes?
2. ¿Cómo puedo comprobar la diferencia entre las varianzas residuales?
3. ¿Cuál es una forma sencilla y eficaz de presentar estas comparaciones?

Una pregunta relacionada, Método para comparar el coeficiente de una variable en dos modelos de regresión El Sr. G. B., sugiere que se vuelva a ejecutar el modelo con una variable ficticia para diferenciar las pendientes, ¿existen opciones que permitan utilizar conjuntos de datos independientes?

Answer 1

Para responder a estas preguntas con el código R, utilice lo siguiente:
1. ¿Cómo puedo comprobar la diferencia entre las pendientes?
Respuesta: Examine el valor p del ANOVA de la interacción de Pétalo.Ancho por Especie, luego compare las pendientes usando lsmeans::lstrends, como sigue.

library(lsmeans)
m.interaction <- lm(Sepal.Length ~ Petal.Width*Species, data = iris)
anova(m.interaction)
 Analysis of Variance Table

 Response: Sepal.Length
                      Df Sum Sq Mean Sq  F value Pr(>F)    
 Petal.Width           1 68.353  68.353 298.0784 <2e-16 ***
 Species               2  0.035   0.017   0.0754 0.9274    
 Petal.Width:Species   2  0.759   0.380   1.6552 0.1947    
 Residuals           144 33.021   0.229                    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Obtain slopes
m.interaction$coefficients
m.lst <- lstrends(m.interaction, "Species", var="Petal.Width")
 Species    Petal.Width.trend        SE  df   lower.CL upper.CL
 setosa             0.9301727 0.6491360 144 -0.3528933 2.213239
 versicolor         1.4263647 0.3459350 144  0.7425981 2.110131
 virginica          0.6508306 0.2490791 144  0.1585071 1.143154

# Compare slopes
pairs(m.lst)
 contrast                 estimate        SE  df t.ratio p.value
 setosa - versicolor    -0.4961919 0.7355601 144  -0.675  0.7786
 setosa - virginica      0.2793421 0.6952826 144   0.402  0.9149
 versicolor - virginica  0.7755341 0.4262762 144   1.819  0.1669

2. ¿Cómo puedo comprobar la diferencia entre las varianzas residuales?
Si entiendo la pregunta, se pueden comparar las correlaciones de Pearson con una transformada de Fisher, también llamada "r-z de Fisher", de la siguiente manera.

library(psych)
library(data.table)
iris <- as.data.table(iris)
# Calculate Pearson's R
m.correlations <- iris[, cor(Sepal.Length, Petal.Width), by = Species]
m.correlations
# Compare R values with Fisher's R to Z
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("setosa", "versicolor"), .N])
paired.r(m.correlations[Species=="setosa", V1], m.correlations[Species=="virginica", V1], 
         n = iris[Species %in% c("setosa", "virginica"), .N])
paired.r(m.correlations[Species=="virginica", V1], m.correlations[Species=="versicolor", V1], 
         n = iris[Species %in% c("virginica", "versicolor"), .N])

3. ¿Cuál es una forma sencilla y eficaz de presentar estas comparaciones?
"Utilizamos una regresión lineal para comparar la relación entre la longitud de los sépalos y la anchura de los pétalos para cada especie. No encontramos una interacción significativa en las relaciones entre la longitud del sépalo y la anchura del pétalo para I. Setosa (B = 0.9), I. Versicolor (B = 1,4), ni I. Virginica (B = 0,6); F (2, 144) = 1,6, p = 0,19. Una comparación de r a z de Fisher indicó que la correlación de Pearson para I. Setosa (r = 0,28) fue significativamente menor (p = 0,02) que I. Versicolor (r = 0.55). Del mismo modo, la correlación para I. Virginica (r = 0,28) fue significativamente más débil (p = 0,02) que la observada para I. Versicolor ."

Por último, visualiza siempre tus resultados.

plotly_interaction <- function(data, x, y, category, colors = col2rgb(viridis(nlevels(as.factor(data[[category]])))), ...) {
  # Create Plotly scatter plot of x vs y, with separate lines for each level of the categorical variable. 
  # In other words, create an interaction scatter plot.
  # The "colors" must be supplied in a RGB triplet, as produced by col2rgb().

  require(plotly)
  require(viridis)
  require(broom)

  groups <- unique(data[[category]])

  p <- plot_ly(...)

  for (i in 1:length(groups)) {
    groupData = data[which(data[[category]]==groups[[i]]), ]
    p <- add_lines(p, data = groupData,
                   y = fitted(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                   x = groupData[[x]],
                   line = list(color = paste('rgb', '(', paste(colors[, i], collapse = ", "), ')')),
                   name = groups[[i]],
                   showlegend = FALSE)
    p <- add_ribbons(p, data = augment(lm(data = groupData, groupData[[y]] ~ groupData[[x]])),
                     y = groupData[[y]],
                     x = groupData[[x]],
                     ymin = ~.fitted - 1.96 * .se.fit,
                     ymax = ~.fitted + 1.96 * .se.fit,
                     line = list(color = paste('rgba','(', paste(colors[, i], collapse = ", "), ', 0.05)')), 
                     fillcolor = paste('rgba', '(', paste(colors[, i], collapse = ", "), ', 0.1)'),
                     showlegend = FALSE)
    p <- add_markers(p, data = groupData, 
                     x = groupData[[x]], 
                     y = groupData[[y]],
                     symbol = groupData[[category]],
                     marker = list(color=paste('rgb','(', paste(colors[, i], collapse = ", "))))
  }
  p <- layout(p, xaxis = list(title = x), yaxis = list(title = y))
  return(p)
}

plotly_interaction(iris, "Sepal.Length", "Petal.Width", "Species")

Answer 2

¿Cómo puedo comprobar la diferencia entre las pendientes?

Incluir un dummy para las especies, dejar que interactúe con $P_i$ y ver si esta variable ficticia es significativa. Sea $L_i$ sea la longitud del sépalo y $P_i$ sea la anchura del pedal y $S_1, S_2, S_3$ son las variables ficticias para las tres especies. La comparación del modelo

$$ E(L_i) = \beta_0 + \beta_1 P_i $$

con el modelo que permite el efecto de $P_i$ ser diferente para cada especie:

$$ E(L_i) = \alpha_0 + \alpha_1 S_2 + \alpha_2 S_3 + \alpha_4 P_i + \alpha_5 P_iS_2 + \alpha_6 P_i S_3 $$

Los estimadores GLS son MLEs y el primer modelo es un submodelo del segundo, por lo que se puede utilizar la prueba de razón de verosimilitudes aquí. Las verosimilitudes pueden extraerse utilizando la función logLik y los grados de libertad de la prueba serán $4$ ya que has borrado $4$ parámetros para llegar al submodelo.

¿Cuál es una forma sencilla y eficaz de presentar la comparación?

Creo que la forma más atractiva sería trazar las líneas de regresión para cada especie en los mismos ejes, quizás con barras de error basadas en los errores estándar. Esto haría que el diferencia (o no diferencia) entre las especies y su relación con $P_i$ muy evidente.

Editar: Me he dado cuenta de que se ha añadido otra pregunta al cuerpo. Por lo tanto, voy a añadir una respuesta a eso:

¿Cómo puedo comprobar la diferencia entre las varianzas residuales?

Para ello, tendrá que estratificar el conjunto de datos y ajustar modelos separados, ya que el modelo basado en la interacción que sugerí restringirá la varianza residual para que sea la misma en cada grupo. Si ajusta modelos separados, esta restricción desaparece. En ese caso, puede seguir utilizando la prueba de razón de verosimilitudes (la verosimilitud del modelo más amplio se calcula ahora sumando las verosimilitudes de los tres modelos separados). El modelo "nulo" depende de con qué se quiera comparar

• si usted sólo quiere probar la varianza, dejando los efectos principales, entonces el modelo "nulo" debería ser el modelo con las interacciones que he escrito arriba. Los grados de libertad para la prueba son entonces $2$ .

• Si quieres probar la varianza conjuntamente con los coeficientes, entonces el modelo nulo debería ser el primer modelo que he escrito arriba. Los grados de libertad para la prueba son entonces $6$ .

Answer 3

Estoy de acuerdo con la sugerencia anterior. Debería ajustar un modelo de regresión múltiple con una variable ficticia para cada conjunto de datos. Esto le permitirá comprobar si los interceptos difieren. Si además quiere saber si los pendientes difieren, entonces hay que incluir también las interacciones entre las variables ficticias y la variable en cuestión. No hay ningún problema con el hecho de que los datos sean independientes. Tenga en cuenta que si son ambos independiente y (por ejemplo) diferentes especies, entonces no se podría saber si la diferencia que se encuentra se debe a las diferentes especies o a los diferentes conjuntos de datos, ya que están perfectamente confundidos. Sin embargo, hay ninguna prueba / tarjeta de salida de la cárcel que le permitirá sortear ese problema sin tener que recoger una nueva muestra y volver a realizar su estudio.

Answer 4

Este enlace ofrece una respuesta sencilla a la pregunta:

https://stat.ethz.ch/pipermail/r-sig-teaching/2011q4/000387.html

Answer 5

Estoy haciendo un tipo de análisis similar, teniendo regiones metiladas de ADN y tengo una pregunta.

En resumen, quiero ver si las pendientes son diferentes en función del fenotipo (variable binaria) y de la edad de los pacientes.

Basado en la respuesta de Kayle Sawyer Tengo una fórmula parecida a esta model=lm(característica~edad*alergia)

Ahora supongamos que quiero utilizar el sexo del paciente como factor de confusión.

¿Puede decirme cómo debo transformar mi modelo? ¿Es correcto algo así? model=lm(característica~edad*alergia + sexo)

Muchas gracias

Leonardos