Modreg.Rmd

# Modelos de regressão

Modelos estocásticos têm sido amplamente utilizados tanto na comunidade científica 
como no mundo dos negócios em geral. Estudos de mercado usando modelos altamente 
estruturados, modelos de predição em séries financeiras, análises de componentes 
de solos para agricultura de precisão, mapeamento de doenças, modelos ecológicos, 
ambientais e climáticos são algumas das áreas de aplicações de tais modelos, 
entre tantas outras.
De todos os métodos estatísticos, os modelos de regressão talvez sejam os mais 
amplamente utilizados na prática. Nestes modelos procura-se explicar, ao menos 
parcialmente, a variabilidade de uma variável considerada como resposta, 
por outras variáveis chamadas de explicativas ou covariáveis.
Procura-se verificar a existência e quantificar o efeito das covariáveis sobre a 
resposta, embora a relação não seja necessariamente de causa e efeito.

A área de modelos de regressão teve seu início com o tradicional modelo de 
regressão linear gaussiano, que, seja por sua aplicabilidade e/ou pela facilidade 
analítica e computacional, foi e continua sendo largamente utilizado para descrever 
a associação entre um conjunto de covariáveis e uma variável resposta. 
Diversas transformações da resposta foram propostas para adaptar a suposição de normalidade. 
Porém tais modelos não são satisfatórios para respostas discretas como no caso de dados binários, 
contagens baixas e respostas restritas a certos intervalos como proporções 
e índices.

A modelagem estatística por modelos de regressão recebeu um importante novo impulso desde o
desenvolvimento dos modelos lineares generalizados no início da década de $70$ 
cujo marco inicial é o trabalho de @Nelder1972 que estendeu de uma forma mais geral
a classe de modelos de regressão, acomodando sob a mesma abordagem diversas 
distribuições agrupadas na forma da família exponencial. 
Variáveis resposta binárias e de contagens foram as que receberam mais atenção, 
talvez pela dificuldade em tratar deste tipo de resposta com transformações do 
modelo linear gaussiano. Desde então modelos de regressão logística e de Poisson, 
passaram a ser utilizados nas mais diversas áreas de pesquisa.
Nas décadas seguintes e até os dias de hoje o modelo linear generalizado (MLG) 
tornou-se uma estrutura de referência inicial, canônica, a partir da qual diversas 
modificações, expansões e alternativas são propostas na literatura.
Classes mais amplas de distribuições para resposta, preditores não lineares ou 
aditivamente lineares, modelagem conjunta de média e dispersão, adição de variáveis 
latentes no preditor são algumas das importantes extensões, às quais adicionam-se métodos 
para inferência e algorítimos para computação. Neste capítulo vamos considerar alguns 
modelos simples de regressão com ilustrações computacionais. No Capítulo seguinte 
consideramos uma extensão obtida pela inclusão de efeitos aleatórios.

Considere $Y_1, \ldots, Y_n$ variáveis aleatórias independentes e identicamente 
distribuídas e $X$ é a matriz de delineamento do modelo conhecida. 
Um modelo de regressão em notação matricial pode ser escrito da seguinte forma:
\begin{align*}
Y|X &\sim f(\boldsymbol{\mu}, \phi) \\
\boldsymbol{\mu} &= g(X;\boldsymbol{\beta}).
\end{align*}

A distribuição de probabilidade  $f(\boldsymbol{\mu}, \phi)$ da variável resposta 
é descrita por dois conjuntos de parâmetros, os de locação (média) e os de 
dispersão (precisão / variância). 
Por simplicidade de apresentação supomos aqui que o parâmetro de dispersão é 
comum a todas as observações e o que muda é apenas a média. 
Esta restrição pode ser flexibilizada, supondo que o parâmetro de dispersão, 
também possa ser descrito por alguma função das covariáveis que é possível 
requerendo uma generalização da implementação computacional. 

Nesta especificação é importante estar atento ao espaço paramétrico de 
$\boldsymbol{\mu}$ e $\phi$. O segundo em geral tem como seu espaço paramétrico 
os reais positivos, já que, o parâmetro de dispersão deve ser positivo. 
Para a parte de média devemos ter mais cuidado, uma vez que a função $g(\cdot)$ 
deve mapear o preditor linear ou não-linear, que pode assumir qualquer valor real, 
para o espaço paramétrico de $\boldsymbol{\mu}$. 
Por exemplo, se estamos trabalhando com dados de contagens onde o vetor $\boldsymbol{\mu}$ 
representa o parâmetro da distribuição de Poisson. A função $g(\cdot)$ deve então 
mapear dos reais para os reais positivos. Em outra situação, por exemplo com dados 
restritos ao intervalo $(0,1)$, a função $g(\cdot)$ deve mapear dos reais para o 
intervalo unitário e assim de forma análoga para outros espaço paramétricos.

Note que para a declaração completa do modelo são necessárias duas suposições. 
A primeira se refere a distribuição de probabilidade atribuída a variável resposta, neste caso $f(\boldsymbol{\mu}, \phi)$. A segunda é a definição de uma função $g(\cdot)$ fazendo com que o preditor, que pode apresentar qualquer valor nos reais, seja mapeado adequadamente para o espaço paramétrico da parte de média do modelo.

A função $g(\cdot)$ deve ser duplamente diferenciável. Uma restrição adicional comum é que seja estritamente monótona, porém em modelos não-lineares, isso nem sempre é necessário. Esta função tem como seus argumentos a matriz de delineamento $X$ conhecida e os parâmetros $\boldsymbol{\beta}$ desconhecidos a serem estimados. Sendo a função $g(\cdot)$ definida e mapeando os valores do preditor, seja ele linear ou não, para espaço paramétrico de $\boldsymbol{\mu}$, o modelo está composto e é passível de ter seus parâmetros estimados pela maximização da função de verossimilhança.

Um fato que pode passar desapercebido é com relação ao suporte da distribuição concordar ou não com os possíveis valores do fenômeno em estudo. Por exemplo, ao estudar o peso de animais é comum atribuir o modelo gaussiano, porém este modelo tem como suporte os reais, enquanto que o fenômeno varia apenas nos reais positivos. Nestes casos a usual justificativa prática é que a probabilidade atribuída a parte negativa é tão pequena que é desprezível na prática.

Uma outra situação é quando trabalhamos com uma variável resposta restrita ao intervalo unitário.
Podemos definir como função $g(\cdot)$ o inverso da função _logit_ e vamos mapear a média do modelo ao intervalo unitário compatível com os dados. 
Porém, na escolha da distribuição nem sempre isso é feito. 
É comum por exemplo, encontrar aplicações de modelos não-lineares atribuindo 
a distribuição gaussiana para $f(\boldsymbol{\mu},\phi)$  que tem suporte nos reais e em 
desacordo com os dados restritos a algum intervalo. 
Apesar deste tipo de construção não ser incomum em aplicações, 
é recomendado atribuir uma distribuição de probabilidade que 
tenha suporte concordante com o campo de variação da variável resposta. 
No caso de variável resposta restrita ao intervalo unitário, a distribuição beta ou simplex são opções razoáveis. Tal fato motiva um dos exemplos deste capítulo.
Quando pensamos na concordância entre dados e modelos é mais natural atribuir ao fenômeno uma distribuição cujo suporte, concorde com o campo de variação do fenômeno.

Dado a relevância de modelos de regressão em aplicações nas mais diversas áreas do conhecimento, 
neste capítulo vamos mostrar como implementar a estimação de alguns modelos de regressão mais usuais 
como regressão de Poisson, Simplex e não-linear gaussiano.
Vamos explorar o algoritmo de Newton-Raphson no modelo de regressão Poisson, para o modelo Simplex, vamos explorar o conceito de verossimilhança concentrada e exemplificar como usar a função escore obtida analiticamente dentro do algoritmo `BFGS` implementado na função `optim()`. No último modelo vamos tratar tudo numericamente apenas como uma implementação simplificada. 
Como material adicional apresentamos um estudo de caso com um modelo de regressão para análise de um experimento com contagens subdispersas.

## Regressão Poisson {#PoisReg}

Em situações onde a resposta é discreta na forma de contagens a escolha da distribuição para a variável resposta recai, ao menos inicialmente, sobre a distribuição de Poisson. 
Conforme mencionado anteriormente, um modelo de regressão qualquer é descrito por duas suposições: 
a distribucional, que se refere a distribuição de probabilidade atribuída a variável resposta 
e a uma função $g(\cdot)$, que liga o preditor à média desta distribuição.
No caso do modelo de regressão de Poisson  a distribuição tem a seguinte forma,
\[  P(Y = y) = \frac{\exp\{- \lambda\} \lambda^y}{y!}, \quad \text{ em que } \quad \lambda \ge 0 \quad 
  \text{ e } y = 0,1,2 \ldots . \]

No modelo Poisson o parâmetro indexador $\lambda$ corresponde à esperança de $Y$. 
O modelo é um caso particular da família de MLG em que temos apenas a componente de média, 
já que a variância é igual a média. 
O parâmetro $\lambda$ deve ser maior que zero e portanto devemos escolher uma função $g(\cdot)$ que mapeie dos reais aos reais positivos e a escolha usual para regressão Poisson é a função exponencial. 
No contexto de modelos lineares generalizados pode-se mostrar que a função de 
ligação canônica para o modelo Poisson é a função logaritmo cuja a inversa é a exponencial, 
por isso sua popularidade. 
Resumindo, supomos que $\mathbf{Y} \sim P(\boldsymbol{\lambda})$ e que $\boldsymbol{\lambda} = \exp\{X \boldsymbol{\beta}\}$ com a suposição adicional que o preditor é linear. 
Esta notação pode ser confusa para leitores familiarizados com distribuições multivariadas.
A princípio parece que estamos atribuindo uma distribuição para o vetor aleatório $\mathbf{Y}$ e
que portanto precisamos de uma distribuição multivariada. Na verdade estamos assumindo que os componentes do vetor $\mathbf{Y}$ são independentes e cada um de seus componentes é Poisson distribuído com o respectivo parâmetro $\lambda$. Porém, tal notação é conveniente para ilustrar as implementações computacionais.

Expressões genéricas  para função de verossimilhança, escore e informação, 
assim como algoritmos para estimação dos parâmetros são disponíveis para toda classe de MLG.
Entretanto como nosso objetivo aqui é ilustrar as computações 
vamos considerar o caso da Poisson individualmente.
A função de log-verossimilhança escrita em formato matricial é dada por,
\[  
l(\boldsymbol{\beta}, \mathbf{y}) = 
  - 1^\top g( X \boldsymbol{\beta}) + \mathbf{y}^\top X \boldsymbol{\beta} + 
  1^\top \log(\mathbf{y} !).
\]
Derivando em relação ao vetor de parâmetros $\mathbf{\beta}$ temos a função escore.
\[
  U(\boldsymbol{\beta}) = \left( \mathbf{y} - \exp\{X \boldsymbol{\beta}\} \right)^\top X 
\]
Derivando novamente obtemos a matriz de informação observada em que $\mathrm{diag}(X \boldsymbol{\beta})$ denota
uma matrix diagonal com elementos $\boldsymbol{\lambda} = X \boldsymbol{\beta}$:
\[
I_o(\boldsymbol{\beta}) = - X^\top [\mathrm{diag}(\exp\{X \boldsymbol{\beta}\})] X . 
\]
Com isso temos todos os elementos para montar o algoritmo de Newton-Raphson e
encontrar as estimativas de máxima verossimilhança para o vetor $\boldsymbol{\beta}$. 

Para exemplificar o processo vamos considerar dados simulados 
supondo que $X \boldsymbol{\beta}$ tem a forma $\beta_0 + \beta_1 x_i$,  em que $x_i$ é uma covariável com valores conhecidos.
Nas simulações a covariável assume $n=10, 50$ e $100$ valores igualmente espaçados entre 0 e 5.  
Assim, o modelo tem apenas dois parâmetros $\boldsymbol{\beta} = (\beta_0, \beta_1)^{\top}$. 
O código abaixo apresenta uma função para simular realizações deste modelo e  
mostramos o comando para gerar o conjunto de dados com $n=10$. 

```{r keep.souce = FALSE}
simula.poisson <- function(formula, beta){
  X <- model.matrix(formula)
  lambda <- exp(X %*% beta)
  y <- rpois(nrow(X), lambda=lambda)
  return(data.frame(y=y, X))}
```

```{r}
set.seed(123)
cov <- seq(0, 5, length=10)
dados10 <- simula.poisson(~cov, c(2,0.5))
```

```{r echo = FALSE, results = 'hide'}
dump("dados10", file="PoisReg.R")
```

A Figura \@ref(fig:modPoisson) apresenta a superfície de log-verossimilhança 
em escala de _deviance_ exata e pela aproximação quadrática, para diferentes tamanhos de amostras.
As linhas de contorno definem regiões de confiança que correspondem aos quantis
99, 95, 90, 70, 50, 30, 10 e  5 \% de distribuição $\chi^2_{2}$.
O ponto indica o valor verdadeiro do parâmetro. Note-se que os gráficos não estão na mesma escala.

```{r echo = FALSE}
## Verossimilhança
ll.poisson <- function(par, formula, dados){
  X <- model.matrix(as.formula(formula), data=dados)
  lambda <- drop(exp(X %*% par))
  ll <- sum(dpois(dados$y, lambda=lambda,log=TRUE))
  return(-ll)
}

hessiano <- function(par, formula, dados){
  X <- model.matrix(as.formula(formula), data=dados)
  mat <- diag(length(dados$y))
  diag(mat) <- - exp(X%*%par)
  H <- crossprod(X,mat)%*%X
  return(H)
}

dev.app <- function(par, par.hat, ...){
  dpar <- par - par.hat
  devpp <- -crossprod(dpar, hessiano(par=par.hat, ...))%*%dpar
  return(devpp)}

## Simulando os conjuntos de dados
set.seed(123)
#dados10 <- simula.poisson(b0=2,b1=0.5,n.simul=10)
cov <- seq(0, 5, length=10)
dados10 <- simula.poisson(~cov, c(2,0.5))
set.seed(123)
cov <- seq(0, 5, length=50)
#dados50 <- simula.poisson(b0=2,b1=0.5,n.simul=50)
dados50 <- simula.poisson(~cov, c(2,0.5))
set.seed(123)
#dados100 <- simula.poisson(b0=2,b1=0.5,n.simul=100)
cov <- seq(0, 5, length=100)
dados100 <- simula.poisson(~cov, c(2,0.5))
#set.seed(123)
##dados200 <- simula.poisson(b0=2,b1=0.5,n.simul=200)
#cov <- seq(0, 5, length=200)
#dados200 <- simula.poisson(~cov, c(2,0.5))

## Ajustando
par10  <- glm(y ~ cov,family=poisson,data=dados10)
par50  <- glm(y ~ cov,family=poisson,data=dados50)
par100 <- glm(y ~ cov,family=poisson,data=dados100)
#par200 <- glm(y ~ cov,family=poisson,data=dados200)

## Gride de parametros
beta0.10 <- seq(1.7,2.7,l=50)
beta1.10 <- seq(0.3,0.56,l=50)
grid.par.10 <- as.matrix(expand.grid(beta0.10,beta1.10))

## Avaliando a vero pra n = 10
ll.10 <- apply(grid.par.10,1,ll.poisson,formula="~cov",dados=dados10)
dev.p10 <- -2*(ll.poisson(par=coef(par10),formula="~cov",dados=dados10) - ll.10)
dev.ap10 <- apply(grid.par.10, 1, dev.app, par.hat = coef(par10), formula="~cov", dados=dados10)

## Avaliando a vero pra n = 50
beta0.50 <- seq(1.65,2.2,l=50)
beta1.50 <- seq(0.43,0.6,l=50)
grid.par.50 <- as.matrix(expand.grid(beta0.50,beta1.50))
ll.50 <- apply(grid.par.50,1,ll.poisson,formula="~cov",dados=dados50)
dev.p50 <- -2*(ll.poisson(par=coef(par50),formula="~cov",dados=dados50) - ll.50)
dev.ap50 <- apply(grid.par.50, 1, dev.app, par.hat = coef(par50), formula="~cov", dados=dados50)

## Avaliando a vero pra n = 100
beta0.100 <- seq(1.8,2.15,l=50)
beta1.100 <- seq(0.45,0.55,l=50)
grid.par.100 <- as.matrix(expand.grid(beta0.100,beta1.100))
ll.100 <- apply(grid.par.100,1,ll.poisson,formula="~cov",dados=dados100)
dev.p100 <- -2*(ll.poisson(par=coef(par100),formula="~cov",dados=dados100) - ll.100)
dev.ap100 <- apply(grid.par.100, 1, dev.app, par.hat = coef(par100), formula="~cov", dados=dados100)

## Avaliando a vero pra n = 200
#beta0.200 <- seq(1.85,2.15,l=50)
#beta1.200 <- seq(0.45,0.55,l=50)
#grid.par.200 <- as.matrix(expand.grid(beta0.200,beta1.200))
#ll.200 <- apply(grid.par.200,1,ll.poisson,formula="~cov",dados=dados200)
#dev.p200 <- -2*(ll.poisson(par=coef(par200),formula="~cov",dados=dados200) - ll.200)
#dev.ap200 <- apply(grid.par.200, 1, dev.app, par.hat = coef(par200), formula="~cov", dados=dados200)
```

```{r modPoisson, echo = FALSE, fig.width = 9, fig.height = 3, fig.cap = 'Superfície de deviance para diferentes tamanhos de amostra - Regressão Poisson, detalhes no texto.'}
par(mfrow=c(1,3),mar=c(2.8, 2.8, 1.7, 1),mgp = c(1.8, 0.7, 0))
LEVELS <- c(0.99,0.95,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.20,0.1,0.05)
contour(beta0.10, beta1.10, matrix(dev.p10,50,50),
        levels=qchisq(LEVELS,df=2),labels=LEVELS,
        xlab=expression(beta[0]),ylab=expression(beta[1]),main="n=10")
par(new=TRUE)
contour(beta0.10,beta1.10,matrix(dev.ap10,50,50),
        levels=qchisq(LEVELS,df=2),labels=LEVELS,
        xlab=expression(beta[0]),ylab=expression(beta[1]),axes=FALSE,
        lty=2,col="red",lwd=2)
points(2,0.5,type="p",pch=19)
## N = 50
contour(beta0.50,beta1.50,matrix(dev.p50,50,50),
        levels=qchisq(LEVELS,df=2),labels=LEVELS,
        xlab=expression(beta[0]),ylab=expression(beta[1]),
        main="n=50")
par(new=TRUE)
contour(beta0.50,beta1.50,matrix(dev.ap50,50,50),
        levels=qchisq(LEVELS,df=2),labels=LEVELS,
        xlab=expression(beta[0]),ylab=expression(beta[1]),axes=FALSE,
        lty=2,col="red",lwd=2)
points(2,0.5,type="p",pch=19)
## N = 100
contour(beta0.100,beta1.100,matrix(dev.p100,50,50),
        levels=qchisq(LEVELS,df=2),labels=LEVELS,
        xlab=expression(beta[0]),ylab=expression(beta[1]),
        main="n=100")
par(new=TRUE)
contour(beta0.100,beta1.100,matrix(dev.ap100,50,50),
        levels=qchisq(LEVELS,df=2),labels=LEVELS,
        xlab=expression(beta[0]),ylab=expression(beta[1]),axes=FALSE,
        lty=2,col="red",lwd=2)
points(2,0.5,type="p",pch=19)
```

Podemos perceber que a aproximação quadrática é excelente para os dois parâmetros, o que era esperado, já que, se tratam de parâmetros de média. 
Entretanto, mesmo esses parâmetros podem ter comportamento assimétrico
em algumas situações, como por exemplo, na presença de dados censurados. 
A relação entre os dois parâmetros também 
é clara nos gráficos mostrando que os parâmetros não são ortogonais. 
Note que com o aumento do tamanho da amostra os contornos vão ficando menores, mais concentrados em torno do verdadeiro valor do parâmetro, 
sugerindo nesta ilustração a propriedade assintótica de consistência e não-viés.

Seguindo com o processo de inferência podemos usar o algoritmo de Newton-Raphson para encontrar as estimativas de máxima verossimilhança de $\beta_0$ e $\beta_1$, para isto precisamos escrever duas funções, 
a `escore()` e a `hessiana()` conforme código \@ref(lem:esc-hess-pois).
Definimos a função hessiana de duas maneiras na liguagem `R`. 
Na primeira programamos diretamente como na expressão. Na sequência 
reescrevemos a função evitando usar a matriz completa de pesos
já que os elementos fora da diagonal são nulos e portanto desnecessários.

```{lemma, esc-hess-pois}
**Função escore e hessiana para a regressão de Poisson.**
```

```{r esc-hess-poiscod}
escore <- function(par, formula, dados){
  mf <- model.frame(formula, dados)
  X <- model.matrix(formula, data=mf)
  esco <- crossprod(model.response(mf) - exp(X %*% par), X)
  return(drop(esco))
}

## Hessiana ("naïve")
hessiano <- function(par, formula, dados){
  X <- model.matrix(formula, data=dados)
  mat <- matrix(0, nrow(X), nrow(X))
  diag(mat) <- -exp(X%*%par)
  H <- t(X) %*% mat %*% X
  return(H)
}

## Hessiana (equivalente a anterior)
hessiano <- function(par, formula, dados){
  X <- model.matrix(formula, data=dados)
  H <- crossprod(X * -exp(drop(X%*%par)),  X)
  return(H)
}
```

Usando o algoritmo de Newton-Raphson já apresentado em \@ref(lem:NR), temos as estimativas 
para o conjunto simulado com $n=10$. 

```{r echo = FALSE}
NewtonRaphson <- function(initial, escore, hessiano, tol=0.0001, 
                          max.iter, n.dim,...){
  solucao <- initial
  for(i in 2:max.iter){
    solucao <- initial-solve(hessiano(initial, ...), escore(initial, ...))
    tolera <- abs(solucao - initial)
    if(all(tolera < tol) == TRUE)break
    initial <- solucao
  }
  return(initial)
}
```

```{r}
(beta10 <- NewtonRaphson(initial=c(0,0), escore=escore, 
                         hessiano=hessiano, max.iter=1000, 
                         formula=y~cov, dados=dados10))
```

Para a construção dos intervalos de confiança assintóticos, 
basta avaliar a inversa da matriz hessiana no ponto encontrado
e os erros padrão das estimativas são dados pelas raízes quadradas dos elementos diagonais.

```{r}
Io <- -hessiano(par=beta10, formula=y~cov, dados=dados10)
Io
(erro.padrao <- sqrt(diag(solve(Io))))
```

Lembrando que as variâncias de $\hat{\beta_0}$ e $\hat{\beta_1}$ são dadas pelos termos da diagonal da inversa da matriz de informação observada, temos pela distribuição assintótica do EMV que um intervalo de $95\%$ de confiança para $\beta_0$ e $\beta_1$ poderia ser obtido por:

```{r}
beta10[1] + c(-1,1)*qnorm(0.975)*erro.padrao[1]
beta10[2] + c(-1,1)*qnorm(0.975)*erro.padrao[2]
```

Nossos resultados coincidem com as estimativas 
retornadas pela função `glm()` do `R`.

```{r}
beta10.glm <- glm(y ~ cov, family=poisson, data=dados10)
coef(beta10.glm)
summary(beta10.glm)$coef
confint(beta10.glm, type="quad")  
```

Os intervalos acima podem ser inadequados quando os parâmetros são não ortogonais 
por considerar isoladamente a curvatura em cada direção. O grau de não ortogonalidade 
também é influenciado pela ordem de magnitude das covariáveis e consequentemente dos parâmetros.
Uma primeira alternativa é a obtenção dos intervalos através de verossimilhança
(ou deviance) perfilhada, que tipicamente  fornece intervalos mais largos e por vezes assimétricos. 
Esta solução pode ser computacionalmente mais cara ou mesmo proibitiva dependendo da quantidade de observações e parâmetros no modelo.

```{r}
confint(beta10.glm)
```

```{r eval = FALSE, echo = FALSE}
llMAX <- -ll.poisson(beta10, formula=y~cov, dados=dados10)

pll.poisson.b0 <- function(b1, b0, ...)
   ll.poisson(par=c(b0, b1), ...)
pl0 <- function(b0, ...)
     optimize(pll.poisson.b0, c(0, 1), b0=b0, ...)$objective
b0 <- seq(1.9, 2.55, len=50)
prof0 <- -sapply(b0, pl0, formula=y~cov, dados=dados10)
plot(b0, -2*(prof0-llMAX), type="l", xlab=expression(beta[0]), ylab=expression(pl(beta[0])))
lines(b0, (b0-beta10[1])^2 * diag(Io)[1], col=2, type="l")
abline(h=c(1.96, 3.84))
abline(v=ab[1,])
abline(v=aa[1,], col=2)

pll.poisson.b1 <- function(b0, b1, ...)
   ll.poisson(par=c(b0, b1), ...)
pl1 <- function(b1, ...)
     optimize(pll.poisson.b1, c(1, 3), b1=b1, ...)$objective
b1 <- seq(0.35, 0.51, len=50)
prof1 <- -sapply(b1, pl1, formula=y~cov, dados=dados10)
plot(b1, -2*(prof1-llMAX), type="l", xlab=expression(beta[1]), ylab=expression(pl(beta[1])))
lines(b1, (b1-beta10[2])^2 * diag(Io)[2], col=2, type="l")
```

Uma outra solução é utilizar covariáveis centradas para obter verossimilhanças (aproximadamente) ortogonal.
Além disto as covariáveis podem ainda ser escalonadas para que os coeficientes tenham ordem de grandeza comparáveis e a função tenha contornos ao menos aproximadamente circulares.
Desta forma, os coeficientes são inicialmente estimados pontualmente e por intervalo para variáveis transformadas. Posteriormente, as estimativas são transformadas para escala original pela equação que define a reparametrização.

Com isso, exemplificamos a construção dos cálculos computacionais do modelo de regressão de Poisson. 
Procedemos a estimação por máxima verossimilhança, obtivemos o vetor escore e a matriz hessiana que possibilitaram o uso do algoritmo de Newton-Raphson. 
Vimos também que a aproximação quadrática no caso do modelo Poisson é bastante acurada, 
sendo que em nossos exemplos é visualmente difícil encontrar diferença entre ela e a exata mesmo para a amostra de tamanho 10. Usando resultados assintóticos construímos intervalos de confiança de duas formas diferentes. Além disso, com os resultados obtidos estamos aptos a proceder os testes de hipóteses, 
por exemplo, $H_0 :\beta_1 = 0$ contra $H_1 : \beta_1 \neq 0$. 
Neste caso, poderíamos usar o tradicional teste de Wald, teste escore 
ou mesmo o teste de razão de verossimilhança.

## Regressão Simplex {#reg-simplex}

Em situações práticas podemos estar diante de variáveis restritas ao intervalo $(0,1)$, como taxas, proporções e índices. O exemplo que motiva essa seção vem da área de ciências sociais na qual é comum a mensuração de variáveis latentes (não observáveis) que muitas vezes buscam medir qualidade, através de indicadores indiretos sobre o fenômeno latente. Um exemplo deste tipo de análise é o IDH - Índice de Desenvolvimento Humano, preconizado pela ONU - Organização das Nações Unidas, em todo o mundo. 

Estes índices que tradicionalmente, por construção, resultam sempre em valores no intervalo unitário, buscam medir indiretamente o nível de desenvolvimento de um país. O índice pode ser adaptado 
como o IDH-M (Índice de Desenvolvimento Humano - Municipal) para municípios, 
estados ou qualquer aglomeração de indivíduos, como por exemplo, bairros de uma cidade. 
É importante ressaltar que a construção leva a um número no intervalo unitário, 
assim o índice é expresso como uma nota em uma escala comum,
permitindo a comparação entre países, estados, municípios ou bairros.

Com foco em tais situações, selecionamos para este exemplo o Índice de Qualidade de Vida de Curitiba (IQVC). Este índice é composto por $19$ indicadores separados em $4$ áreas temáticas: Habitação, Saúde, Educação e Segurança. Estes indicadores buscam de forma indireta medir o nível de vida da população residente em cada um dos $75$ bairros da cidade. A metodologia para sua construção segue as premissas do IDH, e como tal resulta sempre em valores no intervalo unitário. Para sua construção é utilizada a base de microdados do Censo 2000, disponibilizada pelo IBGE (Instituto Brasileiro de Geografia e Estatística). Um interesse comum é relacionar o IQVC com a renda média do bairro medida em salários mínimos vigentes na época da pesquisa, neste caso ano de 2000. Os dados são disponibilizados pelo IPPUC/PR
[www.ippuc.org.br] e uma versão resumida está disponível no arquivo **simplex.txt**
dos complementos _online_.

Para a construção do modelo de acordo com a abordagem que utilizamos neste material, 
é necessário fazer duas suposições, como nos MLG's usuais, a primeira a respeito da distribuição de probabilidade atribuída a variável resposta e a segunda referente a função $g(\cdot)$ que vai mapear o preditor linear ou não ao espaço paramétrico induzido pela distribuição suposta. Para este exemplo vamos supor a distribuição como sendo simplex e a função $g(\cdot)$ como a inversa da função _logit_, comum em modelos lineares generalizados. A apresentação do modelo baseia-se no trabalho de @Miyashiro:2008.

A função densidade probabilidade simplex é indexada por dois tipos de parâmetros locação $\mu \in (0,1)$ e dispersão $\sigma^2 > 0$. Ambos podem ser decompostos em função de covariáveis observadas. Uma variável aleatória $Y$ que segue o modelo simplex tem função densidade dada por
\begin{equation}
(\#eq:simplex)
  f(y; \mu, \sigma^2) = [ 2 \pi \sigma^2 {y(1-y)}^3]^{-1/2} \exp \left\{ - d(y;\mu)/2\sigma^2 \right\}, \quad y \in (0,1),
\end{equation}
em que
\begin{equation*}
  d(y;\mu) = \frac{(y - \mu)^2}{y (1-y) \mu^2 (1- \mu)^2} .
\end{equation*}

Sejam $Y_1, \ldots, Y_n$ variáveis aleatórias independentes, sendo cada $Y_i \sim S(\mu_i, \sigma^2)$, $i = 1,2, \ldots, n.$ O modelo de regressão simplex é definido pela densidade \@ref(eq:simplex), sendo as médias $\mu_i$ dadas por
\begin{equation*}
  \mu_i = g( x_i^\top \underline{\beta}) = g(\eta_i)
\end{equation*}
em que $g(\cdot)$ é uma função que transforma valores dos reais ao intervalo unitário, $\boldsymbol{\beta} = (\beta_1, \beta_2, \ldots, \beta_p)^\top$ é o vetor de parâmetros da regressão $(\boldsymbol{\beta} \in \Re^p)$, $x_i^\top = (x_{i1}, x_{i2}, \ldots, x_{ip})^\top$ são os valores conhecidos de $p$ covariáveis e $\eta_i$ neste caso é o preditor linear. Poderíamos, sem perda de generalidade definir o preditor como

\begin{equation*}
  h(\mu_i) = x_i^\top \boldsymbol{\beta}  = \eta
\end{equation*}
onde agora $h(\cdot)$ é uma função de ligação que transforma do $(0,1)$ em $\Re$, dependendo da situação pode ser mais simples definir de uma forma ou de outra. No caso simplex vamos definir a função $h(\cdot)$ como sendo a função $logit$ que facilita a obtenção do vetor escore. O objetivo deste exemplo é mostrar como usar os otimizadores numéricos implementados na função `optim()`, porém usando o vetor escore obtido analiticamente.

Dada a função densidade de probabilidade, podemos facilmente obter a função de log-verossimilhança,
\[
  l_i(\mu_i,\sigma^2) = - 0.5 [\log(2\pi) + \log(\sigma^2) + 3 \log \{ y_i (1 - y_i)\} + d(y_i;\mu_i)/\sigma^2].
\]

Temos que $h(\mu_i) = x_i^\top \underline{\beta} = \eta_i$, então para obter o vetor escore precisamos derivar em relação a cada $\beta_p$.

\begin{equation*}
  \frac{\partial l(\boldsymbol{\beta}, \sigma^2)}{\partial \beta_p} = 
  \sum_{i=1}^n \frac{\partial l_i(\mu_i, \sigma^2)}{\partial \mu_i} \frac{\partial \mu_i}{ \partial \eta_i} \frac{\partial \eta_i}{\partial \beta_i}.
\end{equation*}

Para obter o vetor escore para a parte de média precisamos de três componentes, o primeiro é:
\begin{eqnarray}
 & \frac{\partial l(\mu_i, \sigma^2)}{\partial \mu_i} = -\frac{1}{2 \sigma^2} \frac{\partial}{\partial \mu_i} d(y_i, \mu_i) = \\ \nonumber
 &  =  \sigma^{-2} \left[\underbrace{ \frac{(y_i - \mu_i)}{(1-\mu)^2 \mu^2 (1-y)y} + \frac{(y_i - \mu_i)^2}{(1-\mu)^2 \mu^3 (1-y)y} - 
      \frac{ (y_i - \mu_i)^2}{(1-\mu_i)^3 \mu^2 (1 - y_i)y_i} }_u \right] ,
\end{eqnarray}
o segundo,
\[
  \frac{\partial \mu_i}{ \partial \eta_i} = \frac{1}{g'(\mu_i)} = \mu_i (1 - \mu_i) ,
\]
e o terceiro,
\[
  \frac{\partial \eta_i}{\partial \beta_i} =  x_{ip} .
\]

Em notação matricial, a função escore para $\boldsymbol{\beta}$ é dada por
\[ U_{\boldsymbol{\beta}}(\boldsymbol{\beta}, \sigma^2) = \sigma^{-2} X^\top T u, \]
onde $X$ é uma matriz $n \times p$, de delineamento do modelo, $T = diag\{ 1/g'(\mu_1), \ldots, 1/g'(\mu_i) \}$ e $u^\top = (u_1, \ldots, u_n)\top.$

A função escore para $\sigma^2$ é facilmente obtida, derivando $l(\boldsymbol{\beta}, \sigma^2)$ em relação a $\sigma^2$,
\begin{equation*}
  U_{\sigma^2}(\boldsymbol{\beta}, \sigma^2) = -\frac{n}{2 \sigma^2} + \frac{1}{2 \sigma^4} \sum_{i=1}^n d(y_i; \mu_i).
\end{equation*}

Note que podemos obter $\hat{\sigma}^2$ analiticamente, sua equação é dada por
\[  \hat{\sigma}^2 =  \frac{1}{n} \sum_{i=1}^n d(y_i, \hat{\mu}_i) . \]

Com isso temos todos os elementos necessários para a implementação computacional do modelo de regressão simplex. Não há uma implementação da densidade simplex entre as funções
básicas do `R`. Começamos então definindo uma função para a densidade simplex.

```{r}
dsimplex <- function(y, mu, sigma, log=TRUE){
  dis <- (y-mu)^2 / (y*(1-y)*mu^2 *(1-mu)^2)
  dsim <- -0.5*log(2*pi) -0.5*log(sigma) - 
                (3/2)*log(y*(1-y)) - (1/(2*sigma))*dis
  if(log == FALSE){dsim <- exp(dsim)}
  return(dsim)
}
```

Dada a aplicação teremos o vetor $\boldsymbol{\beta} = (\beta_0, \beta_1)$ e então o modelo de regressão simplex em `R` é apresentado no código abaixo. 
Definimos a função de ligação _logit_ para $h(\cdot) = g^{-1}(\cdot)$ e sua inversa $g(\cdot)$.
Note que como o parâmetro $\sigma^2$ tem solução analítica e depende apenas de $\mu_i$ não precisamos maximizar numericamente em relação a este parâmetro, estamos novamente usando uma log-verossimilhança concentrada. Basta substituirmos seu estimador na expressão da 
log-verossimilhança.
Em nossas implementações este procedimento foi fundamental para conseguir 
estabilidade nos resultados do algoritmo de maximização da verossimilhança. 

```{r}
inv.logit <- function(x, lambda){ 1/(1+exp(-x)) }
modelo.simplex <- function(b0, b1, formula, data){
  mf <- model.frame(formula, data=data)
  y <- model.response(mf)
  X <- model.matrix(formula, data=mf)
  mu <- inv.logit(X%*%c(b0,b1))
  d0 <- ((y-mu)^2)/(y*(1-y)*mu^2*(1-mu)^2)
  sigma <- sum(d0)/length(y)
  ll <- sum(dsimplex(y, mu=mu, sigma=sigma))
  return(-ll)
}
```

O próximo passo é implementar a função escore já obtida, note novamente que ela é função apenas de dois parâmetros $\beta_0$ e $\beta_1$, uma vez que tenhamos estes substituímos na expressão de $\hat{\sigma}^2$ e o usamos na expressão do vetor escore.

```{r}
escore <- function(b0, b1, formula, data){
  mf <- model.frame(formula, data=data)
  y <- model.response(mf)
  X <- model.matrix(formula, data=mf)
  eta <- X%*%c(b0,b1)
  mu <- inv.logit(eta)
  d0 = ((y - mu)^2) / (y*(1-y)*mu^2 * (1-mu)^2)
  sigma <- sum(d0)/length(y)
  T <- diag(c(mu*(1-mu)))
  u <- (1/(sigma*(1-mu)*mu))*(-((y-mu)^2/((1-mu)^2*mu^2)) +  
         (((y - mu)^2)/((1-mu)^2*mu^2))+((y-mu)/((1-mu)^2*mu^2)))
  es <- (crossprod(X,T) %*% u)/sigma
  return(as.numeric(-es))
}
```

Com isto, estamos aptos a passar a log-verossimilhança e o vetor escore para a `optim()` e obter as estimativas de máxima verossimilhança para $\beta_0$ e $\beta_1$, com estas substituímos na expressão de $\hat{\sigma}^2$ e concluímos o processo de estimação. Antes disso, precisamos do conjunto de dados que será analisado. O arquivo chamado **simplex.txt** traz o conjunto para a análise. O código abaixo lê a base de dados e apresenta um pequeno resumo. A Figura \@ref(fig:histSimplex) apresenta uma análise exploratória gráfica, com um histograma e um diagrama de dispersão relacionando o IQVC com a renda média do bairro em salários mínimos que foi divido por $10$ para evitar problemas numéricos.

```{r}
dados <- read.table("simplex.txt",header=TRUE)
head(dados)
```

O conjunto de dados contêm apenas três colunas, a primeira denominada $ID$ apenas identifica os diferentes bairros, a coluna $y$ apresenta o Índice de Qualidade de Vida de Curitiba (IQVC) e a coluna $MEDIA$ apresenta a renda média do bairro em salários mínimos dividido por $10$.

```{r histSimplex, echo = FALSE, fig.width = 7.5, fig.height = 4, fig.cap = 'Histograma e diagrama de dispersão - IQVC'}
par(mfrow=c(1,2),mar=c(2.8, 2.8, 1.7, 1),mgp = c(1.8, 0.7, 0), cex.lab = 0.8,
    cex.main = 0.8, cex.axis = 0.8)
hist(dados$y,ylab="Frequência",xlab="IQVC",main="")
plot(y ~ MEDIA, ylim=c(0,1), xlim=c(0,4), ylab="IQV", xlab="Renda em S.M", data=dados)
mod <- lm(y ~ MEDIA, data=dados)
abline(mod)
```

Para o ajuste do modelo simplex, vamos novamente usar as facilidades do pacote `bbmle` 
que usa internamente a função `optim()`. O ajuste é dado no seguinte código. 

```{r}
require(bbmle)
simplex.logit.gr <- mle2(modelo.simplex, start=list(b0=0,b1=0),
                         gr=escore, method="BFGS",
                         data=list(formula=y~MEDIA, data=dados))
```

A única diferença é no uso do argumento `gr` onde é passada a função escore com exatamente os mesmos argumentos adicionais usados na função de log-verossimilhança. A função `optim()` tem quatro otimizadores básicos: `Nelder-Mead`, `gradiente conjugado`, `BFGS` e `simulating annealing`. 
O uso de gradiente analítico só é possível nos algoritmos `BFGS` e `gradiente conjugado`.

O resumo tradicional do modelo pela função `summary()`, traz os erros padrões que são usados para a construção de intervalos de confiança assintóticos, que podem ser obtidos diretamente pela função 
`confint(..., method="quad")`.

```{r}
summary(simplex.logit.gr)
confint(simplex.logit.gr, method="quad")
```

Deixamos como exercício para o leitor obter o intervalo de confiança para $\sigma^2$ e os perfis de verossimilhança para $\beta_0$ e $\beta_1$. Para finalizar o exemplo, podemos visualizar o ajuste sobreposto aos dados. Para isto, vamos escrever uma função genérica para gerar os preditos.

```{lemma}
**Função para calcular os valores preditos para um modelo de regressão simplex.**
```

```{r}
my.predict <- function(modelo, formula, newdata){
  X <- model.matrix(formula,data=newdata)
  beta <- coef(modelo)[1:ncol(X)]
  preditos <- inv.logit(X%*%beta)
  return(preditos)
}
```

Usando a função `my.predict()` obtemos os valores preditos pelo modelo para rendas em 
uma sequência de valores.

```{r}
preditos <- my.predict(simplex.logit.gr, formula=~MEDIA, 
                       newdata= data.frame(MEDIA=seq(0,4,l=500)))
```

Por fim, plotamos na Figura \@ref(fig:desSimplex) um diagrama de dispersão com as observações e a média predita pelo modelo para a sequência de valores de renda.

```{r desSimplex, echo = FALSE, fig.width = 6, fig.height = 4.5, fig.cap = 'Diagrama de dispersão e modelo ajustado - IQVC 2000.'}
plot(y ~ MEDIA, ylim=c(0,1), xlim=c(0,4), ylab="IQVC", xlab="Renda em S.M/10", data=dados, cex.lab = 0.8, cex.main = 0.8, cex.axis = 0.8)
lines(preditos ~ seq(0,4,l=500))
```

## Parametrização em modelos não lineares {#naolinear}

Um modelo de regressão é não linear quando pelo menos uma derivada do
modelo em relação aos parâmetros não for livre dos parâmetros.
Em geral, a adoção de modelos de regressão não linear está associado à considerações teóricas sobre a parte mecanística/determinística do fenômeno. Modelos dessa classe são frequentes na física e na química onde os resultados para experimentos simples, como determinar o tempo de queda livre de um objeto e o produto final de uma reação química, podem ser determinados pelo conhecimento das condições experimentais. Por construção, os modelos não lineares possuem parâmetros com significado de forma que o conhecimento desses parâmetros permite uma interpretação fundamentada do fenômeno por meio dos parâmetros que o determinam/governam.

O modelo logístico é um modelo de regressão não linear largamente usado em diversas áreas. Na área biológica é empregado principalmente para descrever o crescimento ou evolução, como crescimento de seres vivos e evolução de doenças. Devido à grande aplicação por diversas áreas, o modelo apresenta diferentes parametrizações, motivadas pela conveniente interpretação que cada uma delas dá ao fenômeno. Entretanto, como já apresentado nos capítulos anteriores, a parametrização de um modelo têm um importante papel com relação a estabilidade numérica do processo de estimação e a qualidade das inferências estatísticas.

De forma simplificada, o modelo logístico tem a seguinte expressão geral
\begin{equation}
  \text{logis}(x, \theta, \beta) = \frac{\theta}{1+f(\exp\{x\},\beta)}
\end{equation}
em que $x$ é a variável independente, $\theta$ é a assíntota do modelo que representa o valor limite da função quando $x$ tende ao infinito, e $f(\cdot)$ é uma função do vetor de parâmetros $\beta$ e do exponencial de $x$. As reparametrizações desse modelo surgem das diferentes formas da função $f(\cdot)$. Vamos aqui considerar quatro parametrizações, a saber
\begin{itemize}
  \item $f_a(x) = \exp\{(a_1-x)/a_2\}$;
  \item $f_b(x) = b_1 \exp\{b_2 x\}$;
  \item $f_c(x) = \exp\{c_1+c_2 x\}$;
  \item $f_d(x) = (-1+1/d_1)\exp\{-d_2 x\}$.
\end{itemize}
Em todas essas parametrizações consideramos $\theta$ conhecido e igual a 1 
simplesmente para facilitar as ilustrações dos resultados e superfícies de verossimilhança.  
O vetor $\beta$ com dois elementos representados por letras arábicas. 
Todas as expressões são reparametrizações de um mesmo modelo. 
Para ilustrar, considere o modelo $f_a(\cdot)$ de onde obtemos que
\begin{center}
\begin{tabular}{lll}
  $b_1=\exp\{a_1/a_2\}$ & $c_1=a_1/a_2 $ &  $d_1=1/(1+\exp\{a_1/a_2\})$\\
  $b_2=-1/a_2$          & $c_2=-1/a_2  $ &  $d_2=1/a_2$. \\
\end{tabular}
\end{center}
Dessa maneira, a partir das estimativas pontuais sob qualquer modelo podemos 
obter as de qualquer outro aplicando as equações de reparametrização. Este é um fato importante pois podemos ajustar aos dados apenas o modelo que seja mais adequado para estimação e inferência e também computacionalmente e a partir deste obter estimativas em outras parametrizações que podem ser desejáveis para interpretações práticas.

Para avaliar a inferência sob diferentes parametrizações, faremos o ajuste de cada uma delas à dados de incidência de ferrugem do pessegueiro, causada pelo fungo _Tranzshcelia discolor_, em função do tempo em plantas da cultivar Chimarrita (figura \@ref(fig:incplot)). O pomar foi implantado em 2004 em Curitiba (25$^\circ$, 55' 10'' S, 49$^\circ$ 57' 26'' W, 945m), com espaçamento de 1m entre plantas de 2,5m entre linhas. Cada planta foi considerada como uma unidade experimental onde foram avaliadas o número de folhas com presença de ferrugem em relação ao total. Isso foi feito em 6 ramos marcados por planta durante 14 avaliações com intervalo de 10 à 16 dias entre avaliações. Para análise considerou-se a proporção observada de folhas atacadas.

Considerou-se o modelo gaussiano para as observações e a função de log-verossimilhança foi implementada de forma a receber qualquer modelo de regressão não linear (código \@ref(lem:logverologis)). 
Esta suposição é usual na prática embora, neste contexto, uma resposta binomial poderia ser considerada.
As funções em `R` para os modelos foram criadas e a estimação foi feita pela `optim()`. Após a estimação foram obtidos os contornos de log-verossimilhança para os pares de parâmetros
mostrados na Figura \@ref(fig:llcontornos).

```{r importacao, echo = FALSE, results = 'hide'}
#------------------------------------------------------------------------------------------
# ler

dados <- structure(list(dia = c(
                          0L, 0L, 0L, 0L, 0L, 0L, 14L, 14L, 14L, 14L, 
                          14L, 14L, 24L, 24L, 24L, 24L, 24L, 24L, 38L, 38L, 
                          38L, 38L, 38L, 38L, 54L, 54L, 54L, 54L, 54L, 54L,
                          69L, 69L, 69L, 69L, 69L, 69L, 81L, 81L, 81L, 81L,
                          81L, 81L, 92L, 92L, 92L, 92L, 92L, 92L, 106L, 106L,
                          106L, 106L, 106L, 106L, 117L, 117L, 117L, 117L, 117L,
                          117L, 128L, 128L, 128L, 128L, 128L, 128L, 142L, 142L,
                          142L, 154L, 154L, 166L, 166L),
                        inc2 = c(
                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                          0.09, 0, 0.13, 0, 0, 0.08, 0.18, 0.05, 0.38, 0.25,
                          0.13, 0.15, 0.5, 0.14, 0.33, 0.4, 0.63, 0.16, 0.5,
                          0.62, 0.56, 0.73, 0.63, 0.3, 0.89, 0.57, 0.92, 0.62,
                          0.75, 0.52, 1, 0.67, 1, 1, 1, 0.53, 1, 0.85, 1, 1, 1,
                          0.93, 1, 1, 1, 1, 1, 0.96, 1, 1, 1, 1, 1, 1, 1, 1, 1)),
                   .Names = c("dia", "inc2"), class = "data.frame")
diff(unique(dados$dia))

#------------------------------------------------------------------------------------------
```

```{lemma, logverologis}
**Função de log-verossimilhança para modelos de regressão não linear.**
```

```{r logverologiscod}
ll <- function(th, y, x, model){
  ex <- do.call(model, list(x=x, th=th))
  sd <- sqrt(crossprod(y-ex)/length(x))
  ll <- sum(dnorm(y, mean=ex, sd=sd, log=TRUE))
  ll
}
```

```{r ajustes}
# parametrizações
f.a <- function(x, th){ 1/(1+exp((th[1]-x)/th[2])) }
f.b <- function(x, th){ 1/(1+th[1]*exp(th[2]*x)) }
f.c <- function(x, th){ 1/(1+exp(th[1]+th[2]*x)) }
f.d <- function(x, th){ 1/(1+(-1+1/th[1])*exp(-th[2]*x)) }

# dados
y <- dados$inc2; x <- dados$dia

# lista com valores iniciais e modelo
init.list <- list(A=list(par=c(80,13), model=f.a),
                  B=list(par=c(120,-0.06), model=f.b),
                  C=list(par=c(5,-0.06), model=f.c),
                  D=list(par=c(0.008, 0.065), model=f.d))

# lista com termos fixos durante otmização
fixed.list <- list(fn=ll, x=x, y=y, method="BFGS",
                   control=list(fnscale=-1))

# otimização em série dos modelos
op.all <-
  lapply(init.list,
         function(i){
           op <- do.call(optim, c(i, fixed.list)); op
         })

# estimativas dos parâmetros
pars <- sapply(op.all, "[[", "par"); pars

# log-verossimilhança dos modelos
ll0 <- sapply(op.all, "[[", "value"); ll0

# contornos de log-verossimilhança
leng <- 100
grid.a <- expand.grid(th1=seq(67,79,l=leng), th2=seq(11,20,l=leng))
grid.a$ll <- apply(grid.a, 1, ll, y=y, x=x, model=f.a)
grid.a$dev <- with(grid.a, -2*(ll-ll0["A"]))

levels <- c(0.05,0.5,0.75,0.9,0.95,0.99)
qchi <- qchisq(levels, df=2)
```

```{r incplot, echo = FALSE, fig.width = 7.5, fig.height = 4.5, fig.cap = 'Valores observados e ajustados de incidência em função do tempo. As parametrizações são iguais em termos de valores ajustados.'}
require(lattice)
plot(inc2~jitter(dia), dados, ylab="Incidência", xlab="Dias", cex.lab = 0.8,
     cex.main = 0.8, cex.axis = 0.8)
curve(f.a(x, th=pars[,1]), 0, 180, add=TRUE, lwd=1.4, cex.lab = 0.8, cex.main = 0.8,
      cex.axis = 0.8)
```

Verifica-se que na parametrização $A$ os contornos de _deviance_ são praticamente elípticos com eixos paralelos aos eixos cartesianos. Do ponto de vista de otimização, os procedimentos baseados em avaliação de gradiente tem melhor taxa de convergência quando à função objetivo apresenta simetria e ortogonalidade. Do ponto de vista de inferência estatística, por exemplo na obtenção de intervalos de confiança, 
a ortogonalidade permite fazer inferência para um parâmetro sem a necessidade de 
correções considerando a curvatura na direção dos demais parâmetros.

A parametrização $C$, em comparação com $A$, apresentou rotação dos eixos do contorno com inclinação negativa e uma leve assimetria dos contornos. A rotação deve-se ao fato de que as funções de reparametrização entre $A$ e $C$ são praticamente lineares e a assimetria devido a uma delas ser um quociente. As parametrizações $B$ e $D$ apresentaram forte desvio da forma elíptica com contornos em formato de "banana". Dentre esses últimos, 
o modelo $B$ parece ser o mais susceptível e pode apresentar problemas numéricos pois apresenta 
assimetria mais acentuada do que o modelo $D$. 
Dessa forma, as parametrizações em termos de conveniência 
para estimação e inferência estatística seguem a ordem de preferência: $A$, $C$, $D$ e $B$. 
Vale relembrar que tanto as estimativas pontuais quanto intervalares do modelo em uma parametrização podem levar aos seus equivalentes em outra parametrização pelo princípio da invariância do estimador de máxima verossimilhança.

```{r echo = FALSE, results = 'hide'}
grid.b <- expand.grid(th1=seq(30,700,l=leng), th2=seq(-0.09,-0.05,l=leng))
grid.b$ll <- apply(grid.b, 1, ll, y=y, x=x, model=f.b)
grid.b$dev <- with(grid.b, -2*(ll-ll0["B"]))

grid.c <- expand.grid(th1=seq(3.5,6.7,l=leng), th2=seq(-0.09,-0.05,l=leng))
grid.c$ll <- apply(grid.c, 1, ll, y=y, x=x, model=f.c)
grid.c$dev <- with(grid.c, -2*(ll-ll0["C"]))

grid.d <- expand.grid(th1=seq(0.0005,0.027,l=leng), th2=seq(0.049,0.09,l=leng))
grid.d$ll <- apply(grid.d, 1, ll, y=y, x=x, model=f.d)
grid.d$dev <- with(grid.d, -2*(ll-ll0["D"]))

grid.full <- rbind(grid.a, grid.b, grid.c, grid.d)
grid.full$para <- rep(c("A","B","C","D"), each=leng^2)

expr <- c(expression(f[a](x) == 1/(1+exp((a[1]-x)/a[2]))),
          expression(f[b](x) == 1/(1+b[1]*exp(b[2]*x))),
          expression(f[c](x) == 1/(1+exp(c[1]+c[2]*x))),
          expression(f[d](x) == 1/(1+(-1+1/d[1])*exp(-d[2]*x))))
```

```{r llcontornos, echo = FALSE, fig.width = 8, fig.height = 7, fig.cap = 'Contornos de \\textit{deviance} para os parâmetros do modelo logístico em cada uma das parametrizações.'}
print(
contourplot(dev~th1+th2|para, data=grid.full, at=qchi, scales="free",
            labels=list(labels=as.character(levels), cex=0.8), #aspect=0.8,
            xlab=expression("(a,b,c,d)"[1]), ylab=expression("(a,b,c,d)"[2]),
            strip=strip.custom(bg="gray90", factor.levels=(expr)),
            par.settings=list(layout.heights=list(strip=2)),
            panel=function(...){
              panel.contourplot(...)
              i <- which.packet()
              panel.abline(v=pars[1,i], h=pars[2,i], lty=2)
            })
  )
```

## Modelo contagem-Gama {#vero-countgama}

Contagens são variáveis aleatórias que assumem valores inteiros não negativos. Esses valores representam o número de vezes que um evento ocorre em um domínio fixo contínuo, como um intervalo de tempo ou espaço, ou discreto, como a avaliação de um indivíduo ou setor censitário.

Em análise de dados de contagem o uso do modelo de regressão Poisson tem ocorrência predominante. Esse modelo tem sido uma escolha natural e imediata para essa classe de variável aleatória devido primeiramente à compatibilidade de suporte. Além disso, dados de contagem apresentam certa assimetria em relação à média e variabilidade que depende desse valor médio, duas características que são inerentes ao modelo de regressão Poisson. Entretanto, também inerente à esse modelo, está a forte suposição de equidispersão, ou seja, a da variância esperada ser igual à média esperada.

Afastamentos dessa suposição são frequentemente observadas em análise de dados, principalmente a superdispersão. Superdispersão corresponde à variância ser maior que a média. Isso pode ser resultado de

- Ausência de covariáveis importantes, sejam ausentes por que não foram observadas ou ausentes porque não foram declaradas, como por exemplo a omissão da interação entre covariáveis;
- Excesso de zeros, que na verdade trata-se de uma mistura de distribuições em que parte dos zeros são gerados pelo processo Poisson e parte é gerado por um processo Bernoulli;
- Diferentes amplitudes de domínio, quando eventos são contados em intervalos de tempo ou espaço em que o tamanho desses intervalos é variável e não considerado na análise, geralmente como _offset_.
- Efeito aleatório à nível de observações ou grupos de observações, quando existe no preditor teórico um componente estocástico, à nível de observações ou grupos de observação (blocos), $\ldots$.


Note que em todas as situações acima foram introduzidos processos que geram mais variabilidade (mistura com Bernoulli, efeito aleatório) ou que não permitem explicar a variabilidade existente (desconhecimento dos domínios, falta de covariáveis ou termos no modelo) de uma variável aleatória Poisson.

O contrário à superdispersão é a subdispersão, relatada com menor frequência na literatura. 
 Processos que reduzam a variabilidade da variável aleatória Poisson são pouco conhecidos. 
Nenhum dos procedimentos descritos acima para superdispersão ou modelos de efeitos aleatórios como os discutidos no Capítulo \@ref(cap:ranef) são empregados para modelar subdispersão.

Uma forma de relaxar a suposição de equidispersão é generalizar o processo Poisson. Para isso exploramos a relação que existe entre a distribuição Exponencial e a distribuição de Poisson. De maneira bem simplificada, suponha que temos eventos que ocorrem sob um domínio unidimensional de forma que o intervalo entre eventos tenha distribuição Exponencial. Se dividirmos o domínio em intervalos de igual tamanho, o número de eventos por intervalo terá distribuição de Poisson. Então, para generalizarmos a distribuição do número de eventos, basta atribuirmos outra distribuição para o intervalo entre eventos. Exatamente empregando esse raciocínio que @Winkelman1995 generalizou a distribuição Poisson empregando a distribuição Gama para o intervalo entre eventos. Entretanto, qualquer outra distribuição de suporte positivo poderia ser empregada. O notável e vantajoso da distribuição Gama é que a distribuição Exponencial é um caso particular dela. Ou seja, a distribuição Gama tem a seguinte função densidade de probabilidade
\[  
\text{Gama}(y; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}\exp\{-\beta x\}\,x^{\alpha-1},
\]
e para $\alpha= 1$ temos a função densidade de probabilidade da distribuição Exponencial como caso particular
\[
  \text{Exp}(y; \beta) = \beta\exp\{-\beta x\} .
\]

Para entender como as distribuições estão de fato ligadas, pense assim: a variância da Gama é $\alpha/\beta^2$, logo, menor valor de $\alpha$ corresponde a menor variância, quer dizer que se concentra uma grande proporção de valores ao redor da média, então intervalos entre eventos serão mais regulares e se dividimos o domínio em intervalos iguais temos um número pouco variável de eventos por intervalo (subdispersão). Por outro lado, maior valor de $\alpha$ corresponde à maior variância, logo, os intervalos entre eventos são mais variáveis/irregulares, de forma que ao dividir o domínio em intervalos temos um número de eventos por intervalo com maior variância (superdispersão).

Na Figura \@ref(fig:countgamma1) foram representados eventos ocorrendo ao longo de um domínio onde a distribuição dos intervalos entre eventos é Gama com média fixa e igual à 1. Em cada um dos casos os parâmetros usados foram sub:$(\alpha=5, \beta=5)$, equi:$(\alpha=1, \beta=1)$, super:$(\alpha=0.2, \beta=0.2)$ que correspondem às variâncias de 0.2, 1 e 5 respectivamente. As funções densidade de probabilidade estão no topo da Figura. No centro da Figura temos o domínio dividido em intervalos unitários e o número de eventos por intervalo representado. Notamos uma distribuição mais uniforme no sub que nos demais o que está de acordo com a argumentação do parágrafo anterior. Para obter os histogramas foram geradas amostras de 50000 elementos e o intervalo teve amplitude de 10 unidades para quantificação do número de eventos (por isso a média passou a ser 10).

```{r simula, echo = FALSE, results = 'hide'}
n <- 100
# alpha=1, beta=1, E(X)=1/1=1, V(X)=1/1²=1, CV=1 (equidisperso)
r1 <- rgamma(n, shape=1, rate=1)
# alpha=5, beta=5, E(X)=5/5=1, V(X)=5/5²=0.2, CV<1 (subdisperso)
r2 <- rgamma(n, shape=5, rate=5)
# alpha=0.2, beta=0.2, E(X)=0.2/0.2=1, V(X)=0.2/0.2²=5, CV>1 (superdisperso)
r3 <- rgamma(n, shape=0.2, rate=0.2)
rr <- c(r2, r1, r3)
ss <- rep(-1:1, each=n)
rr <- do.call(c, tapply(rr, ss, cumsum))
da <- data.frame(rr, ss)
str(da)

#x11(w=12, h=3)
## par(mar=c(3,3.5,3,1), cex.lab = 0.8, cex.main = 0.8, cex.axis = 0.8)
## plot(ss~rr, data=da, ylim=c(-1.25, 2))
## lim <- with(da, c(ceiling(max(tapply(rr, ss, min))), floor(min(tapply(rr, ss, max)))))
## plot(ss~rr, ylim=c(-1.25, 2), data=subset(da, rr>lim[1] & rr<lim[2]),
##      xaxt="n", yaxt="n", xlim=lim, xlab=NA, ylab=NA)
## int <- seq(lim[1], lim[2], 1)
## abline(h=-1:1, v=int, col="gray90")
## mtext(text=rev(c("super","equi","sub")), at=-1:1, side=2)
## mtext(text="Dispersão", side=2, line=1.5)
## mtext(text="Domínio", side=1, line=1)
## mtext(text="Número de eventos por intervalo", side=3, line=1)
## tb <- with(da, tapply(rr, ss, function(x) table(cut(x, int))))
## lapply(tb, function(x) c(sum(x), mean(x), var(x), var(x)/mean(x)))
## lapply(tb, table)
## nn <- cbind(do.call(c, tb), int[-1], rep(-1:1, each=length(int)-1))
## head(nn)
## text(nn[,2]-0.5, nn[,3], label=nn[,1], pos=3)

n <- 50000
r1 <- rgamma(n, shape=1, rate=1)     # equi
r2 <- rgamma(n, shape=5, rate=5)     # sub
r3 <- rgamma(n, shape=0.2, rate=0.2) # super

rr <- c(r2, r1, r3)
ss <- rep(-1:1, each=n)
rr <- do.call(c, tapply(rr, ss, cumsum))
int <- seq(ceiling(min(rr)), floor(max(rr)), by=10)
tb <- tapply(rr, ss, function(x) table(cut(x, int)))

sapply(tb, function(x) c(mean(x), var(x), var(x)/mean(x)))

## par(mfrow=c(1,3), mar=c(4,4,2,1), cex.lab = 0.8, cex.main = 0.8, cex.axis = 0.8)
## lapply(tb, function(x){
##   hist(x, breaks=(min(x)-1):max(x)+0.5, col="gray90")
##   abline(v=mean(x), lwd=2)
##   box()
## })
layout(1)
```

```{r countgamma1, echo = FALSE, results = 'hide', fig.width = 9, fig.height = 8, fig.cap = 'Funções densidade de probabilidade para diferentes distribuições de intervalos entre eventos (acima).  Eventos distribuídos ao longo de um domínio dividido em intervalos com o número de eventos por intervalo representado (centro). Histogramas do número de eventos por intervalo para cada distribuição de intervalo entre eventos (abaixo).'}
lim <- with(da, c(ceiling(max(tapply(rr, ss, min))), floor(min(tapply(rr, ss, max)))))
layout(matrix(c(1,2,3,4,4,4,5,6,7), ncol=3, byrow=TRUE))
# 1 gráfico
m <- 0.001; M <- 3
xseq <- seq(m,M,0.01)
curve(dgamma(x, shape=5, rate=5), 0, 3,
      main=expression(G(alpha==5, beta==5)), ylim=c(0,1), xlab="", ylab="")
polygon(c(xseq,M,m), c(dgamma(xseq,5,5),0,0), col="gray90")
abline(v=1)
curve(dgamma(x, shape=1, rate=1), 0, 3,
      main=expression(G(alpha==1, beta==1)), ylim=c(0,1), xlab="", ylab="")
polygon(c(xseq,M,m), c(dgamma(xseq,1,1),0,0), col="gray90")
abline(v=1)
curve(dgamma(x, shape=0.2, rate=0.2), 0, 3,
      main=expression(G(alpha==0.2, beta==0.2)), ylim=c(0,3), xlab="", ylab="")
polygon(c(xseq,M,m), c(dgamma(xseq,0.2,0.2),0,0), col="gray90")
abline(v=1)
# 2o gráfico
par(mar=c(3,3.5,2,0.5))
plot(ss~rr, ylim=c(-1.25, 2), data=subset(da, rr>lim[1] & rr<lim[2]),
     xaxt="n", yaxt="n", xlim=lim, xlab=NA, ylab=NA)
int <- seq(lim[1], lim[2], 1)
abline(h=-1:1, v=int, col="gray90")
mtext(text=c("super","equi","sub"), at=1:-1, side=2)
mtext(text="Dispersão", side=2, line=1.5)
mtext(text="Domínio", side=1, line=1)
mtext(text="Número de eventos por intervalo", side=3, line=0.5)
tb <- with(da, tapply(rr, ss, function(x) table(cut(x, int))))
nn <- cbind(do.call(c, tb), int[-1], rep(-1:1, each=length(int)-1))
text(nn[,2]-0.5, nn[,3], label=nn[,1], pos=3)
# 3o gráfico
int <- seq(ceiling(min(rr)), floor(max(rr)), by=10)
tb <- tapply(rr, ss, function(x) table(cut(x, int)))
str(tb)
attr(tb[[1]], which="name") <- "subdisperso"
attr(tb[[2]], which="name") <- "equidisperso"
attr(tb[[3]], which="name") <- "superdisperso"
# e gráfico
par(mar=c(4,4,2,1))
lapply(tb, function(x){
  hist(x, breaks=(min(x)-1):max(x)+0.5, col="gray90",
       main=attr(x, which="name"), xlab="", ylab="")
  abline(v=mean(x), lwd=2)
  box()
})
layout(1)
```

@Winkelman1995 obteve que a função de probabilidades de uma variável aleatória de contagem ($N$) pode ser escrita como função da distribuição acumulada dos intervalos entre eventos ($T$). Assim, para intervalos Gama, a distribuição de probabilidades do número de eventos por intervalo unitário é dada por

\begin{equation}
(\#eq:verocountgamma1)
  P(N=n) = G(1, \alpha n, \beta)-G(1, \alpha(n+1), \beta).
\end{equation}
A distribuição de $N$ é superdispersa para $0<\alpha<1$ e subdispersa para $\alpha>1$.

Uma vez conhecida a função de probabilidade de uma variável aleatória, podemos fazer a estimação dos parâmetros pelo método da máxima verossimilhança. Tomando a equação \@ref(eq:verocountgamma1) e estendendo para um modelo de regressão obtemos a seguinte função de log-verossimilhança

\begin{equation}
(\#eq:verocountgamma2)
  l(y_i; x_i, \alpha, \gamma) = \sum_{i=1}^{n} \log\left(
    G(1, \alpha y_i, \alpha\exp(x_i^\top\gamma))-G(1, \alpha (y_i+1), \alpha\exp\{x_i^\top\gamma)\}
  \right).
\end{equation}
em que $y_i$ é o número de eventos observado na observação $i$, $x_i$ um vetor de covariáveis conhecidas, $\gamma$ é o vetor de parâmetros associado as covariáveis $x_i$. No modelo acima estamos descrevendo a média de $T$ por um modelo de regressão uma vez que $E(T|x_i) = \alpha/\beta = \exp(-x_i^\top\gamma)$ que manipulando resulta em $\beta = \alpha\exp(x_i^\top\gamma)$. A média de $N$ por sua vez pode ser obtida pela seguinte expressão

\begin{equation}
(\#eq:en)
  E(N|x_i) = \sum_{i=1}^{\infty} G(\alpha i, \alpha\exp(x_i^\top\gamma)).
\end{equation}

Os dados estão no complemento _online_ do texto.

A subdispersão é uma característica facilmente observada no número de capulhos. Ao dispor a variância amostral como função das médias amostrais percebemos que todos os valores ficaram abaixo da linha 1:1 (Figura \@ref(fig:meansd)), portanto, uma evidência contra a suposição de equidispersão. Nessa situação, o emprego do modelo de regressão Poisson não é adequado e o modelo de regressão contagem-Gama, desenvolvido por @Winkelman1995, será aplicado. A função de log-verossimilhança do modelo está disponível no código \@ref(lem:logveroGammaCount).

```{r dados, echo = FALSE, results = 'hide'}
require(lattice) ; require(latticeExtra)
est <- c("vegetativo","botao","flor","maca","capulho")
cap <- expand.grid(rept=1:5, des=seq(0,100,l=5)/100,
                   est=factor(est, levels=est))
cap <- subset(cap, select=-rept)
cap$nc <- c(10,9,8,8,10,11,9,10,10,10,8,8,10,8,9,9,7,7,8,9,8,6,6,5,
            6,7,8,8,9,10,9,12,7,10,9,8,9,9,10,8,11,10,7,8,8,7,7,7,
            7,8,10,9,8,12,8,7,5,5,7,5,6,5,7,4,7,8,5,7,6,4,5,5,4,4,5,
            8,10,7,8,10,9,6,6,8,6,9,7,11,8,9,6,6,6,6,7,3,3,2,4,3,11,
            7,9,12,11,9,13,8,10,10,9,7,7,9,9,8,8,10,8,10,9,8,10,8,10)
## xyplot(nc~des|est, data=cap)
```

```{r meansd, echo = FALSE, fig.width = 5, fig.height = 5, fig.cap = 'Variância amostral em função da média amostral para o número de capulhos do algodão.'}
require(plyr)
aux <- ddply(cap, .(des, est), summarise, m=mean(nc), v=var(nc))

print(
xyplot(v~m, aux, aspect="xy", type=c("p","g"), col=1,
       jitter.x=TRUE, xlim=c(0,11), ylim=c(0,11),
       xlab=expression("Média amostral"~(bar(x))),
       ylab=expression("Variância amostral"~(s^2)),
       panel=function(...){
         panel.xyplot(...)
         panel.abline(a=0, b=1, lty=2)
       })
  )
```

```{lemma, logveroGammaCount}
**Função de log-verossimilhança para o modelo de regressão contagem-Gama.**
```

```{r logveroGammaCount1}
ll <- function(theta, y, X){
  ## theta: vetor de parâmetros/estimativas
  ## y: vetor de dados observados
  ## X: matrix de covariáveis associadas aos valores observados
  eXb <- exp(X%*%theta[-1]) #*theta[1]
  ## retorna a log verossimilhança para a amostra
  sum(log(pgamma(1, theta[1]*y, theta[1]*eXb)-
          pgamma(1, theta[1]*(y+1), theta[1]*eXb)))
}
```

A análise exploratória dos dados indica que o valor médio do número de capulhos depende do estágio e do nível de desfolha de forma não aditiva. Assim, usamos um modelo de interação entre os fatores com polinômio de segundo grau para efeito do nível de desfolha. Dessa maneira, a matriz do modelo é construida por meio da função `model.matrix()`. Esse modelo considera a estimação de um vetor de 11 parâmetros: intercepto único, efeito linear e efeito quadrático de desfolha por nível de estágio. A estimação usando `optim()` exige valores iniciais para iniciar o algoritmo. Um conjunto de bons valores iniciais são as estimativas dos parâmetros para o modelo de regressão de Poisson, que obtemos no `R` com a função `glm()`.

```{r optim}
X <- model.matrix(~est:(des+I(des^2)), data=cap)
## estimação modelo Poisson
rp <- glm(nc~est:(des+I(des^2)), data=cap, family=poisson)
## estimação modelo Contagem Gamma
op <- optim(c(alpha=1, coef(rp)), ll, y=cap$nc, 
            X=X, hessian=TRUE, method="BFGS", 
            control=list(fnscale=-1))
# 2*diferença da log-verossimilhança
dll <- c(diff.ll=2*abs(op$value-c(logLik(rp)))); dll 
```

Ao comparar as verossimilhanças do modelo Poisson com o contagem-Gama estamos testanto a hipótese do parâmetro $\alpha=1$. Vemos que essa hipótese foi rejeitada. De forma complementar ao teste de hipótese, podemos obter o perfil de log-verossimilhança para o parâmetro $\alpha$ e construir intervalos de confiança. Para isso escrevemos outra função de log-verossimilhança onde $\alpha$ é fixo e a estimação se dá apenas para os demais parâmetros do modelo (código \@ref(lem:logveroalpha)). Para fins de comparação, construimos intervalos de confiança baseados na aproximação quadrática da função de log-verossilhança para avaliar a qualidade da aproximação. Nos códigos abaixo tais procedimentos são ilustrados. A função `uniroot.all()` do pacote `rootSolve` é empregada para obter os intervalos de confiança baseados no perfil da deviance.

```{lemma, logveroalpha}
**Função de perfil de log-verossimilhança para o parâmetro $\alpha$ do modelo de regressão contagem-Gama.**
```

```{r logveroalphacod}
ll.alpha <- function(theta, alpha, y, X){
  ## theta: vetor de parâmetros/estimativas
  ## alpha: escalar do parâmetro alpha fixo
  ## y: vetor de dados observados
  ## X: matrix de covariáveis associadas aos valores observados
  eXb <- exp(X%*%theta) #*theta[1]
  ## retorna a log verossimilhança com um valor FIXO de ALPHA
  sum(log(pgamma(1, alpha*y, alpha*eXb)-
          pgamma(1, alpha*y+alpha, alpha*eXb)))
}
```

```{r perfil}
alpha <- sort(c(seq(3,8,l=30), op$par[1])) # valores para alpha
perfil <- sapply(alpha,
                function(a){
                 op <- optim(coef(rp), ll.alpha, alpha=a, y=cap$nc, X=X,
                             method="BFGS", control=list(fnscale=-1))
                 c(op$value, op$par[1])
                 })

coef <- op$par; vcov <- -solve(op$hessian); llik <- op$value
alp <- coef["alpha"]; sd.alp <- sqrt(vcov["alpha","alpha"])
dev.perf <- 2*(llik-perfil[1,])  # deviance da log-ver perfilhada
dev.quad <- (alp-alpha)^2/sd.alp # deviance da aprox quadrática

require(rootSolve)
qchi <- qchisq(0.95, df=1)
fperf <- approxfun(alpha, dev.perf-qchi)
lim <- uniroot.all(fperf, c(0, 10)) # limites do IC perf
lim2 <- alp+c(-1,1)*1.96*sd.alp     # limites do IC assint
```

O perfil de log-verossimilhança para $\alpha$ está representado na Figura \@ref(fig:plotprofile). A aproximação quadrática apresentou intervalos de mesma amplitude daqueles obtidos pelo perfil de verossimilhança, porém com um leve deslocamento para a esquerda. Nenhum dos intervalos contém o valor sob a hipótese nula e portanto os dados apresentam subdispersão. Assim, o número de capulhos por algodão apresenta uma distribuição mais regular que aquela prevista pela distribuição de Poisson.

```{r plotprofile, echo = FALSE, fig.width = 8, fig.height = 4, fig.cap = 'Perfil de log-verossimilhança representado pelo valor na diferença na deviance (esq.) e a raíz da diferença na \\textit{deviance} (dir.) para o parâmetro $\\alpha$. Setas indicam os limites do intervalo de confiança.'}
keytext <- c("perfil de verossimilhança","aproximação quadrática")
par(mfrow=c(1,2))
# gráfico 1
plot(dev.perf~alpha, type="l", lty=1, ylim=c(-0.5,15),
     xlab=expression(alpha), ylab=expression(2*(L(hat(alpha))-L(alpha))))
lines(dev.quad~alpha, lty=2)
legend("topright", legend=keytext, lty=1:2, bty="n", cex=1, ncol=1)
abline(h=qchi, lty=3)
arrows(lim, qchi, lim, -0.5, lty=1, length=0.1)
arrows(lim2, (alp-lim2)^2/sd.alp, lim2, -0.5, lty=2, length=0.1)
abline(v=alp)
plot(sqrt(dev.perf)~alpha, type="l", lty=1, ylim=c(-0.03,4),
     xlab=expression(alpha), ylab=expression(sqrt(2*(L(hat(alpha))-L(alpha)))))
lines(sqrt(dev.quad)~alpha, lty=2)
legend("topright", legend=keytext, lty=1:2, bty="n", cex=1, ncol=1)
#par(xpd=FALSE)
abline(h=sqrt(qchi), lty=3)
arrows(lim, sqrt(qchi), lim, 0, lty=1, length=0.1)
arrows(lim2, sqrt((alp-lim2)^2/sd.alp), lim2, 0, lty=2, length=0.1)
abline(v=alp)
```

Os valores esperados para o número de capulhos foram obtidos pela equação \@ref(eq:en), bem como os limites do intervalo de predição (Figura \@ref(fig:countpred)).

```{r tabcoef}
tabcoef <- data.frame(Estimativas=coef, ErroPadrão=sqrt(diag(vcov)))
tabcoef$zvalor <- with(tabcoef, Estimativas/ErroPadrão)
tabcoef$pvalor <- with(tabcoef, pnorm(abs(zvalor), lower=FALSE)*2)
tabcoef
```

```{r predict, echo = FALSE}
pred <- expand.grid(est=levels(cap$est), des=seq(0,1,0.05))
Xp <- model.matrix(~est:(des+I(des^2)), pred)
pred$eta <- Xp%*%coef[-1]
pred$y <- exp(pred$eta)

# obter as bandas de confiança
U <- chol(vcov[-1,-1]) # decomposição de Cholesky na vcov
pred$se <- sqrt(apply(Xp%*%t(U), 1, function(x) sum(x^2))) # erros padrões

# valores preditos num grid fino
aux <- c(pred$eta)+outer(pred$se, qnorm(0.975)*c(lwr=-1, fit=0, upr=1), "*")
n <- 1:40
aux <- matrix(sapply(aux, function(m){ sum(pgamma(1, shape=alp*(n), rate=alp*exp(m))) }),
              ncol=3, dimnames=list(NULL, c("lwr","fit","upr")))
pred <- cbind(pred, aux)
#names(pred)
```

```{r countpred, echo = FALSE, fig.width = 9, fig.height = 3.5, fig.cap = 'Valores observados e preditos acompanhados do intervalo de confiança de 95\\% para o número de capulhos do algodão em função do nível de desfolha e estágio fenológico.'}
print(
xyplot(nc~des|est, data=cap, jitter.x=TRUE, col="gray40",
       xlab="Nível artificial de desfolha", ylab="Número de capulhos",
       strip=strip.custom(bg="gray90"), layout=c(5,1))+
  as.layer(xyplot(fit~des|est, data=pred, type="l", col=1, lty=1, lwd=2))+
  as.layer(xyplot(lwr~des|est, data=pred, type="l", col=1, lty=2, lwd=1))+
  as.layer(xyplot(upr~des|est, data=pred, type="l", col=1, lty=2, lwd=1))
)
```