LC for residuals

残差热度图

三种方法的残差

导入数据

#导入美国数据
USA.mortality<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡率1950-2017.csv",header=TRUE)
USA.population<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国人数1950-2017.csv",header=TRUE)
USA.deathpopulation<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡人数1950-2017.csv",header=TRUE)
USA.gdp<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/女性/美国人均gdp.csv",head=TRUE)#1960-2018
m_all<-matrix(rep(0,100*68),100,68)   #死亡率
lnm_all<-matrix(rep(0,100*68),100,68) #ln死亡率
D_all<-matrix(rep(0,100*68),100,68)   #计算的死亡人数
E_all<-matrix(rep(0,100*68),100,68)   #总人数
d_all<-matrix(rep(0,100*68),100,68)   #实际观测的死亡人数

year<-1950:2012#63年
age<-0:99      #100个
cohort<-1851:2012 #100+63-1=162

#女性1950-2017年的所有数据
for (t in 1:68){
  for(x in 1:100){
    E_all[x,t]<-USA.population[x+100*(t-1),3];   #总人数，100行，68列
    m_all[x,t]<-USA.mortality[x+100*(t-1),3];    #女性死亡率
    D_all[x,t]<-USA.deathpopulation[x+100*(t-1),3] #实际观测的死亡人数
    lnm_all[x,t]<-log(m_all[x,t]);
  }
}
#用于实验的数据
E<-E_all[,1:63] #1950-2012
m<-m_all[,1:63] #1950-2012
lnm<-lnm_all[,1:63]
m_original<-lnm_all[,64:68]
D<-D_all[,1:63]
y<-D_all[,1:63]#实际观测到的死亡人数

#估计的参数
USA.guji<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/女性/估计参数.csv",header=TRUE)
USA.kt<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/女性/估计的kt.csv",header=TRUE)

k_svd<-USA.kt$k_svd
k_likelihood<-USA.kt$k_likelihood
k_wls<-USA.kt$k_wls

a_svd<-USA.guji$a_svd
a_likelihood<-USA.guji$a_likelihood
a_wls<-USA.guji$a_wls

b_svd<-USA.guji$b_svd
b_likelihood<-USA.guji$b_likelihood
b_wls<-USA.guji$b_wls

1, svd

#计算残差
e_svd<-matrix(rep(0,100*63),100,63)
for(x in 1:100){
  for(t in 1:63){
    e_svd[x,t]<-lnm[x,t]-(a_svd[x]+b_svd[x]*k_svd[t])
  }
}
#画热度图
library(pheatmap)
library(RColorBrewer)
colnames(e_svd)=c(1950:2012)
rownames(e_svd)=c(0:99)
#colnames(e_svd) = paste("year",1950:2012,sep = "")
#rownames(e_svd) = paste("age", 0:99,sep = "")
pheatmap(e_svd,
         color=brewer.pal(11,"PRGn"),
         cellwidth = 4,#每个单元格的宽度
         cellheight = 4,#每个单元格的高度
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         labels_row=c(0:99),
         fontsize = 6,#热图中字体显示的大小
         #display_numbers = TRUE,
         number_color = "black"
)

2，wls

#计算残差
e_wls<-matrix(rep(0,100*63),100,63)
for(x in 1:100){
  for(t in 1:63){
    e_wls[x,t]<-sqrt(y[x,t])*(lnm[x,t]-(a_wls[x]+b_wls[x]*k_wls[t]))
  }
}
#画热度图
library(pheatmap)
library(RColorBrewer)
colnames(e_wls)=c(1950:2012)
rownames(e_wls)=c(0:99)
pheatmap(e_wls,
         color=brewer.pal(11,"PRGn"),
         cellwidth = 4,#每个单元格的宽度
         cellheight = 4,#每个单元格的高度
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         labels_row=c(0:99),
         fontsize = 6,#热图中字体显示的大小
         #display_numbers = TRUE,
         number_color = "black"
)

3，MLE

#计算残差
D_hat<-matrix(rep(0,100*63),100,63)
for(x in 1:100){
  for(t in 1:63){
    D_hat[x,t]<-E[x,t]*exp(a_likelihood[x]+b_likelihood[x]*k_likelihood[t]) #估计的死亡人数
  } 
}
e_mle<-matrix(rep(0,100*63),100,63)
for(x in 1:100){
  for(t in 1:63){
    e_mle[x,t]<-sign(D[x,t]-D_hat[x,t])*(D[x,t]*log(D[x,t]/D_hat[x,t])-(D[x,t]-D_hat[x,t]))^(1/2)
  }
}
which(is.na(D))
which(D<0)
which(is.na(D_hat))
which(D_hat<0)
#画热度图
library(pheatmap)
library(RColorBrewer)
colnames(e_mle)=c(1950:2012)
rownames(e_mle)=c(0:99)
pheatmap(e_mle,
         color=brewer.pal(11,"PRGn"),
         cellwidth = 4,#每个单元格的宽度
         cellheight = 4,#每个单元格的高度
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         labels_row=c(0:99),
         fontsize = 6,#热图中字体显示的大小
         #display_numbers = TRUE,
         number_color = "black"
)

从图中可以看出，有明显的队列效应未能捕捉到
三种方法估计和预测的死亡率差别不大
但WLS和MLE具有更大的随机性