Curso R - UFAC - Modelos Lineares Generalizados

© Ricardo Solar/UFMG - compartilhamento e utilização não-comercial livres. Não reproduzir sem autorização > DOI: http://doi.org/10.5281/zenodo.7392285

A partir de agora, nossos modelos mudam um pouco de acordo com a aula que tivemos em relação à migração dos modelos lieares (lm) para os modelos lineares generalizados (glm). Como vimos, a escolha inicial da distribuição de erros depende da natureza dos dados. Trabalharemos nessa seção com dados de contagens, dados binários e dados de proporção. De acordo com o que vimos na aula teórica, teremos então:

• dados de contagens - distribuição Poisson (pronuncia-se pô-á-sôn) ou Binomial Negativa, a depender das premissas atendidas. Podemos ainda ter que lançar mão de uma aproximação conhecida como quasiPoisson.
• dados binários - são aqueles dados de 0 ou 1. Seguem normalmente a distribuição Binomial (ou Bernoulli). Em certas situações, podemos ter que usar a aproximação chamada de quasiBinomial
• dados de proporções - seguirão provavelmente a distribuição Binomial (quando possuimos os dados originais) ou regressão Beta (quando possuímos as proporções). Em certas situações, podemos ter que usar a aproximação chamada de quasiBinomial

1 A função de ligação e as suas vantagens

A função de ligação em Modelos Lineares Generalizados (GLMs) é um componente crucial que conecta o preditor linear (uma combinação das variáveis explicativas) à média da distribuição da variável dependente. Em termos mais simples, a função de ligação garante que as previsões feitas pelo modelo sejam válidas para a distribuição assumida da variável de resposta.

1.1 Pontos-Chave sobre a Função de Ligação:

1. Propósito: A função de ligação mapeia o preditor linear (que pode assumir qualquer valor real) para o intervalo da média do valor esperado da variável de resposta, o que depende da distribuição de probabilidade específica assumida para essa variável. Por exemplo, na regressão logística (um tipo de GLM), a função de ligação mapeia o preditor linear para uma probabilidade, que deve estar entre 0 e 1.

2. Funções de Ligação Comuns:

   • Função Identidade: O valor esperado da variável de resposta é o mesmo que o preditor linear. É usada na regressão linear simples.
   • Função Logit: Usada na regressão logística. Ela mapeia o preditor linear para o intervalo [0,1], apropriado para desfechos binários.
   • Função Logarítmica: Usada na regressão de Poisson. Ela garante que o valor esperado da variável de resposta seja positivo, mapeando o preditor linear para valores positivos.
   • Função Inversa: Usada em modelos onde a variável de resposta segue uma distribuição Gamma.

1.2 Diferença entre a Função de Ligação e a Transformação da Variável:

1.2.1 Função de Ligação:

• Parte da Estrutura do Modelo: A função de ligação é uma parte intrínseca do GLM e define como a média da variável de resposta se relaciona com o preditor linear.
• Interpretação: Ela permite o uso de técnicas de modelagem linear enquanto acomoda diferentes tipos de variáveis de resposta (por exemplo, binárias, contagens, contínuas).
• Sem Transformação Direta dos Dados: A função de ligação opera na média da variável de resposta, não na variável de resposta em si.

1.2.2 Transformação da Variável:

• Pré-processamento de Dados: Transformar a variável geralmente se refere a modificar diretamente as variáveis de resposta ou explicativas, como aplicar o logaritmo, a raiz quadrada ou o inverso de uma variável antes de ajustar o modelo.
• Propósito: As transformações são frequentemente usadas para estabilizar a variância, tornar os dados mais normalmente distribuídos ou linearizar relações entre as variáveis.
• Efeito nos Dados: A transformação altera diretamente os dados, potencialmente mudando sua distribuição, escala ou relações com outras variáveis.

1.3 Exemplo para Clarificação:

• Função de Ligação Logit na Regressão Logística: O modelo é \(\text{logit}(p) = \log\left(\frac{p}{1-p}\right) = \mathbf{X}\beta\), onde \(p\) é a probabilidade do evento ocorrer. A função de ligação (logit) mapeia o preditor linear \(\mathbf{X}\beta\) para uma probabilidade \(p\).

• Transformação Logarítmica: Se você aplicar uma transformação logarítmica a uma variável, estará aplicando a transformação diretamente à variável, mudando sua distribuição. Por exemplo, transformar logaritmicamente uma variável de resposta em uma regressão linear levaria a um modelo onde \(\log(Y) = \mathbf{X}\beta\).

2 A função de ligação

A função de ligação é um componente do modelo que relaciona o preditor linear à média da variável de resposta, enquanto a transformação de uma variável é um processo separado que modifica diretamente a variável.

O preditor linear é uma parte fundamental dos Modelos Lineares Generalizados (GLMs) e de outros modelos estatísticos. Ele é, basicamente, uma combinação linear das variáveis explicativas (ou preditoras) do modelo, que são usadas para prever a variável de resposta.

2.1 Explicação Simples:

1. Composição: O preditor linear é formado pela soma dos produtos de cada variável explicativa (também chamada de variável independente ou preditora) por seus respectivos coeficientes. Os coeficientes são números que o modelo estima e que indicam a força e a direção da relação entre cada variável explicativa e a variável de resposta.

2. Fórmula:

   • Imagine que você tem duas variáveis explicativas, \(X_1\) e \(X_2\).
   • O preditor linear, geralmente representado por \(\eta\), seria algo como \(\eta = \beta_0 + \beta_1X_1 + \beta_2X_2\).
   • Aqui, \(\beta_0\) é o intercepto (um valor constante), e \(\beta_1\) e \(\beta_2\) são os coeficientes associados a \(X_1\) e \(X_2\).

3. Função no Modelo: O preditor linear \(\eta\) representa o valor que o modelo “calcula” antes de aplicar a função de ligação. Dependendo da função de ligação escolhida, o preditor linear é então transformado para prever a variável de resposta de maneira que seja adequada ao seu tipo (por exemplo, uma probabilidade, uma contagem, etc.).

2.2 Exemplo Simples:

• Suponha que você queira prever o peso de uma pessoa com base na altura (\(X_1\)) e na idade (\(X_2\)).
• Se o modelo estima que para cada centímetro a mais de altura o peso aumenta em 0,5 kg (\(\beta_1 = 0,5\)) e que para cada ano de idade o peso aumenta em 0,2 kg (\(\beta_2 = 0,2\)), o preditor linear seria algo como: \[ \eta = \beta_0 + 0,5 \times \text{altura} + 0,2 \times \text{idade} \]
• O valor de \(\eta\) obtido dessa forma é a base para a previsão do peso (dependendo de como o modelo é configurado).

Portanto, o preditor linear é a combinação linear das variáveis explicativas e seus coeficientes, que serve como ponto de partida para fazer previsões no modelo.

3 Dados de contagens - distribuição de Poisson

3.1 Carregamento dos dados

Utilizaremos os dados da planilha a6.txt, que contém o número de espécies de macroinvertebrados bentônicos. A pergunta de interesse é se a riqueza de espécies varia em função da Biomassa de vegetação e o pH da água? Vejam que nesse caso é interessante testarmos a interação, uma vez que em diferentes valores de pH, podemos ter respostas diferentes da riqueza à biomassa.

dadosbiom <- read.table("a6.txt", h=T, stringsAsFactors = T)

summary(dadosbiom)

##      pH        Biomassa          Especies    
##  alto :30   Min.   :0.05017   Min.   : 2.00  
##  baixo:30   1st Qu.:1.44132   1st Qu.:12.25  
##  medio:30   Median :3.10836   Median :18.50  
##             Mean   :3.55701   Mean   :19.46  
##             3rd Qu.:5.08570   3rd Qu.:25.75  
##             Max.   :9.98177   Max.   :44.00

head(dadosbiom)

##Construção e teste dos modelos Nossos modelos então ficam assim:

\[Especies \sim Biomassa + pH + Biomassa:pH\] No R, podemos agrupar Biomassa + pH + Biomassa:pH em Biomassa*pH, e nesse caso usando o modelo Riqueza de exemplo, teremos:

\[Especies \sim Biomassa * pH\]

Tanto faz como preferirá escrever isso no R, para modelos com somente duas variáveis explicativas, as duas formas rendem o mesmo modelo. A partir de 3 variáveis, colocar asteriscos faz todas as combinações possíveis, fiquem de olho nisso!!

Agora podemos fazer o modelo. A partir de agora, como mencionei antes, precisamos usar a função glm, que além do modelo exige que declaremos qual a família de erros, com family=Distribuição. Note que podemos usar GLM com distribuição Normal, nesse caso, não é necessário declarar a família. Riqueza de espécies é um dado de contagens, só com números inteiros. Usarei então a princípio GLM com distribuição Poisson.

modeloriq <- glm(Especies~Biomassa*pH, data=dadosbiom, family = poisson)

Devo seguir com a crítica do modelo, através da análise de resíduos. Logo de cara devemos avaliar se há algum problema de sobredispersão. A família Poisson, conforme podemos ver no summary abaixo, assume que (Dispersion parameter for poisson family taken to be 1). Checamos isso também pelo pacote DHARMA.

Primeiro pelo summary do modelo

summary(modeloriq)

## 
## Call:
## glm(formula = Especies ~ Biomassa * pH, family = poisson, data = dadosbiom)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.76812    0.06153  61.240  < 2e-16 ***
## Biomassa         -0.10713    0.01249  -8.577  < 2e-16 ***
## pHbaixo          -0.81557    0.10284  -7.931 2.18e-15 ***
## pHmedio          -0.33146    0.09217  -3.596 0.000323 ***
## Biomassa:pHbaixo -0.15503    0.04003  -3.873 0.000108 ***
## Biomassa:pHmedio -0.03189    0.02308  -1.382 0.166954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 452.346  on 89  degrees of freedom
## Residual deviance:  83.201  on 84  degrees of freedom
## AIC: 514.39
## 
## Number of Fisher Scoring iterations: 4

DHARMa::simulateResiduals(modeloriq, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.03108104 0.6917273 0.9261017 0.9267392 0.3367185 0.8101546 0.6763737 0.5107622 0.4737265 0.1611126 0.2725799 0.2811094 0.3545331 0.5446296 0.8351068 0.1365787 0.8225895 1 0.5494728 0.002006511 ...

O modelo está correto, e vemos que o teste mostra não haver sobredispersão. Podemos conferir isso também pelo summary(modeloriq). Ao final basta dividir a Residual deviance pelo seu degrees of freedom. No nosso caso, essa divisão dá 0.99. Veja abaixo no summary:

summary(modeloriq)

## 
## Call:
## glm(formula = Especies ~ Biomassa * pH, family = poisson, data = dadosbiom)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.76812    0.06153  61.240  < 2e-16 ***
## Biomassa         -0.10713    0.01249  -8.577  < 2e-16 ***
## pHbaixo          -0.81557    0.10284  -7.931 2.18e-15 ***
## pHmedio          -0.33146    0.09217  -3.596 0.000323 ***
## Biomassa:pHbaixo -0.15503    0.04003  -3.873 0.000108 ***
## Biomassa:pHmedio -0.03189    0.02308  -1.382 0.166954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 452.346  on 89  degrees of freedom
## Residual deviance:  83.201  on 84  degrees of freedom
## AIC: 514.39
## 
## Number of Fisher Scoring iterations: 4

Agora sim podemos testar a significância do modelo. Para isso, quando temos distribuição Poisson, não utilizamos mais o teste F, e sim o \(\chi^2\), vamos lá!

anova(modeloriq, test="Chi")

Tudo foi significativo, incluindo a interação. Nesse caso, não há muito mais o que simplificar, mas devemos verificar se os três níveis de pH são distintos entre si. Para isso, não podemos usar os testes par a par clássicos, pois a distribuição não é normal. Ah, observe que essa anova informou que esse modelo é Model: poisson, link: log.

pH2 <- dadosbiom$pH
sort(
  tapply(dadosbiom$Especies, pH2, mean)
  )

##    baixo    medio     alto 
## 11.56667 20.00000 26.80000

Os valores de pH médio e alto têm as médias mais próximas, vou tentar juntá-los

levels(pH2)

## [1] "alto"  "baixo" "medio"

levels(pH2)[1]

## [1] "alto"

levels(pH2)[2]

## [1] "baixo"

levels(pH2)[3]

## [1] "medio"

levels(pH2)[1:2]

## [1] "alto"  "baixo"

levels(pH2)[c(1,3)]

## [1] "alto"  "medio"

levels(pH2)[c(1,3)] <- "alto_medio"

Agora posso tentar fazer um novo modelo, em que os dois níveis de pH alto e médio agora são apenas um.

modeloriq2 <- glm(Especies~Biomassa*pH2, data=dadosbiom, family = poisson)

E agora comparo os modelos modeloriq e modeloriq2

anova(modeloriq, modeloriq2, test="Chi")

Como os modelos são diferentes, tenho evidência para manter separados os três níveis de pH.

Meu modelo então válido é modeloriq.

3.2 Vamos agora ao gráfico

Nosso gráfico agora é bem mais complexo. Temos variação no intercepto e na inclinação das curvas. Felizmente nesse caso, acreditem se quiser, o R permite fazer automático e não preciso construir cada curva! Em sala faremos cada uma individualmente para entender melhor, abaixo deixarei apenas os comandos para fazê-las.

library(ggplot2)

ggplot(dadosbiom, aes(y=Especies, x=Biomassa, colour=pH, fill=pH))+
  geom_point()+
  geom_smooth(method = "glm", method.args = list(family = "poisson"), se=F)

## `geom_smooth()` using formula = 'y ~ x'

3.2.1 Equações das retas

Teremos que fazer três equações. Como havia explicado acima, o R fornece os dados para a primeira curva (em ordem alfabética) e depois as diferenças para as demais curvas. Nesse caso, como a interação foi significativa, temos que observar como modificar intercepto e inclinação. Vamos começar a ver os coeficientes:

coef(modeloriq)

##      (Intercept)         Biomassa          pHbaixo          pHmedio 
##       3.76812359      -0.10712979      -0.81557124      -0.33146257 
## Biomassa:pHbaixo Biomassa:pHmedio 
##      -0.15502818      -0.03189202

Vamos lá, teremos as curvas para [em ordem alfabética]: \(Espécies \sim Biomassa\) para pH alto, baixo e médio. Nos coeficientes, os nomes Intercept e * Biomassa* são a curva para pH alto. Os nomes pHbaixo e pHmedio são a diferença no intercepto. Por fim, os nomes Biomassa:pHbaixo e Biomassa:pHmedio são relativos à diferença na inclinação.

Mas aqui temos uma leve pegadinha, prestem atenção! Quando fizemos o anova(modeloriq, test="Chi"), vimos uma coisa importante lá, lembrando: Model: poisson, link: log, ou seja, a função de ligação desse modelo é log, e para plotar preciso “desfazer” esse log, com a função exp. As curvas então ficam assim:

• pH Alto - \[exp(Intercept + Biomassa \times x)\]
• pH Baixo - \[exp((Intercept+pHbaixo) + (Biomassa+Biomassa:pHbaixo) \times x)\]
• pH Medio - \[exp((Intercept+pHmedio) + (Biomassa+Biomassa:pHmedio) \times x)\]

Calculando as equações da reta

##Criar inverso da ligacao

inversa <- modeloriq$family$linkinv

Coef <- coef(modeloriq)

eqalto <- function(x){exp(Coef[1]+Coef[2]*x)}

eqbaixo <- function(x){exp((Coef[1]+Coef[3])+(Coef[2]+Coef[5])*x)}

eqmedio <- function(x){exp((Coef[1]+Coef[4])+(Coef[2]+Coef[6])*x)}

Agora vou refazer o plot

ggplot(dadosbiom, aes(y=Especies, x=Biomassa, colour=pH))+
  geom_point()+
  stat_function(fun=eqalto, col=1)+
  stat_function(fun=eqbaixo, col=2, xlim=c(0,5))+
  stat_function(fun=eqmedio, col=3)+
  scale_color_manual(values = c(1,2,3))+
  theme_bw()+
  theme(legend.position = "bottom")

4 Dados binários - Distribuição binomial

4.1 Carregamento dos dados

Utilizaremos para essa análise a planilha camundongos.txt, onde é avaliado se há fatores que estão relacionados à maior probabilidade de infecção de camundongos.

A pergunta de interesse é se a infecção por um vírus em camundongos é influenciada pelo peso, idade e sexo biológico do animal. Vamos carregar os dados para ver o que temos:

dadosinfec <- read.table("~/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/GLM/Binomial/camundongos.txt", h=T, stringsAsFactors = T)

Vejam que, a título de exemplo, ao invés de trocar o diretório, eu optei por informar o endereço onde econtra-se a tabela. Agora é o momento de explorar um pouco os dados.

summary(dadosinfec)

##      infec            idade             peso          sexo   
##  Min.   :0.0000   Min.   :  1.00   Min.   : 1.00   femea:28  
##  1st Qu.:0.0000   1st Qu.: 26.00   1st Qu.: 9.00   macho:53  
##  Median :0.0000   Median : 87.00   Median :13.00             
##  Mean   :0.2099   Mean   : 83.65   Mean   :11.49             
##  3rd Qu.:0.0000   3rd Qu.:130.00   3rd Qu.:16.00             
##  Max.   :1.0000   Max.   :206.00   Max.   :18.00

Vejam que a nossa variável resposta, infec, é binária contendo apenas 0 (não infectado) e 1 (infectado). Como vimos em aula, utilizaremos distribuição binomial. Quando os dados contêm apenas 0s ou 1s, não é necessário informar ao R nada além da distribuição. Nesse tipo de modelo, a função de ligação é conhecida por \(logit\). A função logit é descrita assim:

\[logit(p)=log \left( \frac{p}{1-p} \right)\] ## Construindo o modelo

Nesse modelo, temos interesse em entender as interações duplas somente. Então, se eu simplesmente colocar asteriscos, ele fará todas as combinações, incluindo a interação tripla. Discutiremos em sala os problemas de utilizar essas interações triplas.

Para fazer as interações duplas somente, podemos fazer de duas formas, a seguir:

1. infec ~ idade + peso + sexo + idade:peso + idade:sexo + peso:sexo
2. um atalho possível é 1. infec ~ (idade + peso + sexo)^2

modeloinfec <- glm(infec~(peso+sexo+idade)^2, data=dadosinfec, family=binomial)

Como já vimos, tenho que criticar o modelo através de análise de resíduos.

DHARMa::simulateResiduals(modeloinfec, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.330355 0.6106647 0.753566 0.1345021 0.3626166 0.516243 0.4584971 0.6298425 0.04356199 0.7068328 0.8307066 0.2172472 0.193813 0.1631083 0.5370363 0.3113622 0.9344821 0.7449671 0.4973379 0.4664253 ...

Os modelos estão ok, e a análise de resíduos mostra que não temos nenhum desvio das premissas. Estando tudo certo, posso seguir pra testar a significância do modelo.

anova(modeloinfec, test="Chi")

Vejam que agora temos um modelo cheio, que contêm tanto variáveis significativas, quanto variáveis não significativas. Pelo princípio da Navalha de Occam ¹, um modelo deve ser tão mais simples quanto as evidências contrárias à simplificação permitam.

Sendo assim, vamos simplificar nosso modelo. Começaremos pelos termos mais “complexos”, no nosso caso, as interações, seguidos poelos termos simples. Quanto à ordem dentro dos termos, podemos adotar que quanto maiores os valores de P, tiramos essas variáveis primeiro. Vamos à obra.

Uma forma de simplificar os modelos é digitá-los novamente, mas há uma melhor maneira no R para isso, que é o comando update. Nesse comando, primeiro informamos o modelo que será atualizado, em seguida informamos qual a atualização. Nesse momento, lançamos mão de uma possibilidade do R, que é usar um ponto para repetir tudo que estava em determinada situação. No nosso caso, queremos atualizar o modelo, mantendo todo ele menos aqueles termos que quero remover. Feita a atualização passo a passo, devemos testar se o novo modelo pode mesmo substituir o mais complexo. Vejam abaixo:

modeloinfec2 <- update(modeloinfec, .~. -peso:sexo)
anova(modeloinfec, modeloinfec2, test="Chi")

anova(modeloinfec2, test="Chi")

Como os dois modelos, um mais completo e outro simplificado, se apresentaram estatisticamente iguais, mantemos o mais simples. Vou fazer isso com as demais variáveis até que eu tenha somente variáveis significativas retidas.

modeloinfec3 <- update(modeloinfec2, .~. -peso:idade)
anova(modeloinfec2, modeloinfec3, test="Chi")

anova(modeloinfec3, test="Chi")

##

modeloinfec4 <- update(modeloinfec3, .~. -idade:sexo)
anova(modeloinfec3, modeloinfec4, test="Chi")

anova(modeloinfec4, test="Chi")

##

Agora temos apenas variáveis significativas, ou seja, temos em nosso modelo as variáveis peso, sexo e idade como significativos.

Como a interação não é significativa, a melhor forma de representar os dados nesse caso é fazer um gráfico individual para cada variável. Vamos discutir isso em sala.

Para nos auxiliar no plot das variáveis categóricas, vamos fazer aquele resumo do summarySE para calcular os valores de média e erro padrão.

library(Rmisc)

## Loading required package: lattice

## Loading required package: plyr

medias <- summarySE(dadosinfec, measurevar = "infec", groupvars = "sexo")

ggplot(dadosinfec, aes(x=sexo, y=infec))+
  geom_jitter(width=.3)+
  geom_errorbar(data=medias, aes(ymax=infec+se, ymin=infec-se), position = position_dodge(.9), width=.3)+
  geom_point(data = medias, aes(x=sexo, y=infec), size=4, colour="red", position = position_dodge(.9))

ggplot(dadosinfec, aes(x=idade, y=infec))+
  geom_point()+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=F)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(dadosinfec, aes(x=peso, y=infec))+
  geom_point()+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=F)

## `geom_smooth()` using formula = 'y ~ x'

Dá para juntar as coisas?

Podemos começar juntando sexo e idade e também sexo e peso

ggplot(dadosinfec, aes(x=peso, y=infec, colour=sexo))+
  geom_point()+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=F)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(dadosinfec, aes(x=idade, y=infec, colour=sexo))+
  geom_point()+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=F)

## `geom_smooth()` using formula = 'y ~ x'

Por fim, poderíamos juntar idade e peso no mesmo gráfico, mas aí vejam que fica bem mais complexo.

ggplot(dadosinfec, aes(x=idade, y=infec))+
  geom_point(aes(size=peso, colour=peso))+
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=F)

## `geom_smooth()` using formula = 'y ~ x'

Vejam que a representação gráfica é sempre um desafio quando temos muitas variáveis. É importante refletir a respeito de como você quer mostrar os seus dados.

5 Dados de proporções - Distribuição Binomial e Regressão Beta

5.1 Carregamento e exploração dos dados

Para essa análise agora, temos duas espécies de rotíferos (Keratella cochlearis e Polyarthra major) que foram expostos a águas com densidades diferenças. A pergunta é se a densidade da água está relacionada a um maior número de indivíduos suspensos na coluna d’água. A questão aqui é que o número absoluto de indivíduos suspensos é uma medida enviesada, pois algumas réplicas possuíam mais indivíduos como um todo (Expostos). Sendo assim eu preciso ponderar essa resposta, calculando a proporção de indivíduos suspensos, ou seja \(Proporção = \frac{Suspensos}{Expostos}\)

Para termos o exemplo registrado, importaremos os dados diretamente do arquivo Excel, no formato xlsx. Para isso, precisamos do pacote library(readxl).

library(readxl)
dadosexp <- read_excel("~/Dropbox/UFMG/Disciplinas/R UFAC/Aulas/GLM/Binomial/rotiferos.xlsx")
summary(dadosexp)

##    Densidade      Suspensos         Expostos        Especie         
##  Min.   :1019   Min.   :  7.00   Min.   : 14.00   Length:40         
##  1st Qu.:1030   1st Qu.: 12.50   1st Qu.: 47.25   Class :character  
##  Median :1044   Median : 22.00   Median : 75.00   Mode  :character  
##  Mean   :1045   Mean   : 47.02   Mean   :109.12                     
##  3rd Qu.:1060   3rd Qu.: 40.50   3rd Qu.:160.25                     
##  Max.   :1070   Max.   :488.00   Max.   :492.00

Vejam que nesse caso, o R interpretou a variável Espécie como um conjunto de textos simples (Class: character). Precisarei explicar pro programa que não é bem assim. Farei isso contando para o programa que a coluna dadosexp$Especie é na verdade um fator, utilizando o comando as.factor.

dadosexp$Especie <- as.factor(dadosexp$Especie)

summary(dadosexp)

##    Densidade      Suspensos         Expostos              Especie  
##  Min.   :1019   Min.   :  7.00   Min.   : 14.00   K.cochlearis:20  
##  1st Qu.:1030   1st Qu.: 12.50   1st Qu.: 47.25   P.major     :20  
##  Median :1044   Median : 22.00   Median : 75.00                    
##  Mean   :1045   Mean   : 47.02   Mean   :109.12                    
##  3rd Qu.:1060   3rd Qu.: 40.50   3rd Qu.:160.25                    
##  Max.   :1070   Max.   :488.00   Max.   :492.00

Agora tudo certo!

5.1.1 Modelo com distribuição binomial

Primeiramente, como conhecemos os valores que darão origem à proporção, posso fazer o modelo com distribuição binomial. Para isso, preciso informar ao modelo que a variável y é \(Suspensos/Expostos\), que a family=binomial e que o total é mesmo Expostos, usando weights=Expostos.

modeloind <- glm(Suspensos/Expostos~Densidade*Especie, family=binomial, data=dadosexp, weights = Expostos)

Seguindo com a análise de resíduos. Avalie cuidadosamente os resultados abaixo.

DHARMa::simulateResiduals(modeloind, plot=TRUE)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.965655 0.09977238 0.4073607 0.1231166 0.004 0.5253924 0.9700647 1 0.007142732 0.9088863 0.460317 0.996029 0.9721198 0.1324739 0.1161441 0 0.9026699 0 1 0.8395575 ...

summary(modeloind)

## 
## Call:
## glm(formula = Suspensos/Expostos ~ Densidade * Especie, family = binomial, 
##     data = dadosexp, weights = Expostos)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -1.144e+02  4.034e+00 -28.345   <2e-16 ***
## Densidade                 1.087e-01  3.857e-03  28.191   <2e-16 ***
## EspecieP.major            4.629e+00  6.598e+00   0.702    0.483    
## Densidade:EspecieP.major -3.077e-03  6.329e-03  -0.486    0.627    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3180.99  on 39  degrees of freedom
## Residual deviance:  434.02  on 36  degrees of freedom
## AIC: 596.57
## 
## Number of Fisher Scoring iterations: 5

summary(modeloind)$deviance/summary(modeloind)$df.residual

## [1] 12.05605

Vejam que o modelo não está nada bem ajustado! Um problema que já detectamos é que há grande sobredispersão. O valor é 12.06, quando deveria ser 1.

Quando isso ocorre, lançamos mão de um ajuste na distribuição Binomial, conhecida como quasiBinomial. Ela ajusta a dispersão, resolvendo o problema.

modeloindquasi <- glm(Suspensos/Expostos~Densidade*Especie, family=quasibinomial, data=dadosexp, weights = Expostos)

Análise de Resíduos …

DHARMa::simulateResiduals(modeloindquasi, plot=T)

## Error in simulate.lm(object, nsim = nsim, ...): family 'quasibinomial' not implemented

car::qqPlot(resid(modeloindquasi))

## [1] 34 18

car::residualPlot(modeloindquasi)

> o quasibinimial corrige calculando qual a dispersão, ou seja, não considera mais sendo como 1.

summary(modeloindquasi)

## 
## Call:
## glm(formula = Suspensos/Expostos ~ Densidade * Especie, family = quasibinomial, 
##     data = dadosexp, weights = Expostos)
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -1.144e+02  1.495e+01  -7.647 4.74e-09 ***
## Densidade                 1.087e-01  1.430e-02   7.606 5.36e-09 ***
## EspecieP.major            4.629e+00  2.446e+01   0.189    0.851    
## Densidade:EspecieP.major -3.077e-03  2.346e-02  -0.131    0.896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 13.73901)
## 
##     Null deviance: 3180.99  on 39  degrees of freedom
## Residual deviance:  434.02  on 36  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

Veja acima, o pacote DHARMa não funciona, então fiz os resíduos de outra forma. É sempre bom saber fazer a análise na unha.

Agora sim posso testar meu modelo!

anova(modeloindquasi, test="Chi")

Como vimos antes, o termo da interação não é signficativo, preciso tentar removê-lo.

modeloind2 <- update(modeloindquasi, .~.-Densidade:Especie)

anova(modeloind, modeloind2, test="Chi")

anova(modeloind2, test="Chi")

Posso removê-lo!

5.1.1.1 Plot dos resultados

A curva de um modelo binomial possui a função de ligação logit. Sua inversão é feita através do comando plogis. Acompanhe o plot abaixo.

C <- coef(modeloind2)
curva1 <- function(x){plogis(C[1]+C[2]*x)}
curva2 <- function(x){plogis((C[1]+C[3])+C[2]*x)}

ggplot(dadosexp, aes(y=(Suspensos/Expostos), x=Densidade, col=Especie))+
  geom_point()+
  stat_function(fun=curva1, col="purple")+
  stat_function(fun=curva2, col="brown")+
  scale_colour_manual(values = c("purple", "brown"))

5.1.2 Modelo de regressão beta

Outra forma de trabalhar com proporções é a regressão beta. Vamos usar os mesmos dados e analisá-los dessa forma. Para fazer a regressão beta, precisamos de pequenos ajustes nos dados.

1. a variável resposta já deve estar no formato de proporções, então irei calcular isso antes
2. as proporções não podem ser numericamente iguais a 0 e 1. Nesse caso, vou usar uma função para subtrair um valor ínfimo se for igual a 1, e somar um valor ínfimo se for igual a 0. Pra fazer essa operação, usarei a função ifelse. Essa função é bastante útil, e sua sintaxe básica é ifelse(condição, valor se for verdade, valor se for mentira)

proporcao <- dadosexp$Suspensos/dadosexp$Expostos

range(proporcao)

## [1] 0.03533569 1.00000000

proporcao <- ifelse(
  proporcao==0, proporcao+1e-3, 
  ifelse(
    proporcao==1, proporcao-1e-3, 
    proporcao))

Agora posso construir o modelo, fazer a análise de resíduo e fazer os testes de significância

library(betareg)
modelobeta <- betareg(proporcao~Densidade*Especie, data = dadosexp)

car::qqPlot(resid(modelobeta))

## [1] 16 40

plot(modelobeta, which=1)

A análise de resíduo precisa ser feita de maneira manual, sem o pacote DHARMa. A avaliação da significância também é diferente, e preciso usar a função Anova (com letra maiúscula mesmo) do pacote car.

car::Anova(modelobeta)

modelobeta2 <- update(modelobeta, .~.-Densidade:Especie)

car::Anova(modelobeta2)

coef(modelobeta2)

##    (Intercept)      Densidade EspecieP.major          (phi) 
##   -94.07518722     0.08959406     1.00493284     6.28461574

O plot seria feito da mesma forma.

C <- coef(modelobeta2)
curva1 <- function(x){plogis(C[1]+C[2]*x)}
curva2 <- function(x){plogis((C[1]+C[3])+C[2]*x)}

ggplot(dadosexp, aes(x=Densidade, y=proporcao, col=Especie))+
  geom_point()+
  stat_function(fun=curva1, col="purple")+
  stat_function(fun=curva2, col="brown")+
  scale_colour_manual(values = c("purple", "brown"))

6 E como fazer um gráfico fora do GGPLOT2 com qualidade?

#Primeiramente determinar as curvas
C <- coef(modelobeta2)
curva1 <- function(x){plogis(C[1]+C[2]*x)}
curva2 <- function(x){plogis((C[1]+C[3])+C[2]*x)}

plot(proporcao~Densidade, col=as.numeric(Especie), pch=15+as.numeric(Especie), data=dadosexp)
curve(curva1, col=1, add=T)
curve(curva2, col=2, add=T)

## Nem perto, precisarei detalhar muito mais os comandos

# primeiro o canvas
plot(proporcao~Densidade, bty="n", type="n", axes=F, xlab="", ylab="", data=dadosexp)

# agora começo a adicionar os elementos
##pontos
points(proporcao~Densidade, col=as.numeric(Especie), pch=15+as.numeric(Especie), cex=2, data=dadosexp)
##eixos
axis(1, lwd=2)
axis(2, lwd=2, las=1)
## caixa em volta do gráfico
box(bty="l", lwd=3)
## nomes dos eixos
mtext("Densidade", 1, 2)
mtext("Proporção de suspensos", 2, 3)
## curvas
curve(curva1, col=1, add=T, lwd=2)
curve(curva2, col=2, add=T, lwd=2)
##legendas
names <- levels(dadosexp$Especie)
legend("bottomright", names, col=1:2, pch=16:17, lty=1, bty="o")

---

1. Occam’s razor, also spelled Ockham’s razor, also called law of economy or law of parsimony, principle stated by the Scholastic philosopher William of Ockham (1285–1347/49) that pluralitas non est ponenda sine necessitate, “plurality should not be posited without necessity.” The principle gives precedence to simplicity: of two competing theories, the simpler explanation of an entity is to be preferred. The principle is also expressed as “Entities are not to be multiplied beyond necessity.” https://www.britannica.com/topic/Occams-razor↩︎