# Purpose of the program:
# =======================
# This program builds a two wave panel (wave 15 and 16) with variables related to preferred hours of work.
#
# Created by Nicole Watson & Mossamet Nesa
# Date: 12/11/2024 

rm(list=ls(all=TRUE)) # clear out the workspace

install.packages("haven", method="wininet", dependencies=TRUE)
install.packages("tidyverse", method="wininet", dependencies=TRUE)
install.packages("survey", method="wininet", dependencies=TRUE)

library(haven)
library(tidyverse)
library(survey)

startwave <- 15  # Initial wave data files to extract.
endwave <- 16    # Last wave data files to extract.
maxwave <- 24    # Update to the latest wave.
rls <- 240       # Update to the latest release.
origdatdir <- paste0("H:/HILDA/Release ",rls,"/files/STATA ",rls,"c") # Location of original HILDA data files in Stat format
newdatadir <- "H:/HILDA Working/Release 24/Weightingexample/R" # Location of writing new data files


# Section 1: Read in datasets and create "long" responding person longitudinal file 

setwd(origdatdir)
var <- c("xwaveid", "esbrd", "jbhruc", "jbhrcpr", "jbprhr", "hhwtrp", "xhhstrat", "xhhraid", 
         "rwrp1", "rwrp2", "rwrp3", "rwrp4", "rwrp5", "rwrp6", "rwrp7", "rwrp8", "rwrp9", "rwrp10",
         "rwrp11", "rwrp12", "rwrp13", "rwrp14", "rwrp15", "rwrp16", "rwrp17", "rwrp18", "rwrp19", "rwrp20", 
         "rwrp21", "rwrp22", "rwrp23", "rwrp24", "rwrp25", "rwrp26", "rwrp27", "rwrp28", "rwrp29", "rwrp30", 
         "rwrp31", "rwrp32", "rwrp33", "rwrp34", "rwrp35", "rwrp36", "rwrp37", "rwrp38", "rwrp39", "rwrp40", 
         "rwrp41", "rwrp42", "rwrp43", "rwrp44", "rwrp45")

for( i in startwave:endwave) {
  file_list <- paste0("Rperson_", letters[i], rls, "c.dta")
  temp <- read_dta(file_list)
  var_add <- paste0(letters[i], var) # Add wave letter onto the variable names
  #temp <- temp %>% dplyr::select(xwaveid, any_of(var_add))
  temp <- temp %>% dplyr::select(xwaveid, xhhstrat, xhhraid, any_of(var_add))
  # any_of() lets the program avoid selecting the variable not included in a specific wave and set NA to that variable. eg: "hwhmhl" not included in wave 1
  names(temp)[-1] <-substring(names(temp)[-1], 2) # Remove wave letter from variable names except for xwaveid, xhhstrat and xhhraid
  temp$wave <- i
  if (i == startwave ){
    longfile <- temp
  } else {
    longfile <- bind_rows(longfile, temp) # Append the data file from each wave
  }
}

# Save new data set
setwd(newdatadir)
save(longfile, file = "hrpref.Rdata")

longfile<-longfile %>% rename("xhhstrat"="hhstrat")
longfile<-longfile %>% rename("xhhraid"="hhraid")

# add on longitudinal weight for the pair of waves (15&16)
setwd(origdatdir)
longwts_file <- paste0("longitudinal_weights_", letters[maxwave], rls,"c.dta")
longwts <- read_dta(longwts_file)
longwts <- longwts[c("xwaveid", "wlrop")] 

final_data <- merge(longfile, longwts,  by = "xwaveid", all.x = TRUE)
final_data <- final_data[order(final_data$xwaveid, final_data$wave), ] # Sort the dataset by xwaveid and wave

# Save new data set
setwd(newdatadir)
save(final_data, file = "hrpref.Rdata")


# Section 2: prepare data for analysis. Create indicator variables for hours preferences

# want to change hours
# restrict to employed and exclude a few people who did not know or refused

tab1<-table(final_data$jbhrcpr[final_data$wave==15 & final_data$esbrd==1 & final_data$jbhrcpr>=1])
tab2<-table(final_data$jbhrcpr[final_data$wave==16 & final_data$esbrd==1 & final_data$jbhrcpr>=1])
cbind(tab1,tab2)

#hours worked

sum1<-summary(final_data$jbhruc[final_data$wave==15 & final_data$esbrd==1])
sum2<-summary(final_data$jbhruc[final_data$wave==16 & final_data$esbrd==1])
cbind(sum1,sum2)
# note some cases with missing value codes (-6 to -1)

#Create indicator variables

preffew<-ifelse(final_data$jbhrcpr==1, 1, ifelse(final_data$jbhrcpr>0, 0, NA))
prefsame<-ifelse(final_data$jbhrcpr==2, 1, ifelse(final_data$jbhrcpr>0, 0, NA))
prefmore<-ifelse(final_data$jbhrcpr==3, 1, ifelse(final_data$jbhrcpr>0, 0, NA))

# create hours variable
hours<-ifelse(final_data$esbrd==1 & final_data$jbhruc>=0, final_data$jbhruc, NA)


# Section 3: unweighted means
# tab of cases
tab3<-cbind(table(preffew[final_data$wave==15]), table(prefsame[final_data$wave==15]),
            table(prefmore[final_data$wave==15]))

tab4<-cbind(table(preffew[final_data$wave==16]), table(prefsame[final_data$wave==16]),
            table(prefmore[final_data$wave==16]))

cbind(tab3, tab4)

# unweighted proportion of whether want to change hours
cbind(summary(preffew[final_data$wave==15]), summary(prefsame[final_data$wave==15]),summary(prefmore[final_data$wave==15]),
      summary(preffew[final_data$wave==16]), summary(prefsame[final_data$wave==16]),summary(prefmore[final_data$wave==16]))

# unweighted mean of usual hours worked (for those employed)
cbind(summary(hours[final_data$wave==15]),summary(hours[final_data$wave==16]))


# Section 4A: weighted means (weighted by cross-sectional weight). Apply Taylor series linerisation method (option 1)

# need to convert xhhraid to numeric
final_data$xhhraid<-as.numeric(final_data$xhhraid)

# set up the data by specifying the weight, stratification and primary sampling unit 
hildadesign<-svydesign(id=~xhhraid, strata=~xhhstrat, weights=~hhwtrp, nest=TRUE,  data=final_data)

# population estimates of whether want to change hours
preffew1<-ifelse(final_data$wave==15 &final_data$esbrd==1 & final_data$jbhruc>=0, preffew, NA)
prefsame1<-ifelse(final_data$wave==15 &final_data$esbrd==1 & final_data$jbhruc>=0, prefsame, NA)
prefmore1<-ifelse(final_data$wave==15 &final_data$esbrd==1 & final_data$jbhruc>=0, prefmore, NA)
svymean(~preffew1, hildadesign, na.rm=TRUE)
svymean(~prefsame1, hildadesign, na.rm=TRUE)
svymean(~prefmore1, hildadesign, na.rm=TRUE)
confint(svymean(~preffew1, hildadesign, na.rm=TRUE))
confint(svymean(~prefsame1, hildadesign, na.rm=TRUE))
confint(svymean(~prefmore1, hildadesign, na.rm=TRUE))

jbhruc_new<-ifelse(final_data$wave==15 & final_data$esbrd==1 & final_data$jbhruc>=0,final_data$jbhruc, NA) 
svymean(~jbhruc_new, hildadesign, na.rm=TRUE)
confint(svymean(~jbhruc_new, hildadesign, na.rm=TRUE))


# Section 4B: apply replicate weights (option 2)

hildarepwt<-svrepdesign(data=final_data, weight=~hhwtrp, repweights="rwrp[1-45]", type="JK1", scale=44/45, mse=TRUE)

svymean(~preffew1, hildarepwt, na.rm=TRUE)
svymean(~prefsame1, hildarepwt, na.rm=TRUE )
svymean(~prefmore1, hildarepwt, na.rm=TRUE )
confint(svymean(~preffew1, hildarepwt, na.rm=TRUE))
confint(svymean(~prefsame1, hildarepwt, na.rm=TRUE))
confint(svymean(~prefmore1, hildarepwt, na.rm=TRUE))

svymean(~jbhruc_new, hildarepwt, na.rm=TRUE)
confint(svymean(~jbhruc_new, hildarepwt, na.rm=TRUE))


# Section 4C: apply longitudinal weights (in this case, the paired wave longitudinal weight for waves 15 and 16)

# set up the data as a timeseries dataset so can look forward to the values in the next wave
# to get the hours preference in next wave

pref_time<-ts(final_data$jbhrcpr)
final_data$pref_p1<-lead(pref_time,1)

final_data$pref_p1_new<-ifelse(final_data$wave==15 &final_data$esbrd==1 & final_data$jbhrcpr>=0 & final_data$pref_p1>=1, final_data$pref_p1, NA)
final_data$jbhrcpr_new<-ifelse(final_data$wave==15 & final_data$esbrd==1 & final_data$jbhrcpr>=0 & final_data$pref_p1>=1,final_data$jbhrcpr, NA) 

final_data$pref_few_p1few<-ifelse(final_data$jbhrcpr==1 & final_data$pref_p1_new==1, 1, ifelse(final_data$jbhrcpr==1 & final_data$pref_p1_new>=1,0,NA))
final_data$pref_few_p1same<-ifelse(final_data$jbhrcpr==1 & final_data$pref_p1_new==2, 1, ifelse(final_data$jbhrcpr==1 & final_data$pref_p1_new>=1,0,NA))
final_data$pref_few_p1more<-ifelse(final_data$jbhrcpr==1 & final_data$pref_p1_new==3, 1, ifelse(final_data$jbhrcpr==1 & final_data$pref_p1_new>=1,0,NA))

final_data$pref_same_p1few<-ifelse(final_data$jbhrcpr==2 & final_data$pref_p1_new==1, 1, ifelse(final_data$jbhrcpr==2 & final_data$pref_p1_new>=1,0,NA))
final_data$pref_same_p1same<-ifelse(final_data$jbhrcpr==2 & final_data$pref_p1_new==2, 1, ifelse(final_data$jbhrcpr==2 & final_data$pref_p1_new>=1,0,NA))
final_data$pref_same_p1more<-ifelse(final_data$jbhrcpr==2 & final_data$pref_p1_new==3, 1, ifelse(final_data$jbhrcpr==2 & final_data$pref_p1_new>=1,0,NA))

final_data$pref_more_p1few<-ifelse(final_data$jbhrcpr==3 & final_data$pref_p1_new==1, 1, ifelse(final_data$jbhrcpr==3 & final_data$pref_p1_new>=1,0,NA))
final_data$pref_more_p1same<-ifelse(final_data$jbhrcpr==3 & final_data$pref_p1_new==2, 1, ifelse(final_data$jbhrcpr==3 & final_data$pref_p1_new>=1,0,NA))
final_data$pref_more_p1more<-ifelse(final_data$jbhrcpr==3 & final_data$pref_p1_new==3, 1, ifelse(final_data$jbhrcpr==3 & final_data$pref_p1_new>=1,0,NA))

hildadesign_op<-svydesign(id=~xhhraid, strata=~xhhstrat, weights=~wlrop, nest=TRUE, data=final_data)

# row proportions
tb1<-svytable(~final_data$jbhrcpr_new+final_data$pref_p1_new, hildadesign_op) 
rps<-as.matrix(tb1/rowSums(tb1)) 
rps

# row proportions (with standard errors) - there may be a more elegant way to do this but this does the job
## where prefers fewer in current wave
svymean(~pref_few_p1few, hildadesign_op, na.rm=TRUE)
svymean(~pref_few_p1same, hildadesign_op, na.rm=TRUE )
svymean(~pref_few_p1more, hildadesign_op, na.rm=TRUE )
confint(svymean(~pref_few_p1few, hildadesign_op, na.rm=TRUE))
confint(svymean(~pref_few_p1same, hildadesign_op, na.rm=TRUE))
confint(svymean(~pref_few_p1more, hildadesign_op, na.rm=TRUE))

## where prefers same in current wave
svymean(~pref_same_p1few, hildadesign_op, na.rm=TRUE)
svymean(~pref_same_p1same, hildadesign_op, na.rm=TRUE )
svymean(~pref_same_p1more, hildadesign_op, na.rm=TRUE )
confint(svymean(~pref_same_p1few, hildadesign_op, na.rm=TRUE))
confint(svymean(~pref_same_p1same, hildadesign_op, na.rm=TRUE))
confint(svymean(~pref_same_p1more, hildadesign_op, na.rm=TRUE))

## where prefers more in current wave
svymean(~pref_more_p1few, hildadesign_op, na.rm=TRUE)
svymean(~pref_more_p1same, hildadesign_op, na.rm=TRUE )
svymean(~pref_more_p1more, hildadesign_op, na.rm=TRUE )
confint(svymean(~pref_more_p1few, hildadesign_op, na.rm=TRUE))
confint(svymean(~pref_more_p1same, hildadesign_op, na.rm=TRUE))
confint(svymean(~pref_more_p1more, hildadesign_op, na.rm=TRUE))

# proportions (with standard errors) with cell proportions (not row proportions as was done above)
a<-svymean(~interaction(jbhrcpr_new,pref_p1_new), design=hildadesign_op, na.rm=TRUE)
b<-ftable(a, rownames=list(stype=c("Fewer","Same","More"),comp.imp=c("Fewer","Same","More")))
b