hurdle.Rnw

<<setup-child, include=FALSE>>=
set_parent("slides-mrdcr.Rnw")

## Pacotes utilizados nesta seção

@

\begin{frame}[allowframebreaks]{Modelos \textit{Hurdle}}

\begin{itemize}
    \item A variável de interesse é particionada em contagens nulas e 
    não nulas;
    \item Consideram somente os zeros estruturais;
    \item São chamados também de modelos condicionais, hierárquicos ou de 
    duas partes;
    \item Esta abordagem combina um modelo de contagem truncado à esquerda
    do ponto $y = 1$ e um modelo censurado à direita no mesmo ponto
    $y =1$
\end{itemize}

\framebreak

\begin{block}{Distribuição de probabilidades}

\begin{equation*}
    Pr(Y = y) =
    \begin{dcases*}
        f_z(0) & \text{se } y = 0,\\
            (1 - f_z(0)) \frac{f_c(Y = y)}{1 - 
                f_c(Y = 0)} & \text{se } y = 1, 2, \dots
    \end{dcases*}
\end{equation*}

em que $f_z$ é uma função de probabilidades degenerada no ponto 0 e
$f_c$ um função de probabilidades de uma variável $Y^*$, como a Poisson.

\end{block}

\begin{block}{Momentos da distribuição}
\begin{columns}[t]
\column{.48\textwidth}
\begin{center}
{\bf Média}
\end{center}
$$
E(Y) = \frac{E(Y^*)(1-f_z(0))}{1-f_c(Y=0)}
$$
\column{.48\textwidth}
\begin{center}
{\bf Variância}
\end{center}
$$
V(Y) = \frac{1-f_z(0)}{1-f_c(Y=0)} \left [ E(Y^*)
    \frac{(1-f_z(0))}{1-f_c(Y=0)} \right ]
$$

\end{columns}
\end{block}

\end{frame}

\begin{frame}{Distribuição \textit{Hurdle}}

\begin{columns}[c]
\column{.4\textwidth}

\begin{itemize}
    \item $f_z$ é uma função de probabilidades degenerada no ponto $y=0$,
    ou seja, tem toda massa no ponto 0.
    \item $f_c$ é uma função de probabilidades tradicional, que no modelo
    é truncada em $y=1$.
    \item Os modelos de barreira combinam $f_z$ e $f_c$ para descrever $Y$
    \item Para a parte positiva os dados ainda podem apresentar sub,
    superdispersão ou excesso de valores em outro ponto.
\end{itemize}

\column{.6\textwidth}
<<fig.height=3, fig.width=4.5, out.width="1\\textwidth">>=

set.seed(2016)
n <- 1000

lambda <- 3; pi <- 0.9
y <- sapply(rbinom(n, 1, pi), function(x) {
    ifelse(x == 0, 0, rpois(1, lambda))
})

yu <- sort(unique(y))
py_real <- c(prop.table(table(y)))
m0 <- glm(y ~ 1, family = poisson)
py_pois <- dpois(yu, exp(m0$coef))

cols <- c(trellis.par.get("dot.symbol")$col,
          trellis.par.get("superpose.line")$col[2])
key <- list(corner = c(0.95, 0.9), 
            lines = list(lty = 1, col = rev(cols), lwd = 3, size = 3),
            text = list(expression(f[z], f[c])))

ylim <- extendrange(c(0, max(py_real, py_pois)))
xyplot(py_pois ~ yu, type = "h", lwd = 5, grid = TRUE,
       xlab = "",
       ylab = expression(P(Y==y)),
       ylim = ylim, key = key,
       scales = list(x = list(at = yu)),
       panel = function(x, y, ...) {
           panel.xyplot(x, y, ...)
           panel.lines(x = x[1], y = py_real[1], type = "h",
                       col = cols[2], lwd = 5)
           panel.lines(x = x[1], y = py_pois[1], type = "h",
                       col = cols[2], lwd = 5)
       }
)

@

\end{columns}

\end{frame}

\begin{frame}{Combinações comuns}

Pode-se propor diferentes distribuições para $f_z$ e $f_c$. Uma
escolha natural para $f_z$ é a Bernoulli e para $f_c$ a Poisson. Assim

\begin{columns}[c]
\column{.3\textwidth}
\begin{align*}
    &f_z \sim Bernoulli(\pi) \\
    &f_c \sim Poisson(\lambda)
\end{align*}

\column{.1\textwidth}
$$\Rightarrow$$

\column{.6\textwidth}
\begin{flalign*}
    &P(Y = y) = \begin{dcases*}
        1 - \pi & \text{se } y = 0,\\
            \pi \left ( \frac{e^{-\lambda} \lambda^y}{y! 
                (1- e^{-\lambda})} \right ) & \text{se } y = 1, 2, \dots
    \end{dcases*}&
\end{flalign*}
\end{columns}

\vspace{0.8cm}
Embora essa escolha de modelo seja o que tem o maior suporte
computacional, ressalta-se que outras distribuições podem ser escolhidas
para ambas as partes $f_z$ e $f_c$.

\end{frame}

\begin{frame}[allowframebreaks]{Modelos de regressão \textit{Hurdle}}

\begin{itemize}
    \item Incorporando covariáveis em $f_z$ e $f_c$ na forma $h(Z\gamma)$
    e $g(X\beta)$, respectivamente.
    \item As funções $h(.)$ e $g(.)$, são as funções de ligação escolhidas
    conforme modelos $f_z$ e $f_c$.
    \item O modelo de regressão {\it Hurdle} terá, portanto, os vetores de
    parâmetros $\beta$, $\gamma$ e potencialmente $\phi$ (caso um modelo
    com parâmetro de dispersão for considerado)
    \item Se os modelos para $f_z$ e $f_c$ e as respectivas matrizes 
    $Z$ e $X$ forem as mesmas, o teste $H_0: \beta = \gamma$ avalia a
    a necessidade do modelo Hurdle.
\end{itemize}

\framebreak

\begin{columns}[t,onlytextwidth]
\column{.39\textwidth}

\begin{block}{Função de verossimilhança}
    \begin{align*}
        L(\underline{\theta}; &\underline{y}) = 
        \prod_{i=1}^n (1-\mathds{1}) \left ( f_{z_i}(0) \right ) \cdot \\ 
        &\prod_{i=1}^n \mathds{1} \left ( (1-f_{z_i}(0)) \left (
            \frac{f_{c_i}(y_i)}{1 - f_{c_i}(0)}\right ) \right )
    \end{align*}
\end{block}

\column{.58\textwidth}

\begin{block}{Função de log-verossimilhança}
    \begin{align*}
        l(\underline{\theta}; &\underline{y}) = \sum_{i = 1}^n 
        (1-\mathds{1}) \left ( \log(f_{z_i}(0)) \right ) + \\
        &\sum_{i = 1}^n \mathds{1} \left ( \log(1-f_{z_i}(0)) + 
            \log(f_{c_i}(y_i)) - \log(1 - f_{c_i}(0)) \right )
    \end{align*}
\end{block}
\end{columns}

\vspace{0.8cm}
Sendo $\mathds{1}$ a função indicadora que assume o valor 1 se $y > 0$ e 
$\beta$, $\gamma$ e $\phi$, se houver).
0 se $y = 0$ e $\underline{\theta}$ o vetor de parâmetros do modelo.

\end{frame}

\begin{frame}[fragile, allowframebreaks]{Modelos \textit{Hurdle} no R}

Neste minicurso utilizaremos principalmente pacote o {\tt pscl} 
(\textit{Political Science Computational Laboratory, Stanford University})

<<eval=FALSE, echo=TRUE>>=

library(pscl)
hurdle(y ~ fc_preditor | fz_preditor, dist = "poisson", zero.dist = "poisson")

@


\framebreak

Um outro pacote que proporciona diversas funções e podemos adaptar para
o ajuste desses modelos é o {\tt VGAM} ({\it Vector Generalized Linear
and Additive Models})

<<eval=FALSE, echo=TRUE>>=

library(VGAM)
vglm(y ~ preditor, family = zapoisson)

## ou ajustando as partes
vglm(y ~ fc_preditor, family = pospoisson, data = subset(data, y > 0))
vglm(SurvS4(cy, st) ~ fz_preditor, cens.poisson, 
     data = transform(data, cy = pmin(1, y), st = ifelse(y >= 1, 0, 1))

@

\end{frame}