This presents the parameter fitting for my FOI models. This is based on age prevalence data from Human MST scores.

Model Functions

get_catalytic_result = function(lambda, a){
    z = 1 - exp(-lambda*a)
    return(z)
}

get_variable_result = function(lambda,P, a){
  z = P*(1-exp(-lambda*a))
  return(z)
}

get_reverse_result = function(lambda, rho, a){
    z = (lambda / (lambda + rho)) * (1 - exp(-(lambda + rho)*a))
    return(z)
}

Log Likelihood Functions (binomial distribution)

nll.catalytic = function(log.lambda,df){
    lambda <- exp(log.lambda)
    p <- get_catalytic_result(lambda, df$age)
    nll <- - sum(dbinom(x = FOIMST1Alldf$num.positive, size = FOIMST1Alldf$num.tested,
                        prob = p, log = TRUE))
    return(nll)
}

nll.variable = function(par, df){
  lambda <- exp(par[1])
  P <- exp(par[2])
  p <- get_variable_result(lambda, P, df$age)
  nll <- - sum(dbinom(x = FOIMST1Alldf$num.positive, size = FOIMST1Alldf$num.tested,
                      prob = p, log = TRUE))
  return(nll)
}


nll.reverse = function(par, df){
    # params is a vector of two values
    lambda = exp(par[1])
    rho = exp(par[2])
    p = get_reverse_result(lambda, rho, df$age)
    nll <- - sum(dbinom(x = FOIMST1Alldf$num.positive, size = FOIMST1Alldf$num.tested,
                        prob = p, log = TRUE))
    return(nll)
}

Testing Model Parameters

Simple Catalytic Loop

## Set up
lambdavec <- seq(0,0.1,0.01)



## Simple Catalytic Loop

Catnllvec <- data.frame(Lambda = NA, 'NLL Cat' = NA)
i <- 1

for(i in c(1:length(lambdavec))){
  
  lambda <- lambdavec[i]
  
  p.cat <- get_catalytic_result(lambda, FOIMST1Alldf$age)
  nll.cat <- sum(dbinom(x = FOIMST1Alldf$num.positive, size = FOIMST1Alldf$num.tested,
                      prob = p.cat, log = TRUE))
    
  
  Catnllnew <- c(lambda,nll.cat)
    
   Catnllvec <- rbind(Catnllvec, Catnllnew)

  i <- i + 1
  
}

Catnllvec <- Catnllvec[-1,]

Simple Catalytic Model Lambda Plot

Variable Catalytic Model Loop

#Set Up
lambdavec <- seq(0,0.2,0.01)
Pvec <- seq(0.1,1,0.01)



Varnllvec <- data.frame(Lambda = NA, P = NA, 'NLLVar' = NA)
i <- 1
r <- 1

##Variable catalytic loop

for(i in c(1:length(lambdavec))){
  for(r in c(1:length(Pvec))){
    
    lambda <- lambdavec[i]
    P <- Pvec[r]
    
    p.Var <- get_variable_result(lambda, P, FOIMST1Alldf$age)
    nll.Var <- sum(dbinom(x = FOIMST1Alldf$num.positive, size = FOIMST1Alldf$num.tested,
                          prob = p.Var, log = TRUE))
    
    
    Varnllnew <- c(lambda,P, nll.Var)
    
    Varnllvec <- rbind(Varnllvec, Varnllnew)
    
    r <- r + 1}
  i <- i + 1}
Varnllvec <- Varnllvec[-1,]

Variable Catalytic Model Tile Plot

Reverse Model Loop

#Set Up 
lambdavec <- seq(0,0.1,0.005)
rhovec <- seq(0,0.1,0.005)


Revnllvec <- data.frame(Lambda = NA, Rho = NA, 'NLLRev' = NA)
i <- 1
r <- 1

##Reverse catalytic loop

for(i in c(1:length(lambdavec))){
for(r in c(1:length(rhovec))){
  
  lambda <- lambdavec[i]
  rho <- rhovec[r]
  
  p.Rev <- get_reverse_result(lambda,rho, FOIMST1Alldf$age)
  nll.Rev <- sum(dbinom(x = FOIMST1Alldf$num.positive, size = FOIMST1Alldf$num.tested,
                         prob = p.Rev, log = TRUE))
  
  
  Revnllnew <- c(lambda,rho, nll.Rev)
  
  Revnllvec <- rbind(Revnllvec, Revnllnew)
  
  r <- r + 1}
  i <- i + 1}
Revnllvec <- Revnllvec[-1,]

Reverse Model Tile Plot

Optimise Parameters, NLL and AIC

Simple Catalytic Model Optim

cat.log.lambda.guess <- log(0.03)
cat.log.param <- optim(par = cat.log.lambda.guess, fn = nll.catalytic, df = FOIMST1Alldf, 
                   method = "Brent", lower = log(0.01), upper = log(0.04))


cat.lambda <- exp(cat.log.param$par)
cat.nll <- cat.log.param$value
  
FOIMST1Alldf$CatMod <-  get_catalytic_result(cat.lambda,FOIMST1Alldf$age)

##Test Cat Model with AIC


Cat.AIC <- -2*cat.nll+2*1

Variable Catalytic Model Optim

var.log.lambda.guess <- log(0.05)
var.log.P.guess <- log(0.8)
var.log.param <- optim(par = c(var.log.lambda.guess, var.log.P.guess), fn = nll.variable, df = FOIMST1Alldf,
                      method = "Nelder-Mead")


var.lambda <- exp(var.log.param$par)[1]
var.P <- exp(var.log.param$par)[2]
var.nll <- var.log.param$value

FOIMST1Alldf$VarMod <- get_variable_result(var.lambda,var.P,FOIMST1Alldf$age)



## Test Var model with AIC

Var.AIC <- -2*var.nll+2*2

Reverse Model Optim

rev.log.lambda.guess <- log(0.05)
rev.log.rho.guess <- log(0.02)
rev.log.param <- optim(par = c(rev.log.lambda.guess, rev.log.rho.guess), fn = nll.reverse, df = FOIMST1Alldf,
             method = "Nelder-Mead")


rev.lambda <- exp(rev.log.param$par)[1]
rev.rho <- exp(rev.log.param$par)[2]
rev.nll <- rev.log.param$value

FOIMST1Alldf$RevMod <- get_reverse_result(rev.lambda,rev.rho,FOIMST1Alldf$age)


# Test Rev model with AIC


Rev.AIC <- -2*rev.nll+2*2

Plot of All Data Model (plus error bars)

# Tile plots with optimised parameters

Variable Catalytic Model Heatmap with Optim Parameters

  Varnllvec$NLLVar<- factor(Varnllvec$NLLVar)
N <- nlevels(Varnllvec$NLLVar)


colors <- colorRampPalette(c("blue", "green", "yellow", "red"))(length(unique(Varnllvec$NLLVar)))


ggplot(data = Varnllvec)+
  geom_tile(aes(x = Lambda , y = P, fill = NLLVar))+  
  scale_fill_manual(values=colors, breaks=levels(Varnllvec$NLLVar)[seq(1, N, by=100)])+
  geom_point(x = var.lambda, y = var.P)

Reverse Model Heatmap with Optim Parameters

Revnllvec$NLLRev<- factor(Revnllvec$NLLRev)
N <- nlevels(Revnllvec$NLLRev)
colors <- colorRampPalette(c("blue", "green", "yellow", "red"))(length(unique(Revnllvec$NLLRev)))


 ggplot(data = Revnllvec)+
  geom_tile(aes(x = Lambda, y = Rho, fill = NLLRev))+  
   scale_fill_manual(values=colors, breaks=levels(Revnllvec$NLLRev)[seq(1, N, by=50)])+
   geom_point(x = rev.lambda, y = rev.rho)

Warning: Need to debug past this point

Model Fitting loop (Fitting by Group Vs All Data)

This loop begins by fitting model parameters to the data from each group and assess how well the models fit using AIC. Then it sees how well the model fit to all the data fits each group.

#Set Up



FOIMST1dfNew <- data.frame("Group" = NA,"Age.Group" = NA, "num.tested"= NA, "num.positive" = NA,"Proportion.Positive" = NA,
                           "Age.Group.Start" = NA,  "Age.Group.Finish" = NA, "age" = NA,
                           "CatGroupMod" = NA, "CatGroupModLambda" = NA, "CatGroupModNLL" = NA, "CatGroupModAIC" = NA,
                           "VarGroupMod" = NA, "VarGroupModLambda" = NA, "VarGroupModP" = NA,"VarGroupModNLL" = NA, "VarGroupModAIC" = NA,
                           "RevGroupMod" = NA, "RevGroupModLambda" = NA, "RevGroupModRho" = NA,"RevGroupModNLL" = NA,"RevGroupModAIC" = NA,
                           "CatAllDataMod" = NA, "CatAllDataModLambda" = NA, "CatAllDataModNLL" = NA, "CatAllDataModAIC" = NA,
                           "VarAllDataMod" = NA, "VarAllDataModLambda" = NA, "VarAllDataModP" = NA,"VarAllDataModNLL" = NA, "VarAllDataModAIC" = NA,
                           "RevAllDataMod" = NA, "RevAllDataModLambda" = NA, "RevAllDataModRho" = NA,"RevAllDataModNLL" = NA, "RevAllDataModAIC" = NA)

for(G in unique(FOIMST1df$Group)){
  
FOIMST1dfG <-  subset(FOIMST1df, FOIMST1df$Group == G)


  #  Simple Catalytic model 
  #Group Models
  
  cat.log.lambda.guess <- log(0.03)
  cat.log.param <- optim(par = cat.log.lambda.guess, fn = nll.catalytic, df = FOIMST1dfG, 
                         method = "Brent", lower = log(0.01), upper = log(0.09))
  lambda <- exp(cat.log.param$par)
  
  cat.nll <- cat.log.param$value
  cat.aic <-  -2*cat.nll+2*1
  
  
  FOIMST1dfG$CatGroupMod <-  get_catalytic_result(lambda, FOIMST1dfG$age)
  FOIMST1dfG$CatGroupModLambda <- rep(lambda, NROW(FOIMST1dfG))
  FOIMST1dfG$CatGroupModNLL <- rep(cat.nll, NROW(FOIMST1dfG))
  FOIMST1dfG$CatGroupModAIC <- rep(cat.aic, NROW(FOIMST1dfG))
  
  


  
  #  Variable Catalytic model
  
  var.log.lambda.guess <- log(0.05)
  var.log.P.guess <- log(0.8)
  var.log.param <- optim(par = c(var.log.lambda.guess, var.log.P.guess), fn = nll.variable, df = FOIMST1dfG,
                         method = "Nelder-Mead")
  lambda <- exp(var.log.param$par)[1]
  P <-  exp(var.log.param$par)[2]
  var.nll <- var.log.param$value
  var.aic <-  -2*var.nll+2*2
  
  
  FOIMST1dfG$VarGroupMod <-  get_variable_result(lambda,P, FOIMST1dfG$age)
  FOIMST1dfG$VarGroupModLambda <- rep(lambda, NROW(FOIMST1dfG))
  FOIMST1dfG$VarGroupModP <- rep(P, NROW(FOIMST1dfG))
  
  FOIMST1dfG$VarGroupModNLL <- rep(var.nll, NROW(FOIMST1dfG))
  FOIMST1dfG$VarGroupModAIC <- rep(var.aic, NROW(FOIMST1dfG))
  
  
  # Reverse model
  rev.log.lambda.guess <- log(0.05)
  rev.log.rho.guess <- log(0.02)
  rev.log.param <- optim(par = c(rev.log.lambda.guess, rev.log.rho.guess), fn = nll.reverse, df = FOIMST1dfG,
                         method = "Nelder-Mead")
  lambda <- exp(rev.log.param$par)[1]
  Rho <- exp(rev.log.param$par)[2]
  
  rev.nll <- rev.log.param$value
  rev.aic <-  -2*rev.nll+2*2
  
  
  
  FOIMST1dfG$RevGroupMod <-  get_reverse_result(lambda,Rho, FOIMST1dfG$age)
  FOIMST1dfG$RevGroupModLambda <- rep(lambda, NROW(FOIMST1dfG))
  FOIMST1dfG$RevGroupModRho <- rep(Rho, NROW(FOIMST1dfG))
  
  FOIMST1dfG$RevGroupModNLL <- rep(rev.nll, NROW(FOIMST1dfG))
  FOIMST1dfG$RevGroupModAIC <- rep(rev.aic, NROW(FOIMST1dfG))
  
  
  #All data 
  
  
  #Simple catalytic
  
  cat.log.lambda.guess <- log(0.03)
  cat.log.param <- optim(par = cat.log.lambda.guess, fn = nll.catalytic, df = FOIMST1Alldf, 
                         method = "Brent", lower = log(0.01), upper = log(0.04))
  
  
  cat.lambda <- exp(cat.log.param$par)
  
  
  FOIMST1dfG$CatAllDataMod <-  get_catalytic_result(cat.lambda, FOIMST1dfG$age)
  cat.nll <- nll.catalytic(cat.log.param$par,FOIMST1dfG)
  cat.aic <- -2*cat.nll+2*1
  
  
  FOIMST1dfG$CatAllDataModLambda <- rep(cat.lambda, NROW(FOIMST1dfG))
  FOIMST1dfG$CatAllDataModNLL <- rep(cat.nll, NROW(FOIMST1dfG))
  FOIMST1dfG$CatAllDataModAIC <- rep(cat.aic, NROW(FOIMST1dfG))
  
  
  #Variable catalytic
  var.log.lambda.guess <- log(0.05)
  var.log.P.guess <- log(0.8)
  var.log.param <- optim(par = c(var.log.lambda.guess, var.log.P.guess), fn = nll.variable, df = FOIMST1Alldf,
                         method = "Nelder-Mead")
  
  
  var.lambda <- exp(var.log.param$par)[1]
  var.P <- exp(var.log.param$par)[2]

  FOIMST1dfG$VarAllDataMod <- get_variable_result(var.lambda,var.P,FOIMST1dfG$age)
  var.nll <- nll.variable(var.log.param$par,FOIMST1dfG)
  var.aic <- -2*var.nll+2*2
  
  
  
  FOIMST1dfG$VarAllDataModLambda <- rep(var.lambda, NROW(FOIMST1dfG))
  FOIMST1dfG$VarAllDataModP <- rep(var.P, NROW(FOIMST1dfG))
  FOIMST1dfG$VarAllDataModNLL <- rep(var.nll, NROW(FOIMST1dfG))
  FOIMST1dfG$VarAllDataModAIC <- rep(var.aic, NROW(FOIMST1dfG))
  
  
  #Reverse
  
  rev.log.lambda.guess <- log(0.05)
  rev.log.rho.guess <- log(0.02)
  rev.log.param <- optim(par = c(rev.log.lambda.guess, rev.log.rho.guess), fn = nll.reverse, df = FOIMST1Alldf,
                         method = "Nelder-Mead")
  
  
  rev.lambda <- exp(rev.log.param$par)[1]
  rev.rho <- exp(rev.log.param$par)[2]

  FOIMST1dfG$RevAllDataMod <- get_reverse_result(rev.lambda,rev.rho,FOIMST1dfG$age)
  rev.nll <- nll.reverse(rev.log.param$par,FOIMST1dfG)
  rev.aic <-  -2*rev.nll+2*2
  
  
  FOIMST1dfG$RevAllDataModLambda <- rep(rev.lambda, nrow(FOIMST1dfG))
  FOIMST1dfG$RevAllDataModRho <- rep(rev.rho, nrow(FOIMST1dfG))
  FOIMST1dfG$RevAllDataModNLL <- rep(rev.nll, NROW(FOIMST1dfG))
  FOIMST1dfG$RevAllDataModAIC <- rep(rev.aic, NROW(FOIMST1dfG))
  
  
  
  FOIMST1dfNew <- rbind(FOIMST1dfNew, FOIMST1dfG)}

FOIMST1df <- FOIMST1dfNew[-1,]

Comparing the parameters estimated for group models to heatmaps to help debug

FOIMST1dfG <- subset(FOIMST1df, FOIMST1df$Group == 'MGROS') #can change to any group to test


var.log.lambda.guess <- log(0.06)
var.log.P.guess <- log(0.78)
var.log.param <- optim(par = c(var.log.lambda.guess, var.log.P.guess), fn = nll.variable, df = FOIMST1dfG )


var.lambda <- exp(var.log.param$par)[1]
var.P <- exp(var.log.param$par)[2]
var.nll <- var.log.param$value



#Variable catalytic loop

lambdavec <-seq(0.01,1,0.01)
Pvec <- seq(0,1,0.01)

Varnllvec <- data.frame(Lambda = NA, P = NA, 'NLLVar' = NA)



for(i in c(1:length(lambdavec))){
  for(r in c(1:length(Pvec))){
    
    lambda <- lambdavec[i]
    P <- Pvec[r]
    
    p.Var <- get_variable_result(lambda, P, FOIMST1dfG$age)
    nll.Var <- sum(dbinom(x = FOIMST1dfG$num.positive, size = FOIMST1dfG$num.tested,
                          prob = p.Var, log = TRUE))
    
    
    Varnllnew <- c(lambda,P, nll.Var)
    
    Varnllvec <- rbind(Varnllvec, Varnllnew)
    
    r <- r + 1}
  i <- i + 1}
Varnllvec <- Varnllvec[-1,]


x <- subset(Varnllvec, Varnllvec$NLLVar != -Inf)
y <-subset(x, NLLVar == max(x$NLLVar))


#Tile Plot


colors <- colorRampPalette(c("blue", "green", "yellow", "red"))(length(unique(Varnllvec$NLLVar)))
Varnllvec$NLLVar<- factor(Varnllvec$NLLVar)
N <- nlevels(Varnllvec$NLLVar)


ggplot(data = Varnllvec)+
  geom_tile(aes(x = Lambda , y = P, fill = NLLVar))+  
  scale_fill_manual(values=colors, breaks=levels(Varnllvec$NLLVar)[seq(1, N, by=500)])+
  geom_point(x = var.lambda, y = var.P)

FOI Parameter Fitting

Katie Farnes

03/09/2020