inst/method_compair_demo.R

In RussellXu/rtestimate: Estimation Of Reproductive Number

library(xlsx)

options(digits = 2)

#hk_incidenc <- read.csv('./data/hongkong_cases.csv')
# hk_rt <- read.xlsx('./output/True_rt_r2.xlsx', sheetName = 'Sheet1')
# hk_epiestim <- read.xlsx('./output/epiestim_hk2x.xlsx', sheetName = 'Sheet1')
# hk_viral_load <- read.xlsx('./output/viral_load_r2output.xlsx', sheetName = 'Sheet1')
#
# colnames(hk_rt) <- c('date', 'rt')
# colnames(hk_epiestim) <- c('lower_epi95', 'upper_epi95', 'rt_epi95', 'date')
# colnames(hk_viral_load) <- c('time', 'lower_ct95', 'ct_rt','upper_ct95')

# calculate the goodness of fit data
rt2_table <- readxl::read_excel('~/Desktop/prism_use/modified_rt2/new_modified.xlsx', sheet = 1)
rt2_table <-as.data.frame(lapply(rt2_table, as.numeric))
colnames(rt2_table) <- c('time','epiestim','epiestim_u','epiestim_l','viral_load','viral_load_u','viral_load_l','true_rt','true_rt_u','true_rt_l','incidence','incidence_u','incidence_l', 'modified', 'modified_u', 'modified_l')
# get the peak time of incidence value
peak_time <- rt2_table$time[which.max(rt2_table$incidence)]

# before incidence peak
caret::postResample(pred = rt2_table[1:peak_time,]$epiestim, obs = rt2_table[1:peak_time,]$true_rt)
cor(rt2_table[1:peak_time,]$epiestim, rt2_table[1:peak_time,]$true_rt, method = 'pearson')
caret::postResample(pred = rt2_table[1:peak_time,]$viral_load, obs = rt2_table[1:peak_time,]$true_rt)
cor(rt2_table[1:peak_time,]$viral_load, rt2_table[1:peak_time,]$true_rt, method = 'pearson')

# after incidence peak
caret::postResample(pred = rt2_table[peak_time:nrow(rt2_table),]$epiestim, obs = rt2_table[peak_time:nrow(rt2_table),]$true_rt)
cor(rt2_table[peak_time:nrow(rt2_table),]$epiestim, rt2_table[peak_time:nrow(rt2_table),]$true_rt, method = 'pearson')
caret::postResample(pred = rt2_table[peak_time:nrow(rt2_table),]$viral_load, obs = rt2_table[peak_time:nrow(rt2_table),]$true_rt)
cor(rt2_table[peak_time:nrow(rt2_table),]$viral_load, rt2_table[peak_time:nrow(rt2_table),]$true_rt, method = 'pearson')

# Rt>1
caret::postResample(pred = rt2_table[which(rt2_table$true_rt >1),]$epiestim, obs = rt2_table[which(rt2_table$true_rt >1),]$true_rt)
cor(rt2_table[which(rt2_table$true_rt >1),]$epiestim, rt2_table[which(rt2_table$true_rt >1),]$true_rt, method = 'pearson')
caret::postResample(pred = rt2_table[which(rt2_table$true_rt >1),]$viral_load, obs = rt2_table[which(rt2_table$true_rt >1),]$true_rt)
cor(rt2_table[which(rt2_table$true_rt >1),]$viral_load, rt2_table[which(rt2_table$true_rt >1),]$true_rt, method = 'pearson')

# Rt<1
caret::postResample(pred = rt2_table[which(rt2_table$true_rt <1),]$epiestim, obs = rt2_table[which(rt2_table$true_rt <1),]$true_rt)
cor(rt2_table[which(rt2_table$true_rt <1),]$epiestim, rt2_table[which(rt2_table$true_rt <1),]$true_rt, method = 'pearson')
caret::postResample(pred = rt2_table[which(rt2_table$true_rt <1),]$viral_load, obs = rt2_table[which(rt2_table$true_rt <1),]$true_rt)
cor(rt2_table[which(rt2_table$true_rt <1),]$viral_load, rt2_table[which(rt2_table$true_rt <1),]$true_rt, method = 'pearson')

# the whole line
# epiestim method: RMSE & R squared
caret::postResample(pred = rt2_table$epiestim, obs = rt2_table$true_rt)

# pearson
cor(rt2_table$epiestim, rt2_table$true_rt, method = 'pearson')


# epiestim method: modified RMSE & R squared
caret::postResample(pred = rt2_table$modified, obs = rt2_table$true_rt)

# pearson
cor(rt2_table$modified, rt2_table$true_rt, method = 'pearson')

# viral load method: RMSE & R suqared
caret::postResample(pred = rt2_table$viral_load, obs = rt2_table$true_rt)

# pearson
cor(rt2_table$viral_load, rt2_table$true_rt, method = 'pearson')

# ct data figure to prism
# res <- data.frame()
# data <- read.csv('./data/ct_r2_data.csv')
# data2 <- data %>%
#   filter(ct < 40) %>%
#   group_by(t) %>%
#   mutate(row=row_number()) %>%
#   pivot_wider(names_from = t, values_from = ct) %>%
#   select(-row)
# write.csv(data2, './data/ct_rt2_wider.csv')


# calc every single line make a error bar
# R02 : [12:111]
# R05 : [12:35]
# R010: [12:35]
epiestim <- read.csv('~/Desktop/for_modified_data/R010/epiestim.csv')
epiestim_adjust <- read.csv('~/Desktop/for_modified_data/R010/epiestim_adjust.csv')
viral_load <- read.csv('~/Desktop/for_modified_data/R010/viral_load.csv')
ground_truth <- readxl::read_excel('~/Desktop/for_modified_data/R010/true_rt_R10.xlsx')
true_rt <- ground_true$true_rt


assess_rmse <- function(pred) {
  pred <- pred
  obs <- true_rt[12:46]
    rmse<- caret::postResample(pred = pred, obs = obs)[1]
    names(rmse) <- NULL
    return(rmse)
}

assess_rsquare <- function(pred) {
  pred <- pred
  obs <- true_rt[12:46]
  rsquare<- caret::postResample(pred = pred, obs = obs)[2]
  names(rsquare) <- NULL
  return(rsquare)
}
assess_pearson <- function(pred) {
  pred <- pred
  obs <- true_rt[12:46]
  pearson <- cor(pred, obs, method = 'pearson')
  return(pearson)
}


epi_rsquare <- apply(epiestim, 2, FUN = assess_rsquare)
epi_rmse  <- apply(epiestim, 2, FUN = assess_rmse)
epi_pearson <- apply(epiestim, 2, FUN = assess_pearson)
adj_rsquare <- apply(epiestim_adjust, 2, FUN = assess_rsquare)
adj_rmse  <- apply(epiestim_adjust, 2, FUN = assess_rmse)
adj_pearson <- apply(epiestim_adjust, 2, FUN = assess_pearson)
vir_rsquare <- apply(viral_load[12:46,], 2, FUN = assess_rsquare)
vir_rmse  <- apply(viral_load[12:46,], 2, FUN = assess_rmse)
vir_pearson <- apply(viral_load[12:46,], 2, FUN = assess_pearson)

a_r2 <- data.frame(epi_rsquare = epi_rsquare[-1], epi_rmse = epi_rmse[-1], epi_pearson = epi_pearson[-1],
                   adj_rsquare = adj_rsquare[-1], adj_rmse = adj_rmse[-1], adj_pearson = adj_pearson[-1],
                   vir_rsquare = vir_rsquare, vir_rmse = vir_rmse, vir_pearson = vir_pearson)
#write.csv(a_r2, '~/Desktop/for_modified_data/measure_r10.csv')