# Summary of the statistical fitting process for the insecticide property constants 

library(readxl)
#setwd("~/StochasticSimulations/Data")
indata = read_excel("StatisticalFitting.xlsx",sheet="Data")
indata[indata==100] = 99.9 # processing to remove logit Inf 

# variables 
x = seq(0,48,6)
z_phys = matrix(0,sum(indata[,7]=="Physical"),3)
z_chem = matrix(0,sum(indata[,7]=="Chemical"),3)
z_mort = matrix(0,sum(indata[,7]=="Mortality"),3)

# physical constants 

y = log((indata[indata[,7]=="Physical",12:20]/100)/(1-(indata[indata[,7]=="Physical",12:20]/100))) # logit

for(i in 1:nrow(y)){
  xydata = t(rbind(x,as.vector(y[i,])))
  rownames(xydata) = c()
  colnames(xydata) = c("x","y")
  lr = lm(y~x,data=data.frame(xydata))
  z_phys[i,1:2] = lr$coefficients
}

mean(z_phys[,1]) # ph 
mean(z_phys[,2]) # pd

# chemical constants 

y = log(1/indata[indata[,7]=="Chemical",12:20]) # exponential decay

for(i in 1:nrow(y)){
  xydata = t(rbind(x,as.vector(y[i,])))
  rownames(xydata) = c()
  colnames(xydata) = c("x","y")
  lr = lm(y~x,data=data.frame(xydata))
  z_chem[i,1:2] = lr$coefficients
  if(sum(indata[indata[,7]=="Chemical",11][i]=="Pyrethroid")>0){
    z_chem[i,3] = 1
  }
  if(sum(indata[indata[,7]=="Chemical",11][i]=="PBO")>0){
    z_chem[i,3] = 2
  }
  if(sum(indata[indata[,7]=="Chemical",11][i]=="Pyriproxyfen")>0){
    z_chem[i,3] = 3
  }
  if(sum(indata[indata[,7]=="Chemical",11][i]=="Chlorfenapyr")>0){
    z_chem[i,3] = 4
  }
}
dc = rep(0,4)
hc = rep(0,4)
for(j in 1:4){
  dc[j] = mean(z_chem[z_chem[,3]==unique(z_chem[,3])[j],2])
  hc[j] = mean(z_chem[z_chem[,3]==unique(z_chem[,3])[j],1])
}

# small error in 

# mortality 

y = log((indata[indata[,7]=="Mortality",12:20]/100)/(1-(indata[indata[,7]=="Mortality",12:20]/100))) # logit

for(i in 1:nrow(y)){
  xydata = t(rbind(x,as.vector(y[i,])))
  rownames(xydata) = c()
  colnames(xydata) = c("x","y")
  lr = lm(y~x,data=data.frame(xydata))
  z_mort[i,1:2] = lr$coefficients
  for(j in 1:length(unique(indata[,6]))){
    if(indata[indata[,7]=="Mortality",6][i]==unique(indata[,6])[j]){
      z_mort[i,3] = j
    }
  }
}
dm = rep(0,4)
hm = rep(0,4)
for(j in 1:4){
  dm[j] = mean(z_mort[z_mort[,3]==unique(z_mort[,3])[j],2])
  hm[j] = mean(z_mort[z_mort[,3]==unique(z_mort[,3])[j],1])
}

new_grouping = rep(0,nrow(z_mort))
new_grouping[z_mort[,3]==6|z_mort[,3]==7|z_mort[,3]==8] = 1 # pbo
new_grouping[z_mort[,3]==9] = 2 # chlorfenapyr 
new_grouping[z_mort[,3]==10|z_mort[,3]==11] = 3 # pyriproxyfen

intercept = mean(z_mort[new_grouping==0&(exp(z_mort[,1])/(1+exp(z_mort[,1])))>0.7,1])
gradient = mean(z_mort[new_grouping==0&(exp(z_mort[,1])/(1+exp(z_mort[,1])))>0.7,2])

# pyrethroid case above, for other cases see excel sheet