Código
rm(list=ls(all=T))Regressão Beta
Exercício de Regressão Beta
Carregando o banco de dados, para visualizar os dados:
Os dados utilizados neste exemplo foram extraídos de Smithson e Verkuilen (2006). Os dados foram obtidos de Pammer e Kevin (2004). Trata-se de um estudo sobre a habilidade em leitura de um grupo de 44 crianças realizado na escola de psicologia da Universidade Nacional da Austrália. Smithson e Verkuilen (2006) consideram a contribuicão relativa da sensitividade visual (presença de dislexia) e do QI não-verbal na habilidade de leitura das 44 criança. A variável resposta é o escore em um teste de habilidade em leitura (y) e as variáveis independentes são a ocorrência de dislexia (x1) e o escore padronizado de QI não-verbal (x2), chamadas daqui em diante de dislexia e QI, respectivamente. A covariável x1 assume os valores 1 e -1 indicando se a criança é ou não dislexica, respectivamente.
Código
setwd("C:\\Users\\Pessoal\\Downloads")
dislexia <- read.table("dislexia.txt", header = T)
dislexia$dislexia <- ifelse(dislexia$dislexia==1,1,0)
names(dislexia) <- c("identi", "escore", "dislexia", "QI")Observando, o banco de dados:
Código
DT::datatable(dislexia, caption = "Banco de dados")O banco de dados que temos acima, é o que vamos trabalhar.
Estatística Descritiva
Código
plotly::ggplotly()Podemos observar que a variável está distribuida entre 0 e 1.
Código
skimr::skim(dislexia)| Name | dislexia |
| Number of rows | 44 |
| Number of columns | 4 |
| _______________________ | |
| Column type frequency: | |
| numeric | 4 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| identi | 0 | 1 | 22.50 | 12.85 | 1.00 | 11.75 | 22.50 | 33.25 | 44.00 | ▇▇▇▇▇ |
| escore | 0 | 1 | 0.77 | 0.18 | 0.46 | 0.62 | 0.71 | 0.99 | 0.99 | ▂▇▃▁▇ |
| dislexia | 0 | 1 | 0.43 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▆ |
| QI | 0 | 1 | 0.00 | 1.00 | -1.75 | -0.80 | -0.12 | 0.62 | 1.86 | ▆▇▇▇▆ |
Ajustando o modelo de precisão fixa
- Vamos ajustar o modelo de precisão fixa, os autores proporam que existe significativa entre a interação das outras duas variáveis.
Código
library(betareg)Warning: package 'betareg' was built under R version 4.3.1Código
ajuste_logit_pf <- betareg(escore~QI+dislexia+QI*dislexia, data = dislexia)
summary(ajuste_logit_pf)
Call:
betareg(formula = escore ~ QI + dislexia + QI * dislexia, data = dislexia)
Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7445 -0.5157 0.0617 0.5057 2.3407
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.3074 0.2083 11.079 < 2e-16 ***
QI 0.3793 0.2136 1.776 0.0758 .
dislexia -1.9473 0.2669 -7.295 2.98e-13 ***
QI:dislexia -0.4371 0.2691 -1.624 0.1043
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 11.133 2.444 4.556 5.21e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 51.35 on 5 Df
Pseudo R-squared: 0.6068
Number of iterations: 19 (BFGS) + 2 (Fisher scoring) Comentário: De acordo com o modelo ajustado, podemos observar que o pseudo \(R^2\), foi de 60.68%, um valor que não é tão alto. Podemos observar também que temos todos os coeficientes são significativos, ao nível de significância de 10%, exceto a interação entre as duas variáveis independentes, podemos observar também a estimativa do \(\phi\).
Ajustando o modelo de precisão variável
Código
ajuste_logit_pv <- betareg(escore~QI+dislexia+QI*dislexia| QI, link = "logit", data = dislexia)
summary(ajuste_logit_pv)
Call:
betareg(formula = escore ~ QI + dislexia + QI * dislexia | QI, data = dislexia,
link = "logit")
Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7211 -0.5445 0.0585 0.5261 2.4078
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.2919 0.2076 11.042 < 2e-16 ***
QI 0.1553 0.2479 0.626 0.531
dislexia -1.9251 0.2723 -7.071 1.54e-12 ***
QI:dislexia -0.2156 0.2911 -0.741 0.459
Phi coefficients (precision model with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.3887 0.2223 10.746 <2e-16 ***
QI -0.3765 0.2281 -1.651 0.0988 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 51.87 on 6 Df
Pseudo R-squared: 0.5812
Number of iterations: 14 (BFGS) + 7 (Fisher scoring) Comentário: Podemos observar que o modelo de Regressão com variável não foi muito bom, pois podemos observar que o \(R^2\) não foi melhor que o do modelo de regressão com precisão fixa, pois podemos observar que o valor foi de 58.12%. Podemos observar que a estimativa de \(\phi\), foi de -0.37, e não foi significativo.
Podemos observar que os coeficientes foram não significativos, exceto a dislexia que foi significativo.
Observando qual o melhor dos modelos, o do precisão fixa e o de precisão variável
- Realizando o teste da razão de verossimilhança
Código
lmtest::lrtest(ajuste_logit_pf, ajuste_logit_pv)Likelihood ratio test
Model 1: escore ~ QI + dislexia + QI * dislexia
Model 2: escore ~ QI + dislexia + QI * dislexia | QI
#Df LogLik Df Chisq Pr(>Chisq)
1 5 51.350
2 6 51.868 1 1.035 0.309Resposta: Não rejeitamos \(H_0\), ao nível de significância de 5% não rejeitamos \(H_0\), ou seja, o modelo de precisão fixa é o melhor modelo a ser utilizado.
Então, vamos utilizar o modelo de precisão fixa.
- Vamos utilizar o teste
resetpara verificar se o modelo está bem específicado.
Código
lmtest::lrtest(ajuste_logit_pf, .~.+I(predict(ajuste_logit_pf,type="link")^2)) Likelihood ratio test
Model 1: escore ~ QI + dislexia + QI * dislexia
Model 2: escore ~ QI + dislexia + I(predict(ajuste_logit_pf, type = "link")^2) +
QI:dislexia
#Df LogLik Df Chisq Pr(>Chisq)
1 5 51.350
2 6 52.023 1 1.3449 0.2462Resposta: Não rejeitamos \(H_0\), ao nível de significância de 5%, ou seja, ao nível de significância de 5% o modelo está bem especificado. Temos o melhor modelo, com as variáveis.
- Testando as outras funções de ligação
Código
ajuste_probit_pf <- betareg(escore~QI+dislexia+QI*dislexia,data = dislexia, link = "probit")
summary(ajuste_probit_pf)
Call:
betareg(formula = escore ~ QI + dislexia + QI * dislexia, data = dislexia,
link = "probit")
Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7460 -0.5121 0.0618 0.5042 2.3271
Coefficients (mean model with probit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3332 0.1059 12.585 < 2e-16 ***
QI 0.1915 0.1046 1.831 0.0671 .
dislexia -1.1080 0.1487 -7.454 9.08e-14 ***
QI:dislexia -0.2274 0.1458 -1.560 0.1188
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 11.16 2.45 4.556 5.22e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 51.4 on 5 Df
Pseudo R-squared: 0.6311
Number of iterations: 21 (BFGS) + 2 (Fisher scoring) Resposta: Podemos observar que o \(R^2\), foi de 63.11%, um valor maior do que o do modelo com a função logit, observamos também que obtivemos novas estimativas para os coeficientes.
- Vamos testar novas funções:
Código
ajuste_loglog_pf <- betareg(escore~QI+dislexia+QI*dislexia,data = dislexia, link = "loglog")
summary(ajuste_loglog_pf)
Call:
betareg(formula = escore ~ QI + dislexia + QI * dislexia, data = dislexia,
link = "loglog")
Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7441 -0.5164 0.0618 0.5058 2.3447
Coefficients (mean model with loglog link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.3581 0.1981 11.906 < 2e-16 ***
QI 0.3614 0.2049 1.764 0.0777 .
dislexia -1.7218 0.2361 -7.292 3.05e-13 ***
QI:dislexia -0.4067 0.2412 -1.686 0.0918 .
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 11.128 2.442 4.556 5.21e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 51.34 on 5 Df
Pseudo R-squared: 0.5956
Number of iterations: 21 (BFGS) + 2 (Fisher scoring) Código
ajuste_cloglog_pf <- betareg(escore~QI+dislexia+QI*dislexia,data = dislexia, link = "cloglog")
summary(ajuste_probit_pf)
Call:
betareg(formula = escore ~ QI + dislexia + QI * dislexia, data = dislexia,
link = "probit")
Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7460 -0.5121 0.0618 0.5042 2.3271
Coefficients (mean model with probit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3332 0.1059 12.585 < 2e-16 ***
QI 0.1915 0.1046 1.831 0.0671 .
dislexia -1.1080 0.1487 -7.454 9.08e-14 ***
QI:dislexia -0.2274 0.1458 -1.560 0.1188
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 11.16 2.45 4.556 5.22e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 51.4 on 5 Df
Pseudo R-squared: 0.6311
Number of iterations: 21 (BFGS) + 2 (Fisher scoring) Código
ajuste_cauchit_pf <- betareg(escore~QI+dislexia+QI*dislexia,data = dislexia, link = "cauchit")
summary(ajuste_cauchit_pf)
Call:
betareg(formula = escore ~ QI + dislexia + QI * dislexia, data = dislexia,
link = "cauchit")
Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7354 -0.5157 0.0613 0.5049 2.3892
Coefficients (mean model with cauchit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.6003 0.6920 5.203 1.97e-07 ***
QI 1.2000 0.8324 1.442 0.149
dislexia -3.3134 0.7035 -4.710 2.48e-06 ***
QI:dislexia -1.2482 0.8437 -1.480 0.139
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 10.955 2.403 4.559 5.14e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 51.05 on 5 Df
Pseudo R-squared: 0.4642
Number of iterations: 35 (BFGS) + 2 (Fisher scoring) Podemos observar que um dos melhores modelos é o loglog, o cloglog. Vamos observar pelo AIC:
Código
AIC(ajuste_cauchit_pf, ajuste_cloglog_pf, ajuste_logit_pf, ajuste_loglog_pf, ajuste_probit_pf) df AIC
ajuste_cauchit_pf 5 -92.10616
ajuste_cloglog_pf 5 -92.88542
ajuste_logit_pf 5 -92.70079
ajuste_loglog_pf 5 -92.67536
ajuste_probit_pf 5 -92.79519Análise de Resíduos utilizando o modelo logit
- Testaremos a análise de resíduos, com os resíduos de Pearson e depois com os resíduos ponderados
Código
par(mfrow = c(2,3))
plot(ajuste_logit_pf, which = 1:6, type = "pearson",sub.caption = "")
Resposta: Podemos observar que alguns pontos ficaram fora do gráfico do envelope, vamos testar outro tipo de resíduo.
Código
par(mfrow = c(2,3))
plot(ajuste_logit_pf, which = 1:6, type = "sweighted2",sub.caption = "")
Resposta: Podemos observar que com o tipo de resíduo padronizado, foi melhor nesse sentido, temos apenas um ponto fora do gráfico de envelope.