Ferramentas do usuário

Ferramentas do site


Barra lateral

lcf5876:historico-disciplina:2018:programa:07-modelos-lineares:07-05-linear-generalizado
  LCF-5876 Computação no Ambiente R:
Aplicações em Ecologia
e Recursos Florestais
R logo
7. MODELOS LINEARES

7.5. Modelos Lineares Generalizados

7.5.1. Estrutura dos Modelos Lineares Generalizados


Os “Modelos Lineares Generalizados” (GLM) são generalizados no sentido de que a estrutura linear da relação causal é generalizada para outras famílias estocásticas.

A estrutura geral dos modelos lineares clássicos pode ser vista na seguinte forma:

$$y = \mu(x) + \varepsilon$$ sendo
  • $y$ a variável resposta,
  • $\varepsilon_{ij}$ a variação aleatória que segue a família Gaussiana com média nula (zero) e variância constatne ($\sigma^2$), e
  • $\mu(x)$ a média que varia como uma função linear de variáveis explicativas ou preditoras ($x$).

No caso de variáveis explicativas ou preditoras quantitativas a função da média é:

$$\mu(x) = \beta_0 +\beta_1 x_{1} +\beta_2 x_{2} + \cdots +\beta_p x_{p}$$ que tem a forma da regressão linear.

No caso de variáveis explicativas que são fatores a função da média tem a forma da análise de variância:

\begin{eqnarray*} \textrm{Exp. Interiamente Casualizados:}&\qquad& \mu(x) = \mu + X \\ \textrm{Exp. em Blocos Casualizados:} &\qquad& \mu(x) = \mu + X_1 + X_2 \\ \textrm{Exp. com Dois Fatores:} &\qquad& \mu(x) = \mu + X_1 + X_2 + (X_1 X_2) \\ & & \\ \end{eqnarray*}

No caso dos “Modelos Lineares Generalizados” (GLM), quem têm relação lineares com as variáveis explicativas ou preditoras é uma transformação da média:

\begin{eqnarray*} \textrm{Função de Ligação:} &\qquad& \eta = l(\mu) \\ \textrm{Regressão Linear:} &\qquad& \eta(x) = \beta_0 +\beta_1 x_{1} +\beta_2 x_{2} + \cdots +\beta_p x_{p} \\ \textrm{Exp. Interiamente Casualizados:}&\qquad& \eta(x) = \mu + X \\ \textrm{Exp. em Blocos Casualizados:} &\qquad& \eta(x) = \mu + X_1 + X_2 \\ \textrm{Exp. com Dois Fatores:} &\qquad& \eta(x) = \mu + X_1 + X_2 + (X_1 X_2) \\ \end{eqnarray*}

A função que de transformação da média é chamada de função de ligação (link function) e para cada família estocástica ela tem uma forma determinada. A transformação também implica numa forma determinada para a variância das observações.

Funções de Ligação para os Modelos Lineares Generalizados

Família Função de Ligação Nome da Função Forma da Variância Nome
Binomial $$\log(\mu/(1-\mu))$$ logit $$\mu(1-\mu)$$ mu(1-mu)
$$\Phi^{-1}(\mu)$$ probit $$\mu(1-\mu)$$ mu(1-mu)
Gama $$-1/\mu$$ inverse $$\mu^2$$ mu^2
Gaussiana $$\mu$$ identity $$1$$ constant
Inverse Gaussian $$-2/mu^2$$ 1/mu^2 $$\mu^3$$ mu^3
Poisson $$\log \mu$$ log $$\mu$$ mu
  1 ##############################################################################################
  2 ###########  7. MODELOS LINEARES
  3 ###########  7.5. Modelos Lineares Generalizados
  4 ##############################################################################################
  5 # 7.5.1. Estrutura dos Modelos Lineares Generalizados

7.5.2. A Função "glm"


A função glm (generalized linear model) é a função utilizada no R para se ajustar os modelos lineares generalidados. seus argumentos são:

# 7.5.1. Estrutura dos Modelos Lineares Generalizados
  6 # 7.5.2. A Função "glm"
  7 args(glm)
  8


> args(glm)
function (formula, family = gaussian, data, weights, subset,
    na.action, start = NULL, etastart, mustart, offset, control = list(...),
    model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL,
    ...)
NULL
>
  • formula é a fórmula do modelo linear sendo ajustado, que segue o mesmo padrão das fórmulas utilizadas na função lm;
  • family é a família estocástica do modelo, a valor default é a família Gaussiana;
  • data é a tabela de dados que contem as variáveis da fórmula;
  • weights são os pesos caso se deseje um ajuste ponderado do modelo, de modo análogo à regressão ponderada;
  • subset argumento para selecionar uma subamostra dos dados da tabela de dados (data);
  • os demais argumentos se referem a aspectos particulares do procedimento de ajuste do modelo.

O ajuste dos GLM é realizado pelo Método da Máxima Verossimilhança, que consiste em encontrar as estimativas dos parâmetros que maximizam a função de verossimilhança do modelo. No caso dos modelos lineares esse método tem uma solução algébrica explícita. Mas, no caso dos GLM, não há uma solução explícita e a maximização é obtida por métodos de otimização. Assim, a função glm traz uma série de argumentos para se controlar esse processo de otimização.

7.5.3. Modelos Binomiais


Os modelos binomiais segue a Família Estocástica Binomial que modela variáveis binárias, isto é, na forma Sim/Não, Sucesso/Fracasso, Morta/Viva. No modelo binomial tradicional a variável resposta é uma variável binária e o modelo busca explicar o seu comportamento a partir de variáveis explicativas quantitativas ou qualitativas (fatores).

A função de ligação default dos GLM da família Binomial é a função logit:

$$\eta = \log\left(\frac{\mu}{1-\mu}\right)\quad\Rightarrow\quad \mu = \frac{\exp(\eta)}{1+\exp(\eta)}$$

Como a variável resposta é binária, a sua média ($\mu$) é na verdade a probabilidade de sucesso, cuja transformação ($\eta$) varia linearmente segundo as variáveis explicativas.

EXEMPLO: MORTALIDADE DE ÁRVORES

O arquivo parag-con.csv o DAP de árvores de uma floresta nativa medidos no início de um período de medição e a informação se árvore estava morta (1) ou viva (0) no final desse período.

  9 # 7.5.3. Modelos Binomiais
 10 ## EXEMPLO: MORTALIDADE DE ÁRVORES
 11 mortarv <- read.csv("parag-con.csv",as.is=TRUE)
 12 head(mortarv)
 13

Surge a questão se a probabilidade das árvores morrerem durante o período de estudo pode ser explicada pelo DAP das árvores. Esse modelo é um modelo binomial tendo o DAP como variável explicativa.

 14 ### Logit
 15 mortmod <- glm( morte ~ dap, family=binomial, data=mortarv)
 16 class(mortmod)
 17 mortmod
 18

Note que o objeto resultamente é da classe glm e que não foi necessário informar a função de ligação, pois a função logit é a função default para família Binomial.

Como nos modelos lineares (lm), a avaliação do modelo ajustado é realizada utilizando-se a função summary para o caso de interesse em modelo de regressão (variáveis explicativas quantitativas - Regressão Logística) e a função anova para o caso de variáveis de interesse qualitativas (fatores).

 19 summary(mortmod)
 20 anova(mortmod, test="F")
 21 anova(mortmod, test="Chisq")
 22

Note que o teste $F$ não é apropriado para os GLM, sendo substituído pelo teste Qui-quadrado.

Os modelos binomiais também podem ser como função de ligação a função probit. Nesse caso, é necessário informar a função glm que além da família Binomial se deseja a função de ligação.

 23 ### Probit
 24 mortmod.p <- glm( morte ~ dap, family=binomial(link=probit), data=mortarv)
 25 mortmod.p
 26 summary(mortmod.p)
 27 anova(mortmod.p)
 28

As funções de ligação logit e probit são duas abordagens possíveis nos modelos binomiais, que geralmente geram resultados semelhantes. A escolha entre elas, depende de questões teóricas e práticas. Mas, é possível fazer uma comparação empírica entre os dois modelos utilizando o Critério de Informação de Akaike (AIC).

 29 ### Comparação
 30 AIC(mortmod, mortmod.p)
 31

Nesse caso particular, ambas abordagem geram resultados igualmente bons pelo AIC, mas os modelos não são iguais, como mostra a comparação gráfica.

 32 plot(morte ~ dap, data=mortarv,
 33      xlab="DAP (cm)", ylab="Probabilidade de Morte",
 34      cex.lab=1.6, cex.axis=1.4)
 35 abline(h=0:1, lty=2, col="red")
 36 curve( exp(-1.62-0.026*x) / (1 + exp(-1.62-0.026*x)), 10, 220, add=TRUE, lwd=3, col="red")
 37 curve( exp(-0.98-0.013*x) / (1 + exp(-0.98-0.013*x)), 10, 220, add=TRUE, lwd=3, col="blue")
 38 legend(50, 0.8, c("Logit","Probit"), lwd=2, col=c("red","blue"), cex=1.6)
 39


EXEMPLO: REGENERAÇÃO DO PALMITEIRO - AJUSTE DE PROPORÇÕES

Em algumas situações, em que a variável respostas não tem o formato binário, surge a questão de se modelar uma resposta que é em si uma probabilidade de sucesso. O arquivo regen-estudo.csv apresenta os dados de levantamento da regeneração natural do palmiteiro (Euterpe edulis) na região do Vale do Ribeira, Estado de São Paulo.

 40 ## EXEMPLO: REGENERAÇÃO DO PALMITEIRO
 41 regen <- read.csv("regen-estudo.csv",as.is=TRUE)
 42 head(regen)
 43

As variáveis regen1 e regen2 se referem ao número de plântulas (até a altura de 50cm) e de mudas (altura entre 50cm e 1,5m) encontradas nas parcelas, respectivamente. Já a coluna g se refere à área basal da floresta nos mesmos locais das parcelas de regeneração. Podemos nos perguntar se o aumento da área basal, como indicação de floresta mais desenvolvida, pode resultar numa proporção maior de plântulas dentre as regenerantes (plântulas + mudas).

A análise gráfica sugere que esse efeito pode estar acontecendo de fato.

 
 44 regen$regen1.p <- regen$regen1 / (regen$regen1 + regen$regen2)
 45 plot(regen1.p ~ g, data=regen,
 46      xlab = "Área Basal (m2/ha)", ylab="Proporção de Plântulas",
 47      cex.lab=1.6, cex.axis=1.4)
 48 abline(h=0:1, lty=2, col="red")
 49 lines(loess.smooth(regen$g, regen$regen1.p), col="green", lwd=2)
 50

Então, podemos ajustar um modelo binomial, pois a variável regen1.p é a proporção de plântulas no conjunto de regenerantes.

 51 reg.g.logit <- glm( regen1.p ~ g, family=binomial, data=regen)
 52

A função glm ajusta o modelo sem problemas, mas apresenta o aviso de que um número não inteiro de sucessos foi obtido. Aparentemente o modelo é bom e indica o aumento de plântulas com o aumento da área basal da floresta.

 53 anova(reg.g.logit, test="Chisq")
 54

Apesar do bom resultado, o aviso da função glm está indicando que existe uma forma mais adequada de modelar as proporções de plântuals em função da área basal da floresta (g). A forma adequada para proporções ou probabilidades empíricas é fornecer como variável resposta duas colunas com as contagens pareadas de sucessos e fracassos:

cbind(sucessos,fracassos) ~ x

sendo x a variável explicativa.

No caso dos dados de regeneração do palmiteiro, o modelo fica:

 55 reg.g.logit2 <- glm( cbind(regen1,regen2) ~ g, family=binomial, data=regen)
 56 anova(reg.g.logit2, test="Chisq")
 57

Os dois modelos ajustados não podem ser comparados utilizando o AIC porque as variáveis respostas estão em escalas diferentes. Mas, é possível fazer uma comparação gráfica.

 58 plot(regen1.p ~ g, data=regen,
 59      xlab = "Área Basal (m2/ha)", ylab="Proporção de Plântulas",
 60      cex.lab=1.6, cex.axis=1.4)
 61 abline(h=0:1, lty=2, col="red")
 62 lines(loess.smooth(regen$g, regen$regen1.p), col="green", lwd=3)
 63 cf1 <- coef(reg.g.logit)
 64 cf2 <- coef(reg.g.logit2)
 65 curve( exp(cf1[1]+cf1[2]*x) / (1 + exp(cf1[1]+cf1[2]*x)), 0, 45, add=TRUE, lwd=3, col="red")
 66 curve( exp(cf2[1]+cf2[2]*x) / (1 + exp(cf2[1]+cf2[2]*x)), 0, 45, add=TRUE, lwd=3, col="blue")
 67 legend(10, 0.4,
 68        c("Linha de Tendência","Ajuste da Proporção","Ajuste das Contagens Pareadas"),
 69        lwd=3, col=c("green","red","blue"), cex=1.6)
 70

Note que a curva do modelo ajustado tendo as contagens pareadas de sucessos e fracassos é claramente mais próxima da curva de tendência.

7.5.4. Modelos Poisson


Os modelos Poisson seguem a famíliea estocástica Poisson, portanto, eles são apropriados para variáveis respostas que sejam resultantes de contagem. O único parâmetro da família Poisson é a média ($\mu$), mas a relação linear com as variáveis explicativas ocorrem com o logaritmo da média ($\log(\mu)$), assim a função de ligação default (link function) é a função log.

\begin{eqnarray*} \eta &=& \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p \\ \eta = \log(\mu) &\Rightarrow& \mu = \exp(\eta) = \exp\left(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p\right) \end{eqnarray*}

Como consequência da estrutura da família Poisson, nesses modelos a variância é constante e igual à média ($\mu$).


EXEMPLO: REGENERAÇÃO DO PALMITEIRO - NÚMERO DE ARVORETAS

Entre as quatro fases de regeneração do palmiteiro registradas no arquivo regen-estudo.csv, a variável regen4 se refere às arvoretas, isto é, plantas com mais de 2m de altura mas com menos de 5cm de DAP.

Uma questão que pode ser levantada é se a ocorrência das arvoretas do palmiteiro é influência pelo grau de conservação da vegetação (vegetacao) e pela topografia (topo).

 71 # 7.5.4. Modelos Poisson
 72
 73 table(regen$topo, regen$vegetacao)
 74 by(regen$regen4, list(vegetacao=regen$vegetacao, topo=regen$topo), mean)
 75

O número médio de arvoretas é bastante diferente, mas podemos construir um modelo Poisson para testar essas hipóteses.

 76 arvoreta.pois <- glm(regen4 ~ vegetacao + topo, data=regen, family=poisson)
 77 summary(arvoreta.pois)
 78 anova(arvoreta.pois, test="Chisq")
 79

O resultado pode ser explorado graficamente.

 80 par(mar=c(5,12,4,2))
 81 boxplot(regen4 ~ factor(vegetacao) + factor(topo), data=regen,
 82          horizontal=TRUE, las=1, xlab="Número de Arvoretas de Palmiteiro",
 83          cex.axis=1.4, cex.lab=1.6)
 84 par(mar=c(5,4,4,2))
 85

Constata-se uma única observação é bastante distinta de todas as demais. Essa observação estaria influenciando o resultado do teste de hipótese?

 86 arvoreta.pois.2 <- glm(regen4 ~ vegetacao + topo, data=regen, subset=regen4 < 10, family=poisson)
 87 anova(arvoreta.pois.2 , test="Chisq")
 88
lcf5876/historico-disciplina/2018/programa/07-modelos-lineares/07-05-linear-generalizado.txt · Última modificação: 2020/03/06 09:58 por joaoluis