我参照的是——孙佳美的《生命表编制理论与实践》进行经典lee-carter模型死亡率预测。数据如下图我想用1997-2010年死亡率估计2011-2015年的。年龄组有10组。
代码如下：三种参数估计方法——奇异值分解svd，加权最小二乘法WLS，极大似然估计法likelihood。4个地方报错，求解！
library(demography)
library(forecast)
library(xts)
#读取数据
data = read.csv(file = "F://c15total.csv", header = TRUE)
data
break_year<-2010
T<-break_year-1997+1 #用来建模年份
n<-2015-break_year #用来预测的年份
D_all<-matrix(rep(0,10*19),10,19) #年龄组数：10组和年数：14
E_all<-matrix(rep(0,10*19),10,19)
m_all<-matrix(rep(0,10*19),10,19)
lnm_all<-matrix(rep(0,10*19),10,19)
for(t in 1:19){
for(x in 1:10){
D_all[x,t]<-mydata[x+10*(t-1),3];
E_all[x,t]<-mydata[x+10*(t-1),4];
m_all[x,t]<-D_all[x,t]/E_all[x,t];
lnm_all[x,t]<-log(m_all[x,t])
}
}

D<-D_all[,1:T];
E<-E_all[,1:T];
m<-m_all[,1:T];
lnm_all<-lnm_all[,1:T]; #建模数据

m_original<-m_all[,(T+1):19] #预测数据

age<-mydata[1:10,2];
year<-as.numeric(levels(factor(mydata[,1])))
oldyear<-year[1:T];
newyear<-(T+1):(n+T)

#参数估计——加权最小二乘法WLS
a<-c(rep(0,10));
b<-c(rep(1,10));
k<-c(rep(0,T)) #设置初始值
#建模
RSS<-function(a,b,k){
L<-matrix(rep(NA,10*T),10,T)
for(x in 1:10)
{for(t in 1:T)
{L[x,t]<-D[x,t]*((lnm[x,t]-a[x]-b[x]*k[t])^2)
}
}
RSS<-sum(L)
return(RSS)
}
#迭代a(x)年龄死亡率的均值
aa<-function(a,b,k){
for(x in 1:10){
a[x]<-sum(D[X,]*(lnm[x,]-b[x]*k))/sum(D[x,])
}
return(a)
}
#迭代k(t)时间因子
kk<-function(a,b,k){
for(t in 1:15){
k[t]<-sum(D[,t]*k*(lnm[x,]-a))/sum(D[,t]*b^2)
}
return(b)
}
#迭代b(x)年龄因子
bb<-function(a,b,k){
for(x in 1:10){
b[x]<-sum(D[x,]*k*(lnm[x,]-a[x]))/sum(D[x,]*k^2)
}
return(b)
}
#迭代过程
j=1
diff_RSS<-RSS(a,b,k) Error in RSS(a, b, k) : object 'lnm' not found
while(abs(diff_RSS)>10^(-10))
{RSS0<-RSS(a,b,k)
if(j%%3==1){a<-aa(a,b,k)}
if(j%%3==2){k<-kk(a,b,k)}
if(j%%3==0){b<-bb(a,b,k)}
RSS1<-RSS(a,b,k)
diff_RSS<-RSS1-RSS0
j=j+1
} Error in abs(diff_RSS) : object 'diff_RSS' not found
#标化
a1<-a;
b1<-b;
k1<-k
c1<-sum(b1);
c2<-sum(k1)/T
a_wls<-a1+c2*b/c1;
b_wls<-b1/c1;
k_wls<-c1*k1-c2


#参数估计——奇异值分解法SVD
age_svd<-c(1:10)
std.m<-demogdata(m,E,ages = age_svd,years = oldyear,type = "mortality",label = "swden",
name = "c15",lambda = 0)
swden<-lca(std.m,years = std.m$year,ages = std.m$age,adjust = "dt",restype = "logrates")
Warning messages:
1: In stats::nlm(fred, guess, ...) : NA/Inf被换成最大的正值
2: In stats::nlm(fred, guess, ...) : NA/Inf被换成最大的正值
3: In stats::nlm(fred, guess, ...) : NA/Inf被换成最大的正值
4: In newroot(FUN, guess, ...) : No root exists. Returning closest
5: In newroot(FUN, guess, ...) : No root exists. Returning closest

a_svd<-swden$ax;b_svd<-swden$bx;k_svd<-swden$kt

#导出参数估计数据a(x),b(x),k(t)
age<-as.data.frame(age)
result1<-as.data.frame(cbind(age,a_svd,a_wls,a_likelihood,b_svd,b_wls,b_likelihood));
result2<-as.data.frame(cbind(year[1:T],k_svd,k_wls,k_likelihood))

write.table(result1,"parameters.csv",append=T,sep="",col.names=c("age","a_svd","a_wls",
"a_likelihood","b_svd","b_wls","b_likelihood"),row.names=F)
Warning message:
In write.table(result1, "parameters.csv", append = T, sep = "",:
给文件加列名

write.table(result2,"parameters.csv",append=T,sep="",col.names=c("year","k_svd","k_wls",
"k_likelihood"),row.names=F)


Warning message:
In write.table(result1, "parameters.csv", append = T, sep = "",:
给文件加列名

#代码结束
结果如下

后面还有用参数估计出来的k(t)根据ARIMA模型预测k_predict(t)，然后再和a(x)、b(x)算出死亡率。我这部分都没有代码，求代码啊求代码。真是要跪了。没懂是什么意思，我用R软件算出的只有40岁以上10个年龄组的a(x)、b(x)值和2011-2015年k_predict值，最后怎么拟合出死亡率的预测值啊？
求大神解答！！！！有偿，非数理专业，医学统计出身