test.Rmd 2.48 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
---
title: "test"
output: html_document
---

```{r echo=FALSE}

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

# Read data from thrashing analysis - at 60 min time point 
          dat <- read_csv("N2_DR_thrashing.csv") %>% clean_names()
      dat1 <- read_csv("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_spread
```