#LOG
# 2023-9-26 Start of programming
# 2023-10-12  v1
# 2024-3-13   v1.1 
# 2024-4-2    v1.2
# 2024-6-6    v1.3
# 2025-04-22  v2
# 2025-8-16   v2.1

#R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
#Copyright (C) 2024 The R Foundation for Statistical Computing
#Platform: x86_64-w64-mingw32/x64

#Packages used
library(data.table)
library(lubridate)
library(ggplot2)
library(survival)
library(survminer)
library(bit64)
library(dplyr)
library(janitor)
library(openxlsx)
library(tableone)
library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
library(plyr)
library(parallel)
library(bread)
library(haven)
library(table1)
library(kableExtra)

setwd("E:/workdata/704766/personlige_mapper/anders/hrt_all")

#setting constants
options(scipen=99)
aar=365.2
prec=3
prolong=0

###################Chapter 1: Central Persons Register
#setting missing parents to NA
cpr=fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/t_person2022.csv")[,
      n_m:=.N, pnrm][n_m>30, pnrm:=NA][,
       n_f:=.N, pnrf][n_f>30, pnrf:=NA]

#1.1 Migrations
vnds=fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/vnds2022.csv")
first_ind=vnds[,.SD[1], pnr][indud_kode=="I"][,.(pnr, first_ind=haend_dato)]
censor6=vnds[, `:=` (n_haend_dato=shift(haend_dato, n=1, type="lag"),
                     n_indud_kode=shift(indud_kode, n=1, type="lag")), pnr][indud_kode=="I" & n_indud_kode=="U",dif:=haend_dato-n_haend_dato][dif>=180, cens6_dato:=n_haend_dato+180][!is.na(cens6_dato)][,.(pnr, cens6_dato)][cens6_dato>=as.IDate("1995-01-01")][,.SD[1], pnr]

cpr1=merge(x=cpr, y=first_ind, all.x=T, by="pnr")
cpr2=merge(x=cpr1, y=censor6, all.x=T, by="pnr")

#1.2 Country of birth
# bef=fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/bef.csv")
# bef1=unique(bef[,.(pnr, ie_type, fland)])
# bef2=bef1[order(pnr, fland)]
# 
# fwrite(bef2, "data/bef2.csv")

bef2=fread("data/bef2.csv", encoding="Latin-1")[order(pnr, fland)][,.SD[1], pnr][,fland:=gsub("ø", "oe", fland)]

cpr3=merge(x=cpr2, y=bef2, all.x=T, by="pnr")

#1.3 Extracting population
cpr31=cpr3[c_kon=="K" & !(c_status %in% c(5,7))][
  year(fdato)>= 1950 & year(fdato)<=1977][
  c_status %like% "^[2-8]", edato:=d_status_hen_start][
  first_ind<=fdato+(25*aar)| is.na(first_ind)][
    c_status=="90", ddato:=d_status_hen_start][,
      start:=as.IDate(pmax(fdato+(45*aar), as.IDate(1995-01-01), first_ind, na.rm=T))][,
         stop:=as.IDate(pmin(as.IDate("2023-07-31"), edato, cens6_dato, ddato, na.rm=T))][start<stop][,
           .(pnr, pnrm, pnrf, fdato, start, stop, start_age=round((start-fdato)/aar, prec), stop_age=round((stop-fdato)/aar, prec), death_age=round((ddato-fdato)/aar,prec), fland)]

############## Chapter 2: LMDB Prescription register
# lmdb2022=fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/hc2022.csv", encoding="Latin-1")
# 
# #2.1 Extracting MHT from prescription register
# ht=lmdb2022[atc %like% "^G03FA|^G03FB|^G03CA|^G03DA|^G03DB|^G03DC|^G02BA"] #Removing ^G03HB Diane Mite
# 
# #Multiple packs in one day are added together
# ht1=merge(x=cpr31[,.(pnr,fdato)], y= ht, all.x=T, by="pnr")[order(pnr, eksd)][!is.na(eksd)][, id:=.GRP, c("pnr", "eksd", "pname")]
# 
# ht2=dcast(ht1, id+pnr+eksd+pname+atc+streng+packsize+volume+dosform+strnum+fdato~., value.var="apk", fun=sum)
# setnames(ht2, old=".", new="apk")
# 
#Categories
# ht3=ht2[!(pname %like% "Prolutex|Vagifem|Estring|Rewellfem|Oestring|Premarin|Crinone|Cyclogest|Lutinus|Pleyris|Progestan")][!(pname=="Ovestin" & dosform %like% "VAG")][,
#           .(pnr, id, fdato, eksd, apk, pname, atc, streng, packsize, volume, dosform, strnum)][
#             dosform %in% c("TABOTR","TAB","KAPB", "TABFILM"), dosform:="oral"][
#               dosform=="IUTIND", dosform:="intrauterin"][
#                 dosform %in% c("DEPPLA","GEL","KUTSPRO"), dosform:="transderm"][
#                   dosform=="NASSPRO", dosform:= "nosespray"][
#                     dosform=="SUP", dosform:= "rectal"][
#                       dosform %in% c("VAGCREM", "VAGTAB", "VAGIND", "VAG", "VAGGEL"), dosform:= "vaginal"][
#                         dosform=="INJVSK", dosform:= "injection"][,
#                           pname:=gsub("\"","", pname)][,
#                           streng:=gsub("mg","", streng)][,
#                              eksd_stop:=eksd+(apk*volume)+prolong][
#                                atc=="G02BA03" & !(pname %like% "Jaydess"), eksd_stop:=as.IDate(eksd+(5*aar))][
#                                  pname %like% "Jaydess", eksd_stop:=as.IDate(eksd+(3*aar))][
#                                    atc %like% "^G03FA", httype:="combined_monophase"][
#                                      atc %like% "^G03FB|^G03HB01", httype:="combined_cyclic"][
#                                        atc %like% "^G03CA", httype:="estrogen"][
#                                          atc %like% "^G03D[A-C]|^G02BA03", httype:="progestin"]
# 
# 
# estrogen=ht3[httype!="progestin"][,`:=`
#                                    (n=shift(eksd_stop, n=1, type="lag"),
#                                      n1=shift(atc, n=1, type="lag"),
#                                      n2=shift(streng, n=1, type="lag"),
#                                      n3=shift(pname, n=1, type="lag"),
#                                      n4=shift(pnr, n=1, type="lag"))][
#                                        n<eksd |atc != n1 |streng != n2 |pname != n3 | n4!=pnr, split:=1][
#                                          1, split:=1][is.na(split), split:=0][,
#                                               split:=cumsum(split)][,
#                                         min_eksd:=min(eksd,na.rm=T), split][,
#                                           max_eksd:=max(eksd_stop,na.rm=T), split][,
#                                                 ddd:=apk*volume]
# 
# estrogen1=dcast(estrogen, split+pnr+min_eksd+max_eksd+atc+dosform+streng+httype+pname+fdato~., value.var="ddd", fun=sum)[order(pnr, min_eksd)][,split:=1:.N]
# setnames(estrogen1, old=".", new="ddd")
# 
# estrogen2=estrogen1[httype=="estrogen"]
# 
# prog=ht3[httype=="progestin"][,.(pnr, progatc=atc, progstart=eksd, progslut=eksd_stop, progname=pname)]
# 
# ht4=merge(x=estrogen1, y=prog, all.x=T, by="pnr", allow.cartesian = T)[min_eksd<=progstart & max_eksd>=progstart | min_eksd<=progslut & max_eksd>=progslut, prog:=1][order(split, -prog)][prog==1 & httype=="estrogen", httype:="estrogen+progestin"][,.SD[1], split][httype!="estrogen+progestin", `:=` (progatc=NA, progname=NA)][httype=="estrogen+progestin", pname:=paste(pname, progname)][,.(pnr, fdato, htstart=round((as.IDate(min_eksd)-fdato)/aar, prec), htstop=round((as.IDate(max_eksd)-fdato)/aar, prec), htatc=atc, dosform, streng, httype, htname=pname, ddd)]
# 
# fwrite(ht4, "data/ht4.csv")
# 
# ht_pp=ht4[,htstop90:=htstop+(90/aar)][,n:=shift(htstart, n=1, type="lead"), pnr][htstop90<n, keep:=1][order(pnr, -htstop)][, n1:=1:.N, pnr][is.na(keep) & n1==1, keep:=1][order(pnr, htstop)][keep==1][,.SD[1], pnr][,.(pnr, htstop90=round(htstop90, 3))]
# 
# fwrite(ht_pp, "data/ht_pp.csv")

ht4=fread("data/ht4.csv")
ht_pp=fread("data/ht_pp.csv")

#2.2 Diabetes
# dm=lmdb2022[atc %like% "^A10"][pname!= "Wegovy Flextouch"][order(pnr, eksd)][,.SD[1], pnr]
# dm1=merge(x=cpr31[,.(pnr, fdato)], y=dm[,.(pnr, eksd)], all.x=T, by="pnr")[!is.na(eksd)][,.(pnr, diabetes_age=round((eksd-fdato)/aar,prec))]
# fwrite(dm1, "data/dm1.csv")
dm1=fread("data/dm1.csv")

#2.3 Hypertension
# hypertens=lmdb2022[atc %like% "^C"][order(pnr, eksd)]
# 
# hyperten1=hypertens[atc %like% "^C02[A-C]", class:="alfa_adrenerg"][atc %like% "^C02DA|^C02L|^C03A|^C03B|^C03D|^C03E|^C03X|^C07C|^C07D|^C08G|^C09BA|^C09DA|^C09XA52", class:="loop"][atc %like% "^C02DB|^C02DD|^C02DG|^C04|^C05", class:="vasodil"][atc %like% "^C07", class:="beta_block"][atc %like% "^C07F|^C08|^C09BB|^C09DB", class:="calcium_can_block"][atc %like% "^C09", class:="renin_ang_inhib"][!is.na(class)]
# 
# hyperten2=merge(x=cpr31[,.(pnr, fdato)], y=hyperten1, all.x=T, by="pnr")[!is.na(eksd)][,ht_age:=round((eksd-fdato)/aar,prec)]
# 
# hyperten3=hyperten2[order(pnr, eksd)][,n:=1:.N, c("pnr", "class")][n==1][,n:=1:.N, pnr][n==2][,.(pnr, ht_age=round((eksd-fdato)/aar,prec))]
# 
# fwrite(hyperten3, "data/hyperten3.csv")

hypertens=fread("data/hyperten3.csv")

#2.3 Hypercholestrolamia
# statin=lmdb2022[atc %like% "^C10A[A-B]"][order(pnr, eksd)][,.SD[1],pnr]
# 
# statin1=merge(x=cpr31[,.(pnr, fdato)], y=statin[,.(pnr, eksd)], all.x=T, by="pnr")[!is.na(eksd)][,.(pnr, statin_age=round((eksd-fdato)/aar,prec))]
# fwrite(statin1, "data/statin1.csv")

statin1=fread("data/statin1.csv")

#########Chapter 3: LPR Patient register
# 3.1 All surgery 1977-2023
op_all77_23=fread("data/op_all77_23.csv")

op_pop=merge(x=op_all77_23, y=cpr31[,.(pnr, fdato)], all.x=T, by="pnr")[!is.na(fdato)][!is.na(pnr)][,op_age:=round((odto-fdato)/aar,prec)]

hyst=unique(op_pop[okode %like% "^KLCC1[0-1]|^KLCC20|^KLCD|^KLEF13|^KMCA33|^61000|^61020|^61040|^61100|^72110|^72230|^72240|^72650"][,.(pnr, hyst_age=op_age, odto)])[!is.na(odto)][order(pnr, odto)][,.SD[1], pnr]

uso=unique(op_pop[okode %like% "^KLAE1|^KLAF0|^60100|^60101|^60300|^70200|^70300"][,.(pnr, uso_dato=odto, fdato, odto)])[,uso:=1][,uso_sum:=cumsum(uso), pnr][uso_sum>1, bso:=1][bso==1][,bso_age:=round((uso_dato-fdato)/aar, prec)][,.(pnr, bso_age, odto)]#KLCD[3-4]|
bso=op_pop[okode %like% "^KLAF1|^6012[0-1]|^60320|^70210|^70310"][,bso:=1][,.(pnr, bso_age=op_age, odto)]
bso1=unique(rbind(uso,bso)[order(pnr, bso_age)])[!is.na(odto)][order(pnr, odto)][,.SD[1], pnr]

#3.2 Subsetting LPR due to size
# lpr=fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/lpr77_23_ab.csv")
# setkey(lpr1, c_diag)
# 
# lpr1=merge(x=cpr31[,.(pnr, fdato)], y=lpr, all.x=T, by="pnr")
# lpr1=lpr1[!is.na(c_diag)][order(pnr, d_inddto, c_diag)][,diag_age:=round((d_inddto-fdato)/aar, prec)]
# fwrite(lpr1, "data/lpr1.csv")

lpr1=fread("data/lpr1.csv")

#3.3 Incident diagnosis
thromb_inc=lpr1[c_diag %like% "^DD685|^DD686"][,.(pnr, thrombinc_age=diag_age)][,.SD[1],pnr]
liver_inc=lpr1[c_diag %like% "^DK7[0-7]|^57[0-3]"][,.(pnr, liverinc_age=diag_age)][,.SD[1],pnr]
cvd_inc=lpr1[c_diag %like% "^410|^DI2[1-2]|^43[3-4]|^436|^DI63|^DI64"][,.(pnr, cvdinc_age=diag_age)][,.SD[1],pnr]
vte_inc=lpr1[c_diag %like% "^DI26|^DI80[1-3]|^DI80[8-9]|^450|^45100|^4510[8-9]|^45199"][,.(pnr, vteinc_age=diag_age)][,.SD[1],pnr]
breastgyn_inc=lpr1[c_diag %like% "^DC50|DC5[4-7]"][,.(pnr, breastgyninc_age=diag_age)][,.SD[1],pnr]
cancer_inc=lpr1[c_diag %like% "^DC[0-9]"][!(c_diag %like% "^DC44")][,.(pnr, cancerinc_age=diag_age)][,.SD[1],pnr]
afli_inc=lpr1[c_diag %like% "^DI48|^4279[3-4]"][,.(pnr, afli_age=diag_age)][,.SD[1],pnr]
hf_inc=lpr1[c_diag %like% "^DI110|^DI42|^DI50|^DJ819|^425|^4270|^4271"][,.SD[1],pnr]#[,.(pnr, hf_age=diag_age, hfinc_dato=d_inddto)]
#loop=lmdb2022[atc %like% "^C03C"][order(pnr, eksd)][,.SD[1],pnr]
#fwrite(loop, "data/loop.csv")
loop=fread("data/loop.csv")
hf2=merge(x=hf_inc, y=loop[,.(pnr, eksd)], all.x=T, by="pnr")[!is.na(eksd)][, hf_dato:=pmax(d_inddto, eksd)][,.(pnr, hf_age=round((hf_dato-fdato)/aar, prec))]
valv_inc=lpr1[c_diag %like% "^DI05|^DI06|^DI3[4-5]|^39[4-6]|^424[0-1]"][,.(pnr, valv_age=diag_age)][,.SD[1],pnr]

#3.4 All diagnoses
cvd=lpr1[c_diag %like% "^410|^DI2[1-2]|^43[3-4]|^436|^DI63|^DI64"][,.(pnr, cvd_age=diag_age)]
vte=lpr1[c_diag %like% "^DI26|^DI80[1-3]|^DI80[8-9]|^450|^45100|^4510[8-9]|^45199"][,.(pnr, vte_age=diag_age)]
frac=lpr1[c_diag %like% "^DM80|^S72[0-2]^82[0-9]"][,.(pnr, frac_age=diag_age)]

cancer=fread("E:/workdata/704766/personlige_mapper/anders/registre_csv/cancer.csv")[,.(pnr, c_diag=search, d_inddto=d_diagnosedato)]
cancer1=merge(x=cpr31[,.(pnr, fdato)], y=cancer, all.x=T, by="pnr")[!is.na(c_diag)][,diag_age:=round((d_inddto-fdato)/aar, prec)]

breastgyn_cancer=cancer1[c_diag %like% "^DC50|DC5[4-7]"][,.(pnr, breastgyninc_age=diag_age)]
breastgyn_inc1=rbind(breastgyn_inc, breastgyn_cancer)[order(pnr, breastgyninc_age)][,.SD[1], pnr]

cancer_cancer=cancer1[c_diag %like% "^C|^[0-9]"][!(c_diag %like% "^C44|^191[0-9]")][,.(pnr, cancerinc_age=diag_age)]
cancer_inc1=rbind(cancer_cancer, cancer_inc)[order(pnr, cancerinc_age)][,.SD[1], pnr]

#3.5 Number om admissions at age 44 (high use of medical system)
n_diag44=lpr1[c_diagtype=="A" & diag_age>=44 & diag_age<45][,n_diag44:=.N, pnr][,.SD[1], pnr][,.(pnr, n_diag44)][, many_diag:=ifelse(n_diag44>=3, 1, 0)]

######## 4. MFR Medical births register
#Not available yet in updated form, using Steens gravdb
gravdb=fread(file="E:/workdata/704766/personlige_mapper/anders/regstre_csv1/gravdb_steen.csv")
gravdb1=merge(x=cpr31[,.(pnr, fdato)], y=gravdb[,.(pnr, slutdato, gruppe, bmi, ryger)])
gravdb2=gravdb1[gruppe %in% c(1,8)][order(pnr, slutdato)][,parity:=1:.N, pnr][ryger>=1 & ryger <99, ryger:=1][ryger==-1 | is.na(ryger) | ryger==99, ryger:=NA]
gravdb2[,.N, ryger]
parity=gravdb2[,.(pnr, preg_age=round((slutdato-fdato)/aar, prec), parity, bmi=round(bmi,1), ryger)][bmi<10 | bmi>100, bmi:=NA][, bmi:=nafill(bmi, type="locf"), pnr][, ryger:=nafill(ryger, type="locf"), pnr]

# parity_below45=parity[!(is.na(preg_age))][preg_age<45][order(pnr, -preg_age)][,.SD[1],pnr][,.(pnr, parity)]
# parity_over45=parity[!(is.na(preg_age))][preg_age>=45][,.(pnr, tstart=preg_age, parity45=parity)]
# bmi=parity[!is.na(bmi)][,.(pnr, bmi, preg_age)]

######## 5. Demographic registers
#5.1 Education
udda <- fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/udda2022.csv")

udda1=merge(x=cpr31[,.(pnr, fdato)], y=udda, all.x=T, by="pnr")[is.na(uddniv), `:=` (uddniv=99, uddniv1="Unknown")][,udda_age:=round((hf_vfra-fdato)/aar,prec)][order(pnr, udda_age)][uddniv<99, uddniv2:="0_Primary or lower secondary"][uddniv>=30, uddniv2:="1_Upper secondary"][uddniv>=60, uddniv2:="2_Bachelor or equivalent"][uddniv>=70, uddniv2:="3_Master or equivalent"][uddniv==99, uddniv2:="4_Unknown"][,.SD[1], c("pnr", "uddniv2")][,.(pnr, udda_age, uddniv2)][is.na(udda_age), udda_age:=0]

# 5.2 Marriages and divorces
civ_bef=fread("data/civ_bef.csv")

civ_bef1=unique(merge(x=cpr31[,.(pnr, fdato)], y=civ_bef, all.x=T, by="pnr")[is.na(civst), civst:="NA"][is.na(civ_start), civ_start:=fdato][order(pnr, civst)][civst %in% c("G","F", "U", "NA")])[civst=="U", civst1:="0_Unmarried"][civst=="G", civst1:="1_Married"][civst=="F", civst1:="2_Divorced"][civst=="NA", civst1:="3_Unknown"][,.(pnr, civ_age=round((civ_start-fdato)/aar, prec), civst1)]
nrow(civ_bef1[,.SD[1],pnr])

#5.3 Income at age 45 and 55, adjusting for 2% average inflation rate
# ind=fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/ind9521.csv")
# ind1=merge(x=cpr31[,.(pnr, fdato)], y=ind[,.(pnr=as.integer64(pnr), loenmv_13, perindkialt_13, ind_aar, alder_ult_ink, year_ind=year(ind_aar))], by="pnr", all.x=T)
# fwrite(ind1, "data/ind1.csv")
ind1=fread("data/ind1.csv")[!is.na(perindkialt_13)][, adj_income:=perindkialt_13*(1+0.02)^(2023-year_ind)][, income_quantile:=cut(adj_income, breaks=quantile(adj_income, probs=c(0,0.25,0.5,0.75,1), na.rm=T), include.lowest=T, labels=c("q1", "q2", "q3", "q4"))]

ind45=ind1[alder_ult_ink<=45][order(pnr, -alder_ult_ink)][,.SD[1],pnr][,.(pnr, q_ind45=income_quantile, age=45)]
ind55=ind1[alder_ult_ink<=55][order(pnr, -alder_ult_ink)][,.SD[1],pnr][,.(pnr, q_ind55=income_quantile, age=55)]


###### 6 The great merge
cpr4=merge(x=cpr31, y=thromb_inc, all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=liver_inc, all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=cvd_inc, all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=vte_inc, all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=breastgyn_inc1, all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=ht4[,.(pnr, first_ht=htstart)][,.SD[1], pnr], all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=bso1[,.(pnr, bso_age)], all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=parity[preg_age<45][order(pnr, -preg_age)][,.SD[1],pnr][,.(pnr, parity_bl=parity)], all.x=T, by="pnr")[is.na(parity_bl), parity_bl:=0]
cpr4=merge(x=cpr4, y=ind45, all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=n_diag44[,.(pnr, many_diag)], all.x=T, by="pnr")[is.na(many_diag), many_diag:=0]
#cpr4=merge(x=cpr4, y=epikur3[,.(pnr, many_drugs)], all.x=T, by="pnr")[is.na(many_drugs), many_drugs:=0]
cpr4=merge(x=cpr4, y=udda1[udda_age<45][order(pnr,-udda_age)][,n:=1:.N, pnr][n==1][,.(pnr, udda_bl=uddniv2)], all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=civ_bef1[civ_age<45][order(pnr,-civ_age)][,n:=1:.N, pnr][n==1][,.(pnr, civst_bl=civst1)], all.x=T, by="pnr")
cpr4=merge(x=cpr4, y=dm1[diabetes_age<45][,.(pnr, dm_bl=1)], all.x=T, by="pnr")[is.na(dm_bl), dm_bl:=0]
cpr4=merge(x=cpr4, y=hypertens[ht_age<45][,.(pnr, ht_bl=1)], all.x=T, by="pnr")[is.na(ht_bl), ht_bl:=0]
cpr4=merge(x=cpr4, y=afli_inc[afli_age <45][,.(pnr, afli_bl=1)], all.x=T, by="pnr")[is.na(afli_bl), afli_bl:=0]
cpr4=merge(x=cpr4, y=valv_inc[valv_age <45][,.(pnr, valv_bl=1)], all.x=T, by="pnr")[is.na(valv_bl), valv_bl:=0]
cpr4=merge(x=cpr4, y=hf2[hf_age <45][,.(pnr, hf_bl=1)], all.x=T, by="pnr")[is.na(hf_bl), hf_bl:=0]
cpr4=merge(x=cpr4, y=statin1[statin_age <45][,.(pnr, statin_bl=1)], all.x=T, by="pnr")[is.na(statin_bl), statin_bl:=0]

# 6.1 Numbers for fig. 1
a=format(nrow(cpr4),big.mark = ",")
b=format(nrow(cpr4[thrombinc_age<=start_age]),big.mark = ",")
c=format(nrow(cpr4[liverinc_age<=start_age]),big.mark = ",")
d=format(nrow(cpr4[cvdinc_age<=start_age]),big.mark = ",")
e=format(nrow(cpr4[vteinc_age<=start_age]),big.mark = ",")
f=format(nrow(cpr4[breastgyninc_age<=start_age]),big.mark = ",")
g=format(nrow(cpr4[first_ht<=start_age]),big.mark = ",")
h=format(nrow(cpr4[bso_age<=start_age]),big.mark = ",")

i=format(nrow(cpr4[thrombinc_age<=start_age | liverinc_age<=start_age | cvdinc_age<=start_age | vteinc_age<=start_age | breastgyninc_age<=start_age |first_ht<=start_age | bso_age<=start_age]), big.mark=",")

cpr43=cpr4[(thrombinc_age>start_age | is.na(thrombinc_age)) &
               (liverinc_age>start_age | is.na(liverinc_age)) & 
              (cvdinc_age>start_age | is.na(cvdinc_age)) & 
              (vteinc_age>start_age | is.na(vteinc_age)) & 
              (breastgyninc_age>start_age | is.na(breastgyninc_age)) &
              (first_ht>start_age | is.na(first_ht)) &
            (bso_age>start_age | is.na(bso_age))]
               
j=format(nrow(cpr43),big.mark = ",")

#6.2 Years at splits
age2002=cpr43[,.(pnr, age2002=round((as.IDate("2002-01-01")-fdato)/aar,prec))]
age2009=cpr43[,.(pnr, age2009=round((as.IDate("2009-01-01")-fdato)/aar,prec))]
age2016=cpr43[,.(pnr, age2016=round((as.IDate("2016-01-01")-fdato)/aar,prec))]
age2023=cpr43[,.(pnr, age2023=round((as.IDate("2023-01-01")-fdato)/aar,prec))]

######## 7 Tmerge and analysis prep
cpr44=tmerge(cpr43[,.(pnr, start_age, stop_age, year_start=decimal_date(start), q_ind45, fland, parity_bl, many_diag, udda_bl, civst_bl, dm_bl, ht_bl, afli_bl, valv_bl, hf_bl, statin_bl)], cpr43, id=pnr, tstart=start_age, tstop=stop_age, event=event(death_age))
cpr44=tmerge(cpr44, udda1[udda_age>=45], id=pnr, udda=tdc(udda_age))
cpr44=tmerge(cpr44, civ_bef1[civ_age>=45], id=pnr, civst=tdc(civ_age))
cpr44=tmerge(cpr44, parity[preg_age>=45][,.(pnr, preg_age)], id=pnr, parity=cumtdc(preg_age))
cpr44=tmerge(cpr44, hyst, id=pnr, hyst=tdc(hyst_age))
cpr44=tmerge(cpr44, bso1, id=pnr, bso=tdc(bso_age))
cpr44=tmerge(cpr44, thromb_inc, id=pnr, thromb=tdc(thrombinc_age))
cpr44=tmerge(cpr44, liver_inc, id=pnr, liver=tdc(liverinc_age))
cpr44=tmerge(cpr44, cvd_inc, id=pnr, cvd=tdc(cvdinc_age))
cpr44=tmerge(cpr44, vte_inc, id=pnr, vte=tdc(vteinc_age))
cpr44=tmerge(cpr44, breastgyn_inc1, id=pnr, breastgyn=tdc(breastgyninc_age))
cpr44=tmerge(cpr44, cancer_inc1, id=pnr, cancer=tdc(cancerinc_age))
cpr44=tmerge(cpr44, ht4[,.(pnr, htstart)], id=pnr, htstart=cumtdc(htstart))
cpr44=tmerge(cpr44, ht4[,.(pnr, htstop)], id=pnr, htstop=cumtdc(htstop))
cpr44=tmerge(cpr44, age2002, id=pnr, age2002=tdc(age2002))
cpr44=tmerge(cpr44, age2009, id=pnr, age2009=tdc(age2009))
cpr44=tmerge(cpr44, age2016, id=pnr, age2016=tdc(age2016))
cpr44=tmerge(cpr44, age2023, id=pnr, age2023=tdc(age2023))
cpr44=tmerge(cpr44, ht_pp, id=pnr, htstop90=tdc(htstop90))
cpr44=tmerge(cpr44, dm1, id=pnr, dm=tdc(diabetes_age))
cpr44=tmerge(cpr44, hypertens, id=pnr, hypertens=tdc(ht_age))
cpr44=tmerge(cpr44, afli_inc, id=pnr, afli=tdc(afli_age))
cpr44=tmerge(cpr44, hf2, id=pnr, hf=tdc(hf_age))
cpr44=tmerge(cpr44, valv_inc, id=pnr, valv=tdc(valv_age))
cpr44=tmerge(cpr44, statin1, id=pnr, statin=tdc(statin_age))

cpr5=survSplit(Surv(time=tstart, time2=tstop, event=event) ~., cpr44,
               cut=c(50, 55,60, 65,70, 75), episode ="timegroup")

cpr6=as.data.table(cpr5)[,risktime:=tstop-tstart][,
            n_risktime:=shift(risktime, type="lag", n=1), pnr][
  is.na(n_risktime), n_risktime:=0][,
            cum_riskstart:=cumsum(n_risktime),pnr][,
             cum_riskstop:=cumsum(risktime),pnr][,
            year_stop:=year_start+cum_riskstop][,
              year_start:=year_start+cum_riskstart][, 
                year_cat:=cut(year_start, breaks=c(1994,2002,2009,2016,2024), right=F)]

#7.2 time-dep covariates, rolling merge
#parity
cpr61=cpr6[,parity1:=parity_bl+parity][, parity_cat:=cut(parity1, breaks=c(0,1,3,20), right=F, labels=c("0", "1-2",">=3"))]
#civil status
setkey(civ_bef1, pnr, civ_age)
setkey(cpr61, pnr, tstart)
cpr62=civ_bef1[cpr61, roll=T]
setnames(cpr62, old="civ_age", new="tstart")

test=cpr62[tstart==55][,.(pnr, civst1)]
#udda
setkey(udda1, pnr, udda_age, uddniv2)
setkey(cpr62, pnr, tstart)
cpr63=udda1[cpr62, roll=T]
setnames(cpr63, old="udda_age", new="tstart")
#bmi
bmi=parity[!is.na(bmi)][,.(pnr, bmi, preg_age)]
setkey(bmi, pnr, preg_age)
setkey(cpr63, pnr, tstart)
cpr64=bmi[cpr63, roll=T]
setnames(cpr64, old="preg_age", new="tstart")
cpr64[!is.na(bmi)]
#smoking
ryger=parity[!is.na(ryger)][,.(pnr, ryger, preg_age)]
setkey(ryger, pnr, preg_age)
setkey(cpr64, pnr, tstart)
cpr65=ryger[cpr64, roll=T]
setnames(cpr65, old="preg_age", new="tstart")
cpr65[!is.na(ryger)]
#q_ind55
setkey(ind55, pnr, age, q_ind55)
setkey(cpr65, pnr, tstart)
cpr66=ind55[cpr65, roll=T]
setnames(cpr66, old="age", new="tstart")
cpr67=cpr66[,q_ind:=q_ind45][!is.na(q_ind55), q_ind:=q_ind55]

cpr7=cpr67[,.(pnr, tstart, tstop, event, parity1, hyst, bso, thromb, liver, cvd, vte, breastgyn, risktime, cum_riskstart, cum_riskstop, year_cat, q_ind, fland, civst1, htstart=tstart, timegroup, year_start, htstop90, dm, hypertens, afli, hf, valv, statin, many_diag, bmi, ryger, uddniv2, parity_cat)]

cpr71=merge(x=cpr7, y=ht4[,.(pnr, htstart, ddd, dosform, htname, htatc, httype)], all.x=T, by=c("pnr", "htstart"))

#7.3 Adding cumulative doseform
cpr72=cpr71[is.na(ddd), ddd:=0][,
                  ht_sum:=cumsum(ddd)/aar, pnr][,
                ht_cat:=cut(ht_sum, breaks=c(0, 0.001,1,3,5,10,100), right=F)][,
                    ddd_oral:=ddd][dosform!="oral", ddd_oral:=0][, 
                      ht_sum_oral:=cumsum(ddd_oral)/aar, pnr][,
                    ddd_plaster:=ddd][dosform!="transderm", ddd_plaster:=0][,
                      ht_sum_plaster:=cumsum(ddd_plaster)/aar, pnr][,
                        ddd_other:=ddd][dosform %like% "trans|oral", ddd_other:=0][,
                          ht_sum_other:=cumsum(ddd_other)/aar, pnr][
                      ht_sum==0, ht_dosform:=0][
                        ht_sum_oral> ht_sum_plaster & ht_sum_oral> ht_sum_other, ht_dosform:=1][
                        ht_sum_plaster> ht_sum_oral & ht_sum_plaster> ht_sum_other, ht_dosform:=2][
                        ht_sum_other> ht_sum_oral & ht_sum_other> ht_sum_plaster, ht_dosform:=3][,
                        ht_dosform1:=as.numeric(ht_dosform)][,
                          ht_dosform1:=nafill(ht_dosform1, type="locf", fill=NA)][,
                            ht_dosform1:=as.character(ht_dosform1)][
                            ht_dosform1==1, ht_dosform1:="1oral"][
                            ht_dosform1==2, ht_dosform1:="2transder"][
                            ht_dosform1==3, ht_dosform1:="3other"]#[
                            #ht_dosform1=="1oral" & ht_sum_oral>=2, ht_dosform1:="1oral_>=2"][
                            #ht_dosform1=="2transder" & ht_sum_plaster>=2, ht_dosform1:="2transder_>=2"]

#7.4 Adding cumulative progestogen regimen
cpr73=cpr72[httype=="estrogen+progestin" & htname %like% "Jaydess|Kyleena|Mirena|Levosert", httype:= "estrogen+progestin_monophase"][httype=="estrogen+progestin", httype:="estrogen+progestin_cyclic"][,
           ddd_estrogen:=ddd][httype!="estrogen", ddd_estrogen:=0][, 
             ht_sum_estrogen:=cumsum(ddd_estrogen)/aar, pnr][,
               ddd_cyclic:=ddd][!(httype %like% "cyclic"), ddd_cyclic:=0][,
                 ht_sum_cyclic:=cumsum(ddd_cyclic)/aar, pnr][,
                   ddd_monophase:=ddd][!(httype %like% "monophase"), ddd_monophase:=0][,
                      ht_sum_monophase:=cumsum(ddd_monophase)/aar, pnr][
                        ht_sum==0, ht_regimen:=0][
                          ht_sum_estrogen> ht_sum_cyclic & ht_sum_estrogen> ht_sum_monophase, ht_regimen:=1][
                           ht_sum_cyclic> ht_sum_estrogen & ht_sum_cyclic> ht_sum_monophase, ht_regimen:=2][
                              ht_sum_monophase> ht_sum_estrogen & ht_sum_monophase> ht_sum_cyclic, ht_regimen:=3][,
                                ht_regimen1:=as.numeric(ht_regimen)][,
                                 ht_regimen1:=nafill(ht_regimen1, type="locf", fill=NA)][,
                                  ht_regimen1:=as.character(ht_regimen1)][
                                    ht_regimen1==1, ht_regimen1:="1estrogen"][
                                      ht_regimen1==2, ht_regimen1:="2cyclic"][
                                        ht_regimen1==3, ht_regimen1:="3monopha"][,
                                          ht:=ifelse(ht_sum==0, 0,1)][ht==0, htstop90:=0]

#7.5 Overview of MHT, eTable 2
mht1=cpr73[!is.na(htatc)][,.(risktime=sum(risktime)), c("htatc", "dosform", "httype")][,sum:=sum(risktime)][order(htatc)][,prop:=format(round((risktime/sum)*100,1), nsmall=)][,.(htatc, httype, dosform, prop)]

mht2=cpr73[!is.na(htatc)][,.(htname, htatc, httype, dosform)][,n:=.N, c("htname")][order(-n)][,.SD[1], c("htatc", "httype", "dosform")][,.(htname, dosform, httype, htatc)]

mht3=merge(x=mht1, y=mht2, all.x=T, by=c("htatc", "httype", "dosform"))[,.(htatc, httype, dosform, example_name=htname, prop)][order(htatc, httype, -prop, dosform)][prop==" 0.0", prop:="<0.1"]

write.xlsx(mht3, file="output/etabel2.xlsx")

kbl(mht3)%>%
  kable_classic()

cpr74=cpr73[,.(pnr, htstart, tstart, tstop, event, parity1, hyst, bso, thromb,liver, cvd, vte, breastgyn, risktime,  cum_riskstart,  cum_riskstop, year_cat, q_ind, fland, civst1, timegroup, year_start,htstop90,  dm, hypertens, afli, hf,valv,  statin, many_diag, bmi, ryger, ht_cat, ht_dosform1,ht_regimen1, ht, uddniv2, parity_cat, ht_sum)]


# #7.6 Mode imputation
cpr75=cpr74[is.na(civst1), civst1:="3_Unknown"][,civst_mode:=civst1][civst1=="3_Unknown", civst_mode:="1_Married"]
cpr76=cpr75[is.na(q_ind), q_ind:="q5"][,q_ind_mode:=q_ind][q_ind_mode=="q5", q_ind_mode:="q4"][,q_ind_mode:=factor(q_ind_mode)]
cpr77=cpr76[is.na(uddniv2), uddniv2:="4_Unknown"][,udda_mode:=uddniv2][udda_mode=="4_Unknown", udda_mode:="1_Upper secondary"]
cpr8=cpr77[,fland_mode:=ifelse(fland %like% "01 Danmark" | is.na(fland), "1_Denmark","2_Other")][,fland1:=ifelse(fland %like% "01 Danmark", "1_Denmark","2_Other")][is.na(fland), fland1:="3_Unknown"]

####################################### 8 Analyses       ##########################################
#Adjustment set: year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag
test=cpr8[,.N, c("year_cat", "fland_mode", "civst_mode", "udda_mode", "q_ind_mode", "parity_cat", "dm", "hypertens", "afli", "valv", "hf", "statin", "many_diag")]
#8.0 Primary
#summary(cpr8[,.SD[.N],pnr]$tstop)

rt0=cpr8[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m01=coxph(Surv(tstart, tstop, event)~ht, cpr8)
crude01=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m01)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m02=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8)

adj01=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m02)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

shoen1=cox.zph(m02)
ggcoxzph(shoen1)

rt01=merge(x=rt0, y=crude01, all.x=T, by="ht_cat")
rt02=merge(x=rt01, y=adj01, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a0:ever_never"]

#8.1 Primary with time periods
rt1=cpr8[,.(event=sum(event), risktime=sum(risktime)), c("ht_cat")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat][order(ht_cat)][, ht_cat:=paste0("ht_cat",ht_cat)]

m11=coxph(Surv(tstart, tstop, event)~ht_cat, cpr8)
crude1=clean_names(as.data.table(summary(m11)$conf.int[1:5,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=rn)]

m12=coxph(Surv(tstart, tstop, event)~ht_cat+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8)
adj1=clean_names(as.data.table(summary(m12)$conf.int[1:5,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

schoen1=cox.zph(m12)

rt11=merge(x=rt1, y=crude1, all.x=T, by="ht_cat")
rt12=merge(x=rt11, y=adj1, all.x=T, by="ht_cat")[order(ht_cat1)][,ht_cat1:=NULL][,analysis:="a1:main"]

#8.2 censoring at contrindications cancer, vte, cvd, liver
cpr8_cens=cpr8[thromb==0 & vte==0 & cvd==0 & breastgyn==0 & liver==0]

#prop with contraindication before event
round(nrow(cpr8[event==1 & (thromb | vte==1 | cvd==1 | breastgyn==1 | liver==1)])/nrow(cpr8[event==1])*100,1)

rt2=cpr8_cens[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m21=coxph(Surv(tstart, tstop, event)~ht, cpr8_cens)
crude2=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m21)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m22=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8_cens)
adj2=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m22)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt21=merge(x=rt2, y=crude2, all.x=T, by="ht_cat")
rt22=merge(x=rt21, y=adj2, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a2:censored"]

#8.3 Sens. sibling analysis
sib=cpr8[order(pnr, -tstop)][,.SD[1], pnr][,.(pnr, ht_sum)]

sib1=merge(x=sib, y=cpr[,.(pnr, pnrm, pnrf, fdato)], all.x=T, by="pnr")[!is.na(pnrm)][order(pnrm, pnrf)][,n_mor:=.N, pnrm][n_mor>1][,ht_sib:=ifelse(ht_sum==0, 0, 1)][,ht_sib:=max(ht_sib), pnrm][,no_ht_sib:=ifelse(ht_sum==0, 1, 0)][,no_ht_sib:=max(no_ht_sib), pnrm][ht_sib==1 & no_ht_sib==1][,family_id:=.GRP, pnrm][order(pnrm, fdato)][,sib_num:=1:.N, family_id]

sib2=merge(x=cpr8, y=sib1[,.(pnr, family_id, sib_num, fdato)], all.y=T, by="pnr")

sib3=sib2[,.SD[.N],pnr][order(family_id)]

nrow(sib1[,.N, pnr])
round(nrow(sib1[,.N, pnr])/nrow(cpr8[,.N, pnr])*100,1)
rt3=sib2[,.(event=sum(event),risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][, ht_cat:=paste0("ht_cat",ht)]

m31=coxph(Surv(tstart, tstop, event)~ht, sib2)
crude3=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m31)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m32=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag+strata(family_id), sib2)

adj3=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m32)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt31=merge(x=rt3, y=crude3, all.x=T, by="ht_cat")
rt32=merge(x=rt31, y=adj3, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a3:sibling"]

#8.4 Sens censoring at MHT cessation (as treated)
cpr_pp=cpr8[htstop90==0]#[,.(pnr, tstart, tstop, event, ht, ht_cat, htstop90)]

# test0=cpr_pp[event==1 & ht==0][,.(pnr, small_ds=1)]
# test1=cpr8[event==1 & ht==0]
# test2=merge(x=test0, y=test1, all.y=T, by="pnr")
# test2[is.na(small_ds)]

rt4=cpr_pp[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m41=coxph(Surv(tstart, tstop, event)~ht, cpr_pp)
crude4=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m41)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m42=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr_pp)
adj4=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m42)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt41=merge(x=rt4, y=crude4, all.x=T, by="ht_cat")
rt42=merge(x=rt41, y=adj4, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a4:per_protokol"]

#8.5 Strat model dosage type: table, plaster/gel, other
rt5=cpr8[,.(event=sum(event), risktime=sum(risktime)), c("ht_dosform1")][, ir:=round(event*10000/risktime,1)][,ht_cat:=paste0("ht_dosform1", ht_dosform1)][order(ht_dosform1)]

m51=coxph(Surv(tstart, tstop, event)~ht_dosform1, cpr8)
crude5=clean_names(as.data.table(summary(m51)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=rn)]

m52=coxph(Surv(tstart, tstop, event)~ht_dosform1+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8)
adj5=clean_names(as.data.table(summary(m52)$conf.int[1:5,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt51=merge(x=rt5, y=crude5, all.x=T, by="ht_cat")
rt52=merge(x=rt51, y=adj5, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a5:dosform"]

#8.6 stratified by regimen: cyclic, continuous, mono
rt6=cpr8[,.(event=sum(event), risktime=sum(risktime)), c("ht_regimen1")][, ir:=round(event*10000/risktime,1)][,ht_cat:=paste0("ht_regimen1", ht_regimen1)][order(ht_regimen1)]

m61=coxph(Surv(tstart, tstop, event)~ht_regimen1, cpr8)
crude6=clean_names(as.data.table(summary(m61)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=rn)]

m62=coxph(Surv(tstart, tstop, event)~ht_regimen1+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8)
adj6=clean_names(as.data.table(summary(m62)$conf.int[1:3,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt61=merge(x=rt6, y=crude6, all.x=T, by="ht_cat")
rt62=merge(x=rt61, y=adj6, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a6:regimen"]

#8.7 Stratified by hysterectomy
cpr8_hyst=cpr8[breastgyn==0][hyst==1][, ht_cat2:=cut(ht_sum, breaks=c(0, 0.001,5,100), right=F)]
rt7=cpr8_hyst[,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2", "hyst")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(hyst, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]#[hyst==1, ht_cat2:=paste0(ht_cat2, ":hyst")]

m71=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr8_hyst)
crude7=clean_names(as.data.table(summary(m71)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m72=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8_hyst)
adj7=clean_names(as.data.table(summary(m72)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt71=merge(x=rt7, y=crude7, all.x=T, by="ht_cat2")
rt72=merge(x=rt71, y=adj7, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="a7:strat_hysterectomy"][order(hyst, ht_cat1)][,ht_cat1:=NULL][,hyst:=NULL]
setnames(rt72, "ht_cat2", "ht_cat")

#8.8 Restricted by oofrectomy 45-54 years
cpr8_bso=merge(x=cpr8, y=bso1[,.(pnr, bso_age)], all.x=T, by="pnr")[bso_age>=55, bso:=2][,bso:=factor(bso)][, ht_cat2:=cut(ht_sum, breaks=c(0,0.001,5,100), right=F)][breastgyn==0]

rt8=cpr8_bso[bso=="1"][,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2", "bso")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(bso, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]#[bso==1, ht_cat2:=paste0(ht_cat2, ":bso1")][bso==2, ht_cat2:=paste0(ht_cat2, ":bso2")]

m81=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr8_bso[bso==1])
crude8=clean_names(as.data.table(summary(m81)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m82=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8_bso[bso==1])
adj8=clean_names(as.data.table(summary(m82)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt81=merge(x=rt8, y=crude8, all.x=T, by="ht_cat2")
rt82=merge(x=rt81, y=adj8, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="a8:bso<55"][order(bso, ht_cat1)][,ht_cat1:=NULL][,bso:=NULL]
setnames(rt82, "ht_cat2", "ht_cat")

nrow(cpr8_bso[bso==1 & event==1])
summary(cpr8_bso[bso==1 & event==1 & ht==1]$tstop)
summary(cpr8_bso[bso==1 & event==1 & ht==0]$tstop)


#8.9 oophorectomy >=55 years
rt9=cpr8_bso[bso=="2"][,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2", "bso")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(bso, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]#[bso==1, ht_cat2:=paste0(ht_cat2, ":bso1")][bso==2, ht_cat2:=paste0(ht_cat2, ":bso2")]

m91=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr8_bso[bso==2])
crude9=clean_names(as.data.table(summary(m91)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m92=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8_bso[bso==2])
adj9=clean_names(as.data.table(summary(m92)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt91=merge(x=rt9, y=crude9, all.x=T, by="ht_cat2")
rt92=merge(x=rt91, y=adj9, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="a9:bso>=55"][order(bso, ht_cat1)][,ht_cat1:=NULL][,bso:=NULL]
setnames(rt92, "ht_cat2", "ht_cat")

#8.10 Stratified by age of initiation
ht_age=merge(x=cpr8, y=ht4[,.SD[1],pnr][,.(pnr, first_ht=htstart)], by="pnr", all.x=T)[tstart<first_ht |is.na(first_ht), first_ht:=0][, first_ht:=cut(first_ht, breaks=c(0,0.1,52, 57, 150), right=F)]#[,.(pnr, tstart, tstop, event, first_ht, htatc, htname, year_start, year_cat)]

ht_age=cpr8[ht==1, first_ht:=tstart][,first_ht1:=min(first_ht, na.rm=T), pnr][is.na(first_ht), first_ht1:=NA][is.na(first_ht1), first_ht1:=0][, first_ht:=cut(first_ht1, breaks=c(0,0.1,52, 57, 150), right=F)]

rt10=ht_age[,.(event=sum(event), risktime=sum(risktime)), c("first_ht")][, ir:=round(event*10000/risktime,1)][order(first_ht)][, first_ht:=paste0("first_ht",first_ht)]

m101=coxph(Surv(tstart, tstop, event)~first_ht, ht_age)
crude10=clean_names(as.data.table(summary(m101)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), first_ht=rn)]

m102=coxph(Surv(tstart, tstop, event)~first_ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, ht_age)
adj10=clean_names(as.data.table(summary(m102)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), first_ht=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt101=merge(x=rt10, y=crude10, all.x=T, by="first_ht")
rt102=merge(x=rt101, y=adj10, all.x=T, by="first_ht")[order(first_ht)][,analysis:="a10:sub_first_hc"][order(first_ht)]
setnames(rt102, "first_ht", "ht_cat")

#8.11 reviewer sens, stratified by year 2002
ht2002=cpr8[year_start<2003]#[year_cat=="[1994,2002)"][vte==0 & cvd==0 & breastgyn==0 & liver==0]

rt2002a=ht2002[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m111=coxph(Surv(tstart, tstop, event)~ht, ht2002)
crude02a=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m111)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m112=coxph(Surv(tstart, tstop, event)~ht+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, ht2002)#2025check

adj2002a=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m112)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt2002a1=merge(x=rt2002a, y=crude02a, all.x=T, by="ht_cat")
rt2002a2=merge(x=rt2002a1, y=adj2002a, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a11:review_<=2002"]

ht2002_1=cpr8[year_start>=2003]#[year_cat!="(1994,2002]"]

rt2002b=ht2002_1[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m121=coxph(Surv(tstart, tstop, event)~ht, ht2002_1)
crude02b=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m121)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m122=coxph(Surv(tstart, tstop, event)~ht+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, ht2002_1)#2025check

adj2002b=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m122)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt2002b1=merge(x=rt2002b, y=crude02b, all.x=T, by="ht_cat")
rt2002b2=merge(x=rt2002b1, y=adj2002b, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a11:review_>2002"]


#8.10 reviewer sens, adj for hysterectomy and oophorectomy
m132=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag+hyst+bso, cpr8)#2025check

adj13=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m132)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt132=merge(x=rt01, y=adj13, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a12:reviewer_adj_hyst_bso"]

#8.11 reviewer sens, stratfied by hysterectomy 45-54
cpr8_hyst_strat=merge(x=cpr8, y=hyst[,.(pnr, hyst_age)], all.x=T, by="pnr")[hyst_age>=55, hyst:=2][,hyst:=factor(hyst)][, ht_cat2:=cut(ht_sum, breaks=c(0,0.001,5,100), right=F)][breastgyn==0]

rt14=cpr8_hyst_strat[hyst=="1"][,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2", "hyst")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(hyst, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]#[bso==1, ht_cat2:=paste0(ht_cat2, ":bso1")][bso==2, ht_cat2:=paste0(ht_cat2, ":bso2")]

m141=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr8_hyst_strat[hyst==1])
crude14=clean_names(as.data.table(summary(m141)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m142=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8_hyst_strat[hyst==1])
adj14=clean_names(as.data.table(summary(m142)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt141=merge(x=rt14, y=crude14, all.x=T, by="ht_cat2")
rt142=merge(x=rt141, y=adj14, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="a13:hyst<55"][order(bso, ht_cat1)][,ht_cat1:=NULL][,hyst:=NULL]
setnames(rt142, "ht_cat2", "ht_cat")

#8.12 hysterectomy >=55 years
rt15=cpr8_hyst_strat[hyst=="2"][,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2", "hyst")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(hyst, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]#[bso==1, ht_cat2:=paste0(ht_cat2, ":bso1")][bso==2, ht_cat2:=paste0(ht_cat2, ":bso2")]

m151=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr8_hyst_strat[hyst==2])
crude15=clean_names(as.data.table(summary(m151)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m152=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr8_hyst_strat[hyst==2])
adj15=clean_names(as.data.table(summary(m152)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt151=merge(x=rt15, y=crude15, all.x=T, by="ht_cat2")
rt152=merge(x=rt151, y=adj15, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="a14:hyst>=55"][order(bso, ht_cat1)][,ht_cat1:=NULL][,hyst:=NULL]
setnames(rt152, "ht_cat2", "ht_cat")


#8.13 cc bmi
cpr8bmi=cpr8[!is.na(bmi)][,bmi_cat:=cut(bmi, breaks=c(-Inf,30,Inf), right=F)]
rt16=cpr8bmi[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m161=coxph(Surv(tstart, tstop, event)~ht, cpr8bmi)
crude16=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m161)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m162=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag+bmi_cat, cpr8bmi)

adj16=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m162)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt161=merge(x=rt16, y=crude16, all.x=T, by="ht_cat")
rt162=merge(x=rt161, y=adj16, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a15:bmi_cc"]


#8.14 cc smoking
cpr8ryg=cpr8[!is.na(ryger)]
rt17=cpr8ryg[,.(event=sum(event), risktime=sum(risktime)), c("ht")][, ir:=round(event*10000/risktime,1)][order(ht)][,ht_cat:=paste0("ht_cat", ht)]

m171=coxph(Surv(tstart, tstop, event)~ht, cpr8bmi)
crude17=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m171)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat=paste0("ht_cat", rn))]

m172=coxph(Surv(tstart, tstop, event)~ht+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag+ryger, cpr8bmi)

adj17=clean_names(as.data.table(cbind("rn"=c(0,1),rbind(1,summary(m172)$conf.int[,c(1,3:4)])), keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat="ht_cat1", adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)][2]

rt171=merge(x=rt17, y=crude17, all.x=T, by="ht_cat")
rt172=merge(x=rt161, y=adj17, all.x=T, by="ht_cat")[order(ht_cat)][,analysis:="a16:smoking_cc"]


#8.11 combine analyses
rt=rbind(rt02, rt12, rt22, rt32, rt42, rt72, rt82, rt92, rt52, rt62, rt102, rt2002a2, rt2002b2, rt132, rt142, rt152, rt162, rt172, fill=T)[is.na(adj_est), `:=` (adj_est=1, adj=1, crude=1)][event<3, event:=NA][,.(analysis, ht_cat, event, risktime, ir, crude, adj)]

write.xlsx(rt, "output/table2.xlsx")

ptab=rt[,.(analysis,"Years of MHT use"=revalue(ht_cat,c("ht_cat[0,0.001)"="0 years", "ht_cat[0.001,1)"="<1 year", "ht_cat[1,3)"="1-3 years", "ht_cat[3,5)"="3-5 years", "ht_cat[5,10)"="5-10 years","ht_cat[10,100)"=">10 years", "first_ht[0,1)"="No MHT", "first_ht[1,52)" = "45-51 years", "first_ht[52,57)"="52-56 years", "first_ht[57,150)"=">=57 years",  "ht_cat0"="Never use", "ht_cat1"="Past or present use", "ht_cat2[0,0.001)"="0 years", "ht_cat2[0.001,5)"="<5 years", "ht_cat2[5,100)"=">5 years", "ht_dosform10"="No MHT", "ht_dosform11oral"="Oral",  "ht_dosform12transder"="Transdermal", "ht_dosform13other"="Other formulation", "ht_regimen10"="No MHT", "ht_regimen11estrogen"="Estrogen monopherapy", "ht_regimen12cyclic" ="Estrogen+cyclic progestin", "ht_regimen13monopha" ="Estrogen+continuous progestin")), "Mortality events"=format(event, big.mark=","), "Person years"=format(round(risktime),big.mark=","), "Incience rate"=format(ir, small.mark="."), "Crude HR (95% CI)"=crude, "Adjusted HR (95% CI)*"=adj, " "="\t\t\t")]

write.xlsx(ptab, file="output/t2_word.xlsx")


###################### 9 Cause specific mortality
#9.1 import
event_pop=cpr8[year_start < decimal_date(as.IDate("2023-01-01")) & event==1][,.(pnr, death_age=tstop)]
dar=clean_names(fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/dodsaarsag2023.csv"))[,d_dodsdto:=as.IDate(mdy(d_dodsdto))][,.(pnr=k_cpr, tilgrund=c_dod1, dod2=c_dod2, dod3=c_dod3, ddato=d_dodsdto, c_liste_49)]

dar1=clean_names(fread("E:/workdata/704766/personlige_mapper/anders/regstre_csv1/dodsaarsag2022_2.csv"))[,d_dodsdto:=as.IDate(mdy(d_statdato))][,.(pnr=k_cpr, tilgrund=c_dodtilgrundl_acme, ddato=d_dodsdto, c_listea, c_listeb)]

dar2=rbind(dar, dar1, fill=T)[order(pnr, tilgrund)][,.SD[1], pnr][c_liste_49 %in% c(3,4,5,6,7,8,9,10,11,13,14), csm:="cancer"][c_listea=="A-02", csm:="cancer"][c_listea %in% c("A-08", "A-09"), csm:="cvd"][c_liste_49 %in% c(23,24,25,26,27,28), csm:="cvd"][is.na(csm), csm:="other"][,.(pnr, csm)]
event_pop1=unique(merge(x=event_pop, y=dar2, all.x=T, by="pnr")[is.na(csm), csm:="other"])

#10.2 cvd
cpr20=merge(x=cpr8, y=event_pop1, all.x=T, by="pnr")[csm!="cvd", event:=0][, ht_cat2:=cut(ht_sum, breaks=c(0, 0.001,5,100), right=F)][year_start < decimal_date(as.IDate("2023-01-01"))]

rt20=cpr20[,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(hyst, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]

m201=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr20)
crude20=clean_names(as.data.table(summary(m201)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m202=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+many_diag, cpr20)
adj20=clean_names(as.data.table(summary(m202)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt201=merge(x=rt20, y=crude20, all.x=T, by="ht_cat2")
rt20_csm=merge(x=rt201, y=adj20, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="csm1:cvd"]

#10.3 cancer
cpr21=merge(x=cpr8, y=event_pop1, all.x=T, by="pnr")[csm!="cancer", event:=0][, ht_cat2:=cut(ht_sum, breaks=c(0, 0.001,5,100), right=F)][year_start < decimal_date(as.IDate("2023-01-01"))]

rt21=cpr21[,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(hyst, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]

m211=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr21)
crude21=clean_names(as.data.table(summary(m211)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m212=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr21)
adj21=clean_names(as.data.table(summary(m212)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt211=merge(x=rt21, y=crude21, all.x=T, by="ht_cat2")
rt21_csm=merge(x=rt211, y=adj21, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="csm2:cancer"]

#10.4 other
cpr22=merge(x=cpr8, y=dar2, all.x=T, by="pnr")[csm!="other", event:=0][, ht_cat2:=cut(ht_sum, breaks=c(0, 0.001,5,100), right=F)][year_start < decimal_date(as.IDate("2023-01-01"))]

rt22=cpr22[,.(event=sum(event), risktime=sum(risktime)), c("ht_cat2")][, ir:=round(event*10000/risktime,1)][,ht_cat1:=ht_cat2][order(hyst, ht_cat2)][, ht_cat2:=paste0("ht_cat2",ht_cat2)]

m221=coxph(Surv(tstart, tstop, event)~ht_cat2, cpr22)
crude22=clean_names(as.data.table(summary(m221)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(crude=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn)]

m222=coxph(Surv(tstart, tstop, event)~ht_cat2+year_cat+fland_mode+civst_mode+udda_mode+q_ind_mode+parity_cat+dm+hypertens+afli+valv+hf+statin+many_diag, cpr22)
adj22=clean_names(as.data.table(summary(m92)$conf.int[,c(1,3:4)], keep.rownames = T))[,.(adj=paste0(format(round(exp_coef,2),nsmall=2), " (", format(round(lower_95,2),nsmall=2), "-", format(round(upper_95,2),nsmall=2), ")"), ht_cat2=rn, adj_est=exp_coef, adj_low=lower_95, adj_up=upper_95)]

rt221=merge(x=rt22, y=crude22, all.x=T, by="ht_cat2")
rt22_csm=merge(x=rt221, y=adj22, all.x=T, by="ht_cat2")[order(ht_cat1)][,analysis:="csm3:other"]
rt_csm=rbind(rt20_csm, rt21_csm, rt22_csm)[is.na(adj_est), `:=` (adj_est=1, adj=1, crude=1)][event<3, event:=NA]

write.xlsx(rt_csm, "output/table3.xlsx")

ptab_csm=rt_csm[,.(analysis, "Years of MHT use"=revalue(ht_cat2,c("ht_cat2[0,0.001)"="0 years", "ht_cat2[0.001,5)"="<5 years", "ht_cat2[5,100)"="\u22655 years")), "Mortality events"=format(event, big.mark=","), "Person years"=format(round(risktime),big.mark=","), "Incience rate"=ir, "Crude HR\n(95% CI)"=crude, "Adjusted HR\n(95% CI)*"=adj, " "="\t\t")]

write.xlsx(ptab_csm, file="output/t3_word.xlsx")

round(nrow(cpr20[event==1])/(nrow(cpr20[event==1])+nrow(cpr21[event==1])+nrow(cpr22[event==1]))*100,1)
round(nrow(cpr21[event==1])/(nrow(cpr20[event==1])+nrow(cpr21[event==1])+nrow(cpr22[event==1]))*100,1)
round(nrow(cpr22[event==1])/(nrow(cpr20[event==1])+nrow(cpr21[event==1])+nrow(cpr22[event==1]))*100,1)
          

######################    9 Descriptive stats
#9.1 Table 1


last=cpr8[order(pnr, -tstop)][,.SD[1], pnr]
last[,.N, ht_cat]
k=format(nrow(last[ht_cat=="[0,0.001)"]),big.mark = ",")
l=format(nrow(last[ht_cat=="[0.001,1)"]),big.mark = ",")
m=format(nrow(last[ht_cat=="[1,3)"]),big.mark = ",")
n=format(nrow(last[ht_cat=="[3,5)"]),big.mark = ",")
o=format(nrow(last[ht_cat=="[5,10)"]),big.mark = ",")
p=format(nrow(last[ht_cat=="[10,100)"]),big.mark = ",")

#number for result text
#median follow-up
round(quantile(cpr8[,.(sum(risktime)), pnr]$V1),1)

last1=merge(x=last, y=cpr[,.(pnr, fyear=year(fdato))], all.x=T, by="pnr")[,ht:=ifelse(ht_sum==0, 0 ,1)]
nrow(last1[ht==1])
round((nrow(last1[ht==1])/nrow(last1))*100,1)
nrow(last1[event==1])
round((nrow(last1[event==1])/nrow(last1))*100,1)
round(mean(last1[fyear==1950]$ht)*100,1)
round(mean(last1[fyear==1960]$ht)*100,1)
round(mean(last1[fyear==1970]$ht)*100,1)
round(summary(last1[ht_sum!=0]$ht_sum),1)

#prop in sib cohort
nrow(sib3)/nrow(last1)*100

plot(density(last1[ht_sum!=0]$ht_sum))

#number who only pick up 1 prescription for reviewer
# one_pre=ht4[,n_pre:=.N, pnr][n_pre==1][, .(pnr, one_pre=1)]
# one_pre=merge(x=last[ht==1], y=one_pre, all.x=T, by="pnr")[is.na(one_pre), one_pre:=0]
# nrow(one_pre[one_pre==1])/nrow(one_pre)
# 
# statin=lmdb2022[atc %like% "^C10A[A-B]"][,n_pre:=.N,pnr][order(pnr, eksd)][,.SD[1],pnr]
# nrow(statin[n_pre==1])/nrow(statin)

#9.2 Flowchart, fig 1
f1=create_graph()%>%
  add_node(label=paste(a,"Danish women born between 1950-1977 alive at age 45"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=0,
             y=5,
             fixedsize = F)
           )%>%
  add_node(label="",
           node_aes = node_aes(
             height=0.02,
             width=0.02,
             x=0,
             y=4.25))%>%
  add_node(label=paste(j, "Women included in study"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=0,
             y=3.5,
             fixedsize=F))%>%
  add_node(label=paste("Excluded due to\n",
                  b, "due to thrombophilia\n",
                  c, "due to liver disease\n",
                  d, "due to arterial thrombosis\n",
                  e, "due to venous thrombosis\n",
                  f, "due to breast, endometrial, or ovarian cancer\n",
                  g, "due to previous menopausal hormone therapy\n",
                  h, "due to previous bilateral oophorectomy\n",
                  i, "Women excluded"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=4,
             y=4.25,
             fixedsize=F))%>%
  add_node(label="",
           node_aes = node_aes(
             height=0.02,
             width=0.02,
             x=0,
             y=2.75))%>%
  add_node(label=paste(k, "Women with\n no MHT exposure"), 
                                      node_aes = node_aes(
                                        shape = "rectangle",
                                        x=-2,
                                        y=2.25,
                                        fixedsize=F))%>%
  add_node(label=paste(l,"Women with\n <1 years of MHT exposure"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=2.5,
             y=2.75,
             fixedsize=F))%>%
  add_node(label=paste(m, "Women with\n 1-2.9 years of MHT exposure"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=2.5,
             y=2,
             fixedsize=F))%>%
  add_node(label=paste(n, "Women with\n 3-4.9 years of MHT exposure"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=2.5,
             y=1.25,
             fixedsize=F))%>%
  add_node(label=paste(o, "Women with\n 5-9.9 years of MHT exposure"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=2.5,
             y=0.50,
             fixedsize=F))%>%
  add_node(label=paste(p, " Women with\n \u226510 years of MHT exposure"), 
           node_aes = node_aes(
             shape = "rectangle",
             x=2.5,
             y=-0.25,
             fixedsize=F))%>%
  add_node(label="At end of follow-up", 
           node_aes = node_aes(
             height=0.02,
             width=0.02,
             x=-0.7,
             y=2.9))%>%
  add_edge(from=1, to=3,
           edge_aes=edge_aes(headport="n"))%>%
  add_edge(from=2, to=4,
           edge_aes=edge_aes(headport="w"))%>%
  add_edge(from=3, to=5,
           edge_aes=edge_aes(headport="e",
                             arrowhead="none"))%>%
  add_edge(from=5, to=6,
           edge_aes=edge_aes(headport="w"))%>%
  add_edge(from=5, to=7)%>%
  add_edge(from=5, to=8)%>%
  add_edge(from=5, to=9)%>%
  add_edge(from=5, to=10)%>%
  add_edge(from=5, to=11)

render_graph(f1)

export_graph(f1, "output/fig1.pdf")
export_graph(f1, "output/fig1.png")
f1_svg=export_svg(f1)

rsvg(charToRaw(f1_svg), file="output/fig1.jpg")
rsvg::
#cpr8=merge(x=cpr8, y=bef2[,.(pnr, fland)], by="pnr", all.x=T)

ff=cpr8[tstart==55][,ht:=ifelse(ht_sum==0, "No_MHT", "MHT")]


table1=print(CreateTableOne(ff, vars=c("fland1", "parity1", "parity_cat", "civst1", "uddniv2", "q_ind_mode", "year_cat", "dm", "hypertens", "statin", "afli", "hf", "valv", "many_diag",  "hyst","bso"), factorVars=c("fland1", "hyst", "parity_cat", "bso", "civst1", "uddniv2", "q_ind_mode", "year_cat",  "dm", "hypertens", "statin", "afli", "hf", "valv", "many_diag"), strata="ht", test=F))

write.xlsx(as.data.table(table1, keep.rownames = T), file="output/table1.xlsx")

ff1=ff[,year_cat1:=cut(year_start, breaks=c(2004, 2007, 2021, 2024), right=F)]
table(ff1$year_cat1, ff1$ht)
(nrow(ff1[ht=="MHT" & year_cat1=="[2004,2007)"])/nrow(ff1[year_cat1=="[2004,2007)"]))*100
(nrow(ff1[ht=="MHT" & year_cat1=="[2021,2024)"])/nrow(ff1[year_cat1=="[2021,2024)"]))*100
ff1[,.N, year_cat1]

#Missing
(nrow(last1[uddniv2=="4_Unknown"])/nrow(last1))*100
(nrow(last1[civst1=="3_Unknown"])/nrow(last1))*100
(nrow(last1[fland1=="3_Unknown"])/nrow(last1))*100
(nrow(last1[q_ind=="q5"])/nrow(last1))*100

##############################     END     ##############
