R中L-BFGS-B方法使用最优解时的非有限差分误差

Question

我正试图使用我的对数似然函数来获得五个参数估计的置信区间(1π，一个在0和1之间的比例，4贝塔是指数下降率(。我的数据包括接种疫苗后的时间间隔、分类为1和0的免疫状态，以及过去的疫苗接种历史和个别狗的所有权状态。我的数据如下

data4mle <- 
structure(
list(
Interval_18 = c(35, 123, 13, 140, 8, 115, 16, 
131, 132, 8, 13, 146, 44, 15, 138, 3, 143, 149, 12, 153, 17, 
154, 14, 8, 123, 16, 17, 154, 8, 135, 10, 133, 18, 8, 123, 10, 
133, 16, 163, 10, 150, 10, 3, 8, 8, 139, 22, 9, 122, 1, 130, 
16, 140, 16, 140, 20, 142, 42, 122, 16, 8, 126, 8, 8, 139, 20, 
122, 14, 8, 122, 146, 6, 20, 9, 6, 13, 13, 115, 12, 131, 12, 
12, 23, 116, 6, 102, 13, 137, 6, 107, 12, 14, 115, 18, 123, 155, 
158, 163, 35, 155, 16, 163, 16, 140, 44, 128, 9, 21, 36, 46, 
153, 27, 147, 47, 131, 8, 146, 12, 17, 9, 16, 155, 9, 10, 148, 
32, 148, 25, 15, 138, 46, 158, 121, 8, 131, 8, 8, 118, 31, 154, 
32, 144, 15, 8, 16, 8, 101, 15, 144, 12, 20, 10, 15, 142, 13, 
69, 28, 168, 15, 12, 24, 13, 142, 12, 115, 13, 131, 166, 13, 
118, 12, 32, 6, 36, 164, 16, 135, 20, 46, 19, 137, 11, 150, 6, 
107, 10, 148, 13, 118, 19, 164, 10, 16, 127, 19, 161, 156, 14, 
16, 145, 1, 143, 16, 5, 144, 158, 33, 19, 6, 6, 115, 16, 138, 
32, 13, 18, 13, 11, 158, 10, 151, 16, 16, 69, 17, 166, 26, 158, 
11, 11, 26, 16, 173, 16, 32, 155, 12, 16, 149, 16, 149, 16, 163, 
153, 57, 13, 128, 12, 113, 19, 144, 21, 35, 30, 30, 157, 108, 
14, 16, 33, 23, 149, 14, 14, 8, 18, 164, 115, 19, 11, 151, 19, 
144, 21, 13, 142, 54, 150, 16, 163, 22, 13, 11, 155, 155, 16, 
150, 16, 16, 57, 19, 17, 135, 12, 11, 12, 145, 11, 11, 170, 18, 
164, 10, 16, 12, 11, 155, 27, 11, 24, 16, 135, 47, 12, 17, 149, 
11, 12, 18, 137, 10, 131, 11, 151, 12, 131, 11, 158, 57, 12, 
5, 27, 23, 16, 24, 149, 29, 150, 12, 16, 157, 16, 10, 152, 16, 
10, 149, 23, 11, 30, 16, 141, 27, 147, 11, 151, 16, 157, 71, 
172, 32, 151, 11, 151, 10, 131, 27, 147, 22, 23, 20, 139, 10, 
131, 12, 145, 16, 149, 10, 135, 145, 13, 135, 16, 149, 30, 150, 
16, 163, 141, 51, 131, 16, 25, 16, 140, 32, 144, 41, 129, 14, 
8, 123, 16, 149, 101, 14, 94, 10, 139, 8, 10, 130, 96, 8, 24, 
144, 8, 8, 129, 20, 129, 10, 131, 8, 9, 131, 8, 127, 10, 130, 
22, 142, 8, 127, 13, 127, 10, 131, 10, 131, 24, 166, 8, 129, 
10, 131, 9, 20, 146, 8, 127, 10, 134, 10, 131, 10, 132, 10, 15, 
15, 20, 142, 8, 127, 11, 9, 11, 170, 10, 131, 20, 129, 144, 10, 
130, 10, 130, 8, 127, 10, 131, 8, 8, 130, 10, 14, 128, 10, 148, 
13, 127, 10, 131, 8, 127, 9, 125, 14, 129, 13, 13, 129, 27, 147, 
9, 131, 8, 130, 14, 128, 1, 105, 1, 105, 9, 13, 137, 8, 10, 131, 
17, 148, 13, 127, 8, 10, 8, 10, 10, 122, 116, 14, 128, 14, 128, 
11, 159, 10, 131, 20, 122), 
Immune_status = c(1, 1, 0, 0, 1, 
0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 
1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 
1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 
0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 
0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 
1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 
1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 
1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 
0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 
1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 
0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 
1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 
0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 
1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 
0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 
1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 
1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 
1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 
0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1), 
Owned = c("No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "Yes", 
"Yes", "No", "No", "No", "No", "No", "Yes", "Yes", "No", "Yes", 
"Yes", "No", "Yes", "Yes", "Yes", "Yes", "No", "No", "No", "Yes", 
"Yes", "No", "No", "Yes", "Yes", "No", "No", "No", "No", "Yes", 
"Yes", "Yes", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "Yes", "Yes", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "Yes", "Yes", "No", 
"No", "No", "No", "No", "No", "Yes", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "Yes", "Yes", "Yes", 
"Yes", "No", "Yes", "Yes", "No", "No", "No", "No", "No", "No", 
"Yes", "Yes", "No", "No", "Yes", "Yes", "No", "No", "No", "No", 
"No", "Yes", "Yes", "No", "No", "No", "No", "No", "Yes", "No", 
"Yes", "Yes", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "Yes", "Yes", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "Yes", 
"Yes", "Yes", "Yes", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "Yes", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "Yes", "Yes", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "Yes", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "Yes", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "Yes", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "Yes", "No", "Yes", "Yes", "No", "Yes", "Yes", 
"No", "No", "No", "Yes", "No", "No", "No", "No", "Yes", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "Yes", 
"Yes", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "No", "Yes", "Yes", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"Yes", "Yes", "No", "No", "No", "No", "No", "No", "Yes", "Yes", 
"No", "No", "No", "No", "No", "No", "No", "Yes", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "Yes", "Yes", "Yes", "Yes", "No", "No", "No", "Yes", 
"No", "No", "No", "Yes", "Yes", "No", "No", "No", "Yes", "Yes", 
"No", "No", "No", "No", "No", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "No", 
"No", "No", "No", "Yes", "No", "No", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", 
"No", "No"), 
Not_vacc_not_revacc_before_R3 = c("No", "No", "No", 
"No", "No", "No", "No", "No", "Yes", "No", "No", "No", "No", 
"No", "No", "Yes", "Yes", "No", "Yes", "Yes", "No", "No", "Yes", 
"No", "No", "No", "Yes", "Yes", "No", "Yes", "Yes", "Yes", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "Yes", 
"No", "Yes", "Yes", "No", "Yes", "No", "No", "No", "No", "No", 
"Yes", "Yes", "Yes", "Yes", "No", "No", "No", "No", "Yes", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "No", "Yes", "Yes", 
"No", "Yes", "No", "No", "No", "Yes", "Yes", "No", "No", "Yes", 
"Yes", "No", "No", "Yes", "Yes", "No", "No", "No", "No", "No", 
"No", "No", "No", "Yes", "Yes", "No", "No", "No", "Yes", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "Yes", "No", "No", 
"Yes", "Yes", "No", "Yes", "No", "No", "No", "Yes", "Yes", "Yes", 
"No", "No", "No", "No", "Yes", "No", "No", "Yes", "Yes", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "Yes", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "Yes", 
"Yes", "No", "No", "Yes", "Yes", "No", "No", "No", "Yes", "Yes", 
"No", "No", "Yes", "Yes", "No", "No", "No", "No", "Yes", "No", 
"Yes", "Yes", "No", "No", "Yes", "Yes", "No", "No", "No", "No", 
"No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes", "No", "No", 
"No", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes", "No", 
"No", "No", "Yes", "No", "No", "No", "No", "No", "No", "No", 
"Yes", "No", "No", "No", "Yes", "Yes", "No", "No", "No", "No", 
"No", "No", "No", "Yes", "Yes", "Yes", "No", "No", "Yes", "Yes", 
"Yes", "Yes", "Yes", "No", "Yes", "Yes", "No", "No", "No", "No", 
"No", "Yes", "No", "No", "No", "No", "No", "No", "Yes", "No", 
"No", "Yes", "Yes", "No", "No", "No", "Yes", "No", "No", "No", 
"No", "No", "Yes", "Yes", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "Yes", "Yes", "Yes", "Yes", "No", "No", "Yes", 
"Yes", "No", "Yes", "Yes", "Yes", "Yes", "No", "No", "Yes", "No", 
"No", "Yes", "No", "No", "Yes", "Yes", "Yes", "No", "No", "Yes", 
"Yes", "No", "No", "No", "No", "Yes", "Yes", "No", "No", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "No", "No", "No", 
"No", "No", "No", "No", "No", "Yes", "No", "Yes", "No", "No", 
"No", "Yes", "Yes", "No", "No", "No", "Yes", "Yes", "Yes", "No", 
"No", "No", "No", "No", "No", "No", "Yes", "Yes", "Yes", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "No", "No", "No", 
"No", "Yes", "Yes", "Yes", "Yes", "Yes", "No", "No", "No", "Yes", 
"Yes", "No", "No", "Yes", "Yes", "Yes", "Yes", "No", "Yes", "Yes", 
"Yes", "Yes", "No", "No", "No", "No", "No", "No", "No", "No", 
"Yes", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "Yes", "Yes", "No", "Yes", "Yes", 
"No", "No", "No", "No", "Yes", "No", "No", "Yes", "Yes", "Yes", 
"Yes", "Yes", "No", "No", "No", "No", "Yes", "Yes", "Yes", "Yes", 
"Yes", "Yes", "No", "No", "Yes", "Yes", "Yes", "Yes", "Yes", 
"Yes", "No", "No", "Yes", "Yes", "No", "No", "No", "No", "No", 
"Yes", "No", "No", "No", "Yes", "Yes", "No", "No", "No", "No", 
"No", "No", "No", "No", "Yes", "No", "No", "No", "No", "Yes", 
"Yes", "Yes", "No", "No", "No", "No", "No", "No", "No", "No", 
"Yes", "Yes", "Yes", "Yes", "No", "No", "Yes", "Yes", "No", "No", 
"No", "No", "No", "No", "No", "No", "No", "No", "No", "No", "No", 
"No", "No", "No", "No", "No", "No", "Yes", "Yes", "No", "No", 
"No", "No", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "No", "No", 
"No", "No", "Yes", "Yes", "Yes", "Yes", "Yes", "No", "Yes", "Yes", 
"Yes", "No", "No", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")), 
row.names = c(NA,-542L), 
class = c("tbl_df", "tbl", "data.frame"))

我的对数似然函数

nll.18.4betas <- function(v) {
#function(pi,beta) {
pi <- v[1]; beta_o_v <- v[2]; beta_o_nv <- v[3]; beta_no_v <- v[4]; beta_no_nv <- v[5]

data <- data4mle %>% 
filter(Interval_18 > 0 & Interval_18 <= 180 & !is.na(Immune_status))

#Owned, vaccinated/revaccinated
xi_o_v <- data %>% 
filter(Owned == "Yes" & Not_vacc_not_revacc_before_R3 == "Yes") %>% 
.$Immune_status
ti_o_v <- data %>% 
filter(Owned == "Yes" & Not_vacc_not_revacc_before_R3 == "Yes") %>% 
.$Interval_18

#Owned, not vaccinated/revaccinated
xi_o_nv <- data %>% 
filter(Owned == "Yes" & Not_vacc_not_revacc_before_R3 == "No") %>% 
.$Immune_status
ti_o_nv <- data %>% 
filter(Owned == "Yes" & Not_vacc_not_revacc_before_R3 == "No") %>% 
.$Interval_18

#Not owned, vaccinated/revaccinated
xi_no_v <- data %>% 
filter(Owned == "No" & Not_vacc_not_revacc_before_R3 == "Yes") %>% 
.$Immune_status
ti_no_v <- data %>% 
filter(Owned == "No" & Not_vacc_not_revacc_before_R3 == "Yes") %>% 
.$Interval_18

##Not owned, not vaccinated/revaccinated
xi_no_nv <- data %>% 
filter(Owned == "No" & Not_vacc_not_revacc_before_R3 == "No") %>% 
.$Immune_status
ti_no_nv <- data %>% 
filter(Owned == "No" & Not_vacc_not_revacc_before_R3 == "No") %>% 
.$Interval_18

#calculate neg log likelihood
#Owned, vaccinated/revaccinated
-sum(xi_o_v*(log(pi)-beta_o_v*ti_o_v) + (1-xi_o_v)*log(1 - (pi*exp(-beta_o_v*ti_o_v)))) + 
#Owned, not vaccinated/revaccinated
-sum(xi_o_nv*(log(pi)-beta_o_nv*ti_o_nv) + (1-xi_o_nv)*log(1 - (pi*exp(-beta_o_nv*ti_o_nv)))) + 
#Not owned, vaccinated/revaccinated
-sum(xi_no_v*(log(pi)-beta_no_v*ti_no_v) + (1-xi_no_v)*log(1 - (pi*exp(-beta_no_v*ti_no_v)))) +
#Not owned, not vaccinated/revaccinated
-sum(xi_no_nv*(log(pi)-beta_no_nv*ti_no_nv) + (1-xi_no_nv)*log(1 - (pi*exp(-beta_no_nv*ti_no_nv))))

}
trial <- optim(fn = nll.18.4betas, par = c(0.1,  0.003,0.003,0.003, 0.003), 
lower = c(0.1, rep(1e-3, 4)),
upper = c(0.99,rep(1e-1, 4)), 
method = "L-BFGS-B", hessian = T)

当试图估计hessian时，我一直得到非有限差分误差，即使优化方法本身运行良好，并生成五个参数的估计值(即当hessian = FALSE时(。我需要使用L-BFGS-B，因为第一个参数pi是一个必须绑定在0和1之间的比例。

我看了这个问题中的建议(非有限差分值，许多数据在指数后变成inf和NA(，但据我所知，我需要使用Nelder-Mead方法来使用numDeriv包。

类似地，对这个问题的回答(Optim：L-BFGS-B中的非有限有限差分值(似乎不适用，我不确定这里讨论的内容是否与我的问题直接相关(r中的Optim：非有限差分误差(。

我不确定这里发生了什么，特别是因为当我对一个有六个参数(2个pis和4个betas(的嵌套模型进行优化时，"optim"能够生成一个粗麻布，因此我可以计算置信区间。

不确定这是否为某人提供了足够的信息来提供帮助，但如果有任何建议，我们将不胜感激。谢谢

Answer 1

以下代码似乎根据您的限制生成了一个有效的解决方案：

set.seed(100)
trial <- optim(fn = nll.18.4betas, 
par = runif(5), 
method = "Nelder-Mead", hessian = T,
control=list(maxit=5000))
# RESULT
$par
[1] 0.999362670 0.005248378 0.011744520 0.002088864 0.007242104
$value
[1] 258.8472
$counts
function gradient 
1034       NA 
$convergence
[1] 0
$message
NULL
$hessian
[,1]       [,2]       [,3]      [,4]       [,5]
[1,]    9811.427  -4834.189  -9284.515 -106446.3  -43436.32
[2,]   -4834.189 677710.202      0.000       0.0       0.00
[3,]   -9284.515      0.000 354806.186       0.0       0.00
[4,] -106446.252      0.000      0.000 9402389.7       0.00
[5,]  -43436.320      0.000      0.000       0.0 1259387.74