Skip to content

Commit

Permalink
adding new RSCABS, MQJT related functions
Browse files Browse the repository at this point in the history
  • Loading branch information
Zhenglei-BCS committed Mar 1, 2025
1 parent f056762 commit 60dd518
Show file tree
Hide file tree
Showing 8 changed files with 169 additions and 214 deletions.
90 changes: 38 additions & 52 deletions README.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ Some of the functions are adapted from archived packages or single functions of

This package is open source, and any contributions or improvements, especially on the documentation side, are welcome.

*Please note that the documentation website for this package is currently under development. Some articles are still placeholders, and many more are on the way. However, the ongoing development of the website does not impact the usage of this R package. *

## Installation

You can install the development version of drcHelper from [GitHub](https://github.com/) with:
Expand All @@ -51,85 +53,69 @@ pak::pak("Bayer-Group/drcHelper")

## Data Overview

```{r message=FALSE,warning=FALSE}
library(drcHelper)
library(drc)
library(dplyr)
library(purrr)
library(ggplot2)
theme_set(theme_bw())
sum1 <- oecd201 %>% group_by(Time,Treatment) %>% summarise(Yield_mean=mean(Yield),Yield_sd=sd(Yield),GrowthRate_mean=mean(GrowthRate),GrowthRate_sd=sd(GrowthRate))
sum0 <- sum1%>%filter(Treatment=="Control")%>%rename(Yield0=Yield_mean,GrowthRate0=GrowthRate_mean)%>%dplyr::select(c(Time,Yield0,GrowthRate0))
# sum0
sumtab <- left_join(sum1%>%filter(Time>0),sum0) %>% mutate(Yield_Inhibition=(Yield0-Yield_mean)/Yield0*100,GrowthRate_Inhibition=(GrowthRate0-GrowthRate_mean)/GrowthRate0*100) %>% dplyr::select(c(Time,Treatment,Yield_mean,Yield_sd,Yield_Inhibition,GrowthRate_mean,GrowthRate_sd,GrowthRate_Inhibition))
```
## Preliminary Summary


```{r}
sumtab%>%dplyr::select(c(Yield_mean,Yield_sd,Yield_Inhibition))%>%filter(Time==72)%>%knitr::kable(.,digits = 2,caption="<center><strong>Yield Summary at Time 72h<strong><center>",escape = FALSE)##%>% kableExtra::kable_styling(bootstrap_options = "striped")##%>%kableExtra::kable_classic_2()
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
warning = FALSE,
message = FALSE
)
library(tidyverse)
```


```{r}
sumtab%>%dplyr::select(c(GrowthRate_mean,GrowthRate_sd,GrowthRate_Inhibition))%>%filter(Time==72)%>%knitr::kable(.,digits = 2,caption="<center><strong>Growth Rate Summary at Time 72h<strong><center>",escape = FALSE)##%>%kableExtra::kable_classic()
```{r setup}
library(drcHelper)
```

## Model Fitting and Comparison For Yield
```{r}
data("dat_medium")
dat_medium <- dat_medium %>% mutate(Treatment=factor(Dose,levels=unique(Dose)))
dat_medium$Response[dat_medium$Response < 0] <- 0
prelimPlot3(dat_medium)
prelimSummary(dat_medium) %>% knitr::kable(.,digits = 3)
```

```{r warning=FALSE,message=FALSE}
datTn<- subset(oecd201,Time==72)
## Fitting multiple models and rank them.

mod <- drm(Yield~Concentration,data=datTn,fct=LL.3())
fctList <- list(LL2.3(),W2.3(),W1.3(),EXD.3(),EXD.2(),LN.3(),W2.4(),LL.4(),LL2.4())
plot(mod,type="all")
```{r}
mod <- drm(Response~Dose,data=dat_medium,fct=LL.3())
fctList <- list(LN.4(),LL.4(),W1.3(),LL2.2())
# plot(mod,type="all")
res <- mselect.plus(mod,fctList = fctList )
modList <- res$modList
edResTab <- mselect.ED(modList = modList,respLev = c(10,20,50),trend=datTn$Trend_Yield[1])
plot.edList(edResTab)
resComp <- drcCompare(modRes = res,trend="Decrease")
```
Note that by default settings, the fitted models did not converge except for the LL.3 model.
res$Comparison
```{r}
knitr::kable(edResTab[1:3,],caption = "14 day TSL Yield",digits = 3)
drcCompare(modRes=res)
```


```{r}
knitr::kable(resComp[1,],caption = "14 day TSL Yield, Model Comparison",digits = 3)
library(purrr)
edResTab <- mselect.ED(modList = modList,respLev = c(10,20,50),trend="Decrease",CI="inv")
edResTab
```

```{r fig.cap="Yield Model Fits"}
plot.modList(modList,scale="logx")
```

```{r}
plot.modList(modList[c(1,2,3,4)],scale="logx",npts=40)
p <-plot.modList(modList[c(1)],scale="logx",npts=80)+theme(legend.position = "none")+ggtitle("14 day Total Shoot Length, \n3-parameter type II Weibull Model Fit")
addECxCI(p=p,object=modList[[1]],EDres=NULL,trend="Decrease",endpoint="EC", respLev=c(10,20,50),
textAjust.x=0.01,textAjust.y=1,useObsCtr=FALSE,d0=NULL,textsize = 4,lineheight = 1,xmin=0.012)+ylab("Total Shoot Length [cm]") + xlab("Concentration [µg a.s./L]")
## ggsave("TSL_14d_Yield.png")
```
## Plot multiple models together

```{r}
resED <- t(edResTab[1:3, c(2,4,5,6)])
colnames(resED) <- paste("EC", c(10,20,50))
knitr::kable(resED,caption = "Total Shoot Length Growth Yield at 14 day",digits = 3)
p <- plot.modList(modList[1:3])
p
```

```{r}

mod <-modList[[1]]
edres <- ED.plus(mod,c(5,10,20,50),trend="Decrease")
pander::pander(as.data.frame(edres))
## Adding ECx and ECx CI's to the plots

```{r}
## addECxCI(p)
```


```{r}
modsum <- summary(mod)
pander::pander(coef(modsum))
```


## ToDo
Expand Down
Loading

0 comments on commit 60dd518

Please sign in to comment.