MLGExemplo01
Fernando L Pérez
24/05/2021
Dados
Os dados a seguir, na pacote “faraway”, fornecem a proporção da área foliar afetada pela mancha foliar em 10 variedades de cevada em 9 locais diferentes.
O quadro de dados com 90 observações contém as 3 variáveis a seguir:
blotch (mancha): proporção da folha de cevada afetada pela mancha
site (local): a localização física - um fator com níveis 1 2 3 4 5 6 7 8 9
variety (variedade): variedade de cevada - um fator com níveis 1 2 3 4 5 6 7 8 9 10
Fonte:
R.W M. Wedderburn (1974) “Quasilikelihood functions, generalized linear models and the Gauss-Newton method” Biometrika, 61, 439-447.
data(leafblotch, package = "faraway")
head(leafblotch)## blotch site variety
## 1 0.0005 1 1
## 2 0.0000 1 2
## 3 0.0000 1 3
## 4 0.0010 1 4
## 5 0.0025 1 5
## 6 0.0005 1 6Estudo descritivo
library(lattice)
xyplot(blotch ~ site | variety, data = leafblotch)
O gráfico acima mostra o comportamento da resposta, a proporção da área foliar afetada pela mancha foliar (blotch), segundo as covariáveis, localização (site) e variedade (variety), ambas nesta situação fatores. Isto é uma situação na qual o modelo de regressão se reduz à uma ANOVA, neste caso, um modelo ANOVA generalizado.
Percebe-se que, para as variedades 1,2,3 e 4, somente a proporção de área foliar afetada pela mancha foliar acontece em mais do que 20% da área na localização 9. Para as variedades 5 e 6, a área foliar afetada é maior do que 20% nas localizações 8 e 9. Para as restantes variedades (7,8,8,e 10), conforme aumenta o código da localização, aumenta a proporção da área foliar afetada pela mancha foliar.
Fica claro que existe uma interação entre os fatores.
Qual é a distribuição adequada à variável resposta?
min(leafblotch$blotch); max(leafblotch$blotch)## [1] 0## [1] 0.95Temos um problema! a função de densidade Beta aceita como suporte o intervalo (0,1), aberto em ambos extremos. Significa que nesta situação não é possível utilizar esta função de densidade.
Na libraria ou pacote de funções “gamlss” temos disponível a função de densidade BEZI, a qual define a distribuição beta inflacionada, uma distribuição de três parâmetros, para um objeto gamlss.family a ser usado no ajuste GAMLSS usando a função gamlss().
O beta inflado de zero é semelhante à distribuição beta, mas permite zeros como valores y. Esta distribuição é uma extensão da distribuição beta usando uma parametrização da lei beta que é indexada por parâmetros de média e precisão (Ferrari e Cribari-Neto, 2004). O parâmetro extra modela a probabilidade em zero.
Primeiro modelo
library(gamlss)## Warning: package 'gamlss' was built under R version 4.0.4## Loading required package: splines## Loading required package: gamlss.data## Warning: package 'gamlss.data' was built under R version 4.0.4##
## Attaching package: 'gamlss.data'## The following object is masked from 'package:datasets':
##
## sleep## Loading required package: gamlss.dist## Warning: package 'gamlss.dist' was built under R version 4.0.4## Loading required package: MASS## Loading required package: nlme## Loading required package: parallel## ********** GAMLSS Version 5.3-2 **********## For more on GAMLSS look at https://www.gamlss.com/## Type gamlssNews() to see new features/changes/bug fixes.ajuste01 = gamlss(blotch ~ variety + site, data = leafblotch, family = BEZI,
control = gamlss.control(n.cyc = 2000, trace = FALSE))
summary(ajuste01)## ******************************************************************
## Family: c("BEZI", "Zero Inflated Beta")
##
## Call: gamlss(formula = blotch ~ variety + site, family = BEZI,
## data = leafblotch, control = gamlss.control(n.cyc = 2000,
## trace = FALSE))
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: logit
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.2917 0.5404 -9.791 9.53e-15 ***
## variety2 -0.1132 0.4587 -0.247 0.805865
## variety3 0.2483 0.4477 0.555 0.580869
## variety4 0.5402 0.4379 1.234 0.221479
## variety5 0.9828 0.4433 2.217 0.029880 *
## variety6 1.3876 0.4096 3.387 0.001162 **
## variety7 1.9232 0.4451 4.320 5.05e-05 ***
## variety8 2.3219 0.4190 5.541 4.95e-07 ***
## variety9 2.4809 0.4204 5.901 1.17e-07 ***
## variety10 3.1716 0.4330 7.324 3.19e-10 ***
## site2 0.6127 0.4610 1.329 0.188179
## site3 1.8413 0.4472 4.117 0.000103 ***
## site4 1.4020 0.4414 3.176 0.002221 **
## site5 1.9631 0.4331 4.533 2.34e-05 ***
## site6 2.0796 0.4274 4.865 6.80e-06 ***
## site7 2.7007 0.4303 6.276 2.53e-08 ***
## site8 3.4524 0.4343 7.949 2.26e-11 ***
## site9 4.6902 0.4483 10.463 5.89e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3380 0.1776 13.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: logit
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.0681 0.5115 -5.998 7.89e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 90
## Degrees of Freedom for the fit: 20
## Residual Deg. of Freedom: 70
## at cycle: 11
##
## Global Deviance: -242.3657
## AIC: -202.3657
## SBC: -152.3695
## ******************************************************************Inluímos o parâmetro gamlss.control(…) por dois motivos: n.cyc para considerar mais do que 20 iterações do algoritmo de estimação, o valor padrão, enquanto que trace = “FALSE” para que não mostra-se informações das diversas iterações.
Resultados do primeiro modelo
Observemos alguns detalhes: mostramos a matriz de planejamento do modelo (algumas linhas), depois mostramos as primeiras linhas do preditor linear, ou seja, o resultado da multiplicação entre a matriz de planejamento e os coeficientes da regressão. Por último mostramos os primeiros valores estimados da resposta, utilizando para isso a inverss da função de ligação. Esses serão os valores apresentados como ajustados no gráfico final.
model.matrix(ajuste01)[1:10,]## (Intercept) variety2 variety3 variety4 variety5 variety6 variety7 variety8
## 1 1 0 0 0 0 0 0 0
## 2 1 1 0 0 0 0 0 0
## 3 1 0 1 0 0 0 0 0
## 4 1 0 0 1 0 0 0 0
## 5 1 0 0 0 1 0 0 0
## 6 1 0 0 0 0 1 0 0
## 7 1 0 0 0 0 0 1 0
## 8 1 0 0 0 0 0 0 1
## 9 1 0 0 0 0 0 0 0
## 10 1 0 0 0 0 0 0 0
## variety9 variety10 site2 site3 site4 site5 site6 site7 site8 site9
## 1 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0 0
## 9 1 0 0 0 0 0 0 0 0 0
## 10 0 1 0 0 0 0 0 0 0 0model.matrix(ajuste01)[1:10,]%*%coef(ajuste01)## [,1]
## 1 -5.291673
## 2 -5.404830
## 3 -5.043357
## 4 -4.751428
## 5 -4.308889
## 6 -3.904115
## 7 -3.368514
## 8 -2.969793
## 9 -2.810741
## 10 -2.120077exp(model.matrix(ajuste01)[1:10,]%*%coef(ajuste01))/(1+exp(model.matrix(ajuste01)[1:10,]%*%coef(ajuste01)))## [,1]
## 1 0.005008126
## 2 0.004474705
## 3 0.006410692
## 4 0.008565351
## 5 0.013270019
## 6 0.019760433
## 7 0.033294094
## 8 0.048809344
## 9 0.056746492
## 10 0.107160736resposta = c(rep("Observado", length(leafblotch$blotch)), rep("Predito", length(leafblotch$blotch)))
resultados = data.frame(blotch = c(leafblotch$blotch,fitted(ajuste01)), site = rep(leafblotch$site,2), variety = rep(leafblotch$variety,2), resposta = resposta)
xyplot(blotch ~ site | variety, data = resultados, groups = resposta,
auto.key=list(space="top", columns=2, title="Resposta", cex.title=1))
Não foi considerada a interação entre as variedades e as localizações porque implicaria num modelo perfeito; isso devido ao número de níveis que teria o novo fator.