Métodos de Estimación para un MLG de Poisson

1. Modelo de Poisson

Un Modelo MLG de Poisson se puede expresar de la siguiente manera:

\[y_{i} \sim Poisson(\theta_{i})\] \[log(\theta_{i}) = X_{i}^{'}\beta\] El parámetro \(\theta_{i}\) de la distribución de la variable dependiente \(y_{i}\) se puede enlazar con una combinación lineal de variables predictoras \(X_{i}\) a traves de la función \(log(\theta_{i})\).

El objetivo del modelo es estimar el parámetro \(\theta_{i}\) mediante un vector de parámetros \(\beta\) de la combinación lineal \(X_{i}^{'}\beta\) que minimize el error de estimación. Existen varios métodos de estimación y es lo que se tratará en el presente trabajo.

2. Dataset COVID

Se utiliza la dataset “COVID” que puede ser descargado de Casos positivos por COVID-19 - Ministerio de Salud - MINSA.

El dataset “COVID” contine los siguientes atributos:

FECHA_CORTE : Fecha de corte de información
UUID : ID de la persona como caso positivo de covid-19 confirmada con cualquier tipo de prueba.
UBIGEO : Código de Ubicación Geografica que denotan “DDppdd” (Departamento, provincia,distrito)
DEPARTAMENTO : Departamento donde reside la persona confirmada como caso positivo de covid-19
PROVINCIA : Provincia donde reside la persona confirmada como caso positivo de covid-19
DISTRITO : Distrito donde reside la persona confirmada como caso positivo de covid-19
METODODX : Metodo de laboratorio a la que es sometida una prueba de covid-19
EDAD : Edad de la persona confirmada como caso positivo de covid-19
SEXO : Sexo de la persona confirmada como caso positivo de covid-19
FECHA_RESULTADO: Fecha del resultado de la prueba de covid-19

2.1 Importando el dataset

# ==============================================================================
# Liberias
# ==============================================================================
#install.packages("readr")
#install.packages("dplyr") 
#install.packages("tidyverse") 
#install.packages("Metrics") 
library("tidyverse")
library("readr")
library("dplyr")
library("Amelia")
library("Metrics")

#===============================================================================
# Importando los datos 
#===============================================================================
# Eliminar los objetos de Global Environment
rm(list = ls())

ls_ruta = "D:/Doctorado/DCIES/02.Ciclo/DCIES-241 Estadística aplicada/sesion07/Laboratorio2/Dataset"

setwd(ls_ruta)

#df_data= as.data.frame(read_csv2("positivos_covid.csv"))

df_data = data_frame(read_csv2("positivos_covid.csv"))

2.2 Limpieza y depuración de datos

Filtramos los datos correspondientes a las 8 primeras semanas pa partir del 06/03/2020, fecha en que se reportó el primer caso positivo de COVID en el Perú.

#===============================================================================
# Depuracion y Limpieza de datos
#===============================================================================


# Ocho primeras semanas
ls_fecha_ini    = "20200306"
ls_fecha_fin    = "20200502"


df_peru_corte = filter(df_data,FECHA_RESULTADO >=ls_fecha_ini & 
                                  FECHA_RESULTADO <= ls_fecha_fin )

El estudio solo considera a los distritos de Lima Norte: Ancón, Carabayllo, Comas, Independencia, Los Olivos, Puente Piedra, San Martin de Porras y Santa Rosa.

# Lima Norte: 150102 - Ancón, 150106 - Carabayllo, 
#             150110 - Comas, 150112 - Independencia, 
#             150117 - Los Olivos, 150125- Puente Piedra, 
#             150135 - San Martín de Porres y 150139- Santa Rosa.

ls_lima_norte = c('150102','150106','150110','150112',
                  '150117','150125','150135','150139')

df_lima_norte =
          filter(df_peru_corte,UBIGEO %in% ls_lima_norte)

La variable EDAD la transformanos en la variable GRUPO correspondiente a los siguientes grupos étnarios:

Grupos Etnarios

0 Primera Infancia (0-5 años)
1 Infancia (6 - 11 años)
2 Adolescencia (12 - 18 años)
3 Juventud (18 - 26 años)
4 Adultez (27- 59 años)
5 Persona Mayor (60 años o mas) envejecimiento y vejez

df_lima_norte =
          filter(df_peru_corte,UBIGEO %in% ls_lima_norte)

df_lima_norte = mutate(df_lima_norte, GRUPO = case_when(
                                        EDAD>=  0 & EDAD <=  5 ~ 0,
                                        EDAD>=  6 & EDAD <= 11 ~ 1,
                                        EDAD>= 12 & EDAD <= 18 ~ 2,
                                        EDAD>= 18 & EDAD <= 26 ~ 3,
                                        EDAD>= 27 & EDAD <= 59 ~ 4,
                                        EDAD>= 60 ~ 5))
 
# Factor SEXO
df_lima_norte$SEXO = model.matrix(~df_lima_norte$SEXO)[,2]


# Seleccionado: FECHA_RESULTADO,SEXO, GRUPO  
df_lima_norte = select(df_lima_norte,FECHA_RESULTADO, GRUPO, SEXO)

# Ordenando: FECHA_RESULTADO,SEXO, EDAD  
df_lima_norte = arrange(df_lima_norte,FECHA_RESULTADO,GRUPO, SEXO )

## count(df_lima_norte) 6829

# Determinado los casos por grupos etnarios
df_lima_norte_casos = group_by(df_lima_norte,FECHA_RESULTADO,GRUPO, SEXO) %>% 
                                    summarize(CASOS = n())

Por otro lado la variable FECHA_RESULTADO la transformamos en la variable DIA que correponde a la cantidad de días transcurridos desde que se reportó el primer caso COVID en Lima Norte.

# Creando dataset temporal para calcular el día.

ls_fecha = sort(unique(df_lima_norte_casos$FECHA_RESULTADO))

ls_dia   = c()

for ( i in 1:length(ls_fecha)){
  ls_dia = c(ls_dia,i)
}

df_dia_resultado =  data_frame(FECHA_RESULTADO = ls_fecha,
                               DIA             = ls_dia)


# Combinando df_lima_norte_casos y df_dia_resultado
df_lima_norte_casos = merge(x = df_lima_norte_casos, 
                            y = df_dia_resultado, 
                            by = "FECHA_RESULTADO", 
                            all = TRUE)


df_lima_norte_dia_casos = select(df_lima_norte_casos,CASOS, DIA, GRUPO, SEXO)

# sum(df_lima_norte_casos$CASOS) 6829

head(df_lima_norte_dia_casos,20)

##    CASOS DIA GRUPO SEXO
## 1      1   1     0    1
## 2      1   1     2    1
## 3      1   1     4    0
## 4      1   1     4    1
## 5      2   2     3    0
## 6      1   2     4    0
## 7      1   2     4    1
## 8      2   3     4    0
## 9      1   4     3    0
## 10     1   4     4    0
## 11     1   4     4    1
## 12     1   4     5    0
## 13     1   5     4    0
## 14     2   5     5    1
## 15     1   6     3    0
## 16     3   6     4    0
## 17     1   6     4    1
## 18     3   7     4    0
## 19     2   7     4    1
## 20     1   8     3    1

2.3 Analisis Descriptivo

A continuación se muestra algunos gráficos de la evolución de los CASOS COVID en Lima Norte en las primeras 8 semanas.

#===============================================================================
# Análisis Descriptivo
#===============================================================================

df_casos_dia = group_by(df_lima_norte_dia_casos,DIA) %>% 
                               summarize(CASOS = sum(CASOS))

# sum(df_casos_dia$CASOS) 6829

ggplot( data=df_casos_dia, aes(x = DIA, y= CASOS), main = ) +
                ggtitle("Casos totales Lima") + 
                geom_point()  + 
                geom_smooth(method = 'glm',  method.args = list(family = 'poisson'))

df_casos_dia_sexo = group_by(df_lima_norte_dia_casos,DIA,SEXO) %>% 
                    summarize(CASOS = sum(CASOS))
# sum(df_casos_dia_sexo$CASOS) 6829


ggplot( data=df_casos_dia_sexo, aes(x = DIA, y= CASOS)) +
  ggtitle("Casos totales Lima por sexo") + 
  geom_point(aes(colour = factor(SEXO)))+ 
  geom_smooth(method = 'glm',  method.args = list(family = 'poisson'))

df_casos_dia_grupo = group_by(df_lima_norte_dia_casos,DIA,GRUPO) %>% 
  summarize(CASOS = sum(CASOS))

# sum(df_casos_dia_grupo$CASOS) 6829


ggplot( data=df_casos_dia_grupo, aes(x = DIA, y= CASOS) ) +
  ggtitle("Casos totales Lima por grupo etnario") + 
  geom_point(aes(colour = factor(GRUPO)))+ 
  geom_smooth(method = 'glm',  method.args = list(family = 'poisson'))

# 4 Grupo Adulto (27- 59 años)

df_casos_dia_adulto_sexo =
  filter(df_lima_norte_dia_casos, GRUPO == 4)


ggplot( data=df_casos_dia_adulto_sexo, aes(x = DIA, y= CASOS)  ) +
  ggtitle("Casos totales Lima por grupo adulto") + 
  geom_point(aes(colour = factor(SEXO)))+ 
  geom_smooth(method = 'glm',  method.args = list(family = 'poisson'))

3 Selección de modelo

Se propone y estima 2 modelos para el grupo étnario Adultos de Lima Norte. Los datos entrenamiento se encuentran en el dayaset df_casos_dia_adulto_sexo_train y los datos de prueba en el dataset df_casos_dia_adulto_sexo_test.

Consideraremos 80% de los datos para entrenamiento y 20% restante de prueba para evaluar los modelos propuestos. Los datos de entrenamiento y pruebas se selecionan aleatoriamente.

#===============================================================================
# Modelo GLM   - Poisson
#===============================================================================

split = 0.80

set.seed(4321)

Evaluación del modelo 1: CASOS ~ DIA + SEXO

# Modelo 1 : CASOS ~ DIA + SEXO
df_casos_dia_adulto_sexo =
  filter(df_lima_norte_dia_casos, GRUPO == 4)


train = sort(sample(nrow(df_casos_dia_adulto_sexo),floor(nrow(df_casos_dia_adulto_sexo)*split)))

df_casos_dia_adulto_sexo_train = df_casos_dia_adulto_sexo[train,]
df_casos_dia_adulto_sexo_test  = df_casos_dia_adulto_sexo[-train,]

glm_modelo1 = glm(CASOS ~ DIA + SEXO, df_casos_dia_adulto_sexo_train, family = poisson(link = "log"))
summary(glm_modelo1)

## 
## Call:
## glm(formula = CASOS ~ DIA + SEXO, family = poisson(link = "log"), 
##     data = df_casos_dia_adulto_sexo_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.5075  -2.1587  -0.8635   1.0477   6.8402  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.86465    0.07369   11.73   <2e-16 ***
## DIA          0.07604    0.00171   44.47   <2e-16 ***
## SEXO         0.68308    0.03748   18.22   <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: 3306.60  on 73  degrees of freedom
## Residual deviance:  420.55  on 71  degrees of freedom
## AIC: 773.59
## 
## Number of Fisher Scoring iterations: 5

df_casosp = as.data.frame(df_casos_dia_adulto_sexo_test[,c("DIA","SEXO")])

df_casos_dia_adulto_sexo_test$CASOSP= round(predict(glm_modelo1, df_casosp,   type = "response"),0)

rmse(df_casos_dia_adulto_sexo_test$CASOS,df_casos_dia_adulto_sexo_test$CASOSP)

## [1] 16.16526

Evaluación del modelo 2: CASOS ~ DIA

# Modelo 2 : CASOS ~ DIA
df_casos_dia_adulto = group_by(df_casos_dia_adulto_sexo,DIA) %>% 
  summarize(CASOS = sum(CASOS))

train = sort(sample(nrow(df_casos_dia_adulto),floor(nrow(df_casos_dia_adulto)*split)))

df_casos_dia_adulto_train      = df_casos_dia_adulto[train,]
df_casos_dia_adulto_test       = df_casos_dia_adulto[-train,] 

glm_modelo2 = glm(CASOS ~ DIA , df_casos_dia_adulto_train, family = poisson(link = "log"))
summary(glm_modelo2)

## 
## Call:
## glm(formula = CASOS ~ DIA, family = poisson(link = "log"), data = df_casos_dia_adulto_train)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -5.744  -2.905  -1.659   1.439   6.921  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.957009   0.062250   31.44   <2e-16 ***
## DIA         0.076011   0.001611   47.17   <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: 3058.85  on 38  degrees of freedom
## Residual deviance:  422.85  on 37  degrees of freedom
## AIC: 629.64
## 
## Number of Fisher Scoring iterations: 5

df_casosp = as.data.frame(df_casos_dia_adulto_test[,c("DIA")])

df_casos_dia_adulto_test$CASOSP= round(predict(glm_modelo2, df_casosp,   type = "response"),0)

rmse(df_casos_dia_adulto_test$CASOS,df_casos_dia_adulto_test$CASOSP)

## [1] 29.18219

3.3 Selección del modelo

Se seleciona el modelo 1: CASOS ~ DIA + SEXO por los siguientes criterios:

Variables significativas
Menor MASE
Parsimonia

El siguente cuadro muestra el vector de parámetros \(\beta\) y el indicador MASE:

glm_modelo1 = glm(CASOS ~ DIA + SEXO, df_casos_dia_adulto_sexo_train, family = poisson(link = "log"))
summary(glm_modelo1)

## 
## Call:
## glm(formula = CASOS ~ DIA + SEXO, family = poisson(link = "log"), 
##     data = df_casos_dia_adulto_sexo_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.5075  -2.1587  -0.8635   1.0477   6.8402  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.86465    0.07369   11.73   <2e-16 ***
## DIA          0.07604    0.00171   44.47   <2e-16 ***
## SEXO         0.68308    0.03748   18.22   <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: 3306.60  on 73  degrees of freedom
## Residual deviance:  420.55  on 71  degrees of freedom
## AIC: 773.59
## 
## Number of Fisher Scoring iterations: 5

df_casosp = as.data.frame(df_casos_dia_adulto_sexo_test[,c("DIA","SEXO")])

df_casos_dia_adulto_sexo_test$CASOSP= round(predict(glm_modelo1, df_casosp,   type = "response"),0)

rmse(df_casos_dia_adulto_sexo_test$CASOS,df_casos_dia_adulto_sexo_test$CASOSP)

## [1] 16.16526

4 Métodos de estimación

Se procederá a estimar el vector de parámetros \(\beta\) del modelo 1: CASOS ~ DIA + SEXO por los siguientes métodos:

Estimación utilizando la función R glm
Estimación por Maxima Verosimilitud utilizado la función Optim de R
Estimación por el Método Scoring-Fisher
Estimación por el Método Mínimos Cuadros Reponderados Iterativamente - IRLS
Estimación utilizando el algoritmo de Metropolis-Hasting

Finalmente se presentará un cuadro comparativo de los 5 métodos.

Preparado el cuadro comparativo y la matriz de datos:

#===============================================================================
# Modelo Seleccionado
#===============================================================================
# Modelo 1 : CASOS ~ DIA + SEXO
glm_modelo1 = glm(CASOS ~ DIA + SEXO, df_casos_dia_adulto_sexo_train, family = poisson(link = "log"))
summary(glm_modelo1)

## 
## Call:
## glm(formula = CASOS ~ DIA + SEXO, family = poisson(link = "log"), 
##     data = df_casos_dia_adulto_sexo_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.5075  -2.1587  -0.8635   1.0477   6.8402  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.86465    0.07369   11.73   <2e-16 ***
## DIA          0.07604    0.00171   44.47   <2e-16 ***
## SEXO         0.68308    0.03748   18.22   <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: 3306.60  on 73  degrees of freedom
## Residual deviance:  420.55  on 71  degrees of freedom
## AIC: 773.59
## 
## Number of Fisher Scoring iterations: 5

# ==============================================================================
# Cuadro Comparativo y Matriz de Datos
# ==============================================================================

df_cuadro_comparacion = data.frame(Variable = c("(Intercept)","DIA","SEXO"),
                                   "R.glm"               = c(0,0,0),
                                   "R.optim"             = c(0,0,0),
                                   "Scoring.Fisher"      = c(0,0,0),
                                   "IRLS"                = c(0,0,0),
                                   "Metropolis.Hasting"  = c(0,0,0))


# Preparando las matriz de datos
df_casos_dia_adulto_sexo_train$INTERCEPTO = 1

Y = as.matrix(df_casos_dia_adulto_sexo_train$CASOS)

XX = as.matrix(df_casos_dia_adulto_sexo_train[,c("INTERCEPTO","DIA","SEXO")])

colnames(XX) = NULL

4.1 Estimación utilizando la función R glm

La función glm de R estima el vector de parámetros \(\beta\)

# ==============================================================================
# Método 1:  Regresión utilizando librería glm de R
# ==============================================================================

glm_modelo1 = glm(CASOS ~ DIA + SEXO, df_casos_dia_adulto_sexo_train, family = poisson(link = "log"))
summary(glm_modelo1)

## 
## Call:
## glm(formula = CASOS ~ DIA + SEXO, family = poisson(link = "log"), 
##     data = df_casos_dia_adulto_sexo_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.5075  -2.1587  -0.8635   1.0477   6.8402  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.86465    0.07369   11.73   <2e-16 ***
## DIA          0.07604    0.00171   44.47   <2e-16 ***
## SEXO         0.68308    0.03748   18.22   <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: 3306.60  on 73  degrees of freedom
## Residual deviance:  420.55  on 71  degrees of freedom
## AIC: 773.59
## 
## Number of Fisher Scoring iterations: 5

b1 = as.numeric(coefficients(glm_modelo1)[1])
b2 = as.numeric(coefficients(glm_modelo1)[2]) 
b3 = as.numeric(coefficients(glm_modelo1)[3]) 

df_cuadro_comparacion[,"R.glm"] = c(b1,b2,b3)

print(df_cuadro_comparacion)

##      Variable      R.glm R.optim Scoring.Fisher IRLS Metropolis.Hasting
## 1 (Intercept) 0.86464843       0              0    0                  0
## 2         DIA 0.07603686       0              0    0                  0
## 3        SEXO 0.68307873       0              0    0                  0

4.2 Estimación por Maxima Verosimilitud utilizado la función Optim de R

La función Optim de R optimiza (maximiza) la función log verosimilitud del modelo, el resultado es un vector de parámetros \(\beta\)

# ==============================================================================
# Método 2:  Optimización Verosimilitud con función Optim de R
# ==============================================================================
#------------------------------------------------------------------------------
# Funcion de LogVerosimilitud
#------------------------------------------------------------------------------
fLogVerosimilitud = function (B,Y,XX) {
  nrow = nrow(Y)
  prow = nrow(B)
  
  lvs = 0.0
  
  for (i in 1:nrow) {
    Yi = as.numeric(Y[i,1])
    Xi = as.matrix(XX[i,],prow,1)
    
    
    lvs = lvs + Yi*t(Xi)%*%B - exp(t(Xi)%*%B)
    
  }
  
  return(lvs)
}


#------------------------------------------------------------------------------
# Funcion de LogVerosimilitud
#------------------------------------------------------------------------------
fNegativoLogVerosimilitud= function(B,Y,XX) {    
  
  return( - fLogVerosimilitud (B,Y,XX)) 
  
}


Bo = matrix(0, nrow = ncol(XX), ncol = 1)


r = optim( Bo , fNegativoLogVerosimilitud, Y = Y,XX = XX)

b1 = r$par[1]
b2 = r$par[2]
b3 = r$par[3]

df_cuadro_comparacion[,"R.optim"] = c(b1,b2,b3)

print(df_cuadro_comparacion)

##      Variable      R.glm    R.optim Scoring.Fisher IRLS Metropolis.Hasting
## 1 (Intercept) 0.86464843 0.86460627              0    0                  0
## 2         DIA 0.07603686 0.07603777              0    0                  0
## 3        SEXO 0.68307873 0.68308009              0    0                  0

4.3 Estimación por el Método Scoring-Fisher

Se utiliza el algoritmo Scoring-Fisher, se requiere conocer la función log verosimilitud, la Gradiente y el Hessiano. El cuadro de resultados muestra como el vector \(\beta\) itera hasta llegar al optimo de la función Log Verosimilitud.

Se utiliza la norma2 como criterio para evaluar la convergencia del vector que optimiza la función Log Verosimilitud.

# ==============================================================================
# Método 3:  Scoring-Fisher 
# ==============================================================================

#------------------------------------------------------------------------------
# Funcion calcula la norma de la diferencia de vectores
#------------------------------------------------------------------------------
fNorma = function (A,B){
  
  return (sqrt(as.numeric(t(A- B)%*%(A- B))) )
}

#------------------------------------------------------------------------------
# Funcion de Gradiente LogVerosimilitud
#------------------------------------------------------------------------------
fGLogVerosimilitud = function (B,Y,XX) {
  nrow = nrow(Y)
  prow = nrow(B)
  
  G = matrix(0, nrow = prow, ncol = 1)
  
  for (i in 1:nrow) {
    Yi = as.numeric(Y[i,1])
    Xi = as.matrix(XX[i,],prow,1)
    
    K  =  as.numeric( exp(t(Xi)%*%B))  
    
    G =  G  +   Yi*Xi - K*Xi
  }
  
  return(G)
}


#------------------------------------------------------------------------------
# Funcion de Hessiano LogVerosimilitud
#------------------------------------------------------------------------------
fHLogVerosimilitud = function (B,Y,XX) {
  nrow = nrow(Y)
  prow = nrow(B)
  
  H = matrix(0, nrow = prow, ncol = prow)
  
  for (i in 1:nrow) {
    Yi = as.numeric(Y[i,1])
    Xi = as.matrix(XX[i,],prow,1)
    
    K  =  as.numeric(exp(t(Xi)%*%B))
    
    H = H + K*Xi%*%t(Xi) 
  }
  
  return(-H)
}


#------------------------------------------------------------------------------
# Funcion Informacion de Fisher
#------------------------------------------------------------------------------
fInformacionFisher = function (B,Y,XX) {
  
  return(-fHLogVerosimilitud(B,Y,XX))
}




#------------------------------------------------------------------------------
# Método :  Scoring-Fisher
#------------------------------------------------------------------------------

B0 = matrix(0, nrow = ncol(XX), ncol = 1)


B1 = B0

epsilon = 0.000000001

i = 1


df_resultado=   
  
  data.frame (Iteracion        = i,
              Estimacion       = paste(paste("(" ,toString(round(B1,6) )),")"),
              Control          = c(''),
              LogVerosimilitud = fLogVerosimilitud (B1,Y,XX),
              Gradiente        = paste(paste("(" ,toString(round(fGLogVerosimilitud (B1,Y,XX),3) )),")"))


while (( fNorma ( B1,B0) >= epsilon) | (i == 1)){
  
  
  B0 = B1
  
  i = i + 1
  
  
  B1 = B0 + solve(fInformacionFisher (B0,Y,XX))%*%fGLogVerosimilitud (B0,Y,XX)
  
  df_iteracion = 
    data.frame(Iteracion        = i,
               Estimacion       = paste(paste("(" ,toString(round(B1,6) )),")"),
               Control          = fNorma ( B1,B0),
               LogVerosimilitud = fLogVerosimilitud (B1,Y,XX),
               Gradiente        = paste(paste("(" ,toString(round(fGLogVerosimilitud (B1,Y,XX),2) )),")"))
  
  
  df_resultado  = rbind(df_resultado, df_iteracion)
}

# Ultimas 25 iteraciones
tail(df_resultado,25)

##     Iteracion                           Estimacion              Control
## 88         88 ( -118.987087, 2.518457, 21.132936 )    0.994533105876269
## 89         89 ( -118.986656, 2.518448, 20.132943 )    0.999993208411725
## 90         90  ( -118.985445, 2.518424, 19.13292 )     1.00002323508894
## 91         91 ( -118.982151, 2.518358, 18.132859 )     1.00006656461841
## 92         92   ( -118.9732, 2.518179, 17.132693 )     1.00020623874744
## 93         93 ( -118.948878, 2.517691, 16.132241 )     1.00074724250769
## 94         94 ( -118.882834, 2.516368, 15.131015 )     1.00340275640764
## 95         95 ( -118.703825, 2.512782, 14.127693 )     1.01917275561766
## 96         96 ( -118.221033, 2.503109, 13.118734 )     1.11856137647731
## 97         97  ( -116.936161, 2.477367, 12.09491 )      1.6430985789615
## 98         98 ( -113.635165, 2.411227, 11.033832 )     3.46797366385732
## 99         99   ( -105.877905, 2.255777, 9.89111 )     7.84251590050251
## 100       100    ( -90.958001, 1.95667, 8.620722 )     14.9768779297277
## 101       101     ( -70.484442, 1.54578, 7.26312 )      20.522635257448
## 102       102   ( -49.845424, 1.130581, 5.924327 )     20.6865610151056
## 103       103    ( -31.085885, 0.75161, 4.627462 )     18.8081306714128
## 104       104   ( -14.719535, 0.418088, 3.368085 )     16.4181204608888
## 105       105    ( -3.528232, 0.183681, 2.291028 )     11.2454549610108
## 106       106    ( -1.207503, 0.127739, 1.620036 )     2.41643246279755
## 107       107     ( 0.290037, 0.091435, 1.033028 )     1.60888822555317
## 108       108      ( 0.78242, 0.078345, 0.745777 )     0.57019788754276
## 109       109      ( 0.862081, 0.07611, 0.685388 )   0.0999879411779046
## 110       110     ( 0.864645, 0.076037, 0.683082 )  0.00344991517714593
## 111       111     ( 0.864648, 0.076037, 0.683079 ) 4.53859435858739e-06
## 112       112     ( 0.864648, 0.076037, 0.683079 ) 8.09970279354496e-12
##     LogVerosimilitud                                                Gradiente
## 88     -1.257163e+11 ( -125716228433.56, -6158254562022.8, -125716228433.54 )
## 89     -4.624845e+10  ( -46248415755.55, -2265495250840.91, -46248415755.55 )
## 90     -1.701388e+10   ( -17013840860.76, -833429112424.77, -17013840860.76 )
## 91     -6.259080e+09     ( -6259041736.29, -306601419288.42, -6259041736.29 )
## 92     -2.302613e+09     ( -2302572250.98, -112792342261.28, -2302572250.98 )
## 93     -8.471107e+08        ( -847068467.69, -41493967305.47, -847068467.69 )
## 94     -3.116628e+08        ( -311618554.26, -15264761202.62, -311618554.26 )
## 95     -1.146836e+08         ( -114637553.88, -5615576200.24, -114637553.88 )
## 96     -4.222007e+07            ( -42172332.92, -2065841256.4, -42172332.92 )
## 97     -1.556282e+07            ( -15513973.58, -759971692.55, -15513973.55 )
## 98     -5.755728e+06              ( -5707188.19, -279582081.77, -5707187.95 )
## 99     -2.145279e+06              ( -2100149.16, -102888204.16, -2100147.68 )
## 100    -8.112344e+05                 ( -774505.38, -37947337.44, -774498.48 )
## 101    -3.126200e+05                 ( -288136.78, -14117436.76, -288117.22 )
## 102    -1.215783e+05                  ( -109381.75, -5357527.59, -109346.42 )
## 103    -4.466922e+04                    ( -43292.26, -2117477.62, -43227.98 )
## 104    -1.115517e+04                     ( -18616.65, -904126.55, -18462.37 )
## 105     3.303972e+03                       ( -9014.33, -417198.07, -8694.15 )
## 106     8.434195e+03                        ( -3122.81, -143170.07, -3028.8 )
## 107     9.803155e+03                          ( -930.04, -41151.89, -879.61 )
## 108     9.987263e+03                           ( -149.28, -6427.12, -139.17 )
## 109     9.992830e+03                                ( -5.31, -224.71, -4.97 )
## 110     9.992837e+03                                   ( -0.01, -0.3, -0.01 )
## 111     9.992837e+03                                              ( 0, 0, 0 )
## 112     9.992837e+03                                              ( 0, 0, 0 )

b1 = B1[1]
b2 = B1[2]
b3 = B1[3]

df_cuadro_comparacion[,"Scoring.Fisher"] = c(b1,b2,b3)


print(df_cuadro_comparacion)

##      Variable      R.glm    R.optim Scoring.Fisher IRLS Metropolis.Hasting
## 1 (Intercept) 0.86464843 0.86460627     0.86464843    0                  0
## 2         DIA 0.07603686 0.07603777     0.07603686    0                  0
## 3        SEXO 0.68307873 0.68308009     0.68307873    0                  0

Nota: se muestran las últimas 25 iteraciones para donde ser observa la convergencia del vector \(\beta\).

4.4 Estimación por el Método Mínimos Cuadros Reponderados Iterativamente - IRLS

Se utiliza el algoritmo IRLS que se apoya en el método de Mínimo Cuadrados Poderados. Es un método iterativo que busca la convergencia del vector \(\beta\) mediante el ajuste de la variable de respuesta en una nueva variable \(Z\) y el uso de uma matriz de ponderación \(W\).

# ==============================================================================
# Método 4:  IRLS 
# ==============================================================================
fRespuestaAjustadaZ = function(B,Y,XX){
  nrow = nrow(Y)
  prow = nrow(B)
  
  Z = matrix(0, nrow = nrow, ncol = 1)
  
  for (i in 1:nrow) {
    Yi = as.numeric(Y[i,1])
    Xi = as.matrix(XX[i,],prow,1)
    
    Ui = as.numeric((exp(t(Xi)%*%B)))
    
    dgi = 1/Ui
    
    Z[i,1] = t(Xi)%*%B + dgi*(Yi  - Ui)
  }
  return(Z)
  
}


fMatrizPonderacionW = function(B,Y,XX){
  
  nrow = nrow(Y)
  prow = nrow(B)
  
  W = matrix(0, nrow = nrow, ncol = nrow )
  
  diagonal =c()
  
  for (i in 1:nrow) {
    Yi = as.numeric(Y[i,1])
    Xi = as.matrix(XX[i,],prow,1)
    
    Ui = as.numeric((exp(t(Xi)%*%B)))
    
    dgi = 1/Ui
    
    W[i,i] = 1/dgi
    
  }
  
  return(W)
  
}

# Vector inicial Bk = (0,0,0,0)'
B0 = matrix(0, nrow = ncol(XX), ncol = 1)

B1 = B0


epsilon = 0.000000001

i = 1

df_resultado=   
  
  data.frame (Iteracion        = i,
              Estimacion       = paste(paste("(" ,toString(round(B1,6) )),")"),
              Control          = c(''))


while (( fNorma ( B1,B0) >= epsilon) | (i == 1)){
  
  B0 = B1
  
  i = i + 1
  
  
  # Ajustando la respuesta
  Z1 = fRespuestaAjustadaZ(B0,Y,XX)
  
  # Obteniendo la matriz de ponderaciones
  W = fMatrizPonderacionW(B0,Y,XX)
  
  #Mínimos Cuadrados ponderados
  
  B1 = (solve((t(XX)%*%W)%*%XX))%*%((t(XX)%*%W)%*%Z1)
  
  df_iteracion = 
    data.frame(Iteracion        = i,
               Estimacion       = paste(paste("(" ,toString(round(B1,6) )),")"),
               Control          = fNorma ( B1,B0))
  
  
  df_resultado  = rbind(df_resultado, df_iteracion)
}

# Ultimas 25 iteraciones
tail(df_resultado,25)

##     Iteracion                           Estimacion              Control
## 88         88 ( -118.987086, 2.518457, 21.133173 )    0.994572362479622
## 89         89 ( -118.986655, 2.518448, 20.133179 )     0.99999437110685
## 90         90 ( -118.985445, 2.518424, 19.133158 )      1.0000213547841
## 91         91 ( -118.982153, 2.518358, 18.133097 )     1.00006675646682
## 92         92 ( -118.973204, 2.518179, 17.132931 )     1.00020600366013
## 93         93 ( -118.948887, 2.517692, 16.132479 )     1.00074714793928
## 94         94 ( -118.882858, 2.516369, 15.131254 )     1.00340142891027
## 95         95 ( -118.703893, 2.512783, 14.127932 )     1.01916453470187
## 96         96 ( -118.221214, 2.503113, 13.118975 )     1.11851051814157
## 97         97 ( -116.936637, 2.477376, 12.095157 )     1.64286406251313
## 98         98 ( -113.636357, 2.411251, 11.034092 )     3.46728775382232
## 99         99   ( -105.880533, 2.25583, 9.891396 )     7.84109108186723
## 100       100   ( -90.962415, 1.956759, 8.621039 )     14.9750954988769
## 101       101     ( -70.4895, 1.545882, 7.263443 )     20.5219917429712
## 102       102     ( -49.85013, 1.130676, 5.92464 )     20.6869136753883
## 103       103    ( -31.09009, 0.751695, 4.627767 )     18.8086309815324
## 104       104     ( -14.72305, 0.41816, 3.368375 )     16.4188100391537
## 105       105    ( -3.529595, 0.183711, 2.291228 )     11.2476061537707
## 106       106     ( -1.208005, 0.127751, 1.62019 )      2.4172717623774
## 107       107       ( 0.289834, 0.09144, 1.03313 )     1.60918550084859
## 108       108      ( 0.782371, 0.078347, 0.74581 )    0.570366117405952
## 109       109      ( 0.862078, 0.07611, 0.685391 )    0.100043386572328
## 110       110     ( 0.864645, 0.076037, 0.683082 )  0.00345364120766795
## 111       111     ( 0.864648, 0.076037, 0.683079 ) 4.54833873707845e-06
## 112       112     ( 0.864648, 0.076037, 0.683079 ) 8.14981371071344e-12

b1 = B1[1]
b2 = B1[2]
b3 = B1[3]
 

df_cuadro_comparacion[,"IRLS"] = c(b1,b2,b3)

print(df_cuadro_comparacion)

##      Variable      R.glm    R.optim Scoring.Fisher       IRLS
## 1 (Intercept) 0.86464843 0.86460627     0.86464843 0.86464843
## 2         DIA 0.07603686 0.07603777     0.07603686 0.07603686
## 3        SEXO 0.68307873 0.68308009     0.68307873 0.68307873
##   Metropolis.Hasting
## 1                  0
## 2                  0
## 3                  0

Nota: se muestran las últimas 25 iteraciones para donde ser observa la convergencia del vector \(\beta\).

4.5 Estimación utilizando el algoritmo de Metropolis-Hasting

Este método utiliza una Cadena de Markov para estimar el vector de parámetros \(\beta\) a través de muestra cuando luego de lograr la convergencia de la distribución aposteriori.

# ==============================================================================
# Método 5:  Metropolis-Hasting 
# ==============================================================================


#------------------------------------------------------------------------------
# Funcion log de la Distribución Posteriori
#------------------------------------------------------------------------------
flogDistribucionPosteriori = function(Y,XX,B,sigma){
  nrow = nrow(Y)
  pcol = nrow(B)
  
  
  sumaB = t(B)%*%B
  
  
  logapriori = -sumaB/(2*(sigma^2))
  
  
  logVerosimilitud = 0
  
  for (i in 1:nrow){
    Yi = as.numeric(Y[i,1])
    Xi = as.matrix(XX[i,],prow,1)
    
    xb = t(Xi)%*%B
    
    logVerosimilitud = logVerosimilitud + (Yi*xb) - exp(xb)
  }
  
  return(logVerosimilitud + logapriori)
  
}


# tamaño de muestra
M             = 400000
K             = max(c(0.05*M,1000))
S             = max(c(0.1*K,500))

# Hiper parámetros
sigma1         = 4
sigma2         = 0.01
SIGMA2         =  matrix(data = 0, nrow = ncol(XX), ncol = ncol(XX))
diag(SIGMA2)   = sigma2


#------------------------------------------------------------------------------
# Metropolis Hasting
#------------------------------------------------------------------------------

BMetropolis   =  matrix(data = 0, nrow = M, ncol = ncol(XX))

salto = 0

# 1. Valor inicial
Bx = matrix(0,1,ncol = ncol(XX))

cat(sprintf("%.0f",M), " iteraciones a procesar: \n", sep = "" )

## 400000 iteraciones a procesar:

for( i in 1:M){
  #2. Propuesta
  By = mvtnorm::rmvnorm(n = 1, mean = Bx, sigma = SIGMA2)
  
  #3. Tasa de Aceptación
  logPIy = flogDistribucionPosteriori(Y,XX, t(By),sigma1) 
  
  logPIx = flogDistribucionPosteriori(Y,XX, t(Bx),sigma1) 
  
  logalphaxy = min(c( logPIy - logPIx,0))
  
  alphaxy = exp(logalphaxy)
  
  if (is.nan(alphaxy)) { alphaxy = 0.0}
  
  #4. Aceptacion
  if (runif(1) < alphaxy){
    Bx = By
    salto = salto + 1
  }
  
  BMetropolis[i, ] = Bx
  
  #5. Progreso
  
  if (i%%floor(M/10) == 0){
    
    cat(100*round(i/M, 1), "% completado ... \n", sep = "" )
    
  }
}

## 10% completado ... 
## 20% completado ... 
## 30% completado ... 
## 40% completado ... 
## 50% completado ... 
## 60% completado ... 
## 70% completado ... 
## 80% completado ... 
## 90% completado ... 
## 100% completado ...

#Quemado: Descartando las (M - K) iteraciones
UltimosK = BMetropolis[(M-K + 1) : M,]

#Seleccionando MAS de los ultimo K vectores
muestra = UltimosK[sample(1:K,S),]


#------------------------------------------------------------------------------
# Diagnostico de convergencia de parámetros
#------------------------------------------------------------------------------
fDiagnosticoConvergencia = function (pMetropolis,pmuestra){
  
  M = nrow(pMetropolis)
  
  m = nrow(pmuestra)
  
  lmedia = as.numeric(mean(pmuestra))
  
  par(mfrow=c(2,2))
  
  plot( pMetropolis,type="l",xlab = sprintf("Iteraciones:  %.0f",M), ylab = expression(beta),
        main = "Convergencia de parámetro")
  
  
  plot( pmuestra,type="l",xlab = sprintf("Muestra:  %.0f",m), ylab = expression(beta),
        main = "Serie de la muestra")
  
  acf(pmuestra,main = "Autocorrelación de la muestra")
  
  hist(pmuestra,
       prob = TRUE,
       main = "Distribución de parámetro",
       xlab = expression(beta),
       ylab ="Densidad",
       breaks = 100)
  abline(v=lmedia, col='red', lwd=3)
  
  par(mfrow=c(1,1))
  
  return (lmedia)
  
}



b1 = fDiagnosticoConvergencia (BMetropolis[,1],muestra[,1])

b2 = fDiagnosticoConvergencia (BMetropolis[,2],muestra[,2])

b3 = fDiagnosticoConvergencia (BMetropolis[,3],muestra[,3])

df_cuadro_comparacion[,"Metropolis.Hasting"] = c(b1,b2,b3)

print(df_cuadro_comparacion)

##      Variable      R.glm    R.optim Scoring.Fisher       IRLS
## 1 (Intercept) 0.86464843 0.86460627     0.86464843 0.86464843
## 2         DIA 0.07603686 0.07603777     0.07603686 0.07603686
## 3        SEXO 0.68307873 0.68308009     0.68307873 0.68307873
##   Metropolis.Hasting
## 1         0.84106620
## 2         0.07669809
## 3         0.67159181

5. Cuadro Comparativo

El siguiente cuadro muestra los resultados obtenidos por los 5 métodos. Se aprecia que los resultados son similares.

# ==============================================================================
# Comparación de resultados
# ==============================================================================
df_cuadro_comparacion

##      Variable      R.glm    R.optim Scoring.Fisher       IRLS
## 1 (Intercept) 0.86464843 0.86460627     0.86464843 0.86464843
## 2         DIA 0.07603686 0.07603777     0.07603686 0.07603686
## 3        SEXO 0.68307873 0.68308009     0.68307873 0.68307873
##   Metropolis.Hasting
## 1         0.84106620
## 2         0.07669809
## 3         0.67159181