Os modelos lineares generalizados (GLMs) são uma ampliação dos modelos lineares ordinários. Os GLM's são usados quando os resíduos (erro) do modelo apresentam distribuição diferente da normal (gaussiana). A natureza da variável resposta é uma boa indicação do tipo de distribuição de desvios que iremos encontrar nos modelos. Por exemplos, variáveis de contagem são inteiras e apresentam os valores limitados no zero. Esse tipo de variável, em geral, tem uma distribuição de erros assimétrica para valores baixos e uma variância que aumenta com a média dos valores preditos, violando duas premissas dos modelos lineares. Os casos mais comuns de modelos generalizados são de variáveis resposta de contagem, proporção e binária, muito comum nos estudos de ecologia e evolução.

Uma das formas de entendermos os modelos generalizados é separar o modelo em dois componentes: a relação determinística entre as variáveis (resposta e preditora) e o componente aleatório dos resíduos (distribuição dos erros). Em um modelo linear ordinário a relação entre as variáveis é uma proporção constante, o que define uma relação funcional de uma reta. Quando temos uma contagem, essa relação pode ter uma estrutura funcional de uma exponencial. Para esses casos, o glm faz uso de uma função de ligação log, para linearizar a relação determinística entre as variáveis, temos portanto o estrutura deterministica do modelo definida por um preditor linear associado à função de ligação.

O componente aleatório dos resíduos, no caso de uma variável de contagem, segue, em geral, uma distribuição poisson. A distribuição poisson é uma variável aleatória definida por apenas um parâmetro ($\lambda$), equivalente à média da distribuição normal, chamada de lambda. A distribuição poisson tem uma característica interessante, seu desvio padrão é igual à média. Portanto, se a média aumenta, o desvio acompanha esse aumento e a distribuição passa a ter um maior espalhamento.

O preditor linear está associado à estrutura determinística do modelo e está relacionado à linearização da relação, aqui definido como $\eta$:

$$ \eta = \alpha + \beta x$$

A função de ligação é o que relaciona o preditor linear com a esperança do modelo:

$$ \eta = g^{-1}(E_{(y)}) $$

Ou seja, nos modelos generalizados não é a variável resposta que tem uma relação linear com a preditora, e sim o preditor linear que tem uma relação linear com as preditoras.

Para alguns tipos de famílias de variáveis temos funções de ligações padrões. As mais usadas são:

Um exemplo, apresentado no livro do Michael Crawley, The R Book, relata a contagem de espécies de árvores em unidades amostrais de florestas com diferentes biomassa e classificadas em três níveis de ph no solo: baixo, médio e alto. O objetivo desse experimento não manipulativo é verificar se há relação entre riqueza de árvores e as preditora biomassa da floresta e ph do solo.

No menu Graphs, selecione XY conditioningh Plot e selecione as varíáveis, definindo ph como variável de agrupamento, como no gráfico abaixo.

No menu Models>Graphs selecione Predict effect plots… e selecione as variáveis.

Vamos utilizar um exemplo que está presente no livro de W. Venables e B. Ripley, Modern Applied Statistics with S-PLUS²⁾, sobre o número de dias ausentes da escola de crianças na Austrália.

No Rcommader (Rcmdr) vá ao menu Tools > Load package(s) e selecione o pacote MASS.
Caso o pacote não apareça listado, significa que ele já está carregado, então pule esse passo.

Em sequida:

• abra o menu Data > Data in packages > Read data from an attached package…
• selecione o pacote MASS e os dados quine ³⁾

Os dados estão relacionados ao estudo para entender quais variáveis estão relacionados à ausência (falta) do aluno na escola. A observação está relacionada a alunos amostrados aleatoriamente de escolas na Austrália.

Para nosso exercício vamos deixar de lado a variável Lrn por que há dados faltantes nela com relação a outras variáveis. Vamos construir o modelo cheio com todas as outras variáveis (Eth, Sex, Age ) e todas as possibilidades de interações entre elas. Começamos então com um modelo linear simples.

• abra o menu Statistics > Fit model > Linear Model
• construa um modelo cheio com (Age, Eth e Sex) e as suas interações possíveis
• faça a simplificação do modelo para obter o modelo mínimo adequado
• guarde o resultado do modelo selecionado para comparar com o GLM

Já fizemos a simplificação de um modelo de contagem no exemplo anterior de riqueza de espécies. Nesse exemplo a dispersão dos dados foram bem ajustadas pelo modelo. Caso a razão Residual deviance e degrees of freedom seja maior que um, a poisson não conseguiu lidar com a dispersão dos dados. Nesses casos é possível utilizar o modelo quassipoisson que estima mais um parâmetro para lidar com a sobre ou sub dispersão dos dados. Os passos para esse ajuste são descritos abaixo.

Vamos construir o modelo seguindo essa sequência iniciando com a família de erro POISSON e a função de ligação log.

• abra o menu Statistics > Fit model > Generalized Linear Model
• construa um modelo cheio com (Age, Eth e Sex) e as suas interações possíveis
• faça a simplificação do modelo para obter o modelo mínimo adequado

Um dos pressupostos do modelo Poisson é que a variância aumenta linearmente com a esperança (média do modelo). Podemos avaliar isso dividindo a Residual Deviance pelo seu degrees of freedom. Essa razão deve ser próxima a 1. O que não é o caso do nosso modelo. Nesses casos uma das alternativas é:

• ajustar o modelo usando Family: quasipoisson

O gráfico do modelo pode ser obtido no Rcmdr da mesma forma indicada no modelo anterior, no menu: Models>Graphs selecione Predict effect plots… e selecione as variáveis.

Os modelos de proporção de sucessos (sucessos/tentativas), proporção simple (%) ou de resposta binária (presença/ausência, vivo/morto) são modelados, normalmente, com estrutura do erro binomial. Nesses casos os limites dos valores da variável resposta é bem definido: entre 0 e 1. Além disso, a variância não é constante e varia conforme a média. Essas características fazem com que os resíduos apresentem uma estrutura que aumenta e depois diminuí, e normalmente o máximo de desvios é encontrado nos valores intermediários.

A função Bernoulli, que é a base para a binomial é definida pelo parâmetro de probabilidade de sucesso em um evento com duas possibilidades de resultado (binário). O parâmetro da função Bernoulli é a a probabilidade de sucesso. No caso de uma moeda justa seria a probabilidade de 0,50 de sair coroa.

A binomial é uma generalização da Bernoulli, definida pelo número de sucessos em certo número (n) de tentativas (evento Bernoulli).

Conceitos Importantes

• n = número de tentativas
• s = número de sucessos
• f = número de falhas

 Probabilidade de sucesso 

$$p = \frac{s}{n} $$

Probabilidade de falha 

$$q = \frac{f}{n} $$

$$q = 1 - p$$

Chance de sucesso 

$$\text{odds} = \frac{s}{f} $$

$$\text{odds} = \frac{p}{1-p} $$

A estrutura da função de ligação é a mesma para qualquer modelo:

O preditor linear $\eta$ está associado à estrutura determinística do modelo e relacionado à linearização da relação.

$$ \eta = \alpha + \sum\beta_{x_i}$$

A função de ligação é o que relaciona o preditor linear com a esperança do modelo:

$$ \eta = g^{-1}(E_{(y)}) $$

A função de ligação $g()$ canônica ou padrão para modelos com resposta binária ou proporção é chamada de logit ou logaritmo da chance⁵⁾, definida como:

$$ \eta = \log(\frac{p}{1-p})$$

$$ \log(\frac{p}{1-p}) = \alpha + \sum\beta_{x_i} $$

Sendo $\frac{p}{1-p}$ a chance ou odds em inglês.

Para reverter o preditor linear da função logit para a escala de observação usa-se a função inversa:

$$ g^{-1} = \text{logit}^{-1} = \frac{e^{\eta}}{(1+ e^{\eta})} $$

O predito pelo modelo na escala do preditor linear do modelo binário com função de ligação logit está na escala de logaritmo da chance ($\log(\frac{p}{1-p})$). A razão de chance é uma medida muito popular em outras áreas da ciência, como medicina e mede o quanto uma chance é proporcionalmente diferente de outra, geralmente comparando com um nível controle. Ou seja, qual a proporção de mudança na chance do tratamento em relação a chance do controle. Dado que, em variáveis categóricas os coeficientes do modelo são relacionados às diferenças entre o nível do tratamento e o controle:

$$\text{exp}(\log(\text{odds}_{\text{trat}}) - \log(\text{odds}_{\text{control}})) = \frac{\text{odds}_{\text{trat}}}{\text{odds}_{\text{control}}} $$

então, exponenciar os coeficientes do modelo binomial com preditora categórica transforma os coeficientes em razão de chance comparado com o nível basal ⁶⁾.

No caso de variáveis contínuas a razão de chance é relacionada à chance de x+1 comparada com x, ou seja, qual a proporção de mudança na chance com o aumento de uma unidade da variável contínua preditora.

Portanto, uma forma de interpretar os coeficientes do modelo é exponenciar e interpretá-los como razão de chance, sendo o intercepto a chance o nível basal (variável categórica) ou a chance quando a variável contínua é zero.

O conjunto de dados que vamos usar, isolation.txt tem como variável:

O objetivo do estudo que gerou esses dados é saber se a ocorrência da ave está relacionada com o isolamento e tamanho da ilha.

• abra os dados isolation.txt no Rcmdr (a separação de campo é tabulação)
• monte o modelo cheio com todas a variáveis preditoras e interações
• simplifique o modelo para o mínimo adequado

O modelo prevê a ocorrência da ave na escala de logaritmo da chance (log odds-ratio). Para os coeficientes estimados pelo modelo o melhor é usar o a função exp e interpretar a razão de chance entre categorias ou entre x+1 e x. Para interpretar os valores previsto é necessário aplicar a função inversa do logit, ou seja, nosso modelo faz previsões na escala de log(odds-ratio), nosso preditor linear $\eta$, e precisamos retornar para a escala de observação que é a probabilidade de ocorrência ($\hat{y}$):

$$\hat{y} = \frac{e^{\hat{\eta}}}{1+e^{\hat{\eta}}} $$

Mais um exemplo apresentado no livro do Michael Crawley, The R Book. Neste experimento o objetivo foi avaliar a floração de 5 variedades de plantas tratadas com hormônios de crescimento (6 concentrações). Depois de seis semanas as plantas foram classificadas em floridas ou vegetativas.

O objetivo do estudo que gerou esses dados é saber se o evento de floração é influenciado pelo dose de hormônio e a variedade da planta.

• baixe o arquivo flowering.txt
• abra os dados no Rcmdr (a separação de campo é tabulação) com o nome flower
• crie a variável prop pelo menu Data> Manage variables in active data set> Compute new variable…, colocando no campo Expression to compute:

cbind(sucess = flowered, fail = number - flowered)

Esse comando acima cria uma nova variável nos dados flower chamada prop. Essa nova variável tem duas colunas (sucess e fail) contendo o número de plantas floridas e o número de plantas que não floresceram, respectivamente.

• use a variável prop como resposta (sucessos, falhas)
• monte o modelo cheio com todas a variáveis preditoras e interações
• simplifique o modelo para o mínimo adequado

Para interpretar os coeficientes use o mesmo procedimento do exercício anterior, que é aplicar a função exponencial (exp) nos coeficientes previstos e interpretar como chance e razão de chance⁷⁾.

Para interpretar os valores previsto é necessário aplicar a função inversa do logit, ou seja, nosso modelo faz previsões na escala de log(odds-ratio), nosso preditor linear $\eta$, e precisamos retornar para a escala de observação que é a probabilidade de florescer ($\hat{y}$):

$$\hat{y} = \frac{e^{\hat{\eta}}}{1+e^{\hat{\eta}}} $$

(preditoLinear <- predict("nomedomodelo"))
(preditoProp <- exp(preditoLinear)/(1+ exp(preditoLinear)))

predito <- predict("nomedomodelo", type = "response")
predito

Os modelo GLM poisson e binomial apresentam a variância acoplada à média dos valores, diferentemente dos modelos com distribuição normal onde a média e a variância são independentes. Caso haja uma variação maior ou menor nos dados do que o previsto por essas distribuições, o modelo não consegue dar conta. Essa sobre-dispersão ou sub-dispersão dos dados indica que temos mais ou menos variação do que é predito pelos modelos. Isso pode ser decorrência de vários fontes de erro na definição do modelo, alguns exemplos são:

• o resíduo dos dados pode não ter sido gerado por um processo aleatório poisson ou binomial
• há mais variação do que predito pela ausência de preditoras importantes
• muitos zeros, além do predito pelas distribuições, em decorrência de diferentes processos: um que gera a ausência e outro que gera a variação nas ocorrências de sucesso

• Modelos Lineares Generalizados: binomial
• Modelos Lineares Generalizados: contagem