sim_bwselect1 = function(n, blip, truth, truths_h, bw_seq, g0, Q0, kernel)
{
tmledata = gentmledata(n, d = 1, g0, Q0, V=10, formu = NULL)
big_ind = length(bw_seq)
# undebug(gentmle_alt3)
test_bwselect = lapply(bw_seq, FUN = function(h) {
gentmle_alt3(initdata=tmledata, estimate_fun = blipdist_estimate2,
update_fun = blipdist_update, max_iter = 1000, kernel = kernel,
simultaneous.inference = FALSE, blip = blip, h = h)
})
est_vec = unlist(lapply(test_bwselect, FUN = function(x) x$tmleests))
est_derivs = vapply(2:length(bw_seq), FUN = function(i) {
(est_vec[i]-est_vec[i-1])/.01
}, FUN.VALUE = 1)
# est_derivs
SE_vec = unlist(lapply(1:length(bw_seq), FUN = function(i) {
x = test_bwselect[[i]]
sqrt((x$ED2 - x$ED^2)/n)
}))
SE_derivs = vapply(2:length(SE_vec), FUN = function(i) {
(SE_vec[i]-SE_vec[i-1])/.01
}, FUN.VALUE = 1)
# SE_derivs
# start at a sensible spot, then depending on increasing or decr psi in h
# we travel until no longer helpful
for(i in 6:big_ind) {
ans = all(sign(est_vec[i:(i-4)] - est_vec[(i-1):(i-5)])>=0) |
all(sign(est_vec[i:(i-4)] - est_vec[(i-1):(i-5)])<=0)
if (ans) {
start_h = i-5
break
} else start_h = big_ind
}
decre = est_vec[start_h+5] < est_vec[start_h+4]
for (i in start_h:(big_ind-1)){
if (decre) {
end_h = i + 1
if (!(est_vec[i] >= est_vec[i+1])) {
end_h = i
break
}
} else {
end_h = i + 1
if (!(est_vec[i] <= est_vec[i+1])) {
end_h = i
break
}
}
}
if (decre) ind = order((est_vec - 1.96*SE_vec)[start_h:end_h], decreasing = TRUE)[1] else {
ind = order((est_vec + 1.96*SE_vec)[start_h:end_h])[1]
}
ind = ifelse((start_h:end_h)[ind] == big_ind, big_ind, (start_h:end_h)[ind]+1)
CI_lower = (est_vec - 1.96*SE_vec)[1:big_ind]
CI_upper = (est_vec + 1.96*SE_vec)[1:big_ind]
CIs = cbind(est_vec[1:big_ind], CI_lower, CI_upper)
cover = truth >= CIs[,2] & truth <= CIs[,3]
# do linear regresssion and check positivity of coeffcoefficient
x = CIs[,1]
y = bw_seq[1:big_ind]
incre = coef(lm(y~x))[2] >= 0
ests = CIs[,1]
if (incre) {
estsI = pava(y=ests, decreasing = FALSE)
} else {
estsI = pava(y=ests, decreasing = TRUE)
}
SE_vecI = pava(y=SE_vec, decreasing = TRUE)
if (decre) ind_alt = order((est_vec - 1.96*SE_vecI)[start_h:end_h], decreasing = TRUE)[1] else {
ind_alt = order((est_vec + 1.96*SE_vecI)[start_h:end_h])[1]
}
ind_alt = (start_h:end_h)[ind_alt]
CI_lowerI = estsI - 1.96*SE_vecI
CI_upperI = estsI + 1.96*SE_vecI
CIsI = data.frame(ests = estsI, lower = CI_lowerI, upper = CI_upperI)
CI_lowerI_alt = ests - 1.96*SE_vecI
CI_upperI_alt = ests + 1.96*SE_vecI
CIsI_alt = data.frame(ests = ests, lower = CI_lowerI_alt, upper = CI_upperI_alt)
coverI_alt = truth >= CIsI_alt[,2] & truth <= CIsI_alt[,3]
if (incre) {
m = order(CIsI[,3])[1]
indI = max(which(CIsI[,3]==CIsI[m,3]))
indI = ifelse(indI == big_ind, indI, indI + 1)
} else {
m = order(CIsI[,2], decreasing = TRUE)[1]
indI = max(which(CIsI[,2]==CIsI[m,2]))
indI = ifelse(indI == big_ind, indI, indI + 1)
}
CIs_all = cbind(CIs, CIsI, CIsI_alt)
return(list(CIs_all = CIs_all
,CIs = rbind(CIs[ind, ], CIsI_alt[ind_alt,], CIs[indI, ], CIsI_alt[indI,], CIs[big_ind,])
,cover = c(cover[ind], coverI_alt[ind_alt], cover[indI], coverI_alt[indI], cover[big_ind])
,truths = c(truth = truth, truth_h = truths_h)
,incre = incre
,inds = c(ind, ind_alt, indI)
,range = c(start_h, end_h)))
}
sim_bwselect = function(n, blip, truth, truths_h, bw_seq, g0, Q0, kernel)
{
tmledata = gentmledata(n, d = 1, g0, Q0, V=10, formu = NULL)
big_ind = length(bw_seq)
# undebug(gentmle_alt3)
test_bwselect = lapply(bw_seq, FUN = function(h) {
gentmle_alt3(initdata=tmledata, estimate_fun = blipdist_estimate2,
update_fun = blipdist_update, max_iter = 1000, kernel = kernel,
simultaneous.inference = FALSE, blip = blip, h = h)
})
est_vec = unlist(lapply(test_bwselect, FUN = function(x) x$tmleests))
est_derivs = vapply(2:length(bw_seq), FUN = function(i) {
(est_vec[i]-est_vec[i-1])/.01
}, FUN.VALUE = 1)
# est_derivs
SE_vec = unlist(lapply(1:length(bw_seq), FUN = function(i) {
x = test_bwselect[[i]]
sqrt((x$ED2 - x$ED^2)/n)
}))
SE_derivs = vapply(2:length(SE_vec), FUN = function(i) {
(SE_vec[i]-SE_vec[i-1])/.01
}, FUN.VALUE = 1)
# SE_derivs
# start at a sensible spot, then depending on increasing or decr psi in h
# we travel until no longer helpful
for(i in 6:big_ind) {
ans = all(sign(est_vec[i:(i-4)] - est_vec[(i-1):(i-5)])>=0) |
all(sign(est_vec[i:(i-4)] - est_vec[(i-1):(i-5)])<=0)
if (ans) {
start_h = i-5
break
} else start_h = big_ind
}
decre = est_vec[start_h+5] < est_vec[start_h+4]
for (i in start_h:(big_ind-1)){
if (decre) {
end_h = i + 1
if (!(est_vec[i] >= est_vec[i+1])) {
end_h = i
break
}
} else {
end_h = i + 1
if (!(est_vec[i] <= est_vec[i+1])) {
end_h = i
break
}
}
}
if (decre) ind = order((est_vec - 1.96*SE_vec)[start_h:end_h], decreasing = TRUE)[1] else {
ind = order((est_vec + 1.96*SE_vec)[start_h:end_h])[1]
}
ind = ifelse((start_h:end_h)[ind] == big_ind, big_ind, (start_h:end_h)[ind]+1)
CI_lower = (est_vec - 1.96*SE_vec)[1:big_ind]
CI_upper = (est_vec + 1.96*SE_vec)[1:big_ind]
CIs = cbind(est_vec[1:big_ind], CI_lower, CI_upper)
cover = truth >= CIs[,2] & truth <= CIs[,3]
# do linear regresssion and check positivity of coeffcoefficient
x = CIs[,1]
y = bw_seq[1:big_ind]
incre = coef(lm(y~x))[2] >= 0
ests = CIs[,1]
if (incre) {
estsI = pava(y=ests, decreasing = FALSE)
} else {
estsI = pava(y=ests, decreasing = TRUE)
}
SE_vecI = pava(y=SE_vec, decreasing = TRUE)
CI_lowerI = estsI - 1.96*SE_vecI
CI_upperI = estsI + 1.96*SE_vecI
CIsI = data.frame(ests = estsI, lower = CI_lowerI, upper = CI_upperI)
if (incre) {
m = order(CIsI[,3])[1]
indI = max(which(CIsI[,3]==CIsI[m,3]))
indI = ifelse(indI == big_ind, indI, indI + 1)
} else {
m = order(CIsI[,2], decreasing = TRUE)[1]
indI = max(which(CIsI[,2]==CIsI[m,2]))
indI = ifelse(indI == big_ind, indI, indI + 1)
}
coverI = truth >= CIs[indI,2] & truth <= CIs[indI,3]
CIs_all = cbind(CIs, CIsI)
return(list(CIs_all = CIs_all
,CIs = rbind(CIs[ind, ], CIs[indI, ], CIs[big_ind,])
,cover = c(cover[ind], coverI, cover[big_ind])
,truths = c(truth = truth, truth_h = truths_h)
,incre = incre
,inds = c(ind, indI)
,range = c(start_h, end_h)))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.