A função de verossimilhança serve para encontrar os parâmetros de uma distribuição de probabilidade que melhor descreve os dados. Para isso, muitas vezes precisamos permitir que os parâmetros das distribuições variem, em função uma ou mais variáveis preditoras ou independentes.

Nestes tutoriais vamos simular contagens com médias que dependem de algumas outras variáveis que também foram medidas. Em seguida vamos construir passo a passo modelos estatísticos que expressam estas contagens como variáveis Poisson e binomial, cujos valores esperados são funções de variáveis independentes. Estes são os chamados modelos lineares generalizados (glm), para os quais já há rotinas de ajuste prontas para uso no R, e outros pacotes estatísticos.

• Modelos não-gaussianos
• Modelos lineares generalizados
• Função de ligação
• Função logística, probito e logito
• Sobredispersão

Neste tutorial vamos usar modelos Poisson para descrever a distribuição de número de frutos em plantas adultas de Solanum sp.. Mas ao invés de dados reais, criaremos os valores de uma distribuição Poisson. A vantagem da simulação é que podemos avaliar o desempenho dos modelos, pois conhecemos os valores dos parâmetros.

Este é o modelo mais simples, que é o ajuste de uma distribuição Poisson com um único valor para seu parâmetro:

$$Y \ \sim \ P(\lambda = \beta_0)$$

Que lemos:

Iniciamos o R e carregamos os pacote bbmle, car e sads¹⁾. O pacote bbmle foi criado por Ben Bolker e tem algoritmos de optimização para obter a máxima verosimilhança dos parâmetros. Se você ainda não tem algum dos pacotes instalado use o comando install.packages (“pacote”).

library(bbmle)
library(car)
library(sads)

Fixamos o gerador pseudo-aleatório de dados para todos obtermos exatamente o mesmos resultados:

set.seed(1234)

Vamos sortear 1000 valores de uma distribuição Poisson com parâmetro $\lambda=5$, o que simula uma população de mil plantas com uma média de cinco frutos por indivíduo. Ao usarmos uma Poisson simulamos que os frutos ocorrem a uma taxa constante por planta, e que a cada fruto é um evento independente dos demais.

nplan<-1000
lambda=5
frutos<-rpois(nplan,lambda)

Qual é a probabilidade teórica de encontrar uma planta sexualmente madura sem frutos?

dpois(0,lambda)

Compare com a proporção de plantas sem frutos na amostra simulada:

sum(frutos==0)/length(frutos)

Para fazer esta avaliação para todos os valores, comparamos a distribuição de probabilidade teórica com as proporções empíricas. Comece com um gráfico de barras das proporções:

mfrut<-max(frutos)
fa<-factor(frutos, levels=0:mfrut)
prob.obs<-table(fa)/nplan
par(las=1)
plot(0:mfrut,prob.obs, xlab="Numero de frutos", 
     ylab="Probabilidade", type="h", lwd=5)

E adicione os pontos com os valores teóricos esperados:

prob.tr<-dpois(0:mfrut, lambda)
points(0:mfrut,prob.tr, pch=21, col="red")

Definimos a função de máxima verossimilhança do parâmetro lambda da distribuição Poisson:

x<-frutos
poisNLL<-function(lambda){
-sum(dpois(x, lambda, log=TRUE))
}

Há um teorema que garante que essa função de verossimilhança é unimodal e tem um único mínimo. Por isso podemos buscar a estimativa de máxima verossimilhança de $\lambda$ simplesmente experimentando diversos valores até encontrar um mínimo. Para isso criamos uma sequência de valores na vizinhança da média amostral e aplicamos a cada um a função de log-verossimilhança:

xvec<-seq(4.85,5.3, length=1000)
LLest<-sapply(xvec,poisNLL)

Em seguida buscamos o valor mínimo:

minLLest<-min(LLest)
lambdaLL.fb<-xvec[ LLest == min(LLest)]

Vamos comparar as estimativas do lambda obtidas com os métodos de força bruta, numéricos e analíticos²⁾.

Faça a minimização numérica com a função mle2:

lambdaLL.nm<-mle2(poisNLL, start=list(lambda=4))

E agora compare os três valores:

lambdaLL.fb
lambdaLL.nm
mean(frutos)

E vamos fazer a comparação gráfica, plotando a função de verosimilhança e os pontos obtidos:

mfrutos<-mean(frutos)
LLest2<-LLest-min(LLest)
plot(xvec,LLest2, typ="l", xlab="frutos", ylab="loglik")
abline(v=mfrutos, col="blue", lwd = 3)
abline(v=coef(lambdaLL.nm),col ="darkgray")
abline(v=lambdaLL.fb, col="red")

Vamos supor que a fecundidade das plantas poderia incrementar-se com um aumento da concentração de fósforo no solo. Simulamos esta situação fazendo da taxa de frutos por planta uma função de variáveis preditoras, como o fósforo:

$$E[Y] \ = \ \lambda \ = \ e^{\beta_0 + \beta_i X_i}$$

O que é o mesmo que

$$\ln(\lambda) \ = \beta_0 + \beta_iX_i$$

Ou seja, o logarítmo de $\lambda$ é uma função linear das variáveis preditoras $X_i$³⁾. Essa é a funcão de ligação logarítmica. Sua vantagem é preservar a relação linear com as preditoras, e garantir valores positivos de $\lambda$, como deve ser na distribuição Poisson. Nosso modelo agora torna-se:

$$Y \ \sim \ \mathrm{Pois}(\lambda = e^{\beta_0 + \beta_i X_i})$$

Que lemos:

Então vamos simular a dependência de $\lambda$ ao fósforo! Primeiro definimos valores dos parâmetros:

set.seed(1234)
phos<-runif(100,0,10)
a= 1
b= 0.3
x<-phos

E então os valores esperados de frutos em função do nível de fósforo, conforme a função de ligação:

ydet<-exp(a+b*x)
plot(x,ydet)

Agora para simular um processo Poisson sorteamos amostras Poisson com parâmetro lambda igual a estes valores esperados:

fec<-rpois(100,ydet)

E plotamos nossos dados:

par(las=1)
plot(phos, fec, xlab="Fósforo mg/Kg", ylab="Número de frutos")

Definimos a função de verosimilhança para este modelo:

poisglmNLL = function(a,b) {
   ypred= exp(a+b*x)
  -sum(dpois(fec,lambda=ypred, log=TRUE))
}

Agora que simulamos um processo Poisson com valor esperado que é uma função do nível de fósforo, vamos ver o o ajuste de um modelo que descrevesse exatamente isso. Para estimar a máxima verosimilhança temos que fornecer na função de otimização um valor inicial para os parâmetros a e b.

chute <- list(a=2.5, b= 0.33)

<box red center 75% | Pergunta> Definir bons pontos de partida (starting values) para o otimizador é crucial, especialmente para funções de verossimilhança mais complexas, com muitas preditoras. Aqui o caso é simples e pouco sensível aos valores iniciais. Por isso “chutamos” um valor razoável. Mas mesmo neste caso há maneiras de obter bons palpites e ajudar o otimizador. Você tem alguma ideia de como obter uma estimativa inicial dos parâmetros? </box>

E então otimizamos a função de verossimilhança, com a função mle2:

mod.pois<-mle2(poisglmNLL, start= chute)

Conferimos o resumo do modelo:

summary(mod.pois)

E os perfis de verossimilhança das estimativas:

mod.pois.prof <- profile(mod.pois)
par(mfrow=c(1,2))
plotprofmle(mod.pois.prof) ## carregue o pacote sads para ter esta função
par(mfrow=c(1,1))

Numa abordagem frequentista, os perfis são usados para estimar intervalos de confiança dos parâmetros, transformando a log-verossimilhança negativa em uma variável normal-padrão $z$ ⁴⁾. Esta aproximação à normal quando a amostra é grande o suficiente, devido à propriedade de normalidade assintótica dos MLEs. Uma pista se isto é verdade é a forma dos perfis de verossimilhança. Se estes perfis se aproximam de parábolas a aproximação é válida, como parece ser o caso. aqui:

Aplicando a função plot ao objeto da classe perfil você obtem perfis na escala desta transformação:

plot(mod.pois.prof)

As linhas vermelhas neste gráfico marcam no eixo $x$ os quantis correspondentes a diferentes probabilidades acumuladas da normal. Estes quantis são usados como os limites dos intervalos de confiança para estas probabilidades. Compare os valores na figura para probabilidade de 95% com os obtidos com o comando:

confint(mod.pois)

Por fim, plotamos as curvas dos valores esperados de frutos por planta com os parâmetros e suas estimativas:

par(las=1)
plot(phos,fec, xlab="Fósforo mg/Kg", ylab="Número de frutos" )
a.est<-coef(mod.pois)[1]
b.est<-coef(mod.pois)[2]
curve(exp(a+b*x),add=TRUE, col="red")
curve(exp(a.est +b.est*x), add=TRUE, col="blue", lty=2) #estimada
legend("topleft", c("Parâmetro","Estimativa"),col=c("red","blue"), lty=c(1,2))

O que aconteceria se usássemos outro modelo para descrever estes dados? Neste caso sabemos qual é o modelo correto, mas vamos simular esta situação, imaginando que a pesquisadora está experimentando diferentes modelos. Usaremos a função de verosimilhança de uma binomial negativa, que permite agregação. Em nosso novo modelo mantemos a função de ligação logaritimica para a média, e supomos que o parâmetro de dispersão da binomial negativa é contante. Nosso modelo fica assim então:

$$Y \ \sim \ NB(\mu = e^{\beta_0 + \beta_i X_i, \ k = \alpha_0})$$

Que lemos:

Vamos calcular a máxima verosimilhança deste modelo com a binomial negativa e comparar os modelos com o critério de informação de Akaike (AIC).

negbinNLL<- function(a, b, k){
	ypred<-exp(a+b*x)
	-sum(dnbinom(fec, mu=ypred, size=k, log = TRUE))
        }

Tentamos um valor inicial de k usando o método dos momentos. Como sabemos que a variância da binomial negativa é:

$$\sigma^2= \mu + \frac{\mu^2}{k}$$

Então podemos estimar um valor aproximado de k com:

med<-mean(fec)
vari<-var(fec)
k.init <-med^2/(vari-med)
k.init

E então usamos este valor como palpite inicial para o valor de $k$ na otimização:

mod.negbin<- mle2(negbinNLL, start=list(a=2.5, b= 0.33, k=k.init))

E agora conferimos os parâmetros estimados:

summary(mod.negbin)

E podemos comparar os modelos com a função AICtab , que expressa o grau de plausibilidade relativa do modelo Poisson, comparado com a binomial negativa:

AICtab(mod.pois,mod.negbin, delta=T, sort=T, weights = TRUE)

O modelo com mais apoio dos dados, ou mais plausível, é o modelo com distribuição Poisson. Porém, a diferença nos valores do AIC é menor que 2, o que indica que o modelo com a Binomial Negativa tem plausibilidade equivalente, para esta amostra.

Agora pensemos que os frutos das mesmas plantas solanáceas são atacados pelo fungo Fusarium sp., murchando-os. Pensemos que esse fungo habita no solo e portanto só ataca os frutos quando eles amadurecem e caem.

Vamos usar a distribuição binomial para descrever o número de frutos infestados, entre o total de frutos produzidos pelas planta. A probabilidade de infestação dos frutos pode ser constante ou depender de algum fator como por exemplo o grau de umidade relativa do micro-habitat onde se encontra a planta. O número de “tentativas” (parâmetro $N$ da binomial) será variável entre as plantas, pois corresponde ao número de frutos que ela tenha. Vamos simular esta situação e depois verificar como um modelo binomial se sai em descrevê-la. O modelo que vamos construir pode ser representado assim:

$$Y \ \sim \ \mathrm{Bin}(p = f(\beta_0 + \beta_i X_i), \ N = n_i)$$

Que lemos:

Nas seções seguintes explicamos qual será a função de ligação $f(\cdot)$ neste caso.

Já temos um vetor com número de frutos de cada planta em nossa amostra. Primeiro descartamos as plantas sem frutos:

fec1<-fec[fec!=0]
num.plant<-length(fec1)

E simulamos o número de frutos atacados por planta:

Agora vamos simular uma variável ambiental. Por exemplo, o grau de umidade do solo na vizinhança da planta, que varia entre 30 a 100%, com probabilidade uniforme.

set.seed(4444)
ur<-runif(num.plant, 0.3, 1)

E usamos a função de ligação logística para determinar e posteriormente simular o grau de infestação do fungo de acordo a umidade relativa do ambiente. Esta função estabelece seguinte relação do parâmetro $p$ com as variáveis preditoras $X_i$

$$p \ = \frac{e^{\beta_0+\beta_i X_i}}{1+e^{\beta_0+\beta_i X_i}}$$

A logística restringe os resultados esperados ao intervalo entre 0 a 1, mesmo que a variável preditora tenha valores muito grandes ou pequenos, por isso é adequada para descrever o parâmetro $p$ nos modelos de distribuição binomial. Assim como a função de ligação log para processos Poisson, a logística “preserva” uma relação linear das preditoras com a resposta. Mais formalmente, a ligação faz com que uma transformação do parâmetro seja uma função linear das variáveis preditoras. No caso do modelo binomial, aplicar a logística acima é o mesmo que transformar $p$ para seu logito:

$$\ln \frac{p}{1-p} \ = \beta_0+\beta_iX_i $$

E novamente encontramos uma transformação do parâmetro -o logito $\ln(p/(1-p))$- que pode ser expressa como uma função linear das preditoras sem violar pressupostos da distribuição usada no modelo, no caso a binomial⁵⁾.

Vamos então usar tudo isso para construir a simulação de efeito das preditoras sobre a proporção de frutos murchos:

a<- -4
b<- 7
ydet<-exp(a+b*ur)/(1 + exp(a +b*ur))
plot(ur,ydet)
frut.mur<-rbinom(num.plant,size=fec1,prob= ydet)

O número de frutos intactos é o complemento dos murchos.

frut.sau <- fec1-frut.mur 

Para uma avaliação visual, plotamos os fruto sãos e murchos por planta e arranjamos as plantas em ordem crescente exposição à umidade relativa:

datos<-data.frame(frut.sau,frut.mur,ur)
datos.ur<-datos[order(datos$ur),]
conjun<-datos.ur[,1:2]
conjun.m<-as.matrix(conjun)
 
barplot(t(conjun.m),beside=F, col=c("darkgray","darkblue"), ylim=c(0,80), ylab= "Número de Frutos",
xlab= "Plantas", cex.names=0.4)
legend("topright", c("frutos intactos", "frutos murchos"), col = c("darkgray","darkblue"), pch=c(19,19))
arrows(15.8,62,100,62)
text(60,64, "Umidade", cex=1.8)

Para ajustar nosso modelo binomial definimos sua função de verossimilhança:

binomNLL<- function(a,b){
	prob.det=exp(a+b*ur)/(1 +exp(a+b*ur))
	-sum(dbinom(frut.mur, size=N, prob=prob.det, log=TRUE))
	}

E agora vamos encontrar o mínimo da função de log-verossimilhança negativa com otimização computacional com a função mle2.

Otimizamos:

mod.bin<-mle2(binomNLL, start=list(a =-3,b =9), data =list(N= fec1, k = frut.mur))

<box red 75% |Pergunta> Como no tutorial anterior, temos que definir valores dos parâmetros para o otimizador iniciar. O comando acima usou um “chute informado”. Descubra como chegamos a ele. Uma pista: tem a ver com plotar a razão dos frutos murchos/total dos frutos na escala de logitos. Lembre-se que a função logito lineariza a curva logística. </box>

Alternativamente podemos usar a interface de fórmula da função mle2 para fazer os dois passos simultaneamente:

mod.bin2<-mle2(frut.mur~dbinom(prob = exp(a + b*ur)/(1 + exp(a + b*ur)), size = fec1),
start = list(a= -3, b = 9))

Avaliamos os perfis de verossimilhança:

par(mfrow=c(1,2))
plotprofmle(profile(mod.bin)) ## carregue o pacote sads para ter esta função
par(mfrow=c(1,1))

Finalmente comparamos a função logística simulada com a estimada. No gráfico, o tamanho do ponto é proporcional ao número de frutos que a planta tinha.

prop.mur<-frut.mur/fec1 #Proporção de plantas murchas
par(las=1)
plot(ur, prop.mur,cex=fec1/25, xlab="Umidade relativa", ylab="Proporção Frutos Infectados")
a.esti<-coef(mod.bin)[1]
b.esti<-coef(mod.bin)[2]
curve(exp(a+b*x)/(1+exp(a+b*x)),add=TRUE, col="red")
curve(exp(a.esti+b.esti*x)/(1+ exp(a.esti +b.esti*x)), add=TRUE, col="blue", lty=2) #estimada
legend("topleft", c("Parâmetro","Estimativa"),col=c("red","blue"), lty=c(1,2))

Para modelos Poisson e binomiais como os dos tutoriais acima há uma teoria específica e um método de ajuste mais eficiente. São os glm (generalized linear models).

Neste tutorial vamos comparar os resultados usando a função padrão no R para esses modelos, glm, com os resultados obtidos antes. Os principais argumentos da função glm são family, que define a distribuição, e neste o argumento 'link', para definir a função de ligação. O padrão de link para Poisson é a logarítmica e a para binomial a logística, como usamos nos ajustes numéricos anteriores. Para saber mais consulte a ajuda da função e as leituras indicadas.

Uma regressão de uma variável Poisson como a que fizemos equivale a um glm com distribuição Poisson e função de ligação log. Compare os resultados do modelo Poisson que ajustamos com o mle2 com os resultados obtidos com a função glm do R:

glm.pois <- glm(fec~phos,family=poisson(link="log"))
summary(glm.pois)
summary(mod.pois)
AIC(mod.pois)
AIC(glm.pois)

Cada vez que fazemos um glm com distribuição Poisson e função de ligação log, compramos a premissa de uma relação exponencial entre a variável preditora e o valor esperado da variável dependente. Isso é o mesmo que supor uma relação linear entre a variável preditora e o logaritmo do valor esperado da dependente.

Uma regressão de uma variável binomial como a que fizemos equivale a um glm com distribuição binomial e função de ligação logística. Compare os resultados do modelo binomial que ajustamos com o mle2 com os resultados obtidos com a função glm do R. Note que para esta função a variável-resposta, do lado esquerdo da fórmula do modelo, é uma matriz com números de observações poisitivas e negativas.

glm.bin <- glm(cbind(frut.mur,frut.sau)~ur, 
               data=datos,family=binomial(link="logit"))
summary(mod.bin)
summary(glm.bin)               

Em um experimento, girinos em diferentes quantidades ficaram expostos a um predador por um certo tempo. Os dados são o número inicial de girinos expostos e o número predado. Use modelos de regressão binomial para testar se as observações dão suporte à hipótese de que o número final predado é afetado pelo número inicial.

1. Comece inspecionando gráficos do número de (a) girinos predados pelo número inicial e (b) da proporção de predados em relação ao número inicial.
2. Quando está abaixo da saciação o predador não ncessariamente captura 100% das preasas disponíveis. Pode ser que ele capture uma proporção fixa das presas. A função de ligação logística pode ser generalizada para qualquer assíntota $M$:

$$p \ = \ \frac{M e^{a + bX}}{1+ e^{a + bX}}$$

• girinos.csv

Tomou-se uma amostra de 120 peixes, dos quais contaram-se todos os ectoparasitas. De cada peixe tomou-se o comprimento e registrou-se o sexo. Use as funções glm e glm.nb para investigar a relação entre número de parasitas e atributos dos peixes.

Comece inspecionando os gráficos das relações entre a variável resposta e cada uma das preditoras.

• peixes.csv

Use o conjunto de dados a seguir para construir e comparar modelos de regressão Poisson e binomial negativa. Avalie todas as hipóteses de efeito linear da variável preditora com função de ligação log.

Comece enumerando os modelos possíveis. Para a regressão binomial negativa o efeito pode ser sobre a média e/ou o parâmetro de dispersão. Para cada distribuição há os modelos de ausência de efeitos.

glm_q3.csv

Faça os exercícios 205.1 e 205.2 no sistema notaR.

• Bolker, B.M. 2008 Ecological Models and Data in R Princeton: Princeton University Press:
  • Cap.6 - Likelihood and all that, seção 6.3.1
  • Cap.9 - Standard Statistics Revisited, seção 9.4

• Crawley, M.J. 2007. The R Book. (caps.13,14, e 16)

• Dobson, A.J. 1990. An Introduction to Generalized Linear Models. London: Chapman and Hall.
• McCullagh P. & Nelder, J.A. 1989. Generalized Linear Models. London: Chapman and Hall.
• GLMMs FAQs, por Ben Bolker e colaboradores. Excelente guia prático para glms com efeitos mistos
• Rodríguez, G. 2007. Lecture Notes on Generalized Linear Models.

Os resíduos usados em regressão Gaussiana não são os melhores para diagnósticos de glms. Há vários outros tipos de cálculos de resíduos propostos. Um material básico para se aprofundar:

• As revisões clássicas sobre diagnóstico de GLMs são o capítulo 12 do livro de McCullagh & Nelder e seção 2.3.4 do livro da Dobson.

• Um bom resumo sobre diferentes tipos de resíduos disponíveis para objetos glm no R.

• Vinheta do pacote Dharma, explicação muito didática de outro tipo de resíduo, os resíduos quantílicos, que podem ser aplicados a GLMs e também a modelos com efeitos aleatórios.

Quantidades análogas ao $R^2$ dos modelos lineares Gaussianos, propostas para glms e modelos com efeitos mistos. Têm o prefixo “pseudo” porque a definição estrita de $R^2$ como partição de variação devida ao preditor e devida ao erro de observação só é possível adotando-se distribuição Gaussiana.

Ainda assim, estes pseudo-$R^2$, ou índices de replicabilidade (repeatabilities), são úteis para avaliar a variação das observações captadas por diferentes termos dos modelos estatísticos. É um campo de pesquisa recente e que está desenvolvendo muito rápido. Aqui vai um material selecionado para estudo ⁶⁾:

• Ives, A.R. (2019). R2s for Correlated Data: Phylogenetic Models, LMMs, and GLMMs. Syst Biol 68, 234–251.
• Nakagawa, S. & Schielzeth, H. (2013).A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133–142.
• Nakagawa, S., Johnson, P.C.D. & Schielzeth, H. (2017).The coefficient of determination R2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of The Royal Society Interface 14, 20170213.

• Stoffel, M.A., Nakagawa, S. & Schielzeth, H. (2017). rptR: Repeatability estimation and variance decomposition by generalized linear mixed-effects models. Methods in Ecology and Evolution.
• Ives, A. & Li, D. (2018). rr2: An R package to calculate $R^2$s for regression models. JOSS 3, 1028.

GAMs são uma generalização de glms para incluir preditores não lineares.

• Hastie, T.J. and Tibshirani, R.J., 1990. Generalized additive models (Vol. 43). CRC press. A referẽncia clássica
• Capítulo 18 de Crawley, M.J. 2013. The R Book. Wiley & Sons.
• Um tutorial curto e didático