-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path6.asthma_dx.R
124 lines (95 loc) · 3.65 KB
/
6.asthma_dx.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
library(tidyverse)
library(mgcv)
# goal: estimate the "un"diagnosis prob,
# i.e., the prob that an asthma patient "outgrows" asthma and
# become undiagnosed with asthma
#iterative method
chosen_province <- "CA"
starting_year <- 1999
end_year <- 2065
stabilization_year <- 2025
asthma_tuner <- function(chosen_province="CA",
starting_year= 1999,
end_year = 2065,
stabilization_year=2025){
asthma_inc_model <- read_rds("asthma_incidence_model.rds")
asthma_prev_model <- read_rds("asthma_prevalence_model.rds")
asthma_max_age <- 62
asthma_predictor <- function(age,sex,year,type){
age <- pmin(age,asthma_max_age)
year <- pmin(year,stabilization_year)
if(type == "prev"){
return( exp(predict(asthma_prev_model,newdata=data.frame(age,sex,year))) %>%
unlist())
} else{
return( exp(predict(asthma_inc_model,newdata=data.frame(age,sex,year))) %>%
unlist())
}
}
xs <- expand.grid(age=3:110,sex=c(0,1),year=starting_year:end_year) %>%
as.data.frame()
df_asthma <- xs %>%
mutate(inc = asthma_predictor(age,sex,year,type='inc'),
prev = asthma_predictor(age,sex,year,type='prev')) %>%
# let inc be prev for age = 3
mutate(inc = ifelse(age==3,prev,inc))
df_asthma_year <- df_asthma %>%
group_split(year)
results_assessment <- c()
years <- c((starting_year+1):end_year)
ages <- c(4:110)
for( i in 1:(length(years)-1)){
tmp_year <- years[i]
tmp_prev_past <- df_asthma_year[[i]] %>%
select(age,year,sex,prev) %>%
filter(age!=110) %>%
pivot_wider(names_from=sex,values_from=prev) %>%
select(-age,-year)
tmp_prev <- df_asthma_year[[i+1]] %>%
select(age,year,sex,prev) %>%
filter(age!=3) %>%
pivot_wider(names_from=sex,values_from=prev) %>%
select(-age,-year)
tmp_inc <- df_asthma_year[[i+1]] %>%
select(age,year,sex,inc) %>%
filter(age!=3) %>%
pivot_wider(names_from=sex,values_from=inc) %>%
select(-age,-year)
# prev = prev_past*p(reassessment) + inc*(1-prev_past)
# p(reassessment) = (prev-inc*(1-prev_past))/prev_past
# estimate tmp_u by assuming p(correct diagnosis) = 1
tmp_assessment <- (tmp_prev - tmp_inc * (1-tmp_prev_past))/tmp_prev_past
# # key step: bound it!
# tmp_u[tmp_u>1] <- 1
# tmp_u[tmp_u<0] <- 0
results_assessment[[i]] <- cbind(data.frame(age=4:110,year=tmp_year),tmp_assessment)
}
results_assessment <- results_assessment %>%
do.call(rbind,.) %>%
as.data.frame()
colnames(results_assessment)[c(3,4)] <- c("F","M")
results_assessment$province <- chosen_province
return(results_assessment)
}
CA_tuner <- asthma_tuner(chosen_province="CA",end_year=2066,stabilization_year = 2025)
BC_tuner <- asthma_tuner(chosen_province = "BC",end_year=2043,stabilization_year=2025)
# look <- results_assessment %>%
# do.call(rbind,.) %>%
# pivot_longer(3:4,names_to="sex",values_to="reassessment") %>%
# as.data.frame()
#
# tmp <- gam(log(reassessment)~sex + s(year,age),data=look)
#
# ggplot(data = look %>% filter(age==12),
# aes(y=reassessment,x=year)) +
# geom_line()+
# ylim(0,1.1)+
# facet_grid(.~sex)
master_assessment <- rbind(CA_tuner,
BC_tuner)
master_assessment %>%
mutate(M = ifelse(M>1,1,M))
master_assessment <- master_assessment %>% select(year,age,`F`,M,province)
write_csv(master_assessment %>%
mutate(M = ifelse(M>1,1,M),
`F` = ifelse(`F`>1,1,`F`)),"../src/processed_data/master_asthma_reassessment.csv")