1 Modelo de Regressão Poisson

1.1 Aula 04

Considere que a variável de interesse tenha a seguinte distribuição: \[\begin{eqnarray} Y_i|\lambda_i &\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*}\]

  1. 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)\).

  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 \(\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*}\]

  1. 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
}
  1. 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")

  1. 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] 763
c(e_beta1[which(lvero==max(lvero))], beta[1]) #beta1 que resultou no maximo versus valor verdadeiro
## [1] 0.5255255 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)

1.2 Aula 05

1.2.1 Estimativa pontual - Método Escore (nº fixo de iterações)

  1. Obtenha uma estimativa pontual para os coeficientes \(\boldsymbol{\beta}\), usando o método do Escore com número fixo de iterações.
###################################################################
#Suponha agora que beta1 e beta2 sejam desconhecidos.
#Entao lambda eh desconhecido. 
#Podemos estimar beta1 e beta2 usando Estimador de Maxima verossimilhanca.
#Para isso, usaremos Newton-Raphson, que neste caso eh semelhante
#ao metodo do Escore.
###################################################################

#Aplicando o metodo do Newton-Raphson com numero fixo de iteracoes
nite          = 100 #total de iteracoes
est_beta      = matrix(NA, 2, nite) #matriz que guardara os valores de beta
est_beta[,1]  = c(7,5)#rep(0,2) #inicializando beta
epsilon       = matrix(NA, 2, nite-1)
for(ite in 2:nite){
  U = rep(0,ncol(x))
  for(j in 1:ncol(x)){
    for(i in 1:n){
      U[j]          = U[j] - x[i,j] * exp(x[i,]%*%est_beta[,ite-1]) + y[i]*x[i,j]
    }
  }
  Ulinha = matrix(0, ncol(x), ncol(x))
  for(j1 in 1:ncol(x)){
    for(j2 in 1:ncol(x)){
      for(i in 1:n){
        Ulinha[j1, j2]     = Ulinha[j1, j2] - x[i,j1] * x[i, j2] * exp(x[i,]%*%est_beta[,ite-1]) 
      }
    }
  }
  est_beta[,ite]    = est_beta[,ite-1] - solve(Ulinha) %*% U
  epsilon[,ite-1]   = est_beta[,ite] - est_beta[,ite-1]
}

est_beta[,nite]
## [1] 0.5327077 0.6906265
# x11()
 par(mfrow=c(1,2))
 plot(est_beta[1,], bty="n", xlab="iteração", ylab=expression(beta[1]))
 plot(est_beta[2,], bty="n", xlab="iteração", ylab=expression(beta[2]))

# x11()
 par(mfrow=c(1,2))
 plot(epsilon[1,], bty="n", xlab="iteração", ylab="diferença entre estimativas", main=expression(beta[1]))
 plot(epsilon[2,], bty="n", xlab="iteração", ylab="diferença entre estimativas", main=expression(beta[2]))

# x11()
# par(mfrow=c(1,2))
# plot(epsilon[1,400:410], bty="n", xlab="iteração", ylab="diferença entre estimativas", main=expression(beta[1]))
# plot(epsilon[2,400:410], bty="n", xlab="iteração", ylab="diferença entre estimativas", main=expression(beta[2]))
#

1.2.2 Estimativa pontual - Método Escore (tolerância)

  1. Obtenha uma estimativa pontual para os coeficientes \(\boldsymbol{\beta}\), usando o método do Escore de forma que o algoritmo seja executado enquanto a diferença entre o novo valor e o valor anterior seja menor que \(0,01\).
# #Aplicando o metodo do Newton-Raphson ate que epsilon seja menor ou igual a 0.01
 est_beta      = matrix(NA, 2, 0) #matriz que guardara os valores de beta
 est_beta      = cbind(est_beta, rep(0,2)) #inicializando beta
 dif           = matrix(NA, 2, 0)
 epsilon       = 1
 ite           = 1 
while(epsilon>0.01){
  ite               = ite + 1
  U = rep(0,ncol(x))
  for(j in 1:ncol(x)){
    for(i in 1:n){
      U[j]          = U[j] - x[i,j] * exp(x[i,]%*%est_beta[,ite-1]) + y[i]*x[i,j]
    }
  }
  Ulinha = matrix(0, ncol(x), ncol(x))
  for(j1 in 1:ncol(x)){
    for(j2 in 1:ncol(x)){
      for(i in 1:n){
        Ulinha[j1, j2]     = Ulinha[j1, j2] - x[i,j1] * x[i, j2] * exp(x[i,]%*%est_beta[,ite-1])
      }
    }
  }
  aux    = est_beta[,ite-1] - solve(Ulinha) %*% U
  est_beta          = cbind(est_beta, aux)
  dif               = cbind(dif,  abs(est_beta[,ite] - est_beta[,ite-1]))
  epsilon           = max(dif[,ite-1])
}
ite
## [1] 7
dif[,ite-1]
## [1] 0.001741162 0.007075435
epsilon
## [1] 0.007075435
cbind(est_beta[,ite], beta)
##                beta
## [1,] 0.5327219  0.5
## [2,] 0.6906757  0.7
# x11()
par(mfrow=c(1,2))
plot(est_beta[1,], bty="n", xlab="iteração", ylab=expression(beta[1]))
plot(est_beta[2,], bty="n", xlab="iteração", ylab=expression(beta[2]))

#x11()
par(mfrow=c(1,2))
plot(dif[1,], bty="n", xlab="iteração", ylab="diferença entre estimativas", main=expression(beta[1]))
plot(dif[2,], bty="n", xlab="iteração", ylab="diferença entre estimativas", main=expression(beta[2]))

1.2.3 Estimativa pontual - Método Escore (tolerância) (forma 2)

  1. Obtenha uma estimativa pontual para os coeficientes \(\boldsymbol{\beta}\), usando o método do Escore de forma que o algoritmo seja executado enquanto a diferença entre o novo valor e o valor anterior seja menor que \(0,01\) e usando as fórmulas desenvolvidas para o caso geral de MLG (usando z e w).

\[\begin{eqnarray} \hat{\boldsymbol{\beta}}^{(ite)} &=& \left( \boldsymbol{X}^T\boldsymbol{W}\boldsymbol{X}\right)^{-1} \left( \boldsymbol{X}^T\boldsymbol{W}\boldsymbol{Z}\right)\\ W_{ii} &=& e^{\boldsymbol{X}^T_i\boldsymbol{\beta}^{(ite-1)}}\\ z_i &=& \boldsymbol{X}^T_i\boldsymbol{\beta}^{(ite-1)} + y_i e^{\boldsymbol{X}^T_i\boldsymbol{\beta}^{(ite-1)}} - 1 \end{eqnarray}\]

# Estimação por máxima veroximilhança - Método escore
maxite        = 10000   #total maximo de iteracoes
tol           = 10^(-2) #diferenca desejada entre os valores gerados
beta_ant      = matrix(NA,nrow=2,ncol=1)
beta_ant[1,1] = 7
beta_ant[2,1] = 5
w             = NULL
z             = NULL
for(i in 1:n){
  eta     = x[i,] %*% beta_ant[,1]
  w[i]    = exp(eta)
  z[i]    = eta + y[i]*exp(-eta) - 1
}
W             = diag(w)
beta_novo = (solve(t(x)%*%W%*%x))%*%t(x)%*%W%*%z

for(ite in 2:maxite){
  
  #print(ite)
  beta_ant = beta_novo
  for(i in 1:n){
    eta     = x[i,] %*% beta_ant[,1]
    w[i]    = exp(eta)
    z[i]  = eta + y[i]*exp(-eta) - 1
  }
  W         = diag(w)
  beta_novo = solve(t(x)%*%W%*%x)%*%t(x)%*%W%*%z
  
  if(max(abs(beta_novo - beta_ant))<=tol){
    print("Convergiu!")
    print(paste("Iteração:",ite))
    print("----------")
    break()
  }
  
  if (ite==maxite){
    print("Não convergiu!")
  }
}
## [1] "Convergiu!"
## [1] "Iteração: 22"
## [1] "----------"
# Estimativas pontuais dos coeficientes
betachapeu = beta_novo
betachapeu
##           [,1]
## [1,] 0.5327078
## [2,] 0.6906266

1.3 Aula 06

1.3.1 Estimativa Intervalar

Crie um intervalo de confiança para \(\boldsymbol{\beta}\).

alfa        = 0.05
#betachapeu  = est_beta[,ite]
zalfa2      = qnorm(1-alfa/2) 

Ulinha        = matrix(0, ncol(x), ncol(x))
for(j1 in 1:ncol(x)){
  for(j2 in 1:ncol(x)){
    for(i in 1:n){Ulinha[j1, j2] = Ulinha[j1, j2] - x[i,j1] * x[i,j2] * exp(x[i,] %*% betachapeu)}
  }
}
Iobs        = -Ulinha
Iesp        = Iobs
Inv_Iesp    = solve(Iesp) #inversa da matriz de informacao de Fisher Esperada
j=1
print(paste("IC para ",expression(beta[1]),": ", sep=""))
## [1] "IC para beta[1]: "
round(c(betachapeu[j] - zalfa2*sqrt(Inv_Iesp[j,j]), betachapeu[j] + zalfa2*sqrt(Inv_Iesp[j,j])), 2)
## [1] 0.48 0.58
j=2
print(paste("IC para ",expression(beta[2]),": ", sep=""))
## [1] "IC para beta[2]: "
round(c(betachapeu[j] - zalfa2*sqrt(Inv_Iesp[j,j]), betachapeu[j] + zalfa2*sqrt(Inv_Iesp[j,j])),2)
## [1] 0.65 0.73

1.4 Aula 07 - Adequabilidade do modelo

1.4.1 Deviance

  1. Para avaliar a adequabilidade do modelo, plote a variável resposta versus a covariável. Em seguida, acrescente a curva ajustada através do \(\hat{y}\) (y ajustado).
p     = ncol(x)

# Valores ajustados
ychapeu = NULL
ychapeu = exp(betachapeu[1] + x[,2]*betachapeu[2])
#ychapeu
plot(x[,2],y,cex=2,lwd=2,type="p",xlab='x',ylab='y',pch=16)
points(x[,2],ychapeu,lwd=2,col=2,cex=2,cex.lab=1.5,cex.axis=1.5,pch=16)

  1. Encontre o EMV do modelo saturado.

Primeiro vamos entender o modelo saturado. De acordo com este modelo, temos: \[\begin{eqnarray} Y_1 &\sim& Poi(\lambda_1)\\ Y_2 &\sim& Poi(\lambda_2)\\ &\vdots& \\ Y_n &\sim& Poi(\lambda_n) \end{eqnarray}\] e \[\begin{eqnarray} ln(\lambda_1) &=& \beta_1\\ ln(\lambda_2) &=& \beta_2\\ &\vdots& \\ ln(\lambda_n) &=& \beta_n. \end{eqnarray}\]

Esse modelo pode ser reescrito da seguinte forma: \[\begin{eqnarray} Y_i &\sim& Poi(\lambda_i)\\ ln(\lambda_i) &=& \boldsymbol{X}_i^T \boldsymbol{\beta}, \end{eqnarray}\] sendo \(\boldsymbol{X}_i^T = (x_{i1} \ldots x_{in})\), com \(x_{ii}=1\) e \(x_{ij}=0\) para todo \(i \neq j\) e \(\boldsymbol{\beta} = (\beta_1, \ldots, \beta_n)^T\).

Sejam \(\boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_n)^T\) e \(\boldsymbol{y} = (y_1, \ldots, y_n)^T\). Dessa forma, a função de verossimilhança é \[\begin{eqnarray*} L(\boldsymbol{\lambda}; \boldsymbol{y}) &=& \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]. \end{eqnarray*}\]

E então o log da função de verossimilhança é \[\begin{eqnarray*} l(\boldsymbol{\lambda}; \boldsymbol{y}) &=& \sum_{i=1}^n\left[-\lambda_i + y_i \ln(\lambda_i) - \ln(y_i!)\right] \end{eqnarray*}\] e então a derivada do log da função de verossimilhança com respeito a \(\lambda_j\) é \[\begin{eqnarray*} \frac{\partial l(\boldsymbol{\lambda}; \boldsymbol{y})}{\partial \lambda_j} &=& -1 + \frac{y_j}{\lambda_j} = 0 \Leftrightarrow \lambda_j = y_j. \end{eqnarray*}\] e a segunda derivada do log da função de verossimilhança com respeito a \(\lambda_j\) é \[\begin{eqnarray*} \frac{\partial^2 l(\boldsymbol{\lambda}; \boldsymbol{y})}{\partial \lambda_j^2} &=& -\frac{y_j}{\lambda_j^2} \end{eqnarray*}\] Como \(\frac{\partial^2 l(\boldsymbol{\lambda}; \boldsymbol{y})}{\partial \lambda_j^2}< 0\), temos que \(\hat{\lambda}_j = y_j\) é o EMV do modelo saturado.

  1. Calcule a Deviance do modelo ajustado.
#Calculando a Deviance do modelo saturado (tambem chamado de maximo ou completo)
lambdaSaturado = lambdaAjustado = NULL
lveroSaturado = lveroAjustado = 0
for(i in 1:n){
  lambdaSaturado[i] = y[i]
  lambdaAjustado[i] = exp(x[i,] %*% betachapeu)
  lveroSaturado     = lveroSaturado + dpois(y[i], lambdaSaturado[i], log=T)
  lveroAjustado     = lveroAjustado + dpois(y[i], lambdaAjustado[i], log=T)
}
Deviance = 2*(lveroSaturado - lveroAjustado)
Deviance
## [1] 1040.819
  1. Calcule a Deviance do modelo nulo.

Primeiro vamos entender o modelo nulo. De acordo com este modelo, temos: \[\begin{eqnarray} Y_1 &\sim& Poi(\lambda_1)\\ Y_2 &\sim& Poi(\lambda_2)\\ &\vdots& \\ Y_n &\sim& Poi(\lambda_n) \end{eqnarray}\] e \[\begin{eqnarray} ln(\lambda_1) &=& \beta\\ ln(\lambda_2) &=& \beta\\ &\vdots& \\ ln(\lambda_n) &=& \beta. \end{eqnarray}\]

Esse modelo pode ser reescrito da seguinte forma: \[\begin{eqnarray} Y_i &\sim& Poi(\lambda_i)\\ ln(\lambda_i) &=& X_i \beta = \beta, \end{eqnarray}\] sendo \(X_i = 1\), \(\forall i\).

Sejam \(\boldsymbol{\lambda} = (\lambda_1, \ldots, \lambda_n)^T\) e \(\boldsymbol{y} = (y_1, \ldots, y_n)^T\). Dessa forma, o log da função de verossimilhança é \[\begin{eqnarray*} l(\beta; \boldsymbol{y}) &=& \sum_{i=1}^n\left[-e^\beta + y_i \beta - \ln(y_i!)\right] \\ &=& -n e^{\beta} + \beta \sum_{i=1}^n y_i - \sum_{i=1}^n\ln(y_i!). \end{eqnarray*}\] e então a derivada do log da função de verossimilhança com respeito a \(\beta\) é \[\begin{eqnarray*} \frac{\partial l(\beta; \boldsymbol{y})}{\partial \beta} &=& -ne^{\beta} + \sum_{i=1}^n y_i = 0 \Leftrightarrow e^{\beta} = \bar{y} \\ &\Leftrightarrow& \beta = \ln(\bar{y}). \end{eqnarray*}\] e a segunda derivada do log da função de verossimilhança com respeito a \(\beta\) é \[\begin{eqnarray*} \frac{\partial^2 l(\beta; \boldsymbol{y})}{\partial \beta^2} &=& -ne^{\beta} \end{eqnarray*}\] Como \(\frac{\partial^2 l(\beta; \boldsymbol{y})}{\partial \beta^2}< 0\), temos que \(\hat{\beta} = \ln(\bar{y})\) é o EMV do modelo nulo. Logo, \(\hat{\lambda} = \bar{y}\).

# Deviance para o modelo nulo
ybarra      = mean(y)
lambdaNulo  = rep(ybarra, n)
lveroNulo = 0
for(i in 1:n){
  lambdaNulo[i]     = ybarra
  lveroNulo         = lveroNulo + dpois(y[i], lambdaNulo[i], log=T)
}

Deviance_nula = 2*(lveroSaturado - lveroNulo)
Deviance_nula
## [1] 2093.683
  1. Calcule \(\Delta D\), ou seja, a diferença entre a deviance nula e a deviance aplicada no modelo ajustado e interprete esta medida.
delta_D = Deviance_nula - Deviance
delta_D
## [1] 1052.863
alfa  = 0.05 #nível de significancia alfa

##################conferir codigos abaixo
curve(dchisq(x, p-1), from=max(0,p-1-5), to=p-1+5, lwd=2, xlab="",main="Diferença entre as Deviances", bty="n")
segments(x0=Deviance, y0=0, x1=Deviance, 
         y1=dchisq(Deviance, p-1), lty=3, lwd=2, col="red")
segments(x0=n-p, y0=0, x1=p-1, 
         y1=dchisq(p-1, p-1), lty=3, lwd=2, col="purple")
segments(x0=qchisq(1-alfa,p-1), y0=0, x1=qchisq(1-alfa,p-1), 
         y1=dchisq(qchisq(1-alfa,p-1), p-1), lty=3, lwd=2, col="blue")
legend("topright", legend=c("Deviance", "p-1", "qchisq(1-alfa,p-1)"),
       lty=3, lwd=2, col=c("red", "purple", "blue"), bty = "n")
#pintando a regiao critica
u   = seq(qchisq(1-alfa,p-1), p-1+10, length.out=100)
fu  = dchisq(u, p-1)
polygon(c(u, rev(u)), c(fu, rep(0, length(fu))), col="lightblue", border=NA)

#c(Deviance, qchisq(1-alfa,p-1))
if(delta_D < qchisq(1-alfa,p-1)){print("Rejeito a hipotese de bondade do ajuste")}else{print("Não rejeito a hipotese de bondade do ajuste")}
## [1] "Não rejeito a hipotese de bondade do ajuste"
1-pchisq(delta_D, p-1)
## [1] 0
if(1-pchisq(Deviance,p-1)>alfa){print("Rejeito a hipotese de bondade do ajuste")}else{print("Nao rejeito a hipotese de bondade do ajuste")}
## [1] "Nao rejeito a hipotese de bondade do ajuste"

1.4.2 Resíduos

  1. Analise os resíduos.
# Resíduos de Pearson
residuo_pearson = (y - ychapeu)/sqrt(ychapeu)
plot(residuo_pearson,lwd=4,xlab='x',ylab='Resíduos')
abline(h=0,col=2)

#plot(x[,2],residuo_pearson,lwd=4,xlab='x',ylab='Resíduos')
#abline(h=0,col=2)

hist(residuo_pearson,lwd=4,xlab='Resíduos de Pearson', ylab="densidade", freq=F, main="")
curve(dnorm(x), col="red", lwd=2, add=T)

1.4.3 Usando a função glm no R

  1. Faça o ajuste usando a função glm e compare os valores encontrados no ajuste com os calculados anteriormente.
#ajuste_Poi = glm(y~x-1,family=poisson(link="log"))
ajuste_Poi = glm(y~x[,2],family=poisson(link="log"))

#comparando as estimativas pontuais dos coeficientes
cbind(ajuste_Poi$coefficients, betachapeu)                    
##                  [,1]      [,2]
## (Intercept) 0.5327077 0.5327078
## x[, 2]      0.6906265 0.6906266
#comparando as estimativas intervalares dos coeficientes
confint(ajuste_Poi,level=1-alfa)
##                 2.5 %    97.5 %
## (Intercept) 0.4810221 0.5834235
## x[, 2]      0.6491080 0.7321264
rbind(c(betachapeu[1] - zalfa2*sqrt(Inv_Iesp[1,1]), betachapeu[1] + zalfa2*sqrt(Inv_Iesp[1,1])), 
c(betachapeu[2] - zalfa2*sqrt(Inv_Iesp[2,2]), betachapeu[2] + zalfa2*sqrt(Inv_Iesp[2,2])))
##           [,1]      [,2]
## [1,] 0.4815111 0.5839044
## [2,] 0.6491177 0.7321355
# comparando os valores ajustados obtidos pelo glm e calculados anteriormente
#ajuste_Poi$fitted.value     #ou fitted(ajuste_Poi)
plot(ajuste_Poi$fitted.value, ychapeu, bty="n")
lines(ajuste_Poi$fitted.value,ajuste_Poi$fitted.value,col="red",lty=3)

# Resumo da análise: Estimativas pontuais e estatíticas de teste
summary(ajuste_Poi)
## 
## Call:
## glm(formula = y ~ x[, 2], family = poisson(link = "log"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5878  -0.9294  -0.1784   0.5923   3.1284  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.53271    0.02612   20.39   <2e-16 ***
## x[, 2]       0.69063    0.02118   32.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2093.7  on 999  degrees of freedom
## Residual deviance: 1040.8  on 998  degrees of freedom
## AIC: 3193.4
## 
## Number of Fisher Scoring iterations: 5
# Deviance
Dev = deviance(ajuste_Poi) 
c(deviance(ajuste_Poi), ajuste_Poi$deviance, Deviance)
## [1] 1040.819 1040.819 1040.819
#comparando a deviance nula calculada com a ajustada pelo glm
c(ajuste_Poi$null.deviance, Deviance_nula)
## [1] 2093.683 2093.683
# Análise de deviance (comparação de modelos)
anova(ajuste_Poi,test="Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: y
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                     999     2093.7              
## x[, 2]  1   1052.9       998     1040.8 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Se a deviance for significativamente grande, isso sugere que o modelo completo se ajusta significativamente melhor aos dados do que o modelo reduzido, o que pode indicar a importância das variáveis adicionais incluídas no modelo completo. Por outro lado, se a deviance não for significativamente grande, não há evidência para rejeitar a hipótese nula, indicando que o modelo reduzido é adequado.

#Preditor Linear
etachapeu = ajuste_Poi$linear.predictors   #predict(ajuste_Poi)

# Resíduos
residuo = residuals(ajuste_Poi,type="pearson")   #residuals(ajuste_Poi,type="deviance")

#comparando os residuos calculados "manualmente" com os calculados usando a funcao glm
plot(residuo_pearson,residuo,lwd=4,xlab='Residuos calculados manualmente',ylab='Resíduos do GLM')
lines(residuo_pearson,residuo_pearson,col="red", lty=3, lwd=3)

# todas as saídas da função
names(ajuste_Poi)
##  [1] "coefficients"      "residuals"         "fitted.values"    
##  [4] "effects"           "R"                 "rank"             
##  [7] "qr"                "family"            "linear.predictors"
## [10] "deviance"          "aic"               "null.deviance"    
## [13] "iter"              "weights"           "prior.weights"    
## [16] "df.residual"       "df.null"           "y"                
## [19] "converged"         "boundary"          "model"            
## [22] "call"              "formula"           "terms"            
## [25] "data"              "offset"            "control"          
## [28] "method"            "contrasts"         "xlevels"