Como determinar pesos para a regressão WLS em R?

8

Estou tentando prever a idade em função de um conjunto de marcadores de metilação do DNA. Esses preditores são contínuos entre 0 e 100. Ao executar a regressão do OLS, posso ver que a variação aumenta com a idade.

Assim, decidi ajustar um modelo de regressão ponderada. No entanto, estou tendo problemas para decidir como definir os pesos para o meu modelo. Eu usei o método fGLS, assim:

OLSressq <- OLSres^2                 # Square residuals
lnOLSressq <- log(OLSressq)          # Take natural log of squared residuals
aux <- lm(lnOLSressq~X)              # Run auxillary model
ghat <- fitted(aux)                  # Predict g^
hhat <- exp(ghat)                    # Create h^
fGLS <- lm(Y~X, weights = 1/hhat)    # Weight is 1/h^

E estes foram os meus resultados:

Call:
lm(formula = Y ~ X, weights = 1/hhat)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-4.9288 -1.2491 -0.1325  1.2626  5.1452 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 23.1009494  5.2299867   4.417 1.64e-05 ***
XASPA       -0.1441404  0.0474738  -3.036  0.00271 ** 
XPDE4C       0.6421385  0.0812891   7.899 1.83e-13 ***
XELOVL2     -0.2040382  0.0866564  -2.355  0.01951 *  
XELOVL2sq    0.0088532  0.0009381   9.438  < 2e-16 ***
XEDARADD    -0.1965472  0.0348989  -5.632 5.98e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 1.762 on 200 degrees of freedom
Multiple R-squared:  0.9687,    Adjusted R-squared:  0.9679 
F-statistic:  1239 on 5 and 200 DF,  p-value: < 2.2e-16

No entanto, antes de descobrir como executar o método fGLS, eu estava brincando com pesos diferentes apenas para ver o que aconteceria. Usei 1 / (resíduos quadrados do modelo OLS) como pesos e acabei com isso:

Call:
lm(formula = Y ~ X, weights = 1/OLSressq)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.0893 -0.9916 -0.7855  0.9998  2.0238 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 30.8756737  1.1355861   27.19   <2e-16 ***
XASPA       -0.1956188  0.0116329  -16.82   <2e-16 ***
XPDE4C       0.6168490  0.0102149   60.39   <2e-16 ***
XELOVL2     -0.1596969  0.0116723  -13.68   <2e-16 ***
XELOVL2sq    0.0078459  0.0001593   49.26   <2e-16 ***
XEDARADD    -0.2492048  0.0068751  -36.25   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 1 on 200 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.133e+06 on 5 and 200 DF,  p-value: < 2.2e-16

Como o erro padrão residual é menor, R² é igual a 1 (isso é possível?) E a estatística F é muito mais alta, sou tentado a supor que esse modelo é melhor do que o obtido pelo método fGLS. No entanto, parece-me que escolher pesos aleatoriamente por tentativa e erro deve sempre produzir resultados piores do que quando você realmente tenta matematicamente estimar os pesos corretos.

Alguém pode me dar alguns conselhos sobre quais pesos usar no meu modelo? Também li aqui e ali que você não pode interpretar R² da mesma maneira que faria ao executar a regressão OLS. Mas então como deve ser interpretado e ainda posso usá-lo para comparar meu modelo WLS com o meu modelo OLS?

I. Smeers
fonte
2
Eu seria muito cauteloso com isso R-squared = 1. Você tem uma idéia de quais devem ser os pesos potenciais? Parece que 1/(squared residuals of OLS model)foi apenas um palpite. Se você estiver no escuro sobre os pesos, sugiro o uso de GLS ou Mínimos Quadrados Iterativos.
Jon
Devo acrescentar que, ao ajustar o mesmo modelo a um conjunto de treinamento (metade dos meus dados originais), o R-quadrado caiu de 1 para 0,9983. Foi realmente apenas um palpite, e é por isso que finalmente usei o fGLS, conforme descrito acima. É isso que você quer dizer com "sugiro usar o GLS"? Estou apenas confuso sobre o motivo pelo qual parece que o modelo que fiz ao adivinhar os pesos é um ajuste melhor do que o que fiz ao estimar os pesos através do fGLS. Ainda não ouvi falar dos mínimos quadrados ponderados iterativos, mas analisarei a questão. Obrigado.
I. Smeers
@ Jon, o GLS viável exige que você especifique os pesos (enquanto o GLS inviável que usa pesos teoricamente ótimos não é um estimador viável, ou seja, não pode ser usado na prática).
Richard Hardy
Sim esta correto. No entanto, eles poderiam especificar a estrutura de correlação na nlme::glsfunção. Ele nlme::corClassesfornece uma lista de diferentes estruturas de correlação.
19416 Jon

Respostas:

0

Por que você está usando FLGS? Você tem heterocedasticidade e correlação entre os resíduos? E a matriz var-cov da matriz é desconhecida? Tente bptest(your_model)e se o valor de p for menor que o alfa (por exemplo, 0,05), haverá heterocedasticidade. E então você deve tentar entender se há correlação entre os resíduos com um teste de Durbin Watson: dwtest(your_model)se a estatística W estiver entre 1 e 3, não haverá correlação. Portanto, se você tiver apenas heterocedasticidade, deverá usar o WLS, assim:

mod_lin <- lm(Price~Weight+HP+Disp., data=df)
wts     <- 1/fitted( lm(abs(residuals(mod_lin))~fitted(mod_lin)) )^2
mod2    <- lm(Price~Weight+HP+Disp., data=df, weights=wts)

O mesmo mod2acontece com o modelo antigo, agora com o WLS.

R-square = 1, é muito estranho. Talvez haja colinearidade.

Lorenzo Famiglini
fonte
11
Por que um teste DW seria apropriado. Penso nisso como usado apenas para correlação automática e não vejo como isso se aplicaria neste caso.
meh
11
Bem-vindo ao xvalidated! Por favor, especifique a partir do qual as funções do pacote bpteste dwtestvêm de como eles não são parte da distribuição R padrão.
precisa saber é o seguinte
Porque você precisa entender qual estimador é o melhor: como wls, fgls, ols etc.
Lorenzo Famiglini 22/03/19