library(mgcv)
Rdate = as.Date(strptime(RNA$Date, "%m/%d/%Y"))
Rfdate = as.numeric(Rdate)[!is.na(RNA$Ammonia..mg.l.)]
Rfsite = factor(RNA$Site)[!is.na(RNA$Ammonia..mg.l.)]
RfAm = RNA$Ammonia..mg.l.[!is.na(RNA$Ammonia..mg.l.)]
flowobj = gam(log(RfAm + 1) ~ s(Rfdate, k = 90) + Rfsite) #k chosen by preliminary analysis
# flow ammonia model
wwdat$Flow-> FlowMean
Am = wwdat$Ammonia
Pop = rep((sitedata$PercentageNationalPop * sitedata$NationalPop)[isite], length(FlowMean))
# 0 = has flow, 1 = has only ammonia, 2 = has neither
dataType = is.na(FlowMean)*(1+ is.na(Am))
# do we have *any* flow/ammonia data for site?
if ((sum(!is.na(wwdat$Flow)) +sum(!is.na(wwdat$Am)) )>0){
minAm = min(Am, na.rm=T)
#fill in missing ammonias with ammonia from model
if (sum(!is.na(Am)) > 0) Am[is.na(Am)] = exp(predict(flowobj, data.frame(Rfdate = as.numeric(as.Date(strptime(wwdat$Date, "%m/%d/%Y"))), Rfsite = wwdat$Site))[is.na(Am)]) - 1
Am = pmax(Am, minAm) #do not allow ammonia to go below minimum seen
if (sitename %in% ammoniacoef$site){
param = ammoniacoef[ammoniacoef$site == sitename,]
# note adrian is based on log10 values, we reverse the reparameterisation
pred = matrix(rep(10^(param$intercept+param$slope.lgpop*log10( param$Population )+param$slope.lgammonia*log10 (Am[is.na(wwdat$Flow)])), 3), ncol = 3)
} else {
# not in ammonia est, use unknown
param = ammoniacoef[ammoniacoef$site == "Unknown",]
pred = matrix(rep(10^(param$intercept+param$slope.lgpop*log10( Pop[1] )+param$slope.lgammonia*log10 (Am[is.na(wwdat$Flow)])), 3), ncol = 3)
}
FlowMean[is.na(wwdat$Flow)] = pred[,1]
# if prediction fails use mean
FlowMean[is.na(FlowMean)] = mean(FlowMean, na.rm = T)
} else {
# if a site has no flow/am data at all, use a national linear model based on population
logP = log10((sitedata$PercentageNationalPop * sitedata$NationalPop)[match(RNA$Site, sitedata$Simple)])
pred = predict(lm(log10(RNA$Flow) ~logP), newdata = list(logP = log10(Pop)),interval = "prediction", level = WaterPredLev)
FlowMean[is.na(FlowMean)] = 10^pred[,1]
}
var2 = wwdat$N1.Reported*FlowMean/(sitedata$PercentageNationalPop[isite]*sitedata$NationalPop[1]) /1e6 # calc Mgc/p/d