Modelos Lineares 2

1 Modelo de Regressão Poisson

Considere que a variável de interesse tenha a seguinte distribuição: \[\begin{eqnarray} Y_i|\lambda &\stackrel{ind}\sim& Pois(\lambda_i), \quad i=1, \ldots, n. \end{eqnarray}\]

1. Determine a função de verossimilhança aplicada em uma amostra da variável de interesse \(y_1, \ldots, y_n\) com respeito a \(\lambda_1, \ldots, \lambda_n\).

\[\begin{eqnarray*} L(\lambda_1, \ldots, \lambda_n; y_1, \ldots, y_n) &=& \prod_{i=1}^n{p(y_i|\lambda_i)}\\ &=&\prod_{i=1}^n\left[ \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!}\right]=e^{-\sum_{i=1}^n{\lambda_i}} \prod_{i=1}^n\left[ \frac{ \lambda_i^{y_i}}{y_i!}\right] \end{eqnarray*}\]

2. A função acima pertence a família exponencial? Na forma canônica? Se sim, qual a função de ligação canônica?

\[\begin{eqnarray*} p(y_1, \ldots, y_n|\lambda_1, \ldots, \lambda_n) &=& \exp\left\{ \sum_{i=1}^n y_i ln\left(\lambda_i\right) -\sum_{i=1}^n{\lambda_i} - \sum_{i=1}^n ln(y_i!) \right\} \end{eqnarray*}\] Tomando \(a(y_i) = y_i, b(\lambda_i) = ln\left(\lambda_i\right)\), \(c(\lambda_i) = -\lambda_i\) e \(d(y_i) = ln(y_i!)\), temos que a função acima pertence a família exponencial e que está na forma canônica pois \(a(y_i) = y_i\). Logo, a função de ligação canônica é \(g(\lambda_i) = ln\left( \lambda_i\right)\).

3. Determine a função de verossimilhança aplicada em uma amostra da variável de interesse \(y_1, \ldots, y_n\) com respeito a \(\beta_1, \ldots, \beta_K\).

Seja \(g(\lambda_i) = ln\left( \lambda_i\right) = \boldsymbol{X}_i^T \boldsymbol{\beta}\), sendo \(\boldsymbol{X}_i^T = (X_{i1}, \ldots, X_{ik})\) e \(\boldsymbol{\beta} = (\beta_{1}, \ldots, \beta_{k})^T\). \[\begin{eqnarray*} L(\beta_1, \ldots, \beta_k; y_1, \ldots, y_n) &=& \prod_{i=1}^n{p(y_i|\lambda_i = e^{\boldsymbol{X}_i^T \boldsymbol{\beta}})}\\ &=& e^{-\sum_{i=1}^n{e^{\boldsymbol{X}_i^T \boldsymbol{\beta}}}} \prod_{i=1}^n\left[ \frac{ e^{\boldsymbol{X}_i^T \boldsymbol{\beta} y_i}}{y_i!}\right] \end{eqnarray*}\]

4. Simule uma amostra do modelo de regressão Poisson considerando \(n=1000\), \(\boldsymbol{\beta}^T = (0,5 \; ; \; 0,7)\), \(\boldsymbol{X}_i^T = (1, x_{i2})\), \(x_{i2} \stackrel{iid}{\sim} N(0,1)\) e a função de ligação canônica.

###################################################################
# Modelo de regressao Poisson
# Yi ~ Pois(lambda_i), y_i=0,1,2,..., lambda_i>0, i=1,2,...,n
# ln(lambda_i) = x_i^T beta
###################################################################

#Simulando os dados
n       = 1000              #tamanho da amostra
x       = matrix(NA, n, 2)  #matriz de covariaveis
x[,1]   = rep(1,n)          #incluindo intercepto
x[,2]   = rnorm(n)          #vetor de covariavel
beta    = c(0.5, 0.7)       #intercepto + efeito da covariavel
lambda  = NULL              #inicializando o vetor de medias
y       = NULL              #inicializando o vetor da variavel resposta

for(i in 1:n){
  lambda[i]   = exp(x[i,] %*% beta) #calculando a media para cada unidade
  y[i]        = rpois(1,lambda[i]) #gerando a variavel resposta para cada unidade
}

5. Faça um histograma e os seguintes gráficos de dispersão: (a) da covariável simulada com distribuição normal versus a variável resposta; (b) da covariável simulada com distribuição normal versus a média da variável resposta; (c) da covariável simulada com distribuição normal versus o logaritmo neperiano da média da variável resposta.

#Analisando graficamente os dados simulados
#x11()
par(mfrow=c(2,2))
hist(y, freq=F, main="Mod. de reg. Poisson", ylab="densidade")
plot(x[,2], y, bty="n", xlab="x") #grafico de dispersao
plot(x[,2], lambda, bty="n", xlab="x")
plot(x[,2], log(lambda), bty="n", xlab="x")

6. Suponha agora que o intercepto seja desconhecido e que o efeito da covariável seja conhecido. Logo, teremos que as médias \(\lambda_i\), \(i=1, \ldots, n\), sãos desconhecidas. Crie uma grade de possíveis valores para o intercepto e calcule a função de verossimilhança para cada um destes valores. Faço o gráfico com estes valores e verifique qual valor fez com que a função de verossimilhança atingisse o ponto de máximo.

###################################################################
#Suponha agora que beta1 seja desconhecido e que beta2 seja conhecido.
#Entao lambda eh desconhecido. 
#Podemos estimar beta1 usando Estimador de Maxima verossimilhanca.
#Estimando beta1 supondo beta2 conhecido e usando EMV
###################################################################

#e_beta1   = seq(-5,5,l=1000) #criando uma grade de valores para beta1
e_beta1   = seq(-1,1,l=1000) #criando uma grade de valores para beta1
e_lambda  = NULL #criando um
lvero     = NULL
for(r in 1:length(e_beta1)){
  for(i in 1:n){e_lambda[i] = exp(x[i,] %*% c(e_beta1[r], beta[2]))}
  lvero[r] = sum(dpois(y,e_lambda, log=T))
}
#x11()
par(mfrow=c(1,1))
plot(e_beta1, lvero, lwd=2, type="l", xlab=expression(beta[1]), bty="n")
which(lvero==max(lvero)) #indice no qual esta o ponto de maximo da funcao de verossimilhanca

## [1] 741

c(e_beta1[which(lvero==max(lvero))], beta[1]) #beta1 que resultou no maximo versus valor verdadeiro

## [1] 0.4814815 0.5000000

abline(v=e_beta1[which(lvero==max(lvero))], col="red", lwd=2, lty=3)
abline(h=lvero[which(lvero==max(lvero))], col="red", lwd=2, lty=3)

#calculando a funcao de verossimilhanca de outra forma e comparando com a primeira
flvero2 = function(y, beta){
  n           = length(y)
  saida       = NULL
  for(i in 1:n){
    eta_i     = x[i,] %*% beta
    saida[i]  = -exp(eta_i) + eta_i * y[i] - log(factorial(y[i]))}
  return(sum(saida))
}
lvero2=NULL
for(r in 1:length(e_beta1)){
  lvero2[r] = flvero2(y, c(e_beta1[r], beta[2]))
}
plot(lvero, lvero2)

2 Modelo de Regressão Logístico

Considere que a variável de interesse tenha a seguinte distribuição: \[\begin{eqnarray} Y_i|\theta_i &\stackrel{ind}\sim& Bern(\theta_i), \quad i=1, \ldots, n. \end{eqnarray}\]

1. Determine a função de verossimilhança aplicada em uma amostra da variável de interesse \(y_1, \ldots, y_n\) com respeito a \(\theta_1, \ldots, \theta_n\).

2. A função acima pertence a família exponencial? Na forma canônica? Se sim, qual a função de ligação canônica?

3. Determine a função de verossimilhança aplicada em uma amostra da variável de interesse \(y_1, \ldots, y_n\) com respeito a \(\beta_1, \ldots, \beta_K\).

4. Simule uma amostra do modelo de regressão Logístico considerando \(n=1000\), \(\boldsymbol{\beta}^T = (0,5 \; ; \; 0,7)\), \(\boldsymbol{X}_i^T = (1, x_{i2})\), \(x_{i2} \stackrel{iid}{\sim} N(0,1)\) e a função de ligação canônica.

5. Faça um histograma e os seguintes gráficos de dispersão: (a) da covariável simulada com distribuição normal versus a variável resposta; (b) da covariável simulada com distribuição normal versus a média da variável resposta; (c) da covariável simulada com distribuição normal versus o logito da média da variável resposta.

6. Suponha agora que o intercepto seja desconhecido e que o efeito da covariável seja conhecido. Logo, teremos que as médias da variável resposta sãos desconhecidas. Crie uma grade de possíveis valores para o intercepto e calcule a função de verossimilhança para cada um destes valores. Faço o gráfico com estes valores e verifique qual valor fez com que a função de verossimilhança atingisse o ponto de máximo.