Dose Response

April 27, 2018
Dose Response Data Seminar

Load packages

# .libPaths("P:/RLibrary2")
knitr::opts_chunk$set(echo = TRUE)
# install.packages("drc")

library(drc)
library(tidyverse)

Load data

id <- "1UUxeLSz0IAb9LsTmjIRYQ3kSbqb-JWpE"
dose<- read_csv(sprintf("https://docs.google.com/uc?id=%s&export=download", id))
## Parsed with column specification:
## cols(
##   Trt = col_character(),
##   Dose = col_integer(),
##   Rep = col_integer(),
##   control = col_integer(),
##   RUN = col_integer()
## )
head(dose)
## # A tibble: 6 x 5
##   Trt    Dose   Rep control   RUN
##   <chr> <int> <int>   <int> <int>
## 1 0         0     1       0     1
## 2 0         0     2       0     1
## 3 0         0     3       0     1
## 4 0         0     4       0     1
## 5 0.25x   380     1       0     1
## 6 0.25x   380     2       0     1
glimpse(dose)
## Observations: 72
## Variables: 5
## $ Trt     <chr> "0", "0", "0", "0", "0.25x", "0.25x", "0.25x", "0.25x"...
## $ Dose    <int> 0, 0, 0, 0, 380, 380, 380, 380, 660, 660, 660, 660, 12...
## $ Rep     <int> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, ...
## $ control <int> 0, 0, 0, 0, 0, 0, 0, 20, 35, 20, 25, 30, 0, 20, 25, 5,...
## $ RUN     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
ggplot(data = dose) +
  geom_point(aes(x = Dose, y = control, shape = as.factor(RUN), color = as.factor(RUN))) + 
  labs(x = expression(paste("Herbicide (g ha"^"-1",")")),
       y = expression("Visual control (%)"),
       color = "Run",
       shape = "Run") + 
  theme_classic() + 
  theme(legend.position = c(0.9, 0.2),
        legend.background = element_rect(color = "black"))

log_model <- drm(control ~ Dose, curveid = RUN, data = dose, fct = LL.3(names = c("Slope", "Upper Limit", "ED50")))

summary(log_model)
## 
## Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms)
## 
## Parameter estimates:
## 
##                 Estimate Std. Error t-value    p-value    
## Slope:1         -0.85279    0.31428 -2.7135   0.008484 ** 
## Slope:2         -0.73947    0.33005 -2.2405   0.028431 *  
## Upper Limit:1   77.91164   18.45579  4.2215 0.00007583 ***
## Upper Limit:2   78.56784   18.16784  4.3246 0.00005282 ***
## ED50:1        4662.71009 3329.60813  1.4004   0.166084    
## ED50:2        1869.95810 1545.70856  1.2098   0.230681    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  14.9388 (66 degrees of freedom)
modelFit(log_model)
## Lack-of-fit test
## 
##           ModelDf   RSS Df F value p value
## ANOVA          54 10538                   
## DRC model      66 14729 12  1.7900  0.0734
## Comparison of model parameters
compParm(log_model, "Slope", "-")
## 
## Comparison of parameter 'Slope' 
## 
##     Estimate Std. Error t-value p-value
## 1-2 -0.11332    0.45574 -0.2486  0.8044
compParm(log_model, "Upper Limit", "-")
## 
## Comparison of parameter 'Upper Limit' 
## 
##     Estimate Std. Error t-value p-value
## 1-2 -0.65619   25.89762 -0.0253  0.9799
compParm(log_model, "ED50", "-")
## 
## Comparison of parameter 'ED50' 
## 
##     Estimate Std. Error t-value p-value
## 1-2   2792.8     3670.9  0.7608  0.4495
plot(log_model, type = "all")

vsmall <- .Machine$double.eps

newdata <- expand.grid(dose = seq(vsmall, max(dose$Dose), by =  10),
                       RUN = unique(dose$RUN))

pred_data <- cbind(newdata, predict(log_model, newdata = newdata, interval = "confidence"))

head(pred_data)
##           dose RUN   Prediction         Lower        Upper
## 1 2.220446e-16   1 2.593469e-15 -6.646001e-14 7.164695e-14
## 2 1.000000e+01   1 4.107018e-01 -9.423271e-01 1.763731e+00
## 3 2.000000e+01   1 7.385838e-01 -1.371635e+00 2.848803e+00
## 4 3.000000e+01   1 1.039611e+00 -1.664055e+00 3.743277e+00
## 5 4.000000e+01   1 1.323758e+00 -1.877710e+00 4.525226e+00
## 6 5.000000e+01   1 1.595542e+00 -2.037673e+00 5.228758e+00
ggplot(data = pred_data) +
  geom_line(aes(x = dose, y = Prediction, colour = as.factor(RUN))) + 
  geom_ribbon(aes(x = dose, ymax = Upper, ymin = Lower, fill = as.factor(RUN)), alpha = 0.2) + 
  geom_point(data = dose, aes(x = Dose, y = control, colour = as.factor(RUN), shape = as.factor(RUN))) +
  theme_classic() + 
  scale_x_log10() +
  theme(legend.position = "bottom")

Time Series Analysis

April 27, 2018
Time Series Data Seminar

Spatial Autocorrelation

April 18, 2018
Spatial autocorrelation Data Seminar

Partial Least Squares Regression

April 13, 2018
Refelectance Partial Least Squares Regression Data Seminar
comments powered by Disqus