-
Notifications
You must be signed in to change notification settings - Fork 2
/
demographics_kern.R
83 lines (75 loc) · 3.23 KB
/
demographics_kern.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
demographics_kern = function(agesex,
dailycount,
country="ALL",
cohorts,
sev="pts_all",
site.register){
expit = function(mm){mm=as.numeric(mm);return(exp(mm)/(exp(mm)+1))}
if(country=="ALL"){site.sel=site.register$SiteID
} else{site.sel=site.register$SiteID[site.register$Country==country]}
dailycount=filter(dailycount, siteid %in% site.sel, cohort %in% cohorts)
count=NULL
for(ss in site.sel){
junk=dailycount[dailycount$siteid==ss,]
date.keep = max(as.Date(junk$calendar_date))
count = rbind.data.frame(count,
dailycount[dailycount$siteid==ss&dailycount$calendar_date==date.keep,])
}
if(grepl("pts_all",sev,fixed = T)){count=count[,1:4]
}else if(grepl("pts_ever_severe",sev,fixed = T)){count=count[,c(1:3,7)]
}else if(grepl("pts_never_severe",sev,fixed = T)){count=count[,c(1:3,14)]}
count$sc=paste0(count$siteid,count$cohort)
count=count[order(count$sc),]
rownames(count)=c()
agesex=filter(agesex,siteid%in%site.sel,cohort %in% cohorts,mean_age>18)
agesex$sc=paste0(agesex$siteid,agesex$cohort)
agesex=merge(agesex,count,by="sc",all.x=TRUE)
age=filter(agesex,sex=="all",age_group!="all")
sex=filter(agesex,sex!="all",age_group=="all")
### for age groups first
res=NULL
for(aa in unique(age$age_group)){
tryCatch({
junk=filter(age,age_group==aa)
junk.plo=escalc(xi=junk[,sev],
ni=junk[,colnames(count)[4]],
measure="PLO")
rma.res = rma(yi=junk.plo$yi,
vi=junk.plo$vi,
method="DL",
drop00=TRUE)
res=rbind.data.frame(res,
cbind.data.frame("country"=country,
"cohort"=paste(cohorts,collapse = ''),
"setting"=sev,
"subgroup"=aa,
"p"=expit(rma.res$b),
"se"=rma.res$se,
"ci.95L"=expit(rma.res$ci.lb),
"ci.95U"=expit(rma.res$ci.ub)))
},error=function(e){NA})
}
for(ss in unique(sex$sex)){
tryCatch({
junk=filter(sex,sex==ss)
junk.plo=escalc(xi=junk[,sev],
ni=junk[,colnames(count)[4]],
measure="PLO")
rma.res = rma(yi=junk.plo$yi,
vi=junk.plo$vi,
method="DL",
drop00=TRUE)
res=rbind.data.frame(res,
cbind.data.frame("country"=country,
"cohort"=paste(cohorts,collapse = ''),
"setting"=sev,
"subgroup"=ss,
"p"=expit(rma.res$b),
"se"=rma.res$se,
"ci.95L"=expit(rma.res$ci.lb),
"ci.95U"=expit(rma.res$ci.ub)))
},error=function(e){NA})
}
rownames(res)=c()
return(res)
}