#leyendo el número de  hombres vivos por edad y año de exposición
ex<-as.matrix(read.csv2("H91_09_EXPO.csv", header=TRUE, sep=";", dec=","))
#leyendo el número de  hombres difuntos  por edad y año de exposición
dx<-as.matrix(read.csv2("H91_09_DEAD.csv", header=TRUE, sep=";", dec=","))
#pesos
wx=ex*0+1


#Agrupando  edad, año, tasa bruta de mortalidad y pesos
qdata.ES<-structure(list(x=as.integer(c(1:101)), y=as.integer(c(1991:2009)),
etx=structure(c(t(ex)), .Dim = as.integer(c(19,101))), dtx=structure(c(t(dx)),
.Dim = as.integer(c(19,101))), wa=structure(c(t(wx)),  .Dim = as.integer(c(19,101))),
.Names = c("x", "y", "etx", "dtx", "wa")))

 #identificando cada  objeto
 x=qdata.ES$x
 y=qdata.ES$y
 etx=qdata.ES$etx
 dtx=qdata.ES$dtx
 wa=qdata.ES$wa



#logarítmo de la  tasa bruta de mortalidad
ll=log(dtx/etx)

#Ajuste modelo cásico para un año
#con datos del INE

#Calculando logarítmo tasas brutas de mortalidad par el año 2009
Obs=log(dtx[19,]/etx[19,])
AGE=x



#Regresión lineal Simple
Gomp.fit <- lm(Obs ~ AGE)
names(Gomp.fit)
Gomp.fit$coef



#Gráfico de las tasas brutas
par(mfrow=c(1,1))
plot(AGE, Obs, xlab = "Age", ylab = "log(mortality)",
     main = "Gompertz law: CMI data, ages 0 to 101, year 2009")
lines(AGE, Gomp.fit$fit, lwd = 2)
text(10, 0, "log(m) = -9.68 + 0.084 * Age", adj = 0, cex = 1.4)
legend(60, -6.5, legend = c("Unweighted fit", "Weighted fit"),
       lty = c(1,2), lwd = 2, cex = 1.4)

          
          
#ANALISIS DE LOS RESIDUOS
plot(Gomp.fit)
residuals(Gomp.fit)



          
          

#Regresiòn con pesos
Gomp.fit.w <- lm(Obs ~ AGE, weights = dtx[19,])
#Gráfico de las tasas brutas
par(mfrow=c(1,1))
plot(AGE, Obs, xlab = "Age", ylab = "log(mortality)",
     main = "Gompertz law: CMI data, ages 0 to 101, year 2009")
lines(AGE, Gomp.fit$fit, lwd = 2)
text(10, 0, "log(m) = -9.68 + 0.084 * Age", adj = 0, cex = 1.4)
legend(60, -6.5, legend = c("Unweighted fit", "Weighted fit"),
       lty = c(1,2), lwd = 2, cex = 1.4)
lines(AGE, Gomp.fit.w$fit, lwd = 2, lty = 2)
text(10, -1, "log(m) = -9.99 + 0.0891 * Age, weighted",
          adj = 0, cex = 1.4)



plot(Gomp.fit.w)
residuals(Gomp.fit.w)



#Gráfico de las tasas brutas
par(mfrow=c(1,1))
plot(AGE, Obs, xlab = "Age", ylab = "log(mortality)",
     main = "Gompertz law: CMI data, ages 0 to 101, year 2009")
lines(AGE, Gomp.fit$fit, lwd = 2)
text(10, 0, "log(m) = -9.68 + 0.084 * Age", adj = 0, cex = 1.4)
legend(60, -6.5, legend = c("Unweighted fit", "Weighted fit"),
       lty = c(1,2), lwd = 2, cex = 1.4)
lines(AGE, Gomp.fit.w$fit, lwd = 2, lty = 2)
text(10, -1, "log(m) = -9.99 + 0.0891 * Age, weighted",
          adj = 0, cex = 1.4)