R scripts for modeling metal mixture toxicity 2020-01-01 Edited by: Assco Prof. Hao QIU@SJTU
###########################################################################
Type: R Package
Title: Mixture Toxicity Predictions and statistical analysis
Depends on: R (>= 3.0.0)
IA/CA Model according to the following papers:
J. Hazard. Mater. 2019, 121940；Environ. Pollut. 2019, 246: 114-121；Environ. Toxicol. Chem. 2011, 30: 2084-2093
###########################################################################

#Data analysis
library(nortest)
library(drc)
library(lattice)

###Data input###
Data<-­‐read.table("CarbCylMeasured.txt", header=TRUE, na.strings=NA,sep="\t",dec=".")
Datamean<-­‐read.table("CarbCylMeanMeasured.txt", header=FALSE, na.strings=NA,sep="\t",dec=".")

###INDEPENDENT ACTION MODEL=IA####
###CONCENTRATION ADDITION MODEL=CA###

#Control observation are excluded as they will result in NANs by the a-­‐product term 
#E.g. dividing by 0 as for ((Zinc/EI)+(Cadmium/EC))^(-­‐2)
# both Zinc and Cadmium are 0 in control terms, resulting in a final zero in denominator
# Although this is only for the deviation models, we still do is for standard model as well
# Otherwise the deviation and standard model will be  fitted to other datasets
# Then statistical comparison is not so straightforward.

#STANDARD MODEL
#Fit model only to single stressor data, no mixture data included
ModelIAsingle<-­‐nls(Plant~
mean(Data[1:19,3])*1/((1+(Zinc/EI)^BI)*(1+(Cadmium/EC)^BC)),
data=Data[20:49,],start=list(BI=1,BC=1,EI=4,EC=50),
trace=TRUE, na.action=na.omit)

ModelCAsingle<-­‐nls(Plant~
(mean(Data[1:19,3]))/(((Zinc*EC+Cadmium*EI)/(EI*EC))^B+1),
data=Data[20:49,], start=list(B=1,EI=4,EC=50),
trace=TRUE,na.action=na.omit)

summary(ModelIAsingle)
summary(ModelCAsingle)

#Fit standard model to all data including mixtures
ModelIA<-­‐nls(Plant~
mean(Data[1:19,3])*1/((1+(Zinc/EI)^BI)*(1+(Cadmium/EC)^BC)),
data=Data[20:76,], start=list(BI=1,BC=1,EI=4,EC=50),
trace=TRUE, na.action=na.omit)

ModelCA<-­‐nls(Plant~(mean(Data[1:19,3]))/(((Zinc*EC+Cadmium*EI)/(EI*EC))^B+1),
data=Data[20:76,],start=list(B=1,EI=4,EC=50),trace=TRUE,na.action=na.omit)

summary(ModelIA)
summary(ModelCA)

#Synergism-­‐AntagonismModel 

ModelIAS<-­‐nls(Plant~
mean(Data[1:19,3])*pnorm(qnorm(1/((1+(Zinc/EI)^BI)*(1+(Cadmium/EC)^BC)))+(a*(Zinc/EI)*(Cadmium/EC)*((Zinc/EI)+(Cadmium/EC))^(-­‐2))),
data=Data[20:76,],start=list(BI=0.8,BC=2,EI=4.15,EC=65.90,a=0),
trace=TRUE,na.action=na.omit)

ModelCAS<-­‐nls(Plant~
(mean(Data[1:19,3]))/((((Zinc*EC+Cadmium*EI)/(EI*EC))/exp(a*(Zinc/EI)*(Cadmium/EC)*((Zinc/EI)+(Cadmium/EC))^(-­‐2)))^B+1),
data=Data[20:76,],start=list(B=1,EI=4,EC=50,a=0),
trace=TRUE,na.action=na.omit)

summary(ModelIAS)
summary(ModelCAS)

####Comparing Standard model with Deviation model based on F-­‐statistic
anova(ModelIA, ModelIAS)
anova(ModelCA, ModelCAS)

####Comparing IA model with CA model based on ANCOVA
listmodels<-­‐c(rep(0,length(na.omit(Data[5:96,3]))), rep(1,length(na.omit(Data[5:96,3]))))
outputIA<-­‐cbind(na.omit(Data[5:96,3]),fitted(ModelIAS))
outputCA<-­‐cbind(na.omit(Data[5:96,3]),fitted(ModelCAS))
outputtotal<-­‐rbind(outputIA,outputCA)
ancovamodel<-­‐data.frame(listmodels,outputtotal)
test1<-­‐aov(X1~X2*listmodels,data=model)

summary(test1)

##Veryfinng assumptions for F-­‐statistic 
#Assumption of Normality of Residuals
shapiro.test(residuals(ModelIA))
shapiro.test(residuals(ModelIAS))
shapiro.test(residuals(ModelCA))
shapiro.test(residuals(ModelCAS))

#Assumption of Homoscedasticity of Residuals
concrange<-­‐na.omit(Data[20:76,])
concrange<-­‐concrange[,1]+concrange[,2]
leveneTest(residuals(ModelIA),as.factor(concrange))
leveneTest(residuals(ModelIAS),as.factor(concrange))

concrange<-­‐na.omit(Data[20:76,])
concrange<-­‐concrange[,1]+concrange[,2]
leveneTest(residuals(ModelCA),as.factor(concrange))
leveneTest(residuals(ModelCAS),as.factor(concrange))

###Creating data for model predictions for Isobolograms
###Use parameters predicted by model
Zinc<-­‐rep(seq(0,4.5,0.1),100)
Cadmium<-­‐sort(rep(seq(0,45,1),100))

new<-­‐data.frame(cbind(Zinc,Cadmium))
predicteddataIAsingle<-­‐
c(100/((1+(new[,1]/coef(ModelIAsingle)[[3]])^coef(ModelIAsingle)[[1]])*(1+(new[,2]/coef(ModelIAsingle)[[4]])^coef(ModelIAsingle)[[2]])))

contourIAsingle<-­‐data.frame(cbind(predicteddataIAsingle,Zinc,Cadmium))
isobologramIAsingle<-­‐contourIAsingle(predicteddataIAsingle~Zinc*Cadmium, data=contourIAsingle,
cuts=10,region=FALSE,
ylab="Plant root elongation (% of control)",
xlab="Metal dose (μg/L)")

print(isobologramIAsingle)

#Make model predictions for mixture concentrations based on model developed on single stressor data

PredictIAsingle<-­‐predict(ModelIAsingle,Data[50:76,1:2])

#Plot mean model predictions versus observed mean mixture data
plot(Datamean[1:10,],unique(fitted(ModelIAsingle)),
xlab='data', ylab='fitted values')
lines(Datamean[11:19,],unique(PredictIAsingle), type="p", pch=19) 
abline(0,1)

#Plot all model predictions versus all observed mixture data
plot(na.omit(Data[20:49,3]),fitted(ModelIAsingle),
xlab='data', ylab='fitted values')
lines(Data[50:76,3],PredictIAsingle, type="p", pch=19)
abline(0,1)