Data cleaning, management, and analysis of continuity of essential health services

This provides steps and code for data cleaning, management, and analysis of Continuity of Essential Health Services assessment survey. See here for the questionnare. (Standard module version as of 5/12/2021)

This code:
1. Imports and cleans Continuity of Essential Health Services (aka “EHS”) dataset from Lime Survey (i.e., green and blue tabs in the chartbook),
2. Creates field check tables for data quality monitoring, and
3. Creates indicator estimate data for dashboards and chartbook (purple tab in the chartbook).

FOUR parts must be updated, based on the minmum country-specific adaptation.
1. Directories and local macro in A. SETTING per survey implementation information
2. Local macro in E.1. Country speciic code local per section 1
3. Staff variables in E.2. Construct analysis variables per section 2
4. Option for sampling weights in E.3. Merge with sampling weight per sentinel facility selection design

Depending on the extent of additional country-specific adaptation, further revise E.2. Construct analysis variables.

See Section E.2-Addendum (added on August 17, 2021)

(Last updated: 2022-05-31 15:31:32)


A. SETTING

#SET WORKING DIRECTORIES where this markdown file and subfolders are located

setwd("~/Dropbox/0 iSquared/iSquared_WHO/ACTA/3.AnalysisPlan/")
##setwd("C:/Users/YoonJoung Choi/World Health Organization/BANICA, Sorin - HSA unit/2 Global goods & tools/2 HFAs/1 HFAs for COVID-19/4. Implementation support materials/4. Analysis and dashboards/")

chartbookdir<-("~/Dropbox/0 iSquared/iSquared_WHO/ACTA/3.AnalysisPlan/")
##chartbookdir<-("C:/Users/YoonJoung Choi/World Health Organization/BANICA, Sorin - HSA unit/2 Global goods & tools/2 HFAs/1 HFAs for COVID-19/4. Implementation support materials/4. Analysis and dashboards/")

limesurveydir<-("~/Dropbox/0 iSquared/iSquared_WHO/ACTA/3.AnalysisPlan/ExportedCSV_FromLimeSurvey/")
##limesurveydir<-("C:/Users/YoonJoung Choi/World Health Organization/BANICA, Sorin - HSA unit/2 Global goods & tools/2 HFAs/1 HFAs for COVID-19/4. Implementation support materials/4. Analysis and dashboards/DownloadedCSV/")

# Define local macro for the survey 
country<-"COUNTRYNAME" #country name 
round<-1                #round      
year<-2020              #year of the mid point in data collection   
month<-12               #month of the mid point in data collection              
surveyid<-777777        #LimeSurvey survey ID   
# local macro for analysis: no change needed  

B. Import and drop duplicate cases

B.1. Import raw data from LimeSurvey
#Import directly from LimeSurvey - using the correct link, according to the country overview page (edit 6/23/2021)
dtaraw<-read.csv(paste0("https://who.my-survey.host/index.php/plugins/direct?plugin=CountryOverview&docType=1&sid=",surveyid,"&language=en&function=createExport"))

obsraw<-nrow(dtaraw)
cols<-ncol(dtaraw)

Note: For the URL, we need to use part of the country overview page for the data server. For example, suppose the overview page link looks like this for a country named YYY:
https://extranet.who.int/dataformv3/index.php/plugins/direct?plugin=CountryOverview&country=YYY&password=XXXXXXXXX.

Replace part of the link before plugins with the part in the country-specific link. So, in this example, the code should be:

dtaraw<-read.csv(paste0("https://extranet.who.int/dataformv3/index.php/plugins/direct?plugin=CountryOverview&docType=1&sid=",surveyid,"&language=en&function=createExport"))

As opposed to the above:

dtaraw<-read.csv(paste0("https://who.my-survey.host/index.php/plugins/direct?plugin=CountryOverview&docType=1&sid=",surveyid,"&language=en&function=createExport"))
#Use mock data for practice 
dtaraw<-read.csv(paste0(limesurveydir,"LimeSurvey_CEHS_Example_R1.csv"))

obsraw<-nrow(dtaraw)
cols<-ncol(dtaraw)

As of 2022-05-31 15:31:32, the downloaded raw data has 81 observations and 647 variables.


B.2. Export/save the data daily in CSV form with date
write.csv(dtaraw, paste0(limesurveydir,"LimeSurvey_CEHS_", country, "_R", round, "_", date, ".csv"))

B.3. Export the data for chartbook

Then, export the raw data into the chartbook (green tab: “Facility-level raw data”) - as is.

dtaraw<-dtaraw%>%mutate(test=Sys.time())

wb <- loadWorkbook(paste0(chartbookdir, "WHO_CEHS_Chartbook.xlsx"))
writeData(wb, sheet = "Facility-level raw data", dtaraw)
saveWorkbook(wb, paste0(chartbookdir, "WHO_CEHS_Chartbook.xlsx"), overwrite = TRUE)

B.4. Drop duplicate cases

Assess duplicate rows based on facility code: Q101. Among the duplicate rows, keep the latest row, based on the submission date/time.

Data managers must check the submitdate string values. Then modify format when you convert the string value to numeric, using as.POSIXct. See an example below. (Revision on 5/19/2022)

# CHECK submitdate 
str(dtaraw$submitdate)

# Identify duplicates 
dta<-dtaraw %>% 
    mutate(
        submitdate_string=as.character(submitdate),
        
        #use this line if submitedate does not have seconds 
        submitdate=as.POSIXct(as.character(submitdate),format= "%m/%d/%Y %H:%M")
        #use this line if submitedate dos seconds 
        #submitdate=as.POSIXct(as.character(submitdate),format= "%m/%d/%Y %H:%M:%S")
        )%>%
    group_by(Q101)%>% 
    mutate(
        dupe = n()>1, 
        latest = submitdate ==max(submitdate, na.rm = TRUE))%>%
    ungroup()

# Assess duplicates 
    table(dta$dupe, dta$latest) #cross tab
    dtaduplicate<-dta%>%filter(dupe==1)
    obsduplicate<-dtaduplicate%>%nrow() #number of duplicate rows
    obsduplicateunique<-length(unique(dtaduplicate$Q101)) #number of unique facilities among duplicate rows

#keep only unique/latest rows per facility     
dta<-dta%>%
    filter(latest==1)%>%
    select(-dupe, -latest)

#keep only if ID is NOT missing 
dta<-dta%>%
    filter(is.na(Q101)==FALSE)

obs<-nrow(dta)
obsunique<-length(unique(dta$Q101))

A total of 2 duplicate rows from 1 unique facilities were identified based on Q101. Now keeping only one unique row per facility, the dataset has 80 observations.


C. Destring and recoding

C.1. Change variable names to lowercase
names(dta)<-tolower(names(dta))

C.1.a Assess time stamp data
#REVISED 4/20/2021
#   interviewtime is availabl in dataset only when directly downloaded from the server, 
#   not via export plug-in used in this code
#   thus C.1.a is revised - see markdownfile for OLD and NEW code differences  
dta<-dta%>%
    select(-ends_with("time"))%>%
    select(-starts_with("grouptime"))

C.2. Change variable names to drop odd elements “y” “sq” - because of Lime survey’s naming convention
#Assess odd names to change
colnames(dta)

dtanew<-dta%>%
    rename_all(.funs = funs(sub("ysq", "_", .)))%>%
    rename_all(.funs = funs(sub("sqsq", "_", .)))%>%
    rename_at(.vars = vars(ends_with("sq")), .funs = funs(sub("sq", "", .)))%>%
    rename_all(.funs = funs(sub("_sq", "_", .)))%>%
    rename_all(.funs = funs(sub("sq", "_", .)))%>%
    rename_all(.funs = funs(sub("q409b_", "q409_", .)))%>%
    rename_at(.vars = vars(ends_with("_a1")), .funs = funs(sub("_a1", "_001", .)))%>%
    rename_at(.vars = vars(ends_with("_a2")), .funs = funs(sub("_a2", "_002", .)))

colnames(dtanew)    

dta<-dtanew

C.3. Find non-numeric variables and desting

Check variables section by scion, drop prefix “A” in categorical/character variables, and convert to numeric.

#change all factor to numeric
dta<-dta%>%
    mutate_if(is.factor, as.character) 

#####* Section 1
    varlist<-dta%>%select(q104, q105, q106, starts_with("q114"))%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), list(~ sub("-oth-", "88", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) #check variables ready for further processing

#####* Section 2
    varlist<-dta%>%select(starts_with("q2"))%>%colnames()
    str(dta[varlist]) 
    
    varlist<-dta%>%select(q204, starts_with("q205"), starts_with("q207"))%>%colnames()    
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))    

#####* Section 3
    varlist<-dta%>%select(starts_with("q3"))%>%colnames()
    str(dta[varlist])

    varlist<-dta%>%select(q302, q309)%>%colnames() 
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric)) 
    
    str(dta[varlist])     

#####* Section 4
    varlist<-dta%>%select(starts_with("q406"), starts_with("q409"),
                          starts_with("q412"), q414, q415,
                          starts_with("q417"), starts_with("q420"),
                          starts_with("q421") 
                          )%>%colnames()
    str(dta[varlist])

    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) 
    
    dta<-dta%>%
        rename(q410other = q410_other)%>%
        mutate(
            q410_007=NA, 
            q410_007 = ifelse(q410other!="", 1, q410_007),
            q410_007 = ifelse(q410other=="", 0, q410_007),
            q410_007 = ifelse(is.na(q410_001)==TRUE , NA, q410_007)
        )%>%
        rename(q411other = q411_other)%>%
        mutate(
            q411_012=NA, 
            q411_012 = ifelse(q411other!="", 1, q411_012),
            q411_012 = ifelse(q411other=="", 0, q411_012),
            q411_012 = ifelse(is.na(q411_001)==TRUE , NA, q411_012)
        )
    
#####* Section 5
    varlist<-dta%>%select(starts_with("q503"), starts_with("q505"), starts_with("q507"))%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) 

#####* Section 6
    varlist<-dta%>%select(q604, starts_with("q607"), starts_with("q608"))%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) 

#####* Section 7
    varlist<-dta%>%select(starts_with("q701"), starts_with("q702"), starts_with("q703"))%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) 

#####* Section 8
    varlist<-dta%>%select(starts_with("q8"))%>%colnames()
    str(dta[varlist])
    
    varlist<-dta%>%select(starts_with("q802"), starts_with("q803"), starts_with("q805"))%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) 
    
#####* Section 9: General vaccine readiness 
    varlist<-dta%>%select(starts_with("q9"))%>%colnames()
    str(dta[varlist])
    
    varlist<-dta%>%select(q907, q910, q911)%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))
    
    str(dta[varlist]) 

#####* Section 10: COVID vaccine readiness 
    
#####* Section 11: INFRASTRUCTURE
    
#####* Section 12
    varlist<-dta%>%select(starts_with("q120"))%>%colnames()
    str(dta[varlist])

    varlist<-dta%>%select(q1201, q1204)%>%colnames()
    str(dta[varlist])
    
    dta<-dta%>%
        mutate_at(vars(varlist), list(~ sub("^A", "\\1", .)))%>%
        mutate_at(vars(varlist), funs(as.numeric))        

C.4. Recode yes/no & yes/no/NA
# Recode 2=0 for no
varlist<-dta%>%
    select(
        q202, starts_with("q205_"), q206, starts_with("q207"), 
        q301, q304, starts_with("q306_"), q307, q308, q310,
        q402, q403, q404, q405, starts_with("q406"), q407, q408,
            starts_with("q410"), starts_with("q411"), q416, q418, q419, 
        q501, q502, starts_with("q503"), q504, starts_with("q505"), q506,
        q601, q602, q603, q605, q606, q609, q610, q614, q615, 
        starts_with("q701"), starts_with("q702"), starts_with("q703"), q704, 
        q801, q804, 
        q901, q902, q905, q908, q912, q913, q914 , 
        q1003, starts_with("q1005_"), q1006, q1007, q1008, q1009, 
            q1010, q1011, q1012, q1013, q1014
    )%>%
    colnames()

dta[varlist][dta[varlist] == 2 ] <- 0 #no

C.5. Label values

The following is value labels.

    #delimit;   
    
    lab define q104 
        1"1.Urban" 
        2"2.Rural";  
    lab values q104 q104; 
    
    lab define q105 
        1"1.Level2: Dispensary"
        2"2.Level3: Health Centre"
        3"3.Level4: Primary Hospital"
        4"4.Level5: Secondary level hospital and above"
        5"5.Level6: Tertiary level" ; 
    lab values q105 q105; /*corresponding to the model example context*/
    
    lab define q106 
        1"1.Government"
        2"2.Private"
        3"3.NGO"
        4"4.FBO"
        5"5.Other"; 
    lab values q106 q106; /*corresponding to the model example context*/
    
    lab define q302
        1"1. Yes - user fees exempted only for COVID-19 services"
        2"2. Yes - user fees exempted only for other health services"
        3"3. Yes - user fees exempted for both COVID-19 and other health services"
        4"4. No"; 
    lab values q302 q302; 
        
    lab define q305
        1"1. Yes - for COVID-19 case management services"
        2"2. Yes - for other essential health services"
        3"3. Yes - for both COVID-19 and other services"
        4"4. No"
        5"5. Do not know"; 
    lab values q305 q305; 
        
    lab define change   /*KECT - q409a is weird (asks if all increased, all decreased, no changes, mixed), so it cannot be included here, changed from q409* to q409_* to address this issue*/
        1"1.Yes, increased"
        2"2.Yes, decreased"
        3"3.No" 
        4"4.N/A";  
    foreach var of varlist q409_* q412* q414 q415 q417* q420*{;
    lab values `var' change;    
    };
    
    lab define q409a   
        1"1.Yes, increased in all areas"
        2"2.Yes, decreased in all areas"
        3"3.Yes, increased in some but decreased in some" 
        4"4.No changes in any";  
    foreach var of varlist q409a {;
    lab values `var' q409a; 
    };
    
    lab define q417
        1"1.Yes, less frequent"
        2"2.Yes, suspended"
        3"3.No change" 
        4"4.yes, increased" 
        5"5.N/A";  
    foreach var of varlist q417* {;
    lab values `var' q417;  
    };  
    
    lab define q420
        1"1.Yes, planned & implemented"
        2"2.Yes, planned but not yet implemented"
        3"3.No" 
        4"4.N/A";  
    foreach var of varlist q420* {;
    lab values `var' q420;  
    };  
    
    lab define q421
        1"1.Not at all" 
        2"2.Slightly"   
        3"3.Moderately"     
        4"4.Quite a lot"    
        5"5.A great deal";  
    foreach var of varlist q421* {;
    lab values `var' q421;  
    };      

    lab define ppe
        1"1.Currently available for all health workers"
        2"2.Currently available only for some health workers"
        3"3.Currently unavailable for any health workers"
        4"4.Not applicable - never procured or provided" ;
    foreach var of varlist q507* {;
    lab values `var' ppe;   
    };      

    lab define q604
        1"1.Yes, PCR"
        2"2.Yes, RDT"
        3"2.Yes, PCR & RDT"
        4"5.No";
    lab values q604 q604;       

    lab define availfunc 
        1"1.Available, functional"
        2"2.Available, but not functional"
        3"3.Not available";
    foreach var of varlist q802_* q803_* q805_* q903 q904  {;
    lab values `var' availfunc ;    
    };          

    lab define icepack 
        1"1.Yes, a set of ice packs for all cold boxes"
        2"1.Yes, a set of ice packs only for some cold boxes"
        3"3.No";
    foreach var of varlist q907 q910  {;
    lab values `var' icepack ;  
    };      
    
    lab define icepackfreeze 
        1"1.All"
        2"2.Only some"
        3"3.None-no functional freezer" ;
    lab values q911 icepackfreeze ; 
    
    lab define yesno 1"1.Yes" 0"0.No";  
    foreach var of varlist 
        q202  q205_* q206   q207* 
        q301  q304   q306_* q307  q308  q310
        q402 - q405  q406_* q407  q408  q410_* q411_* q416 q418 q419 
        q501  q502  q503_* q504  q505_* q506
        q601  q602  q603  q605  q606 q609  q610 q614-q615
        q701_* q702_* q703_* q704
        q801 q804
        q901 q902  q905 q908 q912 
        {;      
    labe values `var' yesno; 
    };  

    lab define yesnona 
        1"1.Yes" 
        2"2.No" 
        3"3.N/A";   
    foreach var of varlist q204 q309 {;     
    labe values `var' yesnona; 
    };  
    
    lab define alwayssometimesnever 
        1"1.Always" 
        2"2.Sometimes" 
        3"3.Never"; 
    foreach var of varlist q607_* q608_* {;     
    labe values `var' alwayssometimesnever; 
    };  
        
    #delimit cr

D. Create field check tables

As of 2022-05-31, the following are “field check tables.” In Stata program, xls file is created. In R, the results are directly presented in this markdown file.

#REVISED 4/20/2021
#   since interviewlength is not created (see C.1.a) 
#   do not calculate "time" variables 
dtacheck<-dta%>%
    mutate(
        updatedate= as.character(as.Date(Sys.time(  ), format='%d%b%Y')), 
        date=as.POSIXct(as.character(submitdate_string),format= "%m/%d/%Y"), 
        xresult=q1204==1, 
        responserate= xresult==1 #label define responselist 0 "Not complete" 1 "Complete"    
        )

Assess interview characteristics among all interviews.

# Date of field check table update and the total number of interviews
print(date)

[1] “2022-05-31”

# Date of interviews (submission date, final)
table(dtacheck$date)

2020-12-14 2020-12-15 2020-12-16 2020-12-17 2020-12-18 2020-12-19 20 16 12 12 16 4

# Interview response rate (%)
print(paste(as.character(round(mean(dtacheck$responserate, na.rm=TRUE)*100, 1)), "%"))

[1] “94.7 %”

# Number of interviews by facility type
table(dtacheck$q104)

1 2 36 44

#Average interview length (minutes), among completed interviews
#print(paste(as.character(round(mean(dtacheck$time_complete, na.rm=TRUE))), "minutes"))

#Average interview length (minutes), among partly completed interviews
#print(paste(as.character(round(mean(dtacheck$time_incomplete, na.rm=TRUE))), "minutes"))
#####/*the following calcualtes % missing in select questions among completed interviews*/  

#0. Missing survery results (among all interviews)
dtacheck<-dtacheck%>%
    mutate(
        missing0 = is.na(q1204)==1
        )

#1. Missing number of beds when facility provides inpatient services  
dtacheck<-dtacheck%>%
    mutate(
        missing1 = 0,
        missing1 = ifelse(is.na(q111)==1, 1,
                          ifelse(q111==1 & is.na(q112)==1, 1,
                                 ifelse(q111==1 & is.na(q113)==1, 1,
                                        ifelse(q111==0, NA, missing1)))) ) 

#2. Missing number of nurses (either the total number or the number who have been infected) (among completed interviews)
dtacheck<-dtacheck%>%
    mutate(
        missing2 = 0,
        missing2 = ifelse(is.na(q201_002_001)==1, 1, 
                         ifelse(is.na(q201_002_002)==1, 1, missing2)) 
        )

#3. Missing salary
dtacheck<-dtacheck%>%
    mutate(
        missing3 = 0,
        missing3 = ifelse(is.na(q307)==1, 1, 
                          ifelse(is.na(q308)==1, 1, 
                                 ifelse(q308==1 & is.na(q309)==1, 1, missing3))) 
    )

#4. Missing response to strategy change 
varlist<-dtacheck%>%select(starts_with("q406_"))%>%colnames()
str(dtacheck[varlist])
dtacheck[varlist] <- lapply(dtacheck[varlist],
                             function(x){ifelse(is.na(x)==0,0,
                                                ifelse(is.na(x)==1,1,x))})

dtacheck<-dtacheck%>%
    mutate(
        missing_num = rowSums(dtacheck[varlist] ), 
        missing4 = missing_num, 
        missing4 = ifelse(missing_num>=1, 1, missing4) )

#5. Missing response to OPT volume change
varlist<-dtacheck%>%select(starts_with("q409_"))%>%colnames()
str(dtacheck[varlist])
dtacheck[varlist] <- lapply(dtacheck[varlist],
                             function(x){ifelse(is.na(x)==0,0,
                                                ifelse(is.na(x)==1,1,x))})

dtacheck<-dtacheck%>%
    mutate(
        missing_num = rowSums(dtacheck[varlist] ), 
        missing5 = missing_num, 
        missing5 = ifelse(missing_num>=1, 1, missing5) )

#6. Missing response to restoration plan 
varlist<-dtacheck%>%select(starts_with("q420_"))%>%colnames()
str(dtacheck[varlist])
dtacheck[varlist] <- lapply(dtacheck[varlist],
                             function(x){ifelse(is.na(x)==0,0,
                                                ifelse(is.na(x)==1,1,x))})

dtacheck<-dtacheck%>%
    mutate(
        missing_num = rowSums(dtacheck[varlist] ), 
        missing6 = missing_num, 
        missing6 = ifelse(missing_num>=1, 1, missing6) )

#7. Missing response to PPE
varlist<-dtacheck%>%select(starts_with("q507_"))%>%colnames()
str(dtacheck[varlist])
dtacheck[varlist] <- lapply(dtacheck[varlist],
                             function(x){ifelse(is.na(x)==0,0,
                                                ifelse(is.na(x)==1,1,x))})

dtacheck<-dtacheck%>%
    mutate(
        missingppe_num = rowSums(dtacheck[varlist] ), 
        missing7 = missingppe_num, 
        missing7 = ifelse(missingppe_num>=1, 1, missing7) )

#8. Missing lab transportation 
varlist<-dtacheck%>%select(q604, q605)%>%colnames()

dtacheck<-dtacheck%>%
    mutate(
        missing8 = is.na(q604)==1, 
        missing8 = ifelse(q604==4 & is.na(q605)==1, 1, missing8) )

#9. Missing response to medicines questions  
varlist<-dtacheck%>%select(starts_with("q701_"))%>%colnames()
str(dtacheck[varlist])
dtacheck[varlist] <- lapply(dtacheck[varlist],
                             function(x){ifelse(is.na(x)==0,0,
                                                ifelse(is.na(x)==1,1,x))})

dtacheck<-dtacheck%>%
    mutate(
        missingmed_num = rowSums(dtacheck[varlist] ), 
        missing9 = missingmed_num, 
        missing9 = ifelse(missingmed_num>=1, 1, missing9) )

Assess level of missing responses among COMPLETED interviews.

dtacheck<-dtacheck%>%
    filter(xresult==1)

#1. Missing number of beds when facility provides inpatient services  
#print("Percent of interviews with missing responses for inpatient or ICU beds")   
print(paste(as.character(round(mean(dtacheck$missing1, na.rm=TRUE)*100, 1)), "%"))

[1] “0 %”

#3. Missing overtime and overtime payment responses 
#print("Percent of interviews with missing number of nurses")   
print(paste(as.character(round(mean(dtacheck$missing2, na.rm=TRUE)*100, 1)), "%"))

[1] “0 %”

#3. Missing overtime and overtime payment responses 
#print("Percent of interviews with missing responses for overtime and overtime payment")
print(paste(as.character(round(mean(dtacheck$missing3, na.rm=TRUE)*100, 1)), "%"))

[1] “0 %”

#4. Missing response to strategy changes
#print("Percent of interviews with one or missing responses in strategy changes")
print(paste(as.character(round(mean(dtacheck$missing4, na.rm=TRUE)*100, 1)), "%")) 

[1] “0 %”

#5. Missing response to OPT volume changes
#print("Percent of interviews with one or missing responses in OPT volume changes")
print(paste(as.character(round(mean(dtacheck$missing5, na.rm=TRUE)*100, 1)), "%")) 

[1] “22.2 %”

#6. Missing response to restoration plan questions 
#print("Percent of interviews with one or missing responses in restoration plan questions ")
print(paste(as.character(round(mean(dtacheck$missing6, na.rm=TRUE)*100, 1)), "%")) 

[1] “38.9 %”

#7. Missing response to PPE (%)
#print("Percent of interviews with one or missing responses in PPE questions")
print(paste(as.character(round(mean(dtacheck$missing7, na.rm=TRUE)*100, 1)), "%")) 

[1] “0 %”

#print("Of those, average number of missing responses")
temp<-dtacheck%>%filter(missing7==1)
round(mean(temp$missingppe_num, na.rm=TRUE), 1)

[1] NaN

#8. Missing response to lab referral transportation  
#print("Percent of interviews with missing responses in lab referral transportatin questions")
print(paste(as.character(round(mean(dtacheck$missing8, na.rm=TRUE)*100, 1)), "%")) 

[1] “0 %”

#9. Missing response to medicines questins 
#print("9. Percent of interviews with one or missing responses in medicines questions")   
print(paste(as.character(round(mean(dtacheck$missing9, na.rm=TRUE)*100, 1)), "%")) 

[1] “0 %”

#print("9.1. Of those, average number of missing responses")
temp<-dtacheck%>%filter(missing9==1)
round(mean(temp$missingmed_num, na.rm=TRUE), 1)

[1] NaN


E. Create analytical variables

E.1. Country speciic code local

Update code here based on the questionnaire in the country

urbanmin<-1
urbanmax<-1

minlow  <-1 #/*lowest code for lower-level facilities in Q105*/
maxlow  <-3 #/*highest code for lower-level facilities in Q105*/
minhigh <-4 #/*lowest code for hospital/high-level facilities in Q105*/
maxhigh <-6 #/*highest code for hospital/high-level facilities in Q105*/
primaryhospital <-4 #/*district hospital or equivalent */   

pubmin<-1
pubmax<-1

maxtraining <-7 
maxtrainingsupport <-9 

maxdrug <-17 #/*total medicines asked in q401*/

E.2. Construct analysis variables

Created analysis variables section by section.
* give prefix z for background characteristics, which can be used as analysis strata
* give prefix x for binary variables, which will be used to calculate percentage
* give prefix y for integer/continuous variables, which will be used to calculate total number

#####* Section 1 
dta<-dta%>%
    mutate(
        country = country, 
        round = round,      
        facilitycode = q101, 
        month = month, 
        year = year, 
        
        zurban  = q104>=urbanmin & q104<=urbanmax,  
        zlevel  =q105,
        zlevel  ="", 
            zlevel  =ifelse(q105==1, "Level2", zlevel),
            zlevel  =ifelse(q105==2, "Level3", zlevel),
            zlevel  =ifelse(q105==3, "Level4", zlevel),
            zlevel  =ifelse(q105==4, "Level5", zlevel),
            zlevel  =ifelse(q105==5, "Level6", zlevel),
        
        zlevel_hospital     =q105>=minhigh & q105<=maxhigh,
        zlevel_primhospital =q105==primaryhospital, 
        zlevel_low          =q105>=minlow  & q105<=maxlow, 
        
        zpub    =q106>=pubmin & q106<=pubmax
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("z")), 
              ~replace(., is.na(.)==TRUE, 0))
    

    #lab define zurban 0"Rural" 1"Urban"
    #lab define zlevel_hospital 0"Non-hospital" 1"Hospital"
    #lab define zpub 0"Non-public" 1"Public"

    #lab values  zurban zurban
    #lab values  zlevel_hospital zlevel_hospital
    #lab values  zpub zpub
    
    #lab var id "ID generated from Lime Survey"
    #lab var facilitycode "facility ID from sample list" /*this will be used to merge with sampling weight, if relevant*/
#####* Section 2: staff
dta<-dta%>%
    mutate(
        staff_num_total_md=q201_001_001 , 
        staff_num_covid_md=q201_001_002 , 
        
        staff_num_total_nr=q201_002_001 , 
        staff_num_covid_nr=q201_002_002 , 
        
        staff_num_total_othclinical=rowSums(dta[ c("q201_003_001", "q201_004_001", "q201_005_001", "q201_006_001", "q201_007_001")], na.rm = TRUE), 
        staff_num_covid_othclinical=rowSums(dta[ c("q201_003_002", "q201_004_002", "q201_005_002", "q201_006_002", "q201_007_002")], na.rm = TRUE))
        
dta<-dta%>%
    mutate(        
        staff_num_total_clinical=rowSums(dta[ c("staff_num_total_md", "staff_num_total_nr", "staff_num_total_othclinical")], na.rm = TRUE) , 
        staff_num_covid_clinical=rowSums(dta[ c("staff_num_covid_md", "staff_num_covid_nr", "staff_num_covid_othclinical")], na.rm = TRUE) , 
        
        staff_num_total_nonclinical=rowSums(dta[ c("q201_008_001", "q201_009_001", "q201_010_001")], na.rm = TRUE) , 
        staff_num_covid_nonclinical=rowSums(dta[ c("q201_008_002", "q201_009_002", "q201_010_002")], na.rm = TRUE) )

dta<-dta%>%
    mutate(
        staff_num_total_all=rowSums(dta[ c("staff_num_total_clinical", "staff_num_total_nonclinical")], na.rm = TRUE),  
        staff_num_covid_all=rowSums(dta[ c("staff_num_covid_clinical", "staff_num_covid_nonclinical")], na.rm = TRUE), 
        
        xstaff_covax = q201a==1, 
        staff_num_covaxany  = q201b, 
        staff_num_covaxfull = q201c, 
        staff_num_covaxany  =replace(staff_num_covaxany, xstaff_covax!=1, NA), 
        staff_num_covaxfull  =replace(staff_num_covaxfull, xstaff_covax!=1, NA),    
        
        xabsence=q202==1,
        xabsence_medical=   q203_003==1 | q203_004==1 ,
        xabsence_structure= q203_005==1 | q203_006==1 | q203_007==1 ,
        xabsence_social=    q203_008==1 | q203_009==1 | q203_010==1 ,
                
        xhr=q204==1,
        xhr_shift=      q205_001,
        xhr_increase=   q205_002==1 | q205_003==1 | q205_004==1 | q205_005==1 | q205_006==1 ,
        xhr_increase_exp=   q205_002==1 | q205_003==1 ,
        xhr_increase_new=   q205_004==1 | q205_005==1 | q205_006==1 ,
        xhr_secondment= q205_007==1,
        xhr_decrease=   q205_008==1,
    
        xtraining=q206==1
        )%>%
    mutate_at(vars(xabsence, starts_with("xabsence_"), xhr, starts_with("xhr_"), starts_with("xtraining")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    rename_at(.vars = vars(starts_with("q207_")), .funs = funs(sub("q207_", "xtraining__", .)))%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xtraining__")), 
              ~replace(., is.na(.)==TRUE, 0))

dta<-dta%>% 
    mutate(
            max=maxtraining,
            temp=rowSums(dta[ c("xtraining__001", "xtraining__002", "xtraining__003", "xtraining__004", "xtraining__005")], na.rm=TRUE) , 
        xtraining_score =100*(temp/max) , 
        xtraining_100   =xtraining_score>=100,
        xtraining_50    =xtraining_score>=50,
        
            max=maxtrainingsupport,
            temp=rowSums(dta[colnames(select(dta, starts_with("xtraining__")))], na.rm=TRUE ) , 
        xtrainingsupport_score  =100*(temp/max) , 
        xtrainingsupport_100    =xtrainingsupport_score>=100,
        xtrainingsupport_50     =xtrainingsupport_score>=50
        
    )
#####* Section 3: finance
dta<-dta%>%
    mutate(
        xuserfee=       q301==1 ,
        xexempt=        q302==1 | q302==2 | q302==3 | q303==1 , 
        xexempt_covid=  q302==1 | q302==3 ,
        xexempt_other=  q302==2 | q302==3 , 
        xexempt_vulnp=  q303==1  ,
        xfeeincrease=   q304==1 )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xuserfee"), starts_with("xexempt"), starts_with("xfeeincrease")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate_at(vars(starts_with("xexempt"), xfeeincrease), 
              funs(ifelse(xuserfee!=1, NA, .)) )%>%
    mutate(
        xaddfund = q305==1 | q305==2 | q305==3,
        xaddfund_covid = q305==1 | q305==3,
        xaddfund_ehs = q305==2 | q305==3,
        xaddfund_gov    = q306_001==1 , 
        xaddfund_other  = q306_002==1 | q306_003==1 | q306_004==1 | q306_005==1 , 
            
        xfinance_salaryontime   = q307==1 ,
        xfinance_ot             = q308==1  )%>%
    mutate_at(vars(starts_with("xaddfund"), starts_with("xfinance")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate(
        xfinance_otontime   = q309==1 | q309==3 ,
            xfinance_otontime = ifelse(xfinance_ot==0, NA, xfinance_otontime) ,
            
        xfinance_ontime = xfinance_salaryontime ==1 & (xfinance_otontime==1 | is.na(xfinance_otontime)==TRUE), 
        
        xfinance_PBF = q310==1
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(xfinance_ontime), 
              ~replace(., is.na(.)==TRUE, 0))
#####* Section 4: service delivery & utlization  
    
        varlist=dtacheck%>%select(q402, q403, starts_with("q406"))%>%colnames()

dta<-dta%>%
    mutate(
        xguideline_EHS      =q404==1 , 
        xguideline_EHScovid =q405==1 , 
        
        ###*** Service delivery strategy 
        xstrategy= rowSums(dta[varlist], na.rm = TRUE)>=1, 
        
        xstrategy_reduce=   q402==1  | q403==1  | q406_001==1 | q406_002==1 | q406_003==1 | q406_004==1 | q406_005==1,
        xstrategy_reduce_closure=   q402==1  ,
        #REVISION on 10/29/2021. Correct xstrategy_reduce_hrchange
        #xstrategy_reduce_hrchange=     q403==1  | q406_001==1 | q406_002==1 | q406_003==1 | q406_004==1 | q406_005==1,
        xstrategy_reduce_hrchange=  q403==1, 
        #END OF REVISION
        xstrategy_reduce_reduce=    q406_001==1 | q406_002==1 | q406_003==1,
        xstrategy_reduce_redirect=  q406_004==1,
        xstrategy_reduce_priority=  q406_005==1,
        xstrategy_reduce_combine =  q406_006==1,
        xstrategy_self=     q406_007==1,
        xstrategy_home=     q406_008==1,
        xstrategy_remote=   q406_009==1 ,
        xstrategy_prescription=     q406_010==1 | q406_011==1 | q406_012==1 ,
        
        ###*** Referral/transportation for COVID patients 
        xcvd_ref        = q407==1   , 
        xcvd_reftrans   = q407==1 & q408==1      
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xstrategy"), starts_with("xcvd_")), 
              ~replace(., is.na(.)==TRUE, 0))

        ###***** OPT

        temp1<-dta%>%select(starts_with("q409_"))%>%
            rename_all(.funs = funs(sub("q409_", "xopt_increase__", .)))%>%
            mutate_all( funs(ifelse(.==4, NA,
                                ifelse(. ==1 , 1, 
                                  ifelse(. !=1, 0, .)))) )
        
        
        temp2<-dta%>%select(starts_with("q409_"))%>%
            rename_all(.funs = funs(sub("q409_", "xopt_decrease__", .)))%>%
            mutate_all( funs(ifelse(.==4, NA,
                                ifelse(. ==2 , 1, 
                                  ifelse(. !=2, 0, .)))) )
        
        temp3<-dta%>%select(starts_with("q409_"))%>%
            rename_all(.funs = funs(sub("q409_", "xopt_nochange__", .)))%>%
            mutate_all( funs(ifelse(.==4, NA,
                                ifelse(. ==3 , 1, 
                                  ifelse(. !=3, 0, .)))) )
        
        temp4<-dta%>%select(starts_with("q409_"))%>%
            rename_all(.funs = funs(sub("q409_", "xopt_changeNA__", .)))%>%
            mutate_all( funs(ifelse(. ==4 , 1, 
                                  ifelse(. !=4, 0, .))) )%>%
            mutate_all( funs(ifelse(is.na(.)==TRUE , 0, .))) 
        
        temp5<-dta%>%select(starts_with("q409_"))%>%
            rename_all(.funs = funs(sub("q409_", "xopt_changeMissing__", .)))%>%
            mutate_all( funs(ifelse(is.na(.)==TRUE , 1, 
                                      ifelse(is.na(.)==FALSE, 0, .))) ) 

dta<-cbind(dta, temp1, temp2, temp3, temp4, temp5)

dta<-dta%>%
    mutate(
        xopt_increase_num=rowSums(dta[colnames(select(dta, starts_with("xopt_increase__")))], na.rm=TRUE) , 
        xopt_decrease_num=rowSums(dta[colnames(select(dta, starts_with("xopt_decrease__")))], na.rm=TRUE) , 
        xopt_nochange_num=rowSums(dta[colnames(select(dta, starts_with("xopt_nochange__")))], na.rm=TRUE) , 
        xopt_increase = xopt_increase_num>=1, 
        xopt_decrease = xopt_decrease_num>=1    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(xopt_increase, xopt_decrease, 
                   starts_with("xopt_increase_num"), starts_with("xopt_decrease_num")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate(
        xopt_change=NA, 
            xopt_change=ifelse(xopt_decrease==0 & xopt_increase==0, 1, xopt_change), #/*no change in any*/
            xopt_change=ifelse(xopt_decrease==1 & xopt_increase==1, 2, xopt_change), #/*mixed*/
            xopt_change=ifelse(xopt_decrease==1 & xopt_increase==0, 3, xopt_change), #/*no service volume increased + at least one decreased*/
            xopt_change=ifelse(xopt_decrease==0 & xopt_increase==1, 4, xopt_change) #/*no service volume decreased + at least one increased*/
    )%>%
    
    #REVISION 11/1/2021: generate overall service volumne change var 
    mutate(
        xoptoverall_increaseall = q409a ==1,
        xoptoverall_decreaseall = q409a ==2,
        xoptoverall_mixed = q409a ==3,
        xoptoverall_nochange = q409a ==4,        
        xoptoverall_increaseall = ifelse(is.na(q409a)==TRUE, NA, xoptoverall_increaseall), 
        xoptoverall_decreaseall = ifelse(is.na(q409a)==TRUE, NA, xoptoverall_decreaseall), 
        xoptoverall_mixed = ifelse(is.na(q409a)==TRUE, NA, xoptoverall_mixed), 
        xoptoverall_nochange = ifelse(is.na(q409a)==TRUE, NA, xoptoverall_nochange) 
    )
    #ENF OF REVISION

        ###***** REASONS for OPT volume changes

        temp1<-dta%>%select(xopt_increase, starts_with("q410_"))%>%
            rename_all(.funs = funs(sub("q410_", "xopt_increase_reason__", .)))%>%
            mutate_at(vars(starts_with("xopt_increase_reason__")), 
                      funs(ifelse(is.na(.)==TRUE, 0, .)) )%>%
            select( - xopt_increase)
        
        temp2<-dta%>%select(xopt_decrease, starts_with("q411_"))%>%
            rename_all(.funs = funs(sub("q411_", "xopt_decrease_reason__", .)))%>%
            mutate_at(vars(starts_with("xopt_decrease_reason__")), 
                      funs(ifelse(is.na(.)==TRUE, 0, .)) )%>%
            select( - xopt_decrease)
        
dta<-cbind(dta, temp1, temp2)

dta<-dta%>%
    mutate(
        xopt_increase_reason_covidnow   = xopt_increase==1 & (xopt_increase_reason__001==1 | xopt_increase_reason__002==1 ) , 
        xopt_increase_reason_covidafter = xopt_increase==1 & (xopt_increase_reason__003==1 | xopt_increase_reason__004==1 | xopt_increase_reason__005==1), 
        xopt_increase_reason_gbv        = xopt_increase==1 & (xopt_increase_reason__006==1 ), 
    
        xopt_decrease_reason_comdemand  = xopt_decrease==1 & (xopt_decrease_reason__001==1 | xopt_decrease_reason__002==1 | xopt_decrease_reason__005==1  ),
        xopt_decrease_reason_enviro     = xopt_decrease==1 & (xopt_decrease_reason__003==1 | xopt_decrease_reason__004==1 ),
        xopt_decrease_reason_intention  = xopt_decrease==1 & (xopt_decrease_reason__007==1 | xopt_decrease_reason__008==1 | xopt_decrease_reason__009==1 | xopt_decrease_reason__010==1 ),
        xopt_decrease_reason_disruption = xopt_decrease==1 & (xopt_decrease_reason__011==1 | xopt_decrease_reason__012==1) 
        
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xopt_increase_reason"), starts_with("xopt_decrease_reason")), 
              ~replace(., is.na(.)==TRUE, 0))

        ###***** ER
dta<-dta%>%
    mutate(
        xer = q114_001==1,
        xer_increase= q412_001 ==1,
        xer_decrease= q412_001 ==2,
        xer_nochange= q412_001 ==3) %>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xer")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate_at(vars(starts_with("xer_")), 
              funs(ifelse(xer==0, NA, .)) )

        temp1<-dta%>%select(starts_with("q412_"))%>%
            select(-q412_001)%>%
            rename_all(.funs = funs(sub("q412_", "xer_increase__", .)))%>%
            mutate_at(vars(starts_with("xer_increase__")), 
                      funs(ifelse(. !=1, 0, 
                                  ifelse(. ==1, 1, .))))%>%
            mutate_all(~replace(., is.na(.), 0)) 
        
        temp2<-dta%>%select(starts_with("q412_"))%>%
            select(-q412_001)%>%
            rename_all(.funs = funs(sub("q412_", "xer_decrease__", .)))%>%
            mutate_at(vars(starts_with("xer_decrease__")), 
                      funs(ifelse(. !=2, 0, 
                                  ifelse(. ==2, 1, .))))%>%
            mutate_all(~replace(., is.na(.), 0)) 

dta<-cbind(dta, temp1, temp2)%>%
    mutate_at(vars(starts_with("xer_")), 
              funs(ifelse(xer==0, NA, .)) )

        ###***** IPT 
dta<-dta%>%
    mutate(
        xipt = q111==1, 
        xbed = q112,
        xbedrate = round(100*(q413/q112), 1),
        xipt_increase = q414==1,
        xipt_decrease = q414==2,
        xipt_nochange = q414==3
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xipt")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate_at(vars(starts_with("xipt_")), 
              funs(ifelse(xipt==0, NA, .)) )
            
        ###***** Pre hospital ER
dta<-dta%>%    
    mutate(
        xpreer = q415<=3,
        
        xpreer_increase = q415==1,
        xpreer_decrease = q415==2,
        xpreer_nochange = q415==3
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xpreer")), 
              ~replace(., is.na(.)==TRUE, 0))%>%    
    mutate_at(vars(starts_with("xpreer_")), 
              funs(ifelse(xpreer==0, NA, .)) )
            
        ###***** community outreach
dta<-dta%>%
    mutate(
        xout = q416==1)%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(xout), 
              ~replace(., is.na(.)==TRUE, 0))

        temp1<-dta%>%select(starts_with("q417_"))%>%
            rename_all(.funs = funs(sub("q417_", "xout_decrease__", .)))%>%
            mutate_all( funs(ifelse(. <=2 , 1, 
                                  ifelse(. >=3, 0, .))) )%>%
            mutate_all(~replace(., is.na(.), 0)) 

dta<-cbind(dta, temp1)

dta<-dta%>%
    mutate(
        xout_decrease_num=rowSums(dta[colnames(select(dta, starts_with("xout_decrease__")))], na.rm=TRUE ) , 
        xout_decrease= xout_decrease_num>=1         
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xout_decrease")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate_at(vars(starts_with("xout_decrease")), 
              funs(ifelse(xout==0, NA, .)) )

        ###***** missed appointment 
dta<-dta%>%
    mutate(     
        xresto = q418==1 & q419==1)%>%
    mutate_at(vars(xresto), 
              ~replace(., is.na(.)==TRUE, 0)) %>%
    mutate(
        xresto_imp_preg        = q420_001==1 ,
        xresto_imp_immunization= q420_002==1 ,
        xresto_imp_chronic     = q420_003==1 ,
        xresto_imp_tb          = q420_004==1 ,
        xresto_imp_hiv         = q420_005==1 ,
        
        xresto_imppln_preg          = q420_001==1 | q420_001==2,
        xresto_imppln_immunization  = q420_002==1 | q420_002==2,
        xresto_imppln_chronic       = q420_003==1 | q420_003==2,
        xresto_imppln_tb            = q420_004==1 | q420_004==2,
        xresto_imppln_hiv           = q420_005==1 | q420_005==2,
        
        xresto_imp_preg = ifelse(q420_001>=4 | is.na(q420_001)==TRUE, NA, xresto_imp_preg), 
        xresto_imppln_preg = ifelse(q420_001>=4 | is.na(q420_001)==TRUE, NA, xresto_imppln_preg), 

        xresto_imp_immunization = ifelse(q420_002>=4 | is.na(q420_002)==TRUE, NA, xresto_imp_immunization), 
        xresto_imppln_immunization = ifelse(q420_002>=4 | is.na(q420_002)==TRUE, NA, xresto_imppln_immunization),
        
        xresto_imp_chronic = ifelse(q420_003>=4 | is.na(q420_003)==TRUE, NA, xresto_imp_chronic),
        xresto_imppln_chronic = ifelse(q420_003>=4 | is.na(q420_003)==TRUE, NA, xresto_imppln_chronic),
        
        xresto_imp_tb = ifelse(q420_004>=4 | is.na(q420_004)==TRUE, NA, xresto_imp_tb),
        xresto_imppln_tb = ifelse(q420_004>=4 | is.na(q420_004)==TRUE, NA, xresto_imppln_tb),
        
        xresto_imp_hiv = ifelse(q420_005>=4 | is.na(q420_005)==TRUE, NA, xresto_imp_hiv),
        xresto_imppln_hiv = ifelse(q420_005>=4 | is.na(q420_005)==TRUE, NA, xresto_imppln_hiv)
        )

        ###***** disruption  
        temp1<-dta%>%select(starts_with("q421_"))%>%
            rename_all(.funs = funs(sub("q421_", "xdisrupt__", .)))%>%
            mutate_all(~replace(., is.na(.), 0))%>% 
            mutate_at(vars(starts_with("xdisrupt__")),
                      funs(ifelse( . <=3, 0,
                                  ifelse( . ==4 | . ==5, 1, .))))
            
        
dta<-cbind(dta, temp1)%>%
    rename(
        xdisrupt__hr = xdisrupt__001, 
        xdisrupt__finance = xdisrupt__002, 
        xdisrupt__ipc = xdisrupt__003, 
        xdisrupt__medsupp = xdisrupt__004
    )%>%
    mutate_if(is.logical, as.numeric)
#####* Section 5: IPC 
dta<-dta%>%
    mutate( 
        xipcpp= q501==1, 
        xsafe= q502==1, 
        xguideline= q504==1, 
        xppe=q506==1)%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(xipcpp, xsafe, xguideline, xppe), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    rename_all(.funs = funs(sub("q503_", "xsafe__", .)))%>%   
    rename_all(.funs = funs(sub("q505_", "xguideline__", .)))%>%        
    mutate_at(vars(starts_with("xsafe__"), starts_with("xguideline__")), 
              ~replace(., is.na(.)==TRUE, 0))

    varlist<-dta%>%select(starts_with("xsafe__"))%>%colnames()

dta<-dta%>%
    mutate(
            max=11,
            temp=rowSums(dta[ , varlist], na.rm=TRUE),         
        xsafe_score =100*(temp/max),
        xsafe_100       =xsafe_score>=100,
        xsafe_50        =xsafe_score>=50
    )

    varlist<-dta%>%select(starts_with("xguideline__"))%>%colnames()
    
dta<-dta%>%
    mutate(
        
            max=6,
            temp=rowSums(dta[colnames(select(dta, starts_with("xguideline__")))],
                         na.rm=TRUE ),  
        xguideline_score    =100*(temp/max),
        xguideline_100      =xguideline_score>=100,
        xguideline_50       =xguideline_score>=50
    )

        temp1<-dta%>%select(starts_with("q507_"))%>%
            rename_all(.funs = funs(sub("q507_", "xppe_allsome__", .)))%>%
            mutate_all( funs(ifelse(. <=2 , 1, 
                                  ifelse(. >2, 0,  .))) )%>%
            mutate_all(~replace(., is.na(.)==TRUE, 0)) 
        
        temp2<-dta%>%select(starts_with("q507_"))%>%
            rename_all(.funs = funs(sub("q507_", "xppe_all__", .)))%>%
            mutate_all( funs(ifelse(. ==1 , 1, 
                                  ifelse(. >1, 0, .))) )%>%
            mutate_all(~replace(., is.na(.)==TRUE, 0)) 

dta<-cbind(dta, temp1, temp2)

dta<-dta%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate(
            max=6,
            temp=rowSums(dta[colnames(select(dta, starts_with("xppe_allsome__")))], na.rm=TRUE ) , 
        xppe_allsome_score  =100*(temp/max),
        xppe_allsome_100        =xppe_allsome_score>=100,
        xppe_allsome_50         =xppe_allsome_score>=50, 
        
            max=6,
            temp=rowSums(dta[colnames(select(dta, starts_with("xppe_all__")))], na.rm=TRUE ) , 
        xppe_all_score  =100*(temp/max),
        xppe_all_100        =xppe_all_score>=100,
        xppe_all_50         =xppe_all_score>=50, 
        
        xppedispose = q508==1 
    )
#####* Section 6: COVID management in primary care setting 
dta<-dta%>%
    mutate(
        xcvd_team       =q601==1,
        xcvd_sop        =q602==1,
            
        xcvd_spcm       =q603==1,
        xcvd_test       =q603==1 & q604<=3,
        xcvd_test_pcr       =q603==1 & (q604==1 | q604==3),
        xcvd_test_rdt       =q603==1 & (q604==2 | q604==3),
        xcvd_spcmtrans      =q603==1 & q604==4 & q605==1,
            
        xcvd_pt = q606==1        
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xcvd_")), 
              ~replace(., is.na(.)==TRUE, 0))

            temp1<-dta%>%select(xcvd_pt, starts_with("q607_"))%>%
                rename_all(.funs = funs(sub("q607_", "xcvd_pt__", .)))%>%
                mutate_at(vars(starts_with("xcvd_pt__")), 
                          funs(ifelse(xcvd_pt==1 & .>=2, 0,
                                    ifelse(xcvd_pt==1 & .==1, 1,.))) )%>%
                mutate_all(~replace(., is.na(.), 0))%>%
                select(-xcvd_pt)

dta<-cbind(dta, temp1)

dta<-dta%>%
    mutate(
            max=7,
            temp=rowSums(dta[colnames(select(dta, starts_with("xcvd_pt__")))], na.rm=TRUE ) , 
        xcvd_pt_score   =100*(temp/max),
        xcvd_pt_100     =xcvd_pt_score>=100,
        xcvd_pt_50      =xcvd_pt_score>=50 
    )%>%
    mutate_at(vars(xcvd_pt_score, xcvd_pt_100, xcvd_pt_50, starts_with("xcvd_pt__")),
              funs(ifelse(xcvd_pt==0, NA, .)) )%>%
    mutate_if(is.logical, as.numeric)

dta<-dta%>%
    mutate(xcvd_pthbsi = (q607_006==1 | q607_006==2))%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(xcvd_pthbsi), 
              funs(ifelse(is.na(.)==TRUE, 0, .))) 

            temp1<-dta%>%select(xcvd_pthbsi, starts_with("q608_"))%>%
                rename_all(.funs = funs(sub("q608_", "xcvd_pthbsi__", .)))%>%
                mutate_at(vars(starts_with("xcvd_pthbsi__")), 
                          funs(ifelse(xcvd_pthbsi==1 & .>=2, 0,
                                    ifelse(xcvd_pthbsi==1 & .==1, 1,.))) )%>%
                mutate_all(~replace(., is.na(.), 0))%>%
                select(-xcvd_pthbsi)            

dta<-cbind(dta, temp1)%>%
    #Revision 2021/10/29. leave it as missing if xcvd_pthbsi==0
    mutate_at(vars(starts_with("xcvd_pthbsi__")), 
              funs(ifelse(xcvd_pthbsi==0, NA,.))
    )
    #ENF OF REVISION

dta<-dta%>%
    mutate(
        xcvd_guide_casemanage = q609==1 ,           
        
        xcvd_info = q610==1 ,
        
        xcvd_info_moh   = q610==1   & q611_001==1,
        xcvd_info_localgov  = q610==1   & q611_002==1,
        xcvd_info_who   = q610==1   & q611_003==1,
        xcvd_info_prof  = q610==1   & q611_004==1,
        xcvd_info_other = q610==1   & q611_005==1
    )%>%
    mutate_if(is.logical, as.numeric)
#####* Section 7: Therapeutics
dta<-dta%>%
    mutate( 
            max=maxdrug,
            temp=rowSums(dta[colnames(select(dta, starts_with("q701_")))], na.rm=TRUE ) , 
        xdrug_score =100*(temp/max),
        xdrug_100   =xdrug_score>=100,
        xdrug_50    =xdrug_score>=50, 
    
            max=3,
            temp=rowSums(dta[colnames(select(dta, starts_with("q702_")))], na.rm=TRUE ) , 
        xsupp_score =100*(temp/max),
        xsupp_100   =xsupp_score>=100,
        xsupp_50    =xsupp_score>=50,

            max=6,
            temp=rowSums(dta[colnames(select(dta, starts_with("q703_")))], na.rm=TRUE ) , 
        xvac_score  =100*(temp/max),
        xvac_100    =xvac_score>=100,
        xvac_50     =xvac_score>=50
        )%>%
    rename_all(.funs = funs(sub("q701_", "xdrug__", .)))%>%     
    rename_all(.funs = funs(sub("q702_", "xsupply__", .)))%>%   
    rename_all(.funs = funs(sub("q703_", "xvaccine__", .)))%>%   
    mutate_at(vars(starts_with("xdrug__"), starts_with("xsupply__"), starts_with("xvaccine__")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate(
        xvaccine_child=q409_005>=1 & q409_005<=3
    )%>%
    mutate_at(vars(xvaccine_child), 
              ~replace(., is.na(.)==TRUE, 0))%>%    
    mutate_at(vars(starts_with("xvac_"), starts_with("xvaccine__")), 
                          funs(ifelse(xvaccine_child!=1, NA, .)) )%>%
    mutate(
        xdisrupt_supply = q704==1
    )
            temp1<-dta%>%select(starts_with("q802_"))%>%
                rename_all(.funs = funs(sub("q802_", "xdiag_av_a", .)))%>%
                mutate_all(~replace(., . ==1 |.==2, 1))%>%
                mutate_all(~replace(., . ==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))
    
            temp2<-dta%>%select(starts_with("q803_"))%>%
                rename_all(.funs = funs(sub("q803_", "xdiag_av_h", .)))%>%
                mutate_all(~replace(., . ==1 |.==2, 1))%>%
                mutate_all(~replace(., . ==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))
            
            temp3<-dta%>%select(starts_with("q805_"))%>%
                rename_all(.funs = funs(sub("q805_", "ximage_av_", .)))%>%
                mutate_all(~replace(., . ==1 |.==2, 1))%>%
                mutate_all(~replace(., . ==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))        
            
            temp4<-dta%>%select(starts_with("q802_"))%>%
                rename_all(.funs = funs(sub("q802_", "xdiag_avfun_a", .)))%>%
                mutate_all(~replace(., . ==1, 1))%>%
                mutate_all(~replace(., . ==2 |.==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))
            
            temp5<-dta%>%select(starts_with("q803_"))%>%
                rename_all(.funs = funs(sub("q803_", "xdiag_avfun_h", .)))%>%
                mutate_all(~replace(., . ==1, 1))%>%
                mutate_all(~replace(., . ==2 |.==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))         
            
            temp6<-dta%>%select(starts_with("q805_"))%>%
                rename_all(.funs = funs(sub("q805_", "ximage_avfun_", .)))%>%
                mutate_all(~replace(., . ==1, 1))%>%
                mutate_all(~replace(., . ==2 |.==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))          

dta<-cbind(dta, temp1, temp2, temp3, temp4, temp5, temp6)

dta<-dta%>%
    mutate(
        xdiag=q801==1, 
        ximage=q804==1)%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(xdiag, ximage), 
              ~replace(., is.na(.)==TRUE, 0))%>%

    mutate(
            max=5,
            tempbasic=rowSums(dta[colnames(select(dta, starts_with("xdiag_avfun_a")))], na.rm=TRUE ) ,     
        xdiagbasic_score    =100*(tempbasic/max),
        xdiagbasic_100  =xdiagbasic_score>=100,
        xdiagbasic_50   =xdiagbasic_score>=50,
            
            max=5, 
            max=ifelse(zlevel_hospital==1 , 10, max), 
            temphosp=rowSums(dta[colnames(select(dta, starts_with("xdiag_avfun_h")))], na.rm=TRUE ) ,   
            temprelevant=tempbasic,   
            temprelevant=ifelse(zlevel_hospital==1 , tempbasic+ temphosp, temprelevant), 
        xdiag_score =100*(temprelevant/max),
        xdiag_100   =xdiag_score>=100,
        xdiag_50    =xdiag_score>=50,
    
            max=4,
            tempimage=rowSums(dta[colnames(select(dta, starts_with("ximage_avfun_")))], na.rm=TRUE ) ,     
        ximage_score    =100*(tempimage/max),
        ximage_100  =ximage_score>=100,
        ximage_50   =ximage_score>=50
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xdiag_av_h"), starts_with("xdiag_avfun_h"),
                   starts_with("ximage")), 
              funs(ifelse(zlevel_hospital!=1, NA, .) ) )
#####* Section 9: general vaccine readiness
dta<-dta%>%
    mutate( 
        xvac= q901==1 | q902==1, 
        
        xvac_av_fridge        = q903==1 | q903==2, 
        xvac_avfun_fridge     = q903==1 ,
        xvac_avfun_fridgetemp = q903==1 & q904==1,
        
        xvac_av_coldbox = q905==1,
        
        xvac_avfun_coldbox_all      = q905==1 & (q906>=1 & is.na(q906)==FALSE) & q907==1,
        xvac_avfun_coldbox_all_full = q905==1 & (q906>=1 & is.na(q906)==FALSE) & q907==1 & q911==1,
        
        yvac_avfun_coldbox_all      =NA, 
        yvac_avfun_coldbox_all      = ifelse(xvac_avfun_coldbox_all==1, q906, 
                                             yvac_avfun_coldbox_all), 
        yvac_avfun_coldbox_all_full =NA, 
        yvac_avfun_coldbox_all_full = ifelse(xvac_avfun_coldbox_all==1 & q911==1, q906, 
                                             yvac_avfun_coldbox_all_full),  
        
        xvac_av_carrier = q908==1,
        
        xvac_avfun_carrier_all      = q908==1 & (q909>=1 & is.na(q909)==FALSE) & q910==1,
        xvac_avfun_carrier_all_full = q908==1 & (q909>=1 & is.na(q909)==FALSE) & q910==1 & q911==1, 
        
        yvac_avfun_carrier_all      =NA, 
        yvac_avfun_carrier_all      = ifelse(xvac_avfun_carrier_all==1, q909, 
                                             yvac_avfun_carrier_all),   
        yvac_avfun_carrier_all_full =NA, 
        yvac_avfun_carrier_all_full = ifelse(xvac_avfun_carrier_all==1 & q911==1, q909, 
                                             yvac_avfun_carrier_all_full), 
            
        xvac_av_outreach = xvac_av_coldbox ==1 | xvac_av_carrier ==1,  
        xvac_avfun_outreach_all_full = xvac_avfun_coldbox_all_full ==1 | xvac_avfun_carrier_all_full==1,
        
        xvac_sharp = q912==1,
        
        xvac_aefikit = q913==1, 
        xvac_aefireport = q914==1 
    ) %>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xvac_av"), xvac_sharp, starts_with("xvac_aefi")), 
              ~replace(., is.na(.)==TRUE, 0))%>%
    mutate_at(vars(starts_with("xvac_av"), starts_with("yvac_av"), 
                   xvac_sharp, starts_with("xvac_aefi")), 
              funs(ifelse(xvac!=1, NA, .)) )

    #lab var xvac_av_fridge "has fridge"
    #lab var xvac_avfun_fridge "has functioning fridge"
    #lab var xvac_avfun_fridgetemp "has functioning fridge with temp log"
    
    #lab var xvac_av_coldbox "has coldbox"
    #lab var xvac_avfun_coldbox_all "has functioning coldbox, all"
    #lab var xvac_avfun_coldbox_all_full "has functioning coldbox with icepacks, all"
    
    #lab var xvac_av_carrier "has carrier"
    #lab var xvac_avfun_carrier_all "has functioning carrier, all"
    #lab var xvac_avfun_carrier_all_full "has functioning carrier with icepacks, all"
#####* Section 10: COVID vaccine readiness
dta<-dta%>%
    mutate( 
        #*This part is HIDDEN cause these were created in the above section for the mock data 
        #*In real, however, only section 9 OR 10 will be implemented. 
        #*Open up this part, if section 10 is used. 
        #*Delete this whole Section 10 if section 9 is used. 
        #   xvac_av_fridge        = q1001==1 | q1001==2, 
        #   xvac_avfun_fridge     = q1001==1 ,
        #   xvac_avfun_fridgetemp = q1001==1 & q1002==1,     
        
        xcovax=q1003==1)

            temp1<-dta%>%select(starts_with("q1004_"))%>%
                rename_all(.funs = funs(sub("q1004_", "xcovax_offer__", .)))%>%
                mutate_all(~replace(., . ==1 |.==2, 1))%>%
                mutate_all(~replace(., . ==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))        
            
            temp2<-dta%>%select(starts_with("q1004_"))%>%
                rename_all(.funs = funs(sub("q1004_", "xcovax_offerav__", .)))%>%
                mutate_all(~replace(., . ==1, 1))%>%
                mutate_all(~replace(., . ==2 |.==3, 0))%>%
                mutate_all(~replace(., is.na(.), 0))     
            
            temp3<-dta%>%select(starts_with("q1005_"))%>%
                rename_all(.funs = funs(sub("q1005_", "xcovax_train__", .)))%>%
                mutate_all(~replace(., . ==1, 1))%>%
                mutate_all(~replace(., . !=1, 0))%>%
                mutate_all(~replace(., is.na(.), 0))       

dta<-cbind(dta, temp1, temp2, temp3)%>%
    mutate(
    xcovax_syr      =q1006==1,
    xcovax_sharp    =q1007==1,
    xcovax_strtemp  =q1008==1,
    xcovax_strtemp_w=q1008==1 & q1009==1,
    
    xcovax_infrtrn      =q1010==1,
    xcovax_infrtrn = replace(xcovax_infrtrn, 
                             (q1004_001==3 & q1004_002==3 & q1004_003==3),
                             NA),  
    xcovax_infside      =q1011==1,
    xcovax_infaewhat    =q1012==1,
    
    xcovax_aefikit      = q1013==1,
    xcovax_aefireport   = q1014==1          
    )%>%
    mutate_if(is.logical, as.numeric)%>%
    mutate_at(vars(starts_with("xcovax_")), 
              funs(ifelse(xcovax!=1, NA, .) ) )
            temp<-dta%>%select(starts_with("qa1_002_"))
            temp1<-temp%>%select(1:4)%>%
                rename_all(.funs = funs(sub("qa1_002_", "vol_opt_now_", .)))
            temp2<-temp%>%select(5:8)%>%
                rename_all(.funs = funs(sub("qa1_002_", "vol_opt_last_", .)))

dta<-cbind(dta, temp1, temp2)

            temp<-dta%>%select(starts_with("qa1_003_"))
            temp1<-temp%>%select(1:4)%>%
                rename_all(.funs = funs(sub("qa1_003_", "vol_ipt_now_", .)))
            temp2<-temp%>%select(5:8)%>%
                rename_all(.funs = funs(sub("qa1_003_", "vol_ipt_last_", .)))

dta<-cbind(dta, temp1, temp2)

            temp<-dta%>%select(starts_with("qa1_004_"))
            temp1<-temp%>%select(1:4)%>%
                rename_all(.funs = funs(sub("qa1_004_", "vol_del_now_", .)))
            temp2<-temp%>%select(5:8)%>%
                rename_all(.funs = funs(sub("qa1_004_", "vol_del_last_", .)))

dta<-cbind(dta, temp1, temp2)

            temp<-dta%>%select(starts_with("qa1_005_"))
            temp1<-temp%>%select(1:4)%>%
                rename_all(.funs = funs(sub("qa1_005_", "vol_dpt_now_", .)))
            temp2<-temp%>%select(5:8)%>%
                rename_all(.funs = funs(sub("qa1_005_", "vol_dpt_last_", .)))

dta<-cbind(dta, temp1, temp2)

E.2-Addendum. Rename indicators for easier data use

Rename indicators ending with sub-question numbers with more friendly names. These names are used in the dash board. Thus, it is important to ensure the indicator names are correct, if questionnaire is adapted beyond minimum requirements. (Addendum on August 17, 2021)

#dta%>%select(grep("__", names(dta)))%>%colnames() 
dta<-dta%>%
    rename(  
        xtraining__ipc   =  xtraining__001   , 
        xtraining__ppe   =  xtraining__002   , 
        xtraining__triage    =  xtraining__003   , 
        xtraining__emerg     =  xtraining__004   , 
        xtraining__remote    =  xtraining__005   , 
        xtraining__mental    =  xtraining__006   , 
        xtraining__ss_ipc    =  xtraining__007   , 
        xtraining__ss_ppe    =  xtraining__008   , 
        xtraining__ss_c19cm  =  xtraining__009   , 
                    
        xopt_increase__undiff    =  xopt_increase__001   , 
        xopt_increase__fp    =  xopt_increase__002   , 
        xopt_increase__anc   =  xopt_increase__003   , 
        xopt_increase__pnc   =  xopt_increase__004   , 
        xopt_increase__immun     =  xopt_increase__005   , 
        xopt_increase__sickchild     =  xopt_increase__006   , 
        xopt_increase__hiv   =  xopt_increase__007   , 
        xopt_increase__tb    =  xopt_increase__008   , 
        xopt_increase__std   =  xopt_increase__009   , 
        xopt_increase__malaria   =  xopt_increase__010   , 
        xopt_increase__cvd   =  xopt_increase__011   , 
        xopt_increase__crd   =  xopt_increase__012   , 
        xopt_increase__diabetes  =  xopt_increase__013   , 
        xopt_increase__cancer    =  xopt_increase__014   , 
        xopt_increase__mental    =  xopt_increase__015   , 
        xopt_increase__violence  =  xopt_increase__016   , 
        xopt_increase__ntd   =  xopt_increase__017   , 
        xopt_increase__rehab     =  xopt_increase__018   , 
        xopt_decrease__undiff    =  xopt_decrease__001   , 
        xopt_decrease__fp    =  xopt_decrease__002   , 
        xopt_decrease__anc   =  xopt_decrease__003   , 
        xopt_decrease__pnc   =  xopt_decrease__004   , 
        xopt_decrease__immun     =  xopt_decrease__005   , 
        xopt_decrease__sickchild     =  xopt_decrease__006   , 
        xopt_decrease__hiv   =  xopt_decrease__007   , 
        xopt_decrease__tb    =  xopt_decrease__008   , 
        xopt_decrease__std   =  xopt_decrease__009   , 
        xopt_decrease__malaria   =  xopt_decrease__010   , 
        xopt_decrease__cvd   =  xopt_decrease__011   , 
        xopt_decrease__crd   =  xopt_decrease__012   , 
        xopt_decrease__diabetes  =  xopt_decrease__013   , 
        xopt_decrease__cancer    =  xopt_decrease__014   , 
        xopt_decrease__mental    =  xopt_decrease__015   , 
        xopt_decrease__violence  =  xopt_decrease__016   , 
        xopt_decrease__ntd   =  xopt_decrease__017   , 
        xopt_decrease__rehab     =  xopt_decrease__018   , 
        xopt_nochange__undiff    =  xopt_nochange__001   , 
        xopt_nochange__fp    =  xopt_nochange__002   , 
        xopt_nochange__anc   =  xopt_nochange__003   , 
        xopt_nochange__pnc   =  xopt_nochange__004   , 
        xopt_nochange__immun     =  xopt_nochange__005   , 
        xopt_nochange__sickchild     =  xopt_nochange__006   , 
        xopt_nochange__hiv   =  xopt_nochange__007   , 
        xopt_nochange__tb    =  xopt_nochange__008   , 
        xopt_nochange__std   =  xopt_nochange__009   , 
        xopt_nochange__malaria   =  xopt_nochange__010   , 
        xopt_nochange__cvd   =  xopt_nochange__011   , 
        xopt_nochange__crd   =  xopt_nochange__012   , 
        xopt_nochange__diabetes  =  xopt_nochange__013   , 
        xopt_nochange__cancer    =  xopt_nochange__014   , 
        xopt_nochange__mental    =  xopt_nochange__015   , 
        xopt_nochange__violence  =  xopt_nochange__016   , 
        xopt_nochange__ntd   =  xopt_nochange__017   , 
        xopt_nochange__rehab     =  xopt_nochange__018   , 
                    
        xopt_increase_reason__more_ari   =  xopt_increase_reason__001    , 
        xopt_increase_reason__redirect   =  xopt_increase_reason__002    , 
        xopt_increase_reason__backlog    =  xopt_increase_reason__003    , 
        xopt_increase_reason__react  =  xopt_increase_reason__004    , 
        xopt_increase_reason__comms  =  xopt_increase_reason__005    , 
        xopt_increase_reason__gbv    =  xopt_increase_reason__006    , 
        xopt_increase_reason__other  =  xopt_increase_reason__007    , 
        xopt_decrease_reason__changerecs     =  xopt_decrease_reason__001    , 
        xopt_decrease_reason__fear   =  xopt_decrease_reason__002    , 
        xopt_decrease_reason__lockdown   =  xopt_decrease_reason__003    , 
        xopt_decrease_reason__transport  =  xopt_decrease_reason__004    , 
                    
        xopt_decrease_reason__othercom   =  xopt_decrease_reason__005    , 
        xopt_decrease_reason__servreduc  =  xopt_decrease_reason__006    , 
        xopt_decrease_reason__disrupt    =  xopt_decrease_reason__007    , 
        xopt_decrease_reason__hours  =  xopt_decrease_reason__008    , 
        xopt_decrease_reason__closure    =  xopt_decrease_reason__009    , 
        xopt_decrease_reason__drugs  =  xopt_decrease_reason__010    , 
        xopt_decrease_reason__staff  =  xopt_decrease_reason__011    , 
        xopt_decrease_reason__otherfac   =  xopt_decrease_reason__012    , 
                    
        xer_increase__injuries   =  xer_increase__002    , 
        xer_increase__surgeries  =  xer_increase__003    , 
        xer_increase__ncd    =  xer_increase__004    , 
        xer_increase__blood  =  xer_increase__005    , 
        xer_decrease__injuries   =  xer_decrease__002    , 
        xer_decrease__surgeries  =  xer_decrease__003    , 
        xer_decrease__ncd    =  xer_decrease__004    , 
        xer_decrease__blood  =  xer_decrease__005    , 
                    
        xout_decrease__immun     =  xout_decrease__001   , 
        xout_decrease__malaria   =  xout_decrease__002   , 
        xout_decrease__ntd   =  xout_decrease__003   , 
        xout_decrease__cbc   =  xout_decrease__004   , 
        xout_decrease__home      =  xout_decrease__005   , 
                    
        xsafe__entrance_screening    =  xsafe__001   , 
        xsafe__staff_entrance    =  xsafe__002   , 
        xsafe__sep_room  =  xsafe__003   , 
        xsafe__triage_c19    =  xsafe__004   , 
        xsafe__isolatareas   =  xsafe__005   , 
        xsafe__triage_guidelines     =  xsafe__006   , 
        xsafe__distancing    =  xsafe__007   , 
        xsafe__hygiene_instructions  =  xsafe__008   , 
        xsafe__hygiene_stations  =  xsafe__009   , 
        xsafe__ppe   =  xsafe__010   , 
        xsafe__cleaning  =  xsafe__011   , 
                    
        xguideline__screening    =  xguideline__001  , 
        xguideline__c19_manage   =  xguideline__002  , 
        xguideline__ppe  =  xguideline__003  , 
        xguideline__c19_surveillance     =  xguideline__004  , 
        xguideline__deadbody     =  xguideline__005  , 
        xguideline__waste    =  xguideline__006  , 
                    
        xppe_allsome__gown   =  xppe_allsome__001    , 
        xppe_allsome__gloves     =  xppe_allsome__002    , 
        xppe_allsome__goggles    =  xppe_allsome__003    , 
        xppe_allsome__faceshield     =  xppe_allsome__004    , 
        xppe_allsome__respirator     =  xppe_allsome__005    , 
        xppe_allsome__mask   =  xppe_allsome__006    , 
                    
        xppe_all__gown   =  xppe_all__001    , 
        xppe_all__gloves     =  xppe_all__002    , 
        xppe_all__goggles    =  xppe_all__003    , 
        xppe_all__faceshield     =  xppe_all__004    , 
        xppe_all__respirator     =  xppe_all__005    , 
        xppe_all__mask   =  xppe_all__006    , 
                    
        xcvd_pt__sep_room    =  xcvd_pt__001     , 
        xcvd_pt__c19_check   =  xcvd_pt__002     , 
        xcvd_pt__o2_measure  =  xcvd_pt__003     , 
        xcvd_pt__refer   =  xcvd_pt__004     , 
        xcvd_pt__diagtest    =  xcvd_pt__005     , 
        xcvd_pt__home_isolate    =  xcvd_pt__006     , 
        xcvd_pt__tele    =  xcvd_pt__007     , 
                    
        xcvd_pthbsi__tele    =  xcvd_pthbsi__001     , 
        xcvd_pthbsi__homevisit   =  xcvd_pthbsi__002     , 
        xcvd_pthbsi__followup_fac    =  xcvd_pthbsi__003     , 
        xcvd_pthbsi__safetyinstruct  =  xcvd_pthbsi__004     , 
        xcvd_pthbsi__complianceassess    =  xcvd_pthbsi__005     , 
        xcvd_pthbsi__report  =  xcvd_pthbsi__006     , 
                    
        xdrug__salbutamol    =  xdrug__001   , 
        xdrug__metformin     =  xdrug__002   , 
        xdrug__hydrochlorothiazide   =  xdrug__003   , 
        xdrug__paracetamol_susp  =  xdrug__004   , 
        xdrug__carbamazapine     =  xdrug__005   , 
        xdrug__amoxicillin_tabs  =  xdrug__006   , 
        xdrug__oralcontracept    =  xdrug__007   , 
        xdrug__oxytocin  =  xdrug__008   , 
        xdrug__magnesiumsulphate     =  xdrug__009   , 
        xdrug__heparin   =  xdrug__010   , 
        xdrug__hydrocortisone    =  xdrug__011   , 
        xdrug__epinephrine   =  xdrug__012   , 
        xdrug__artemether    =  xdrug__013   , 
        xdrug__efavirenz     =  xdrug__014   , 
        xdrug__isoniazid     =  xdrug__015   , 
        xdrug__ivfluids  =  xdrug__016   , 
        xdrug__oxygen    =  xdrug__017   , 
                    
        xsupply__syringes    =  xsupply__001     , 
        xsupply__ivsets  =  xsupply__002     , 
        xsupply__gauze   =  xsupply__003     , 
                        
        xvaccine__mcv    =  xvaccine__001    , 
        xvaccine__dtp    =  xvaccine__002    , 
        xvaccine__polio  =  xvaccine__003    , 
        xvaccine__bcg    =  xvaccine__004    , 
        xvaccine__pneumo     =  xvaccine__005    , 
                        
        xdiag_av__malaria    =  xdiag_av_a001    , 
        xdiag_av__bloodglucose   =  xdiag_av_a002    , 
        xdiag_av__urineglucose   =  xdiag_av_a003    , 
        xdiag_av__urineprotein   =  xdiag_av_a004    , 
        xdiag_av__pregnancy  =  xdiag_av_a005    , 
        xdiag_avfun__malaria     =  xdiag_avfun_a001     , 
        xdiag_avfun__bloodglucose    =  xdiag_avfun_a002     , 
        xdiag_avfun__urineglucose    =  xdiag_avfun_a003     , 
        xdiag_avfun__urineprotein    =  xdiag_avfun_a004     , 
        xdiag_avfun__pregnancy   =  xdiag_avfun_a005     , 
        xdiag_av__h_hiv  =  xdiag_av_h001    , 
        xidag_av__h_tb   =  xdiag_av_h002    , 
        xdiag_av__h_hbg  =  xdiag_av_h003    , 
        xdiag_av__h_bloodtype    =  xdiag_av_h004    , 
        xdiag_av__h_bloodcreatine    =  xdiag_av_h005    , 
        xdiag_avfun__h_hiv   =  xdiag_avfun_h001     , 
        xidag_avfun__h_tb    =  xdiag_avfun_h002     , 
        xdiag_avfun__h_hbg   =  xdiag_avfun_h003     , 
        xdiag_avfun__h_bloodtype     =  xdiag_avfun_h004     , 
        xdiag_avfun__h_bloodcreatine     =  xdiag_avfun_h005     , 
                    
        ximage_av__xray  =  ximage_av_001    , 
        ximage_av__mri   =  ximage_av_002    , 
        ximage_av__ultrasound    =  ximage_av_003    , 
        ximage_av__ct    =  ximage_av_004    , 
        ximage_avfun__xray   =  ximage_avfun_001     , 
        ximage_avfun__mri    =  ximage_avfun_002     , 
        ximage_avfun__ultrasound     =  ximage_avfun_003     , 
        ximage_avfun__ct     =  ximage_avfun_004     , 
                    
        xcovax_offer__pfizer     =   xcovax_offer__001   , 
        xcovax_offerav__pfizer   =   xcovax_offerav__001     , 
        xcovax_offer__moderna    =   xcovax_offer__002   , 
        xcovax_offerav__moderna  =   xcovax_offerav__002     , 
        xcovax_offer__astra  =   xcovax_offer__003   , 
        xcovax_offerav__astra    =   xcovax_offerav__003     , 
        xcovax_offer__jj     =   xcovax_offer__004   , 
        xcovax_offerav__jj   =   xcovax_offerav__004     , 
                    
        xcovax_train__storage    =   xcovax_train__001   , 
        xcovax_train__admin  =   xcovax_train__002   , 
        xcovax_train__manage_adverse     =   xcovax_train__003   , 
        xcovax_train__report_adverse     =   xcovax_train__004   , 
                    
        vol_opt_now_m1_opt   =  vol_opt_now_001  , 
        vol_opt_now_m2_opt   =  vol_opt_now_002  , 
        vol_opt_now_m3_opt   =  vol_opt_now_003  , 
        vol_opt_now_m4_opt   =  vol_opt_now_004  , 
        vol_opt_last_m1_opt  =  vol_opt_last_005     , 
        vol_opt_last_m2_opt  =  vol_opt_last_006     , 
        vol_opt_last_m3_opt  =  vol_opt_last_007     , 
        vol_opt_last_m4_opt  =  vol_opt_last_008     , 
        vol_ipt_now_m1_ipt   =  vol_ipt_now_001  , 
        vol_ipt_now_m2_ipt   =  vol_ipt_now_002  , 
        vol_ipt_now_m3_ipt   =  vol_ipt_now_003  , 
        vol_ipt_now_m4_ipt   =  vol_ipt_now_004  , 
        vol_ipt_last_m1_ipt  =  vol_ipt_last_005     , 
        vol_ipt_last_m2_ipt  =  vol_ipt_last_006     , 
        vol_ipt_last_m3_ipt  =  vol_ipt_last_007     , 
        vol_ipt_last_m4_ipt  =  vol_ipt_last_008     , 
        vol_del_now_m1_del   =  vol_del_now_001  , 
        vol_del_now_m2_del   =  vol_del_now_002  , 
        vol_del_now_m3_del   =  vol_del_now_003  , 
        vol_del_now_m4_del   =  vol_del_now_004  , 
        vol_del_last_m1_del  =  vol_del_last_005     , 
        vol_del_last_m2_del  =  vol_del_last_006     , 
        vol_del_last_m3_del  =  vol_del_last_007     , 
        vol_del_last_m4_del  =  vol_del_last_008     , 
        vol_dpt_now_m1_dpt   =  vol_dpt_now_001  , 
        vol_dpt_now_m2_dpt   =  vol_dpt_now_002  , 
        vol_dpt_now_m3_dpt   =  vol_dpt_now_003  , 
        vol_dpt_now_m4_dpt   =  vol_dpt_now_004  , 
        vol_dpt_last_m1_dpt  =  vol_dpt_last_005     , 
        vol_dpt_last_m2_dpt  =  vol_dpt_last_006     , 
        vol_dpt_last_m3_dpt  =  vol_dpt_last_007     , 
        vol_dpt_last_m4_dpt  =  vol_dpt_last_008     
        
    )   

#dta%>%select(grep("__", names(dta)))%>%colnames() 

E.3. Merge with sampling weight

Select Option A (if threre are sampling weights) or Option B (if there is no sampling weight).
Option A

# read sampling weight in the chartbook provided by the country. Makesure there are no duplicates 
dtaweight<-read_excel(paste0(chartbookdir, "WHO_CEHS_Chartbook_10.21.xlsx"), sheet = "Weight")

    names(dtaweight)<-tolower(names(dtaweight))

    dtaweight<-dtaweight%>%
        rename_all(.funs = funs(sub(" ", "", .)))%>%
        select(facilitycode, weight)%>%
        distinct(facilitycode, .keep_all=TRUE)
    
    #normalize weight 
    dtaweight<-dtaweight%>%
        mutate(weight = weight / mean(dtaweight$weight))    

# check datasets     
dim(dta)
dim(dtaweight)
str(dta$facilitycode)
str(dtaweight$facilitycode)

# check datasets 
dta<-left_join(dta, dtaweight, by = c("facilitycode"))%>%
    arrange(id)#/*this is generated from Lime survey*/

# confirm dimension 
dim(dta)

Option B

# Give "sampling weights" 1 to all observations (edit 6/14/2021)
dta<-dta%>%
    mutate(weight=1)

E.4. Export clean Respondent-level data for chartbook
write.csv(dta, paste0("CEHS_", country, "_R", round, ".csv"))
dta<-dta%>%mutate(test=Sys.time())

wb <- loadWorkbook(paste0(chartbookdir, "WHO_CEHS_Chartbook.xlsx"))
writeData(wb, sheet = "Facility-level cleaned data", dta)
saveWorkbook(wb,paste0(chartbookdir, "WHO_CEHS_Chartbook.xlsx"),overwrite = TRUE)

F. Create and export indicator estimate data

F.1. Calculate estimates
dtatemp<-dta%>%
    mutate(
        obs=1,  
        
        obs_userfee=NA,
        obs_ipt=NA,
        obs_er=NA,  
        obs_vac=NA, 
        obs_covax=NA, 
        obs_primary=NA, 
        obs_hospital=NA, 
        
        obs_userfee=ifelse( xuserfee==1, 1, obs_userfee),
        obs_ipt=ifelse( xipt==1, 1, obs_ipt),
        obs_er=ifelse( xer==1, 1, obs_er),
        obs_vac=ifelse( xvac==1, 1, obs_vac),
        obs_covax=ifelse( xcovax==1, 1, obs_covax),
        obs_primary=ifelse( zlevel_low==1, 1, obs_primary), 
        obs_hospital=ifelse( zlevel_low==0, 1, obs_hospital), 
        
        #REVISION 2021/10/29 : CREATE ADDITIONAL "OBS" VARIABLES TO HAVE CORRECT DENOMINATORS IN THE CHARTBOOK 
        obs_cvd_pt=NA,          
        obs_cvd_pthbsi=NA,      
        obs_cvd_pt=ifelse( xcvd_pt==1,1,obs_cvd_pt),
        obs_cvd_pthbsi=ifelse( xcvd_pthbsi==1,1,obs_cvd_pthbsi),
        # END OF REVISION 
        
        obshmis_opt=rowSums(dta[colnames(select(dta, starts_with("vol_opt_")))], na.rm=TRUE )>0 , 
        obshmis_ipt=rowSums(dta[colnames(select(dta, starts_with("vol_ipt_")))], na.rm=TRUE )>0 , 
        obshmis_del=rowSums(dta[colnames(select(dta, starts_with("vol_del_")))], na.rm=TRUE )>0 , 
        obshmis_dpt=rowSums(dta[colnames(select(dta, starts_with("vol_dpt_")))], na.rm=TRUE )>0 
    )

dtatempx<-dtatemp%>%select(country, round, month, year, weight, starts_with("z"),
                           starts_with("x"))

dtatempy<-dtatemp%>%
    select(country, round, month, year, weight, starts_with("z"), 
           starts_with("obs"), starts_with("y"), starts_with("staff"),
           starts_with("vol"))%>%
    mutate_if(is.logical, as.numeric)
# Among all facilities 
dtasummaryx<-dtatempx%>%
    group_by(country, round, month, year)%>%
    summarize_at(vars(starts_with("x")),
                 funs(weighted.mean(., weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(group="All", grouplabel="")

dtasummaryy<-dtatempy%>%
    group_by(country, round, month, year)%>%
    summarize_at(vars(starts_with("obs"), starts_with("yvac"), 
                      starts_with("staff"), starts_with("vol")),
                 funs(sum(.*weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(group="All", grouplabel="")

dtasummaryall<-left_join(dtasummaryx, dtasummaryy, 
                      by = c("country", "round", "month", "year", "group", "grouplabel"))

# By residential area
dtasummaryx<-dtatempx%>%
    group_by(country, round, month, year, zurban)%>%
    summarize_at(vars(starts_with("x")),
                 funs(weighted.mean(., weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(
        group="Location", 
        grouplabel="",
        grouplabel= ifelse(zurban==0, "1.1 Rural" , grouplabel),
        grouplabel= ifelse(zurban==1, "1.2 Urban" , grouplabel) )

dtasummaryy<-dtatempy%>%
    group_by(country, round, month, year, zurban)%>%
    summarize_at(vars(starts_with("obs"), starts_with("yvac"), 
                      starts_with("staff"), starts_with("vol")),
                 funs(sum(.*weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(
        group="Location", 
        grouplabel="",
        grouplabel= ifelse(zurban==0, "1.1 Rural" , grouplabel),
        grouplabel= ifelse(zurban==1, "1.2 Urban" , grouplabel) )

dtasummarylocation<-left_join(dtasummaryx, dtasummaryy, 
                      by = c("country", "round", "month", "year", "group", "grouplabel"))%>%
    select(-starts_with("z"))

# By facility type
dtasummaryx<-dtatempx%>%
    group_by(country, round, month, year, zlevel_hospital)%>%
    summarize_at(vars(starts_with("x")),
                 funs(weighted.mean(., weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(
        group="Level", 
        grouplabel="",
        grouplabel= ifelse(zlevel_hospital==0, "2.1 Non-hospitals" , grouplabel),
        grouplabel= ifelse(zlevel_hospital==1, "2.2 Hospitals" , grouplabel) )

dtasummaryy<-dtatempy%>%
    group_by(country, round, month, year, zlevel_hospital)%>%
    summarize_at(vars(starts_with("obs"), starts_with("yvac"), 
                      starts_with("staff"), starts_with("vol")),
                 funs(sum(.*weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(
        group="Level", 
        grouplabel="",
        grouplabel= ifelse(zlevel_hospital==0, "2.1 Non-hospitals" , grouplabel),
        grouplabel= ifelse(zlevel_hospital==1, "2.2 Hospitals" , grouplabel) )

dtasummarylevel<-left_join(dtasummaryx, dtasummaryy, 
                      by = c("country", "round", "month", "year", "group", "grouplabel"))%>%
    select(-starts_with("z"))

# By facility managing authority
dtasummaryx<-dtatempx%>%
    group_by(country, round, month, year, zpub)%>%
    summarize_at(vars(starts_with("x")),
                 funs(weighted.mean(., weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(
        group="Sector", 
        grouplabel="",
        grouplabel= ifelse(zpub==0, "3.2 Non-public" , grouplabel),
        grouplabel= ifelse(zpub==1, "3.1 Public" , grouplabel) )

dtasummaryy<-dtatempy%>%
    group_by(country, round, month, year, zpub)%>%
    summarize_at(vars(starts_with("obs"), starts_with("yvac"), 
                      starts_with("staff"), starts_with("vol")),
                 funs(sum(.*weight, na.rm = TRUE)))%>%
    ungroup()%>%
    mutate(
        group="Sector", 
        grouplabel="",
        grouplabel= ifelse(zpub==0, "3.2 Non-public" , grouplabel),
        grouplabel= ifelse(zpub==1, "3.1 Public" , grouplabel) )

dtasummarysector<-left_join(dtasummaryx, dtasummaryy, 
                      by = c("country", "round", "month", "year", "group", "grouplabel"))%>%
    select(-starts_with("z"))

# Append all 
dim(dtasummaryall)
dim(dtasummarylocation)
dim(dtasummarylevel)
dim(dtasummarysector)

dtasummary<-rbind(dtasummaryall, dtasummarylocation, dtasummarylevel, dtasummarysector)
dim(dtasummary)

dtasummary<-dtasummary%>%
    select(-starts_with("z"))%>%
    mutate_at(vars(starts_with("x")), 
              funs(round((. * 100), 0)))%>%
    mutate_at(vars(xbedrate, ends_with("_score"), ends_with("_num")), 
              funs(round((. / 100), 0)))%>%
    mutate(
        #***** generate staff infection rates useing the pooled data    
        staff_pct_covid_md = round(100* (staff_num_covid_md / staff_num_total_md ), 1),
        staff_pct_covid_nr = round(100* (staff_num_covid_nr / staff_num_total_nr ), 1),
        staff_pct_covid_othclinical = round(100* (staff_num_covid_othclinical / staff_num_total_othclinical ), 1),
        staff_pct_covid_clinical = round(100* (staff_num_covid_clinical / staff_num_total_clinical ), 1),
        staff_pct_covid_nonclinical = round(100* (staff_num_covid_nonclinical / staff_num_total_nonclinical ), 1),
        staff_pct_covid_all = round(100* (staff_num_covid_all / staff_num_total_all ), 1), 
        
        #***** generate COVID-19 vaccine among staff using the pooled data  
        staff_pct_covaxany  = round(100* (staff_num_covaxany  / staff_num_total_all), 0),
        staff_pct_covaxfull = round(100* (staff_num_covaxfull / staff_num_total_all), 0)
        
    )%>%
    #round number of observations, in case sampling weight was used (edit 5/22/2021)
    mutate_at(vars(starts_with("obs")), 
              funs(round((.), 0)))%>%
    #drop categorial var
    select(-xopt_change)%>% 
    arrange(country, round, group, grouplabel)

#move up identification variables 
setcolorder(dtasummary, 
            c("country", "round", "month", "year", "group", "grouplabel", "obs"))
colnames(dtasummary[ , 1:10])

#deal with NaN
dtasummary[sapply(dtasummary, is.nan)] <- NA

F.2. Export indicator estimate data for chartbook
dtasummary<-dtasummary%>%
    mutate(
        updatedate = as.Date(Sys.time(  ), format='%d%b%Y'), 
        updatetime = Sys.time()
    )

write.csv(dtasummary, 
          paste0("summary_CEHS_", country, "_R", round, ".csv"))
wb <- loadWorkbook(paste0(chartbookdir, "WHO_CEHS_Chartbook.xlsx"))
writeData(wb, sheet = "Indicator estimate data", dtasummary)
saveWorkbook(wb,paste0(chartbookdir, "WHO_CEHS_Chartbook.xlsx"),overwrite = TRUE)
# for cross check against Stata results 
write.csv(dtasummary, 
          paste0("summary_CEHS_", country, "_R", round, "_R.csv"), 
          na="")

END OF MARKDOWN FILE