drc_test.R

# Dose response curve
# A.Bailey 12th June 2018



library(tidyverse)
library(drc)
library(janitor)

# Read data from thrashing analysis - at 60 min time point 
      dat <- read_csv("Learning/N2_DR_thrashing.csv") %>% clean_names()
      dat1 <- read_csv("Learning/bus17_DR_thrashing.csv") %>% 
        clean_names()

# Drop dose and convert to matrix
dat_m <- dat %>% dplyr::select(-dose) %>% as.matrix()
dat_m1 <- dat1 %>%  dplyr::select(-dose) %>%  as.matrix()
# Find NAs
k <- which(is.na(dat_m), arr.ind=TRUE)
k1 <- which(is.na(dat_m1), arr.ind=TRUE)
# Copied the matrix
m <- dat_m
m1 <- dat_m1
# Fill missing values with row mean, mean for that dose
m[k] <- rowMeans(dat_m, na.rm=TRUE)[k[,1]]
m1[k1] <- rowMeans(dat_m, na.rm=TRUE)[k1[,1]]
# Convert back to tibble
d_m <- as.tibble(m)
d_m1 <- as.tibble(m1)
# Re-bind the dose column to complete data
dat_clean <- bind_cols(dose = dat$dose,d_m)

dat_clean1 <- bind_cols(dose = dat$dose,d_m1)
#combine data - look into this here 
dat_comb <- bind_rows(dat_clean, dat_clean1)
# Gather according to dose
dat_melt <- dat_clean %>% gather(key = vars,value = response, -dose) #%>% dplyr::select(-vars)

dat_melt2 <- dat %>% gather(key = vars,value = response,-dose) %>% dplyr::select(-vars)
# Mutate log transformation
dat_melt <- dat_melt %>% mutate(log_dose = log(dose))
# Fit model with imputed data
model <- drc::drm(data = dat_melt, response ~ dose, fct = L.3())
# Fit original data
model2 <- drc::drm(data = dat_melt2, response ~ dose, fct = L.3(), na.action = na.omit)

summary(model2)
plot(model2)
plot(model)

preds <- as.tibble(model$predres) %>% clean_names()

dat_all <- dat_melt %>% mutate(predict = preds$predicted_values)

dat_spread <- dat_all %>% 
  dplyr::select(-log_dose) %>% 
  group_by(dose) %>% 
  mutate(mean_response = mean(response),
         std_dev = sd(response),
         se = std_dev/sqrt(length(response))) %>% 
  dplyr::select(dose,mean_response,vars,se,predict) %>% 
  spread(vars,mean_response) %>% 
  dplyr::select(dose, se, mean_response = x10, predict) %>% 
  ungroup() %>% 
  mutate(norm_rep = (mean_response - min(mean_response)) / 
           (max(mean_response) - min(mean_response)),
         norm_pred = (predict - min(predict)) / (max(predict) - min(predict)),
         norm_se = se/(mean_response))


dat_spread %>% 
  ggplot(aes(dose,norm_rep)) + 
  geom_point() +
  geom_errorbar(aes(ymin = norm_rep - norm_se,ymax = norm_rep + norm_se)) +
  geom_line(aes(dose,norm_pred)) +
  #geom_smooth(aes(dose,norm_pred)) +
  scale_x_log10()
dat_sprea