Introdução

O presente trabalho irá abordar o ajuste de um Modelo linear generalizado(MLG), aplicado ao banco de dados disponivel no proprio software Rstudio, de nome “CO2”, para estes daddos seram aplicados duas distribuições para modelagem, sendo elas a Gama e a Normal Inversa.

1. Definição de um MLG

Os modelos lineares generalizados são ajustados usando a função glm ( ). A forma da função glm é dada por:

glm(Formula, family = Tipo da Distribuição,(link =“Função de link”), data =Banco de Dados).

1.1 Tabela das Distribuições do MLG

Familia Função de Link Padrão
Binomial (link = “logit”)
Gaussiana (link = “identy”)
Gamma (link = “inverse”,“log”,“Identy”)
Inversa (link = “1/mu^2”)
Poisson (link = “log”)
Inverse.gaussian (link = “inverse”,“log”,“Identy”)

2. Quais distribuições posso escolher?

Gaussiano: uma distribuição gaussiana (normal).

Binomial: uma distribuição binomial para proporções.

Poisson: uma distribuição de Poisson para contas.

Gama: uma distribuição gama para dados contínuos positivos.

Inverse.gaussian: uma distribuição Gaussiana inversa para dados contínuos positivos.

3. Pacotes Necessários Para Análise

library(corrplot) 
library(ggplot2)
library(hnp)
library(DescTools)
library(fBasics)
library(nortest)

3.1 Função de cada pacote

corrplot: Tem como sua função exibir gráficamente uma matriz de correlação, intervalo de confiança.

ggplot2: É um sistema para criação de gráficos “declarativamente”, baseado em “The Grammar of Graphics”. Onde você fornece os dados, informa ao ‘ggplot2’ como mapear as variáveis melhorando a estética, sabendo quais primitivas gráficas usar, e cuida atentamente dos detalhes.

hnp: Tem como função produzir um gráfico (meio) normal de um objeto de modelo ajustado para uma variedade de modelos diferentes. Extensível para classes de modelo não implementadas.

DescTools: È uma função para gerar a estimativa do pseudo R2.

fBasics: É uma função para análise descritiva dos dados.

nortest: É uma função para o teste de Normalidade de Anderson-Darling.

4. Banco de Dados

O banco de dados contém 84 obsevações e 5 variáveis,de um experimento sobre a tolerância ao frio das espécies de gramíneas Echinochloa crus-galli,(Tipo de grama), como não ha problemas nos dados (NA´s), mas há algumas variáveis que não apresntam relevancia para entrarem no ajuste do modelo, como previamente já fiz a análise com essas variáveis e não foram siginficativas, irei retirar as variáveis Pant e Treatment para seguirmos com a análise.

data("CO2")
head(CO2)
##   Plant   Type  Treatment conc uptake
## 1   Qn1 Quebec nonchilled   95   16.0
## 2   Qn1 Quebec nonchilled  175   30.4
## 3   Qn1 Quebec nonchilled  250   34.8
## 4   Qn1 Quebec nonchilled  350   37.2
## 5   Qn1 Quebec nonchilled  500   35.3
## 6   Qn1 Quebec nonchilled  675   39.2

4.1 Retirando As Variáveis Plant e Treatment, Gerando um Novo Banco de Dados

CO2<-CO2[,-c(1,3)]

4.2 Definição Das Variáveis

** Variáveis Explicativas**

Type: um fator com níveis Quebec e Mississippi dando a origem da planta.

Conc: um vetor numérico de concentrações de dióxido de carbono no ambiente (mL / L).

** Variável Resposta**

Uptake: um vetor numérico de taxas de absorção de dióxido de carbono (umol / m ^ 2 seg).

4.3 Visualizando as primeiras 6 observações

head(CO2)
##     Type conc uptake
## 1 Quebec   95   16.0
## 2 Quebec  175   30.4
## 3 Quebec  250   34.8
## 4 Quebec  350   37.2
## 5 Quebec  500   35.3
## 6 Quebec  675   39.2

4.4 Verificando a Dimensão dos Dados

dim(CO2)
## [1] 84  3

4.5 Identificando as variáveis Pelo Nome

attach(CO2)

5. Atribuindo as Variável (Type) como Fator

CO2$Type <- as.factor(CO2$Type)

6. Análise Descritiva Dos Dados Qualitativos

Tabela 6.1 Descritiva da variável Tipo (Type)

# table(Type)
Quebec Mississippi
42 42

Comentário: Como esses dados são de um planejamento experimental, contém a mesma quantidade para cada tipo de Grama.

7. Análise Descritiva Dos Dados Quantitativos

 # round(fBasics::basicStats(CO2[4:5]),3)

** Definição das Variáveis**.

Conc: Concentração de CO2 na Plata.

uptake: Absorção de CO2 pela Planta.

conc uptake
n 84 84
Mín 95 7.7
Máx 1000 45.5
1.Quart 175 17.9
3.Quart 675 37.1
Mean 435 27.2
Median 350 28.3
Se Mean 32.2 1.2

Comentário : Com uma amostra de tamanho 84, obteve-se um valor mínimo tanto para concentração de CO2 95, quanto para absorção 7.7, variando até um valor máximo para concentração e absorção 1000 e 45.5 respectivamente. A média de concentração variou em torno de 435 e absorção variou entre 27.2.

8. Histograma da variável uptake (Absorção)

ggplot(data=CO2) + geom_histogram(mapping=aes(uptake),fill="black",color="white")+theme_bw(8)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Comentário : Com base no histograma, há indicios de que a variavél resposta não segue normalidade.

9. Teste de Normalidade da variável resposta uptake (Absorção)

shapiro.test(uptake)
## 
##  Shapiro-Wilk normality test
## 
## data:  uptake
## W = 0.94105, p-value = 0.0007908
ad.test(uptake)
## 
##  Anderson-Darling normality test
## 
## data:  uptake
## A = 1.5793, p-value = 0.0004271

Comentário : Com base nos testes de normalidade, como p P-valor < 0.05, então rejeita a hipótese nula de que a variavel tem distribuição normal. Então caimos na hipótese alternativa que a variavel resposta não segue normalidade, ou seja, seguiremos a análise com um Modelo Linear Generalizado (MLG).

10. Boxplot da variável uptake (Absorção)

qplot(conc,uptake, data=CO2, geom="boxplot",color=conc)+theme_bw()

Comentário : Com base no boxplot, observa-se que a variabilidade dos dados é alta e está em torno da mediana, portanto não há prensça de outliers nesta variável conc.

12. Correlação entre as Variáveis conc e uptake

# cor((CO2[4:5]))
Correlações Conc Uptake
Conc 1 0.49
Uptake 0.49 1

Comentário: as variáveis conc e uptake tem uma correlação que varia de moderada pra fraca.

13. Modelagem Linear Generalizada

Para Ser Realizado a modelagem de um MLG, a variável resposta deve pertencer a familia exponencial, desta forma, como a característica dos dados é continua e assimétrica, seram testadas as distribuições Gamma e Normal Inversa para modelagem.

14. Modelagem com a Distribuição Gamma

Para a modelagem de um ajuste de MLG com a distribuição Gama vou utilizar as funções de ligação “inverse”,“log”,“identity”, para o ajuste do modelo.

14.1 Modelo 1

glmod1 <- glm(uptake ~ Type + conc,data = CO2,family = Gamma(link="inverse"))

14.1.1 Resumo do Modelo 1

summary(glmod1)
## 
## Call:
## glm(formula = uptake ~ Type + conc, family = Gamma(link = "inverse"), 
##     data = CO2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.88842  -0.23414  -0.00797   0.20570   0.51168  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.983e-02  2.487e-03  16.014  < 2e-16 ***
## TypeMississippi  1.753e-02  2.635e-03   6.652 3.14e-09 ***
## conc            -1.996e-05  3.592e-06  -5.558 3.39e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.09354562)
## 
##     Null deviance: 15.8858  on 83  degrees of freedom
## Residual deviance:  8.5349  on 81  degrees of freedom
## AIC: 594.86
## 
## Number of Fisher Scoring iterations: 5

14.2 Modelo 2

glmod2 <- glm(uptake ~ Type + conc, data = CO2, family = Gamma(link="log"))

14.2.2 Resumo do Modelo 2

summary(glmod2)
## 
## Call:
## glm(formula = uptake ~ Type + conc, family = Gamma(link = "log"), 
##     data = CO2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.87631  -0.24400  -0.00133   0.19066   0.57684  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.1893615  0.0669744  47.621  < 2e-16 ***
## TypeMississippi -0.4681381  0.0654629  -7.151 3.44e-10 ***
## conc             0.0006938  0.0001113   6.235 1.93e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.08999319)
## 
##     Null deviance: 15.8858  on 83  degrees of freedom
## Residual deviance:  8.0573  on 81  degrees of freedom
## AIC: 589.94
## 
## Number of Fisher Scoring iterations: 5

14.3 Modelo 3

glmod3 <- glm(uptake ~ Type + conc, data = CO2, family = Gamma(link="identity"))

14.3.3 Resumo do Modelo 3

summary(glmod3)
## 
## Call:
## glm(formula = uptake ~ Type + conc, family = Gamma(link = "identity"), 
##     data = CO2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.86748  -0.22733   0.01168   0.18922   0.60932  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      23.558844   1.776150  13.264  < 2e-16 ***
## TypeMississippi -11.375372   1.667753  -6.821 1.49e-09 ***
## conc              0.021657   0.003156   6.863 1.24e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.08552522)
## 
##     Null deviance: 15.8858  on 83  degrees of freedom
## Residual deviance:  7.7388  on 81  degrees of freedom
## AIC: 586.5
## 
## Number of Fisher Scoring iterations: 8

15. Comparação Entre os Modelos Gamma

#ajuste = c('glmod1', 'glmod2','glmod3')
#aic    = c(AIC(glmod1),AIC(glmod2),AIC(glmod3))
#deviance    = c(deviance(glmod1),deviance(glmod2),deviance(glmod3))
#verossimilhanca = c(logLik(glmod1),logLik(glmod2),logLik(glmod3))
#data.frame(ajuste, aic,verossimilhanca,deviance)
Ajuste AIC Deviance
glmod 1 594.86 8.53
glmod 2 589.94 8.05
glmod 3 586.50 7.73

Comentário: por meio dos ajustes do modelo gama,com os tipos de funções de ligação, o melhor modelo que obdece os pré-requisitos de menor AIC e menor deviance, foi o modelo (glmod3) com a função de ligação (link = “identity”).

15.1: pseudoR2 para modelagem da Gama.

# PseudoR2(glmod1)
# PseudoR2(glmod2)
# PseudoR2(glmod3)
Modelo McFadem
glmod1 0.083
glmod2 0.091
glmod3 0.096

Comentário: Com base na estimativa do PseudoR2, no qual,aquele modelo que apresentar maior valor dessa estimativa é candidato a ser escolhido como melhor modelo, neste caso, apenas confirmou o que ja tinha sido constatado atravéz do AIC e dad Deviance, que atendeu todas essas exigencias foi o glmod3.

16. Gráfico (Envelope simulado) Normalmente Distribuido do Modelo Escolhido.

hnp(glmod3$residuals, sim = 99,resid.type ='deviance',how.many.out=T ,
    conf = 0.95,scale = T)
## Half-normal plot with simulated envelope generated assuming the residuals are 
##         normally distributed under the null hypothesis. 
## Estimated mean: -0.07878169 
## Estimated variance: 50.64003

## Total points: 84 
## Points out of envelope: 0 ( 0 %)

Comentário: Com base no gráfico acima, observamos que a modelagem gama se ajustou bem para modelar este tipo de dados continuos, e apresentou um bom ajuste, englobando todos os pontos dentro do envelope . Mas podemos testar uma outra opção de distribuição para modelar um novo ajuste, neste caso a distribuição Normal Inversa, e comparar os resultados.

17. Modelagem com a Distribuição Normal Inversa

Para a modelagem de um ajuste de MLG com a distribuição Normal Inversa vou utilizar as funções de ligação “inverse”,“log”,“identity”, para o ajuste do modelo.

17.1 Modelo 1

 NImod1 <- glm(uptake ~ Type + conc, data = CO2, family= inverse.gaussian(link="inverse"))

17.1.1 Resumo do Modelo 1

 summary(NImod1)
## 
## Call:
## glm(formula = uptake ~ Type + conc, family = inverse.gaussian(link = "inverse"), 
##     data = CO2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.20785  -0.04707  -0.00052   0.04138   0.11538  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      4.164e-02  2.816e-03  14.790  < 2e-16 ***
## TypeMississippi  1.788e-02  2.806e-03   6.371 1.07e-08 ***
## conc            -2.392e-05  4.288e-06  -5.579 3.11e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 0.004270211)
## 
##     Null deviance: 0.71256  on 83  degrees of freedom
## Residual deviance: 0.42464  on 81  degrees of freedom
## AIC: 610.95
## 
## Number of Fisher Scoring iterations: 2

17.2 Modelo 2

 NImod2 <- glm(uptake ~ Type + conc , data = CO2, family= inverse.gaussian(link="log"))

17.2.2 Resumo do Modelo 2}

 summary(NImod2)
## 
## Call:
## glm(formula = uptake ~ Type + conc, family = inverse.gaussian(link = "log"), 
##     data = CO2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.20450  -0.04820  -0.00014   0.03773   0.13073  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.1277424  0.0744169  42.030  < 2e-16 ***
## TypeMississippi -0.4601046  0.0717683  -6.411 9.01e-09 ***
## conc             0.0008368  0.0001321   6.336 1.25e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 0.004154242)
## 
##     Null deviance: 0.71256  on 83  degrees of freedom
## Residual deviance: 0.39938  on 81  degrees of freedom
## AIC: 605.8
## 
## Number of Fisher Scoring iterations: 8

17.3 Modelo 3

 NImod3 <- glm(uptake ~ Type + conc, data = CO2, family= inverse.gaussian(link="identity"))

17.3.3 Resumo do Modelo 3

 summary(NImod3)
## 
## Call:
## glm(formula = uptake ~ Type + conc, family = inverse.gaussian(link = "identity"), 
##     data = CO2)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.197528  -0.053560  -0.005741   0.036152   0.134413  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      20.683190   1.851119  11.173  < 2e-16 ***
## TypeMississippi -10.253313   1.722774  -5.952 6.50e-08 ***
## conc              0.028483   0.003769   7.558 5.55e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 0.003797027)
## 
##     Null deviance: 0.71256  on 83  degrees of freedom
## Residual deviance: 0.37157  on 81  degrees of freedom
## AIC: 599.74
## 
## Number of Fisher Scoring iterations: 14

18. Comparação Entre os Modelos Normal Inversa

# ajuste = c('NImod1', 'NImod2','gNImod3')
# aic    = c(AIC(NImod1), AIC(NImod2),AIC(NImod3))
# deviance    = c(deviance(NImod1), 
# deviance(NImod2),deviance(NImod3))
# verossimilhanca = 
# c(logLik(NImod1),logLik(NImod2),logLik(NImod3))
# data.frame(ajuste, aic, verossimilhanca,deviance)
Ajuste AIC Deviance
NImod1 610.94 0.42
NImod2 605.79 0.39
NImod3 599.73 0.37

Comentário: Por meio dos ajustes do modelo Normal Inverssa,com os tipos de funções de ligação, o melhor modelo que obdece os pré-requisitos de menor AIC e menor deviance, foi o modelo (glmod3) com funções de ligação (link = “identity”).

18.1 Pseudo R2 para modelagem da Normal Inversa.

# PseudoR2(NImod1)
# PseudoR2(NImod2)
# PseudoR2(NImod3)
Modelo McFadem
NImod1 0.067
NImod2 0.075
NImod3 0.084

Comentário: Com base na estimativa do PseudoR2, no qual,aquele modelo que apresentar maior valor dessa estimativa é candidato a ser escolhido como melhor modelo, neste caso, apenas confirmou o que ja tinha sido constatado atravéz do AIC e dad Deviance, que atendeu todas essas exigencias foi o NImod3.

19. Envelope Simulado Gerado Assumindo que os Resíduos são normalmente distribuído sob a hipótese nula.

hnp(NImod3$residuals, sim = 99,resid.type ='deviance',how.many.out=T ,
    conf = 0.95,scale = T)
## Half-normal plot with simulated envelope generated assuming the residuals are 
##         normally distributed under the null hypothesis. 
## Estimated mean: -0.7334257 
## Estimated variance: 60.46173

## Total points: 84 
## Points out of envelope: 0 ( 0 %)

Comentário: Com base no gráfico acima, observamos que a modelagem Normal Inversa serve para modelar dados continuos, neste caso proporcionou o englobamento de todos os pontos dentro do envelope como se espera de um bom ajuste.

20. Escolha da distribuição

Para escolher qual a melhor distribuição para o ajuste deste banco de dados será montada uma tabela com os valores de AIC, Deviance e Pseudo R2, respectivamente. Como tanto a Gama quanto a Normal Inversa englobaram todos os pontos dentro do envelope, o método de seleção da distribuição será basiado no menor valor de AIC e Deviance, e maior Pseudo R2.

Tabela 20.1 : Tabela com os valores de AIC, Deviance e Pseudo R2, para as distribuições Gama e Normal Inversa.

Distribuição Modelo AIC Deviance Mcfadem
Gama glmod3 586.50 7.73 0.096
Normal Inversa NImod3 599.73 0.37 0.084

Comentário: Como podemos observar nos valores da tabela 20.1, o modelo (glmod3) da distribuição gama, apresentou menor AIC e maior valor na estatística de McFadem, porém apresentou maior Deviance comparando com os valores obtidos na Normal Inversa. Desta forma, a distribuição gama ganha em dois aspectos importantes para selecionar o melhor modelo, assim para este tipo de dados a distribuição que mehor se ajusta aos dados de um experimento sobre a tolerância ao frio das espécies de gramíneas Echinochloa crus-galli,(Tipo de grama), será adistribuição Gama.

21. Interpretação dos Coeficientes do melhor modelo da distribuição Gama.

# glmod3$coefficients
Intercept TypeMississippi conc
23.55884375 -11.37537230 0.02165683

Comentário: Quando é plantado a espécie graminea do tipo Mississippi ela tem uma redução na absorção de CO2 11 vezes menor em relação a do tipo Quebec, como também a concentração de CO2 encontrada na espécie de grama Mississippi é 2% maior comparada a do tipo Quebec.

Conclusão

Como vimos no decorrer desta análise, a Modelagem Linear Generalizada (MLG), proporcionou uma boa estimativa no ajuste dos modelos tanto para distribuição gama quanto para Normal Inversa, por meio de critérios de seleção, a gama foi a que proporcionou melhor ajuste para este banco de dados com a função de Ligação (link = “identity”).