Nothing
# modify input arguments so not order dependent; a1 = a1, etc...
#' Negative log-likelihood function: recovery model
#'
#' Function returning negative log-likelihood (nll) for patterns of survival in
#' infected and uninfected treatments, when infected hosts can recover from
#' infection.
#'
#' This model assumes all the hosts in an infected treatment are all initially
#' infected, and they can all potentially recover from infection. Recovered
#' hosts are assumed to only experience background mortality equivalent to that
#' experienced by matching uninfected or control individuals; no assumptions are
#' made as to whether they are still infected or infectious. It is also assumed
#' that the timing of recovery from infection is not directly observed, but an
#' individual's infected/recovery status can be determined after they have died
#' or been censored.
#'
#' The probability that an infection 'survives' over time, i.e., the host does
#' not recover from infection, is assumed to follow a probability distribution
#' which acts independently of the probability distributions determining
#' background mortality or mortality due to infection.
#'
#' This function only estimates location and scale parameters as constants, it
#' is not designed to estimate them as functions.
#'
#'
#' @param a1,b1 location & scale parameters for background mortality
#' @param a2,b2 location & scale parameters for mortality due to infection
#' @param a3,b3 location & scale parameters for how long infection 'survives'
#' @param data a data.frame with the data
#' @param d1,d2,d3 probability distributions for background mortality, mortality
#' due to infection & how long infection 'survives' ("Weibull", "Gumbel",
#' "Frechet")
#' @return numeric
#' @section Warning: requires the data to be specified in a specific format;
#' see vignette 'data format' for details
#' @examples
#' \donttest{
#' # NB the data to analyse needs to be in a data frame of a specific form
#' head(recovery_data)
#'
#' # step #1: prepare nll function for analysis
#' m01_prep_function <- function(a1, b1, a2, b2, a3, b3){
#' nll_recovery(a1, b1, a2, b2, a3, b3,
#' data = recovery_data,
#' d1 = "Weibull", d2 = "Weibull", d3 = "Weibull"
#' )}
#'
#' # step #2: send 'prep_function' to mle2 for maximum likelihood estimation,
#' # specifying starting values
#' m01 <- mle2(m01_prep_function,
#' start = list(a1 = 2, b1 = 0.5, a2 = 2, b2 = 0.5, a3 = 2, b3 = 0.5)
#' )
#'
#' summary(m01)
#'
#' # values used to simulate data were for the Weibull distribution;
#' # a1 = 2.8, b1 = 0.5, a2 = 2.2, b2 = 0.35, a3 = 2.35, b3 = 0.35
#' }
nll_recovery <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""
){
nll <- P_recovery_calc_likelihood_matrix(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
}
### subfunctions used by above ###
P_recovery_calc_survival_functions <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""){
# function returns matrix with values for various survival functions
# from t[1] -> t[max]
# given values of a1, b1, etc.
tmax <- max(data$t)
calc.matrix <- matrix(0,tmax,13)
colnames(calc.matrix) <- c("t", "f1", "f2", "f3", "S1", "S2", "S3",
"h1", "h2", "h3", "f1S2S3", "f2S1S3", "f3S1S2")
for (t in 1:tmax){
z1 <- P_get_zx(t, a1, b1, d1)
z2 <- P_get_zx(t, a2, b2, d2)
z3 <- P_get_zx(t, a3, b3, d3)
f1 <- P_get_fx(t, z1, b1, d1)
f2 <- P_get_fx(t, z2, b2, d2)
f3 <- P_get_fx(t, z3, b3, d3)
S1 <- P_get_Sx(t, z1, d1)
S2 <- P_get_Sx(t, z2, d2)
S3 <- P_get_Sx(t, z3, d3)
h1 <- P_get_hx(t, z1, b1, d1)
h2 <- P_get_hx(t, z2, b2, d2)
h3 <- P_get_hx(t, z3, b3, d3)
f1S2S3 <- f1 * S2 * S3
f2S1S3 <- f2 * S1 * S3
f3S1S2 <- f3 * S1 * S2
calc.matrix[t,1] <- t
calc.matrix[t,2] <- f1
calc.matrix[t,3] <- f2
calc.matrix[t,4] <- f3
calc.matrix[t,5] <- S1
calc.matrix[t,6] <- S2
calc.matrix[t,7] <- S3
calc.matrix[t,8] <- h1
calc.matrix[t,9] <- h2
calc.matrix[t,10] <- h3
calc.matrix[t,11] <- f1S2S3
calc.matrix[t,12] <- f2S1S3
calc.matrix[t,13] <- f3S1S2
}
calc.matrix
}
###
P_recovery_calc_f3S1S2 <- function(a1 = a1, b1 = b1,
a2 = a2, b2 = b2,
a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""){
# returns matrix with calculations for probabilities of recovery
# while host still alive at different times given values of a1, b1, etc...
tmax <- max(data$t)
matrix01 <- matrix(0,tmax,tmax)
matrix02 <- matrix(0,tmax,tmax)
for (t in 1:tmax){
for (u in 1:t){
matrix01[t,u] <- u
zu1 <- P_get_zu(u, a1, b1, d1)
zu2 <- P_get_zu(u, a2, b2, d2)
zu3 <- P_get_zu(u, a3, b3, d3)
Su1 <- P_get_Su(u, zu1, d1)
Su2 <- P_get_Su(u, zu2, d2)
fu3 <- P_get_fu(u, zu3, b3, d3)
matrix02[t,u] <- fu3 * Su1 * Su2
}
}
matrix02
}
###
P_recovery_calc_St_Sr <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""){
# returns matrix with calculations for probability of surviving
# from time of recovery, S[r], to time t, S[t]: S[t]/S[r]
# given values of a1, b1, etc.
tmax <- max(data$t)
matrix03 <- matrix(0,tmax,tmax)
matrix04 <- matrix(0,tmax,tmax)
for (t in 1:tmax){
for (u in 1:t){
matrix03[t,u] <- t
z1t <- P_get_zx(t, a1, b1, d1)
z1u <- P_get_zu(u, a1, b1, d1)
S1t <- P_get_Sx(t, z1t, d1)
S1u <- P_get_Su(u, z1u, d1)
matrix04[t,u] <- S1t / S1u
}
}
matrix04
}
###
P_recovery_calc_f3S1S2_St_Sr <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""){
# calculates probabilities of recovering at different times
# and surviving to time t
# from product of recovery_calc_f3S1S2 * recovery_St.Sr
# given values of a1, b1, etc.
f3S1S2 <- P_recovery_calc_f3S1S2(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
St.Sr <- P_recovery_calc_St_Sr(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
f3S1S2_St_Sr <- f3S1S2 * St.Sr
}
###
P_recovery_calc_sum_f3S1S2_St_Sr <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""){
# calculates sum of probabilities of recovering at different times
# and surviving to time t given values of a1, b1, etc.
tmax <- max(data$t)
f3S1S2.St.Sr <- P_recovery_calc_f3S1S2_St_Sr(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
vector.sum.f3S1S2.St.Sr <- rowSums(f3S1S2.St.Sr)
sum_f3S1S2_St_Sr <- matrix(vector.sum.f3S1S2.St.Sr, , 1)
}
###
P_recovery_calc_likelihood_matrix <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = ""){
# calculates likelihoods using 'recovery' sub-functions
# values a1, b1, etc...
# define maximum time
tmax <- max(data$t)
# convert data.frame(data) to matrix
matrix01 <- as.matrix(data)
# calculate various survival functions (defined elsewhere) based on values a1,a2,...
survival.functions <- P_recovery_calc_survival_functions(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
f3S1S2 <- P_recovery_calc_f3S1S2(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
St.Sr <- P_recovery_calc_St_Sr(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
f3S1S2.St.Sr <- P_recovery_calc_f3S1S2_St_Sr(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
sum.f3S1S2.St.Sr <- P_recovery_calc_sum_f3S1S2_St_Sr(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, a3 = a3, b3 = b3,
data = data, d1 = "", d2 = "", d3 = "")
h1 <- matrix(survival.functions[, 'h1'], , 1)
# create & fill matrix for likelihood calculations
likelihood.matrix <- matrix(0,nrow(matrix01), 15)
likelihood.matrix[, 01] <- matrix01[, "t"]
likelihood.matrix[, 02] <- matrix01[, 'control.d'] *
survival.functions[, 'f1']
likelihood.matrix[, 03] <- matrix01[, 'control.c'] *
survival.functions[, 'S1']
likelihood.matrix[, 04] <- matrix01[, 'infected.d'] *
(survival.functions[, 'f1S2S3'] + survival.functions[, 'f2S1S3'])
likelihood.matrix[, 05] <- matrix01[, 'infected.c'] *
(survival.functions[, 'S1'] *
survival.functions[, 'S2'] *
survival.functions[, 'S3'])
likelihood.matrix[, 06] <- sum.f3S1S2.St.Sr * h1 # recovered, died
likelihood.matrix[, 07] <- matrix01[, 'recovered.d']
likelihood.matrix[, 08] <- likelihood.matrix[, 6] * likelihood.matrix[, 7]
likelihood.matrix[, 09] <- sum.f3S1S2.St.Sr # recovered, censored
likelihood.matrix[, 10] <- matrix01[, 'recovered.c']
likelihood.matrix[, 11] <- likelihood.matrix[, 9] * likelihood.matrix[, 10]
likelihood.matrix[, 12] <- likelihood.matrix[, 2] +
likelihood.matrix[, 3] +
likelihood.matrix[, 4] +
likelihood.matrix[, 5] +
likelihood.matrix[, 8] +
likelihood.matrix[, 11]
likelihood.matrix[, 13] <- log(likelihood.matrix[, 12])
likelihood.matrix[, 14] <- likelihood.matrix[, 13] * matrix01[, 'fq']
likelihood.matrix[, 15] <- sum(likelihood.matrix[, 14])
# return negative log-likelihood calculated in likelihood matrix
neg.log.likelihood <- -(likelihood.matrix[1, 15])
# separately uncheck the lines below & re-run function to see the
# calculations for each component contributing to likelihood calculation
# likelihood.matrix
# matrix01
# survival.functions
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.