Chapter 10 simulation
The HR of HCC did not attenuated by after control liver cirrhosis in current multivariable Cox-Proportional hazard model. That result suggest the effect of aspirin use for HCC prevention are independent to status of liver cirrhosis or not. However, the stratification analysis show different pattern according to HCC status. The HR for HCC by using aspirant in liver cirrhosis patient is 0.81 with statistical significance, but that in without liver cirrhosis patient are not. The HR and 95%CI for HCC in without liver cirrhosis patient is 0.88 and 0.74-1.04. That might be the relative small number of aspirin usage in current study. 9837 of 161673 (00.00%) are used aspirin in current study, while 00.00% in western country (ref:nhanes.. havard boston Nursing study). Hence, we simulate the data as increased number of aspirin user with fixing the ratio of prevalence of HCC and HR of each group.
= c('non', 'asp')
asp = c(2068, 168)
hcc = c(118317, 8251)
nhcc= data.frame(asp, hcc, nhcc)
test library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.1 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
%>%
test mutate(inc = hcc/nhcc) %>%
mutate(hr = inc /0.02036117 )
## asp hcc nhcc inc hr
## 1 non 2068 118317 0.01747847 0.8584216
## 2 asp 168 8251 0.02036117 0.9999999
= function(x1, n1, x2, n2){
ci95 = log((x1/n1)/(x2/n2)) - 1.96*sqrt(((n1 - x1)/x1)/n1 + ((n2-x2)/x2)/n2)
ll = log((x1/n1)/(x2/n2)) + 1.96*sqrt(((n1 - x1)/x1)/n1 + ((n2-x2)/x2)/n2)
ul = list(exp(ll), exp(ul))
ci
}
= 1
i = ci95(2068*i, 118317*i, 168*i, 8251*i)
gg
= lapply(c(60:70)/100,
gg2 function(i){
ci95(2068*i, 118317*i, 168*i, 8251*i)
})
gg2
## [[1]]
## [[1]][[1]]
## [1] 0.7021556
##
## [[1]][[2]]
## [1] 1.049465
##
##
## [[2]]
## [[2]][[1]]
## [1] 0.7033178
##
## [[2]][[2]]
## [1] 1.047731
##
##
## [[3]]
## [[3]][[1]]
## [1] 0.7044537
##
## [[3]][[2]]
## [1] 1.046042
##
##
## [[4]]
## [[4]][[1]]
## [1] 0.7055641
##
## [[4]][[2]]
## [1] 1.044395
##
##
## [[5]]
## [[5]][[1]]
## [1] 0.7066501
##
## [[5]][[2]]
## [1] 1.04279
##
##
## [[6]]
## [[6]][[1]]
## [1] 0.7077126
##
## [[6]][[2]]
## [1] 1.041225
##
##
## [[7]]
## [[7]][[1]]
## [1] 0.7087524
##
## [[7]][[2]]
## [1] 1.039697
##
##
## [[8]]
## [[8]][[1]]
## [1] 0.7097703
##
## [[8]][[2]]
## [1] 1.038206
##
##
## [[9]]
## [[9]][[1]]
## [1] 0.7107671
##
## [[9]][[2]]
## [1] 1.03675
##
##
## [[10]]
## [[10]][[1]]
## [1] 0.7117434
##
## [[10]][[2]]
## [1] 1.035328
##
##
## [[11]]
## [[11]][[1]]
## [1] 0.7127001
##
## [[11]][[2]]
## [1] 1.033938
= function(non, asp){
sim ci95(non*0.01747847 , non, asp*0.02036117, asp)
}
=lapply(9837:20000, function(x)list(sim(126568, x))[[1]][[2]])
upperplot(unlist(upper))
%>% unlist() %>% tibble() %>%
upper setNames('upper') %>%
mutate(ranks = row_number()) %>%
filter(upper <=1)
## # A tibble: 10,164 × 2
## upper ranks
## <dbl> <int>
## 1 0.991 1
## 2 0.991 2
## 3 0.991 3
## 4 0.991 4
## 5 0.991 5
## 6 0.991 6
## 7 0.991 7
## 8 0.991 8
## 9 0.990 9
## 10 0.990 10
## # … with 10,154 more rows