R: Cálculo da média e erro padrão da média para fatores com lm ​​() vs.

8

Ao lidar com dados com fatores R, pode ser usado para calcular as médias de cada grupo com a função lm (). Isso também fornece os erros padrão para as médias estimadas. Mas esse erro padrão difere do que recebo de um cálculo manualmente.

Aqui está um exemplo (retirado daqui Prevendo a diferença entre dois grupos em R )

Primeiro calcule a média com lm ​​():

    mtcars$cyl <- factor(mtcars$cyl)
    mylm <- lm(mpg ~ cyl, data = mtcars)
    summary(mylm)$coef

                Estimate Std. Error   t value     Pr(>|t|)
  (Intercept)  26.663636  0.9718008 27.437347 2.688358e-22
  cyl6         -6.920779  1.5583482 -4.441099 1.194696e-04
  cyl8        -11.563636  1.2986235 -8.904534 8.568209e-10

A interceptação é a média para o primeiro grupo, os 4 carros com cilindros. Para obter os meios por cálculo direto, uso o seguinte:

  with(mtcars, tapply(mpg, cyl, mean))

         4        6        8 
    26.66364 19.74286 15.10000 

Para obter os erros padrão das médias, calculo a variação padrão da amostra e divido pelo número de observações em cada grupo:

 with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )

         4         6         8 
   1.3597642 0.5493967 0.6842016 

O cálculo direto dá a mesma média, mas o erro padrão é diferente para as duas abordagens, eu esperava obter o mesmo erro padrão. O que está acontecendo aqui? Está relacionado ao lm () ajustar a média para cada grupo e um termo de erro?

Editado: Após a resposta de Svens (abaixo), posso formular minha pergunta de forma mais concisa e clara.

Para dados categóricos, podemos calcular a média de uma variável para diferentes grupos usando lm () sem interceptação.

  mtcars$cyl <- factor(mtcars$cyl)
  mylm <- lm(mpg ~ cyl, data = mtcars)
  summary(mylm)$coef

      Estimate Std. Error
  cyl4 26.66364  0.9718008
  cyl6 19.74286  1.2182168
  cyl8 15.10000  0.8614094

Podemos comparar isso com um cálculo direto das médias e seus erros padrão:

  with(mtcars, tapply(mpg, cyl, mean))

         4        6        8 
    26.66364 19.74286 15.10000 

  with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )

         4         6         8 
   1.3597642 0.5493967 0.6842016 

Os meios são exatamente os mesmos, mas os erros padrão são diferentes para esses 2 métodos (como Sven também observa). Minha pergunta é por que eles são diferentes e não são os mesmos?

(ao editar minha pergunta, devo excluir o texto original ou adicionar minha edição como fiz)

SRJ
fonte

Respostas:

7

A diferença nos erros padrão ocorre porque na regressão você calcula uma estimativa combinada da variação, enquanto no outro cálculo você calcula estimativas separadas da variação.

Glen_b -Reinstate Monica
fonte
2
Obrigado por esclarecer isso. Acabei de encontrar uma resposta muito boa para uma pergunta semelhante aqui, com um bom exemplo: stats.stackexchange.com/questions/29479/…
SRJ:
Sim, isso parece relevante. Bem manchado.
Glen_b -Reinstala Monica
5

A lmfunção não estima médias e erros padrão dos níveis de fator, mas dos controles associados aos níveis de fator.

Se nenhum contraste for especificado manualmente, os contrastes de tratamento serão usados ​​em R. Esse é o padrão para dados categóricos.

O fator mtcars$cylpossui três níveis (4,6 e 8). Por padrão, o primeiro nível, 4, é usado como categoria de referência. A interceptação do modelo linear corresponde à média da variável dependente na categoria de referência. Mas os outros efeitos resultam da comparação de um nível de fator com a categoria de referência. Portanto, a estimativa e o erro padrão para cyl6estão relacionados à diferença entre cyl == 6e cyl == 4. O efeito cyl8está relacionado à diferença entre cyl == 8e cyl == 4.

Se você deseja que a lmfunção calcule as médias dos níveis de fator, você deve excluir o termo de interceptação ( 0 + ...):

summary(lm(mpg ~ 0 + as.factor(cyl), mtcars))

Call:
lm(formula = mpg ~ 0 + as.factor(cyl), data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2636 -1.8357  0.0286  1.3893  7.2364 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
as.factor(cyl)4  26.6636     0.9718   27.44  < 2e-16 ***
as.factor(cyl)6  19.7429     1.2182   16.21 4.49e-16 ***
as.factor(cyl)8  15.1000     0.8614   17.53  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Residual standard error: 3.223 on 29 degrees of freedom
Multiple R-squared: 0.9785, Adjusted R-squared: 0.9763 
F-statistic: 440.9 on 3 and 29 DF,  p-value: < 2.2e-16 

Como você pode ver, essas estimativas são idênticas às médias dos níveis dos fatores. Mas observe que os erros padrão das estimativas não são idênticos aos erros padrão dos dados.

A propósito: Os dados podem ser facilmente agregados com a aggregatefunção:

aggregate(mpg ~ cyl, mtcars, function(x) c(M = mean(x), SE = sd(x)/sqrt(length(x))))

  cyl      mpg.M     mpg.SE
1   4 26.6636364  1.3597642
2   6 19.7428571  0.5493967
3   8 15.1000000  0.6842016
Sven Hohenstein
fonte
Obrigado pela resposta. Eu já sei que coeficientes não são os meios, como escrevi, a interceptação é a média do primeiro nível, os outros coeficientes são a diferença na média dos outros níveis para esse nível. Você também notou que, com sua observação, "os erros padrão das estimativas não são idênticos aos erros padrão dos dados". Isso significa que lm () estima as médias e calcula erros padrão para essas estimativas?
SRJ
Ops, eu queria editar esse comentário para fins de clareza, mas não sabia que só podia editar por 5 minutos. Posso excluir um comentário? Eu não sabia que poderia obter estimativas médias diretamente, omitindo a interceptação, obrigado por essa dica. Se o entendi corretamente, os erros padrão das médias estimadas não são os mesmos que os erros padrão calculados diretamente a partir dos dados. É um conjunto diferente de equações usado em cada caso? E quais são essas equações? Eu gostaria de ter mais alguns detalhes para u nderstand a diferença melhor
SRJ
1

Além do que Sven Hohenstein disse, os mtcarsdados não são equilibrados . Normalmente, é usado aovpara lm com dados categóricos (que são apenas um invólucro lm) que diz especificamente sobre ?aov:

o aov é projetado para projetos balanceados e os resultados podem ser difíceis de interpretar sem equilíbrio: cuidado com o fato de que valores ausentes na (s) resposta (s) provavelmente perderão o equilíbrio.

Eu acho que você também pode ver isso nas correlações estranhas da matriz do modelo:

mf <- model.matrix(mpg ~ cyl, data = mtcars)
cor(mf)
            (Intercept)       cyl6       cyl8
(Intercept)           1         NA         NA
cyl6                 NA  1.0000000 -0.4666667
cyl8                 NA -0.4666667  1.0000000
Warning message:
In cor(mf) : the standard deviation is zero

Portanto, os erros padrão obtidos de aov(ou lm) provavelmente serão falsos (você pode verificar isso se comparar lmeou com lmererros padrão.

Henrik
fonte
Como você aplicaria o lme aqui?
SRJ
As correlações dos valores da matriz do modelo não são estranhas. Como a constante (interceptação) é inerentemente igual a um, não há variação entre seus valores. Por esse motivo, não é possível calcular um coeficiente de correlação entre uma variável e a constante.
Sven Hohenstein
-1
Y = matrix(0,5,6)
Y[1,] = c(1250, 980, 1800, 2040, 1000, 1180)
Y[2,] = c(1700, 3080,1700,2820,5760,3480)
Y[3,] = c(2050,3560,2800,1600,4200,2650)
Y[4,] = c(4690,4370,4800,9070,3770,5250)
Y[5,] = c(7150,3480,5010,4810,8740,7260)

n = ncol(Y)
R = rowMeans(Y)
M = mean(R)

s = mean(apply(Y,1,var))

v = var(R)  -s/n


#z = n/(n+(E(s2)/var(m)))
Q = 6/(6+(s/v))
t = Q*R[1] + (1-Z)*M
user257426
fonte
Isso é ilegível e carece de qualquer tipo de comentário sobre o que significa ou faz. Você pode editá-lo para dar mais clareza?
Mdewey 27/08/19
Não é possível entender a resposta. Requer algum comentário para explicar o que você está fazendo.
Michael R. Chernick