##Probability Unit (Probit) Analysis R Code from the following:

#Journal name: Journal of Visualized Experiments

#Article title: Topical Application Assay to Quantify Insecticide Toxicity for Mosquitoes and Fruit Flies

#Article authors: Brook Jensen(1), Rachel Althoff(1), Sarah Rydberg(1), Emma Royster(1), Alden Estep(2), Silvie Huijben(1)

#Author affiliations: (1) Center for Evolution and Medicine, School of Life Sciences, Arizona State University, Life Sciences C, 427 East Tyler Mall, Tempe, AZ 85281, USA; (2) United States Department of Agriculture, Agricultural Research Service, Center for Medical, Agricultural and Veterinary Entomology, 1600 SW 23rd Drive, Gainesville, FL 32608, USA 

#R version used: 4.1.1 

#Code below used with "Supplemental_Example.Data" 
    #Description of column headings within example data file provided at the end of script
 
#Upload necessary libraries 
library(aod)
library(ggplot2)
library(MASS)

#Select data file
entryf<-file.choose()

#Read data file into R
dataf<-(read.table(entryf,header=T,stringsAsFactors=T)) 

#Database calculations 
dataf$mort=dataf$dead/dataf$total #Calculates mortality
dataf$dose_mosq=(dataf$syringe*dataf$concentration*1000)/(dataf$weight/dataf$total.mosq) #Calculates dose in ng per average weight (mg) of mosquito

dataf$exp=as.factor(paste(as.character(as.integer(dataf$date)),dataf$sex)) #Create unique name per experiment (for Abbott's formula)
dataf$mort_abbot=NA #Creates vector to store Abbott-corrected mortality rates

for (i in 1:length(dataf$id)) { #Performs Abbott's mortality correction
c=dataf$exp[i] #Identifies correct control for correction
dataf$mort_abbot[i]<-(dataf$mort[i]-dataf$mort[dataf$dose==0 & dataf$exp==c])/(1-dataf$mort[dataf$dose==0 & dataf$exp==c])
i+1}

dataf$mort_abbot=pmax(dataf$mort_abbot,0) #Changes negative values following Abbott's following correction to 0
dataf$mort_f=dataf$mort_abbot #Creates new vector for final mortality manipulations
dataf$mort_f[dataf$mort_abbot==1]<-NA #Replaces all 1's with NA 
dataf$mort_f[dataf$mort_abbot==0]<-NA #Replaces all 0's with NA
dataf$mort_probit=5+sapply(dataf$mort_f,qnorm) #Applies probit transformation, add 5 to have positive values


#Check which strains, sexes, and ages of insects are being tested and on what days testing took place
strains<-levels(dataf$strain); strains 
sexes<-levels(dataf$sex); sexes
dates<-levels(dataf$date); dates
ages<-levels(dataf$age); ages

#Subset data by strain and sex and age
ROCK.all <- subset(dataf, strain=="ROCK"); ROCK.all
ROCK.females <- subset(ROCK.all, sex=="F"); ROCK.females
ROCK.males <- subset(ROCK.all, sex=="M"); ROCK.males
IICC.all <- subset(dataf, strain=="IICC"); IICC.all
IICC.females <- subset(IICC.all, sex=="F"); IICC.females
IICC.males <- subset(IICC.all, sex=="M"); IICC.males
IICC.females.young <- subset(IICC.females, age=="young"); IICC.females.young
IICC.females.old <- subset(IICC.females, age=="old"); IICC.females.old
IICC.males.young <- subset(IICC.males, age=="young"); IICC.males.young
IICC.males.old <- subset(IICC.males, age=="old"); IICC.males.old

#Linear regression and chi squared test for each strain
ROCKprobit <- glm((mort_f) ~ log10(dose_mosq)  + sex, family =quasibinomial(link = probit), data = ROCK.all)
summary.ROCK<-summary(ROCKprobit); summary.ROCK
confint(ROCKprobit) #Calculates confidence intervals 
wald.test(b = coef(ROCKprobit), Sigma = vcov(ROCKprobit), Terms = 1:3) #Chi-Square test to assess glm

IICCprobit <- glm((mort_f) ~ log10(dose_mosq)  + sex, family =quasibinomial(link = probit), data = IICC.all)
summary(IICCprobit)
confint(IICCprobit) #Calculates confidence intervals 
wald.test(b = coef(IICCprobit), Sigma = vcov(IICCprobit), Terms = 1:3) #Chi-Square test to assess glm

#Linear regression for trend line
ROCK_fprobit=glm((mort_f) ~ log10(dose_mosq), family =quasibinomial(link = probit), data = ROCK.females)
IICC_fprobit=glm((mort_f) ~ log10(dose_mosq), family =quasibinomial(link = probit), data = IICC.females.young)
ROCK_mprobit=glm((mort_f) ~ log10(dose_mosq), family =quasibinomial(link = probit), data = ROCK.males)
IICC_mprobit=glm((mort_f) ~ log10(dose_mosq), family =quasibinomial(link = probit), data = IICC.males.young)

quartz(w=8,h=10) #Opens window for plotting for mac users
#windows(w=8,h=10) #Opens window for plotting for windows users
par(mfrow=c(2,1)) #Starts plotting routine, 2 plots in 2 rows and 1 column

#Graph Probit values and log(Dose) - females
options(scipen=999) #Changes the 'scientific number' format to just 0.0005 etc. (for remainder of script, to revert: options(scipen=0))
plot((ROCK.females$mort_probit)-5 ~ ROCK.females$dose_mosq, log='x',xaxt='n',yaxt='n',type="p",col="dodgerblue4", xlim=c(0.0001,3), ylim=c(-2.3, 2.33),ylab = "Mortality", xlab = "Deltamethrin dose (ng/mg mosquito)",pch=16) #Creates starting plot
#Manually create desired tickmarks for logarithmic scale on x-axis
x1=seq(log10(0.0005)-1,log(0.5)+1) #Adds range of doses for main tickmarks
ticksat <- as.vector(sapply(x1, function(p) (1:10)*10^p)) #Calculates places for small tickmarks
axis(1, 10^x1); axis(1, ticksat, labels=NA, tcl=-0.25, lwd=0, lwd.ticks=1) #Adds small and large tickmarks to plot
#Manually create desired tickmarks for probit scale on y-axis
y1<-c(0.01,0.05,seq(0.1,0.9,0.1),0.95,0.99)
axis(2, at=qnorm(y1),labels=y1,las=2, adj=0) #Converts percent mortality into probit
points((IICC.females.young$mort_probit)-5 ~ IICC.females.young$dose_mosq, type="p", pch=17, col="firebrick3")#Adds next data points
abline(ROCK_fprobit,lwd=2)
abline(IICC_fprobit,lwd=2)
title(main= substitute(paste(italic("Aedes aegypti")," - females")))#Creates title with scientific name in italics
legend(0.0001, 2.3, c("Rockefeller","IICC"), col=c("dodgerblue4","firebrick3"), pch=c(16,17))#Creates a legend

#Graph Probit values and log(Dose) - males
plot((ROCK.males$mort_probit)-5 ~ ROCK.males$dose_mosq, log='x',xaxt='n',yaxt='n',type="p",col="dodgerblue4", xlim=c(0.0001,3), ylim=c(-2.3, 2.33), ylab = "Mortality", xlab = "Deltamethrin dose (ng/mg mosquito)",pch=16) #Creates starting plot
#Manually create desired tickmarks for logarithmic scale on x-axis
x1=seq(log10(0.0005)-1,log(0.5)+1) #Adds range of doses for main tickmarks
ticksat <- as.vector(sapply(x1, function(p) (1:10)*10^p)) #Calculates places for small tickmarks
axis(1, 10^x1); axis(1, ticksat, labels=NA, tcl=-0.25, lwd=0, lwd.ticks=1) #Adds small and large tickmarks to plot
#Manually create desired tickmarks for probit scale on y-axis
y1<-c(0.01,0.05,seq(0.1,0.9,0.1),0.95,0.99) 
axis(2, at=qnorm(y1),labels=y1,las=2, adj=0) #Converts percent mortality into probit
points((IICC.males.young$mort_probit)-5 ~ IICC.males.young$dose_mosq, type="p", pch=17, col="firebrick3")#Adds next data points
abline(ROCK_mprobit,lwd=2)
abline(IICC_mprobit,lwd=2)
title(main= substitute(paste(italic("Aedes aegypti")," - males")))#Creates title with scientific name in italics
legend(0.0001, 2.3, c("Rockefeller","IICC"), col=c("dodgerblue4","firebrick3"), pch=c(16,17))#Creates a legend

text(0.00005, 12.9,"A.",cex=2,xpd=NA) #Adds figure panel letter "A"
text(0.00005, 3.2,"B.",cex=2,xpd=NA) #Adds figure panel letter "B"
#End of graphing routine

#LD50 Calculation (95%CI = 1.96 * SE)
ROCK.LD50<-10^(dose.p(ROCKprobit,p=c(0.5))); ROCK.LD50
IICC.LD50<-10^(dose.p(IICCprobit,p=c(0.5))); IICC.LD50
#Compare 95%CIs to determine if the strains have significantly different dose-responses (if the CIs don't overlap, they are significantly different)

#Resistance Ratio Calculation
RR_ROCK <- ROCK.LD50/ROCK.LD50; RR_ROCK
RR_IICC <- IICC.LD50/ROCK.LD50; RR_IICC #ROCK is used as the reference strain

##Descriptions of column headings within example data file: 

# "id" = identification code of each data point
# "species" = species name (e.g., Aedes aegypti)
# "insecticide" = name of insecticide topically applied (e.g., Deltamethrin)
# "strain" = name of mosquito strain (e.g., ROCK)
# "date" = start date topical application 
# "sex" = sex of mosquitoes  
# "age" = age of mosquitoes (young = 3-5 day-old; old = 4 weeks old)
# "total.mosq" = total number of mosquitoes weighed in batch
# "weight" = weight (mg) of all mosquitoes within batch
# "concentration" = concentration of insecticide (ug/ml)  
# "syringe" = droplet volume (mL) of syringe  
# "dose" = amount of insecticide active ingredient applied to each mosquito (ng) 
# "total" = number of mosquitoes in each cup 
# "dead" = number of dead mosquitoes in each cup