run_MCMC <- function(start.value, N.iter, proposal.sd, posterior.fun, npar){ posterior.fun <- match.fun(posterior.fun) stopifnot(length(proposal.sd) == npar) chain = array(dim = c(N.iter + 1, npar)) all.lik <- rep(NA, N.iter + 1) chain[1,] = start.value all.lik[1] <- posterior.fun(start.value) for (i in 1:N.iter){ proposal <- chain[i,] + rnorm(npar, mean = 0, sd = proposal.sd) all.lik[i + 1] <- posterior.fun(proposal) probab = exp(all.lik[i + 1] - all.lik[i]) if(runif(1) < probab){ chain[i+1,] = proposal }else{ chain[i+1,] = chain[i,] all.lik[i + 1] <- all.lik[i] } } return(list(chain, all.lik)) } get.fourier <- function(data, cons_name, R2threshold){ # this returns the fourier series coefficients for a given antibiotic, returning minimum number of coefficients such that R2 > R2threshold this.cons <- data[, cons_name] time <- data$year.cat TT <- length(this.cons) n <- TT/2 stopifnot(floor(TT/2) == n) # if the data has odd number will be a bit different # compute coefficient of Fourier series: alphas <- rep(NA, n+1) betas <- rep(NA, n+1) alphas[1] <- mean(this.cons) betas[1] <- 0 alphas[2:n] <- sapply(1:(n-1), function(j) 2 * mean(this.cons * cos(2 * pi * j * (0:(TT-1))/ TT))) betas[2:n] <- sapply(1:(n-1), function(j) 2 * mean(this.cons * sin(2 * pi * j * (0:(TT-1))/ TT))) alphas[n+1] <- mean(this.cons * cos(2 * pi * n * (0:(TT-1))/ TT)) betas[n+1] <- 0 # check that we recover this.cons A <- sqrt(alphas^2 + betas^2) P <- atan2(betas, alphas) y2 <- sapply(0:(TT-1), function(t){ sum(A * cos(2 * pi * (0:n) * t / TT - P)) }) stopifnot(max(abs(y2-this.cons)) < 1e-10) TSS <- sum((this.cons - mean(this.cons))^2) # total sum of square for(ii in 1:length(A)){ sub <- order(A, decreasing = T)[1:ii] pred.cons <- sapply(0:(TT-1), function(t){ sum(A[sub] * cos(2 * pi * (0:n)[sub] * t / TT - P[sub])) }) resid <- pred.cons - this.cons # residuals RSS <- sum(resid^2) # Residual sum of squares R2 <- 1- RSS / TSS #print(R2) if(R2 > R2threshold) break } return(cbind((0:n)[sub], A[sub], P[sub])) } get.cons <- function(t, fourier_coeff){ # a vectorised function to get the antibiotic consumption at any time point from the fourier coefficients TT <- 60 # number of time points frequencies <- fourier_coeff[,1] A <- fourier_coeff[,2] P <- fourier_coeff[,3] pred.conso <- sapply(t, function(t){ sum(A * cos(2 * pi * frequencies * t / TT - P)) }) return(pred.conso) } get.integral <- function(c, p0, fourier_coeff){ # get the integral, that is the sum in equation 7 of the paper timeseq <- seq(0, 59, 1) TT <- length(timeseq) n.harmonics <- nrow(fourier_coeff) frequencies <- fourier_coeff[2:n.harmonics, 1] A <- fourier_coeff[2:n.harmonics, 2] P <- fourier_coeff[2:n.harmonics, 3] omegas <- 2 * pi * fourier_coeff[2:n.harmonics, 1] / TT phi <- atan2(omegas, c) # phase difference amp <- 1/sqrt(c^2 + omegas^2) # amplitude pred.integral <- sapply(timeseq, function(t){ sum(amp * A * cos(omegas * t - P - phi)) }) pred.transient <- sapply(timeseq, function(t){ # transient effects sum(A * exp(- c * t) * (omegas * sin(P) - c * cos(P)) / (c^2 + omegas^2)) }) pred.integral <- pred.integral + pred.transient # add the term corresponding to transient effects return(pred.integral) } get.pred.multi2 <- function(c, b.vec, pbar, p0, cons_name.vec, fourier_coeff.list){ # predicted frequency of resistance under the full model (equation 7 of the paper) timeseq <- seq(0, 59, 1) stopifnot(length(b.vec) == length(cons_name.vec)) stopifnot(length(b.vec) == length(fourier_coeff.list)) integral <- list() for(ii in 1:length(b.vec)){ integral[[ii]] <- b.vec[ii] * get.integral(c, p0, fourier_coeff.list[[ii]]) } integral <- do.call(rbind, integral) integral <- colSums(integral) pred.res <- pbar + exp(-c * timeseq) * (p0 - pbar) + integral # add the other form of transient effect return(pred.res) } get.pred.multi.approx2 <- function(ratio.vec, pbar, cons_name.vec, fourier_coeff.list){ # predicted frequency of resistance under the approximate model (equation 5 of the paper) stopifnot(length(ratio.vec) == length(cons_name.vec)) stopifnot(length(ratio.vec) == length(fourier_coeff.list)) timeseq <- seq(0, 59,1) ### UPDATE BELOW ##### integral <- list() for(ii in 1:length(ratio.vec)){ n.harmonics <- nrow(fourier_coeff.list[[ii]]) integral[[ii]] <- ratio.vec[ii] * sapply(timeseq, get.cons, fourier_coeff = fourier_coeff.list[[ii]][2:n.harmonics,]) # in that case integral is simply ratio * the consumption } integral <- do.call(rbind, integral) integral <- colSums(integral) pred.res <- pbar + integral return(pred.res) } lik.multi2 <- function(par, data, res_name, cons_name.vec, Ni, fourier_coeff.list){ # full likelihood function # average number of samples Ni = 37 for Bedouin kids, Ni = 23 for Jewish kids stopifnot(length(par)==length(cons_name.vec)+2) c <- par[1] b.vec <- par[2:(2 + length(cons_name.vec) - 1)] pbar <- par[2 + length(cons_name.vec)] pred <- get.pred.multi2(c, b.vec, pbar, data[1, res_name], cons_name.vec, fourier_coeff.list) if(any(is.na(pred))) return(-10^9) if(any(pred < 0 | pred > 1)) return(-10^9) data <- data[, res_name] stopifnot(length(pred)==length(data)) likelihood <- choose(Ni, round(data * Ni)) * pred^(Ni * data) * (1-pred)^(Ni * (1-data)) return(sum(log(likelihood))) } lik.multi.approx2 <- function(par, data, res_name, cons_name.vec, Ni, fourier_coeff.list){ # simplified likelihood function when c is high stopifnot(length(par) == length(cons_name.vec) + 1) ratio.vec <- par[1:length(cons_name.vec)] pbar <- par[length(cons_name.vec) + 1] pred <- get.pred.multi.approx2(ratio.vec, pbar, cons_name.vec, fourier_coeff.list) if(any(pred < 0 | pred > 1)) return(-10^9) data <- data[, res_name] #stopifnot(length(pred)==length(data)) likelihood <- choose(Ni, round(data * Ni)) * pred^(Ni * data) * (1-pred)^(Ni * (1-data)) return(sum(log(likelihood))) } # THIS COMPUTES THE LIKELIHOOD OF DYNAMICS OF THE FREQUENCY OF RESISTANCE GIVEN ANTIBIOTIC PRESCRIPTION RATES # ALL FUNCTIONS HAVE BEEN TESTED AND CHECKED EXTENSIVELY