Que teste posso usar para comparar inclinações de dois ou mais modelos de regressão?

Eu gostaria de testar a diferença na resposta de duas variáveis para um preditor. Aqui está um exemplo mínimo reproduzível.

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")

Eu posso ver que os coeficientes da inclinação são 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

Eu tenho três perguntas:

Como posso testar a diferença entre pistas?
Como posso testar a diferença entre as variações residuais?
Qual é uma maneira simples e eficaz de apresentar essas comparações?

Uma questão relacionada, Método para comparar o coeficiente de variável em dois modelos de regressão , sugere re-executar o modelo com uma variável dummy para diferenciar as inclinações; existem opções que permitiriam o uso de conjuntos de dados independentes?

r data-visualization multivariate-analysis hypothesis-testing Abe
fonte

Em relação à primeira pergunta, consulte stats.stackexchange.com/questions/55501/… .

russellpierce

Respostas:

Como posso testar a diferença entre pistas?

Incluir um manequim para a espécie, deixá-lo interagir com , e ver se este manequim é significativo. Vamos ser o comprimento sépalas e seja a largura do pedal e ser as variáveis indicadoras para as três espécies. A comparar o modelo $P_i$ $L_i$ $P_i$ $S_1, S_2, S_3$

E (L_{i}) = β_{0} + β_{1} P_{i}

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

$P_i$

E (L_{i}) = α_{0} + α_{1} S_{2} + α_{2} S_{3} + α_{4} P_{i} + α_{5} P_{i} S_{2} + α_{6} P_{i} S_{3}

$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$

logLik $4$ $4$

Qual é uma maneira simples e eficaz de apresentar a comparação?

$P_i$

Editar: notei que outra pergunta foi adicionada ao corpo. Então, eu estou adicionando uma resposta para isso:

Como posso testar a diferença entre as variações residuais?

Para isso, você precisará estratificar o conjunto de dados e ajustar modelos separados, pois o modelo baseado em interação que eu sugeri restringirá a variação residual para ser a mesma em todos os grupos. Se você ajustar modelos separados, essa restrição desaparecerá. Nesse caso, você ainda pode usar o teste da razão de verossimilhança (a probabilidade do modelo maior agora é calculada somando as probabilidades dos três modelos separados). O modelo "nulo" depende do que você deseja compará-lo

Macro
fonte

S_{1}

$S_1$ gls(Sepal.Length ~ species:Petal.Width, data = iris)

S_{1}

$S_1$

α_{0} + α_{4} P_{i}

$\alpha_0 + \alpha_4 P_i$ . Se speciesfor uma variável categórica, acho gls(Sepal.Length ~ species*Petal.Width, data=iris)que seria a sintaxe.

Macro

@ Macro Boa resposta (+1)! Gostaria de saber se você poderia ajustar o glsmodelo, mas permitindo diferentes variações residuais para cada espécie com a opção weights=varIdent(form=~1|Species)(em relação à segunda pergunta)?

COOLSerdash

Para responder a essas perguntas com o código R, use o seguinte:
1. Como posso testar a diferença entre inclinações?
Resposta: Examine o valor de p ANOVA da interação Petal.Width by Species e compare as inclinações usando lsmeans :: lstrends, como a seguir.

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. Como posso testar a diferença entre as variações residuais?
Se eu entendi a pergunta, você pode comparar as correlações de Pearson com uma transformação de Fisher, também chamada de "r-to-z de Fisher", como segue.

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. Qual é uma maneira simples e eficaz de apresentar essas comparações?
"Utilizamos regressão linear para comparar a relação entre o comprimento do septo e a largura da pétala para cada espécie. Não encontramos uma interação significativa nas relações entre o comprimento do sepal e a largura da pétala para I. Setosa (B = 0,9), I. Versicolor (B = 1,4), nem I. Virginica (B = 0,6); F (2, 144) = 1,6, p = 0,19 Uma comparação de Fisher para z indicou que a correlação de Pearson para I. Setosa (r = 0,28) foi significativamente menor (p = 0,02) do que I. Versicolor (r = 0,55). Da mesma forma, a correlação para I. Virginica (r = 0,28) foi significativamente mais fraca (p = 0,02) do que a observada para I. Versicolor. "

Por fim, sempre visualize seus 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")

Kayle Sawyer
fonte

Eu concordo com a sugestão anterior. Você deve ajustar um modelo de regressão múltipla com uma variável dummy para cada conjunto de dados. Isso permitirá que você teste se as interceptações diferem. Se você também deseja saber se o cartão / sai da cadeia sem solução que o levará a solucionar esse problema sem reunir uma nova amostra e executar o estudo novamente. inclinações diferem, precisará incluir também interações entre os manequins e a variável em questão. Não há problema com o fato de os dados serem independentes. Nota que, se eles são ambos independente e (por exemplo) de espécies diferentes, então você não seria capaz de dizer se a diferença que você encontrará é devido às espécies diferentes ou os conjuntos de dados diferentes, como eles são perfeitamente confundido. No entanto, não há teste

- Reinstate Monica
fonte

Parece que postamos respostas bastante semelhantes quase ao mesmo tempo. +1

Macro

@ Macro, sim, mas o seu é principalmente melhor (+1 antes); você abordou todas as três perguntas que eu perdi na minha primeira (não completa) leitura da questão. Minha contribuição aqui é a parte sobre confusão.

gung - Restabelece Monica

Sim, esse é um bom argumento. Suponho que, se você estivesse fazendo essa investigação, teria que estar operando sob a suposição de que os conjuntos de dados estavam medindo a mesma coisa, etc ... com a única diferença de que as espécies eram diferentes.

Macro

Do meu modo de pensar, vocês dois deveriam receber votos positivos, que é o que estou fazendo.

22812 Michael R. Chernick

A sugestão da variável dummy é boa, desde que a variação do erro não seja significativamente diferente entre os modelos. Caso contrário, você poderá aplicar um teste t Satterthwaite-Welch (que tem a vantagem singular de estar disponível quando apenas estatísticas resumidas são conhecidas, como costuma ser o caso ao ler artigos publicados) ou usar mínimos quadrados ponderados para se ajustar ao modelo combinado.

whuber