library(readxl)
library(drc)
## Loading required package: MASS
## 
## 'drc' has been loaded.
## Please cite R and 'drc' if used for a publication,
## for references type 'citation()' and 'citation('drc')'.
## 
## Attaching package: 'drc'
## The following objects are masked from 'package:stats':
## 
##     gaussian, getInitial
library(magrittr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(broom)
library(ggplot2)
datosT <- read_excel("C:/Users/Heymdall/Documents/Paula/Echinochloa/Ensayo 6.1/DRC 6.1.xlsx")
datosF <- read_excel("C:/Users/Heymdall/Documents/Paula/Echinochloa/Ensayo 6.1/DRC 6.1.xlsx", 
    sheet = "Hoja3", range = "A1:F70")
Datos <-datosT %>% 
  mutate_at(vars(ID, Dosis), as.factor) %>%
  as.data.frame()
head(Datos)
##   ID Dosis Altura (cm) % Daño % Sobrevivencia Peso seco (g) Sobrevivencia
## 1 UM     0    13.16667      0            66.6         0.086            NA
## 2 UM     0    18.83333      0           100.0         0.174            NA
## 3 UM     0    10.66667      0           100.0         0.055            NA
## 4 UM     0    15.00000      0           100.0         0.106            NA
## 5 UM   140    15.83333     15           100.0         0.124            NA
## 6 UM   140    22.16667     10           100.0         0.218            NA
DatosDRC <-Datos %>% 
  as.data.frame()

Altura

## Altura MA
p1 <- Datos %>% 
  filter(ID == "MA") %>%
  ggplot(aes(x = Dosis, y = `Altura (cm)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 50) +
  labs(title = "Altura MA",
       x = "Dosis (gr i.a)",
       y = "Altura (cm)")

## Altura UM
p2 <- Datos %>% 
  filter(ID == "UM") %>%
  ggplot(aes(x = Dosis, y = `Altura (cm)`)) + 
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 50) +
  labs(title = "Altura UM",
       x = "Dosis (gr i.a)",
       y = "Altura (cm)")

## Altura LP
p3 <- Datos %>% 
  filter(ID == "LP") %>%
  ggplot(aes(x = Dosis, y = `Altura (cm)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 50) +
  labs(title = "Altura LP",
       x = "Dosis (gr i.a)",
       y = "Altura (cm)")

## Altura VE
p4 <- Datos %>% 
  filter(ID == "VE") %>%
  ggplot(aes(x = Dosis, y = `Altura (cm)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 50) +
  labs(title = "Altura VE",
       x = "Dosis (gr i.a)",
       y = "Altura (cm)")

## Altura AR
p5 <- Datos %>% 
  filter(ID == "AR") %>%
  ggplot(aes(x = Dosis, y = `Altura (cm)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 50) +
  labs(title = "Altura AR",
       x = "Dosis (gr i.a)",
       y = "Altura (cm)")

grid.arrange(p1, p2,p3, p4, p5, ncol = 2)

fct = LL.4(names=c(“Hill slope”,“Min”,“Max”,“EC50”)))

modelAlt <- drm(data=datosF, 
           formula = `Altura (cm)`~Dosis, 
           curveid = ID,
           fct = LL.4())
tidy(modelAlt)
## Warning in sqrt(diag(varMat)): Se han producido NaNs
## # A tibble: 12 × 6
##    term  curve    estimate std.error statistic     p.value
##    <chr> <chr>       <dbl>     <dbl>     <dbl>       <dbl>
##  1 b     UM       0.853     NaN       NaN      NaN        
##  2 b     LP       2.98        2.68      1.11     2.71e-  1
##  3 b     AR      -0.000186  NaN       NaN      NaN        
##  4 c     UM      13.2         0.698    18.9      1.42e- 26
##  5 c     LP      20.5         1.57     13.1      4.29e- 19
##  6 c     AR      20.9         1.35     15.5      2.00e- 22
##  7 d     UM      14.8         1.29     11.5      2.06e- 16
##  8 d     LP      27.4         1.34     20.4      3.21e- 28
##  9 d     AR      31.7         0.0337  943.       1.69e-121
## 10 e     UM     263.        NaN       NaN      NaN        
## 11 e     LP    1663.        904.        1.84     7.09e-  2
## 12 e     AR      36.6      3479.        0.0105   9.92e-  1
plot(modelAlt, legendPos = c(12000, 40))

ED (modelAlt, respLev = 50, interval = "delta")
## Warning in sqrt(dEDval %*% varCov %*% dEDval): Se han producido NaNs

## Warning in sqrt(dEDval %*% varCov %*% dEDval): Se han producido NaNs
## 
## Estimated effective doses
## 
##          Estimate Std. Error     Lower     Upper
## e:AR:50    36.615   3479.140 -6930.243  7003.472
## e:LP:50  1663.355    903.894  -146.662  3473.372
## e:UM:50   263.206        NaN       NaN       NaN
compParm(modelAlt,"e","-")
## 
## Comparison of parameter 'e' 
## 
##       Estimate Std. Error t-value p-value  
## UM-LP -1400.15     694.81 -2.0152 0.04861 *
## UM-AR   226.59    3430.76  0.0660 0.94757  
## LP-AR  1626.74    3594.64  0.4525 0.65259  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelFit(modelAlt)
## Lack-of-fit test
## 
##           ModelDf    RSS Df F value p value
## ANOVA          46 591.33                   
## DRC model      57 830.17 11  1.6890  0.1062

Daño (%)

p1 <- Datos %>% 
  filter(ID == "MA") %>%
  ggplot(aes(x = Dosis, y = `% Daño`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Daño MA",
       x = "Dosis (gr i.a)",
       y = "% Daño")

## Altura UM
p2 <- Datos %>% 
  filter(ID == "UM") %>%
  ggplot(aes(x = Dosis, y = `% Daño`)) + 
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Daño UM",
       x = "Dosis (gr i.a)",
       y = "% Daño")

## Altura LP
p3 <- Datos %>% 
  filter(ID == "LP") %>%
  ggplot(aes(x = Dosis, y = `% Daño`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Daño LP",
       x = "Dosis (gr i.a)",
       y = "% Daño")

## Altura VE
p4 <- Datos %>% 
  filter(ID == "VE") %>%
  ggplot(aes(x = Dosis, y = `% Daño`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Daño VE",
       x = "Dosis (gr i.a)",
       y = "% Daño")

## Altura AR
p5 <- Datos %>% 
  filter(ID == "AR") %>%
  ggplot(aes(x = Dosis, y = `% Daño`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Daño AR",
       x = "Dosis (gr i.a)",
       y = "% Daño")

grid.arrange(p1, p2,p3, p4, p5, ncol = 2)

modelDano <- drm(data=datosF, 
           formula = `% Daño`~Dosis, 
           curveid = ID,
           fct = LL.4())
tidy(modelDano)
## # A tibble: 12 × 6
##    term  curve  estimate std.error statistic  p.value
##    <chr> <chr>     <dbl>     <dbl>     <dbl>    <dbl>
##  1 b     UM    -8.81e- 1     0.179 -4.91e+ 0 7.95e- 6
##  2 b     LP    -2.12e+ 0     1.04  -2.04e+ 0 4.61e- 2
##  3 b     AR    -6.98e- 3     0.110 -6.32e- 2 9.50e- 1
##  4 c     UM    -1.19e+ 0     5.59  -2.14e- 1 8.32e- 1
##  5 c     LP    -3.84e- 1     4.60  -8.36e- 2 9.34e- 1
##  6 c     AR    -2.34e- 4     6.35  -3.69e- 5 1.00e+ 0
##  7 d     UM     1.15e+ 2     6.51   1.77e+ 1 4.14e-25
##  8 d     LP     5.77e+ 1     5.45   1.06e+ 1 4.44e-15
##  9 d     AR     7.31e+ 0     9.96   7.34e- 1 4.66e- 1
## 10 e     UM     5.28e+ 2     9.98   5.29e+ 1 1.99e-50
## 11 e     LP     2.35e+ 3    10.0    2.35e+ 2 4.10e-87
## 12 e     AR     2.56e+17    10      2.56e+16 0
plot(modelDano, legendPos = c(50, 100))

ED (modelDano, respLev = 50, interval = "delta")
## 
## Estimated effective doses
## 
##           Estimate Std. Error      Lower      Upper
## e:AR:50 2.5603e+17 1.0000e+01 2.5603e+17 2.5603e+17
## e:LP:50 2.3484e+03 9.9987e+00 2.3284e+03 2.3685e+03
## e:UM:50 5.2765e+02 9.9813e+00 5.0766e+02 5.4764e+02
compParm(modelDano,"e","-")
## 
## Comparison of parameter 'e' 
## 
##          Estimate  Std. Error     t-value   p-value    
## UM-LP -1.8208e+03  1.4128e+01 -1.2888e+02 < 2.2e-16 ***
## UM-AR -2.5603e+17  1.4129e+01 -1.8121e+16 < 2.2e-16 ***
## LP-AR -2.5603e+17  1.4141e+01 -1.8105e+16 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelFit(modelDano)
## Lack-of-fit test
## 
##           ModelDf   RSS Df F value p value
## ANOVA          46  9244                   
## DRC model      57 11360 11  0.9574  0.4970

Sobrevivencia (%)

p1 <- Datos %>% 
  filter(ID == "MA") %>%
  ggplot(aes(x = Dosis, y = `% Sobrevivencia`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Sobrevivencia MA",
       x = "Dosis (gr i.a)",
       y = "% Sobrevivencia")

## Altura UM
p2 <- Datos %>% 
  filter(ID == "UM") %>%
  ggplot(aes(x = Dosis, y = `% Sobrevivencia`)) + 
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Sobrevivencia UM",
       x = "Dosis (gr i.a)",
       y = "% Sobrevivencia")

## Altura LP
p3 <- Datos %>% 
  filter(ID == "LP") %>%
  ggplot(aes(x = Dosis, y = `% Sobrevivencia`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Sobrevivencia LP",
       x = "Dosis (gr i.a)",
       y = "% Sobrevivencia")

## Altura VE
p4 <- Datos %>% 
  filter(ID == "VE") %>%
  ggplot(aes(x = Dosis, y = `% Sobrevivencia`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Sobrevivencia VE",
       x = "Dosis (gr i.a)",
       y = "% Sobrevivencia")

## Altura AR
p5 <- Datos %>% 
  filter(ID == "AR") %>%
  ggplot(aes(x = Dosis, y = `% Sobrevivencia`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 100) +
  labs(title = "% Sobrevivencia AR",
       x = "Dosis (gr i.a)",
       y = "% Sobrevivencia")

grid.arrange(p1, p2,p3, p4, p5, ncol = 2)

modelSobr <- drm(data=datosF[-(47:70),], 
           formula = `% Sobrevivencia`~Dosis, 
           curveid = ID,
           fct = LL.4())
## Control measurements detected for level: AR
tidy(modelSobr)
## # A tibble: 8 × 6
##   term  curve  estimate std.error statistic  p.value
##   <chr> <chr>     <dbl>     <dbl>     <dbl>    <dbl>
## 1 b     UM        2.01       1.43    1.40   1.68e- 1
## 2 b     LP        1.79       1.71    1.04   3.03e- 1
## 3 c     UM       -0.294     13.8    -0.0213 9.83e- 1
## 4 c     LP     -160.      1288.     -0.125  9.01e- 1
## 5 d     UM       90.5       10.7     8.42   3.22e-10
## 6 d     LP       99.6        5.63   17.7    3.00e-20
## 7 e     UM      797.       217.      3.67   7.34e- 4
## 8 e     LP    43957.    170333.      0.258  7.98e- 1
plot(modelSobr, legendPos = c(50, 40))

ED (modelSobr, respLev = 50, interval = "delta")
## 
## Estimated effective doses
## 
##           Estimate Std. Error      Lower      Upper
## e:LP:50   43956.82  170333.20 -300864.71  388778.35
## e:UM:50     796.74     216.89     357.67    1235.82
compParm(modelSobr,"e","-")
## 
## Comparison of parameter 'e' 
## 
##       Estimate Std. Error t-value p-value
## UM-LP   -43160     170333 -0.2534  0.8013
modelFit(modelSobr)
## Lack-of-fit test
## 
##           ModelDf   RSS Df F value p value
## ANOVA          30 10376                   
## DRC model      38 14084  8  1.3399  0.2626

Peso seco (g)

p1 <- Datos %>% 
  filter(ID == "MA") %>%
  ggplot(aes(x = Dosis, y = `Peso seco (g)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 1.5) +
  labs(title = "Peso seco (g) MA",
       x = "Dosis (gr i.a)",
       y = "Peso seco (g)")

## Altura UM
p2 <- Datos %>% 
  filter(ID == "UM") %>%
  ggplot(aes(x = Dosis, y = `Peso seco (g)`)) + 
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 1.5) +
  labs(title = "Peso seco (g) UM",
       x = "Dosis (gr i.a)",
       y = "Peso seco (g)")

## Altura LP
p3 <- Datos %>% 
  filter(ID == "LP") %>%
  ggplot(aes(x = Dosis, y = `Peso seco (g)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 1.5) +
  labs(title = "Peso seco (g) LP",
       x = "Dosis (gr i.a)",
       y = "Peso seco (g)")

## Altura VE
p4 <- Datos %>% 
  filter(ID == "VE") %>%
  ggplot(aes(x = Dosis, y = `Peso seco (g)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 1.5) +
  labs(title = "Peso seco (g) VE",
       x = "Dosis (gr i.a)",
       y = "Peso seco (g)")

## Altura AR
p5 <- Datos %>% 
  filter(ID == "AR") %>%
  ggplot(aes(x = Dosis, y = `Peso seco (g)`)) +
  geom_boxplot() +
  scale_fill_manual(values = "black") + ylim(0, 1.5) +
  labs(title = "Peso seco (g) AR",
       x = "Dosis (gr i.a)",
       y = "Peso seco (g)")

grid.arrange(p1, p2,p3, p4, p5, ncol = 2)

modelpeso <- drm(data=datosF, 
           formula = `Peso seco (g)`~Dosis, 
           curveid = ID,
           fct = LL.4())
tidy(modelpeso)
## # A tibble: 12 × 6
##    term  curve   estimate   std.error statistic  p.value
##    <chr> <chr>      <dbl>       <dbl>     <dbl>    <dbl>
##  1 b     UM       21.0       226.        0.0932 9.26e- 1
##  2 b     LP        0.942       0.362     2.60   1.18e- 2
##  3 b     AR      -14.7        42.5      -0.347  7.30e- 1
##  4 c     UM        0.0422      0.0160    2.63   1.09e- 2
##  5 c     LP       -0.466       1.86     -0.250  8.03e- 1
##  6 c     AR        0.227       0.0113   20.0    8.76e-28
##  7 d     UM        0.0943      0.0138    6.82   6.30e- 9
##  8 d     LP        0.286       0.0178   16.0    4.17e-23
##  9 d     AR        0.288       0.0199   14.5    4.88e-21
## 10 e     UM      800.       2976.        0.269  7.89e- 1
## 11 e     LP    67613.     241249.        0.280  7.80e- 1
## 12 e     AR     6345.       6376.        0.995  3.24e- 1
plot(modelpeso, legendPos = c(100, 0.22))

ED (modelpeso, respLev = 50, interval = "delta")
## 
## Estimated effective doses
## 
##           Estimate Std. Error      Lower      Upper
## e:AR:50    6345.12    6375.91   -6422.42   19112.66
## e:LP:50   67612.74  241248.74 -415479.52  550705.01
## e:UM:50     799.56    2976.49   -5160.77    6759.88
compParm(modelpeso,"e","-")
## 
## Comparison of parameter 'e' 
## 
##       Estimate Std. Error t-value p-value
## UM-LP -66813.2   241267.1 -0.2769  0.7828
## UM-AR  -5545.6     7036.5 -0.7881  0.4339
## LP-AR  61267.6   241333.0  0.2539  0.8005
modelFit(modelpeso)
## Lack-of-fit test
## 
##           ModelDf      RSS Df F value p value
## ANOVA          46 0.062779                   
## DRC model      57 0.130396 11  4.5040  0.0001