library("testthat")
# The code below is taken verbatim (with permission) from the SCCS package by
# Yonas Ghebremichael-Weldeselassie, Heather Whitaker, and Paddy Farrington.
### Functions used to calculate the weights to be used as offset in the new model ###
#--------------------------------------------------------#
# Weight function for Exponential- Weibull (Age) #
# mixture Model #
#--------------------------------------------------------#
# p<-p_ewad2
wsmall_ewad2 <- function(t, p, present, astart, aend, Dmatrix) {
thetaA <- p[which.max(Dmatrix)]
thetaB <- p[(length(Dmatrix)) + (which.max(Dmatrix))] + p[2 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(astart))
eta <- p[3 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[4 * (length(Dmatrix)) + (which.max(Dmatrix))] * t
gamma0 <- p[5 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[6 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(astart))
lamA <- (exp(-thetaA)) # 1/rho in the paper
lamB <- (exp(-thetaB)) # 1/mu
pi0 <- (exp(eta) / (1 + exp(eta))) # pi
nu0 <- (exp(gamma0)) # nu
val <- ((1 - present) * log(pi0 * lamA * exp(-lamA * (aend - t)) +
(1 - pi0) * nu0 * lamB * ((aend * lamB)^(nu0 - 1)) * exp(-((aend * lamB)^nu0 - (t * lamB)^nu0))) +
present * log(pi0 * exp(-lamA * (aend - t)) +
(1 - pi0) * exp(-((aend * lamB)^nu0 - (t * lamB)^nu0))))
# print(paste(t, exp(val)))
exp(val)
}
#--------------------------------------------------------------#
# Weight function for Exponential- Weibull (Interval) #
# mixture Model #
#--------------------------------------------------------------#
# p<-p_ewid2
wsmall_ewid2 <- function(t, p, present, aend, Dmatrix) {
thetaA <- p[which.max(Dmatrix)]
thetaB <- p[(length(Dmatrix)) + (which.max(Dmatrix))] + p[2 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(t))
eta <- p[3 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[4 * (length(Dmatrix)) + (which.max(Dmatrix))] * t
gamma0 <- p[5 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[6 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(t))
lamA <- exp(-thetaA) # 1/rho in the paper
lamB <- exp(-thetaB) # 1/mu
pi0 <- exp(eta) / (1 + exp(eta)) # pi
nu0 <- exp(gamma0) # nu
int <- aend - t
val <- ((1 - present) * log(pi0 * lamA * exp(-lamA * int) +
(1 - pi0) * nu0 * lamB * ((int * lamB)^(nu0 - 1)) * exp(-((int * lamB)^nu0))) +
present * log(pi0 * exp(-lamA * int) +
(1 - pi0) * exp(-((int * lamB)^nu0))))
exp(val)
}
#--------------------------------------------------------#
# Weight function for Exponential- Gamma (Age) #
# mixture Model #
#--------------------------------------------------------#
# p<-p_egad2
wsmall_egad2 <- function(t, p, present, astart, aend, Dmatrix) {
thetaA <- p[which.max(Dmatrix)]
thetaB <- p[(length(Dmatrix)) + (which.max(Dmatrix))] + p[2 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(astart))
eta <- p[3 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[4 * (length(Dmatrix)) + (which.max(Dmatrix))] * t
gamma0 <- p[5 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[6 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(astart))
lamA <- exp(-thetaA) # 1/rho in the paper
lamB <- exp(-thetaB) # 1/mu
pi0 <- exp(eta) / (1 + exp(eta)) # pi
nu0 <- exp(gamma0) # nu
rate0 <- nu0 * lamB
# val<- ((1-present)*log(pi0*lamA*exp(-lamA*(aend-t))+
# (1-pi0)*dgamma(aend,shape=nu0,rate=rate0)/pgamma(t,shape=nu0,rate=rate0,lower.tail=F)) +
# present*log(pi0*exp(-lamA*(aend-t))+
# (1-pi0)*pgamma(aend,shape=nu0,rate=rate0,lower.tail=F)/pgamma(t,shape=nu0,rate=rate0,lower.tail=F)))
val <- ((1 - present) * log(pi0 * lamA * exp(-lamA * (aend - t)) +
(1 - pi0) * dgamma(aend, shape = nu0, rate = rate0) / ifelse(pgamma(t, shape = nu0, rate = rate0, lower.tail = F) == 0, 0.000000001, pgamma(t, shape = nu0, rate = rate0, lower.tail = F))) +
present * log(pi0 * exp(-lamA * (aend - t)) +
(1 - pi0) * pgamma(aend, shape = nu0, rate = rate0, lower.tail = F) / ifelse(pgamma(t, shape = nu0, rate = rate0, lower.tail = F) == 0, 0.000000001, pgamma(t, shape = nu0, rate = rate0, lower.tail = F))))
exp(val)
}
#--------------------------------------------------------#
# Weight function for Exponential- Gamma (Interval)#
# mixture Model #
#--------------------------------------------------------#
# p<-p_egid2
wsmall_egid2 <- function(t, p, present, astart, aend, Dmatrix) {
thetaA <- p[which.max(Dmatrix)]
thetaB <- p[(length(Dmatrix)) + (which.max(Dmatrix))] + p[2 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(t))
eta <- p[3 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[4 * (length(Dmatrix)) + (which.max(Dmatrix))] * t
gamma0 <- p[5 * (length(Dmatrix)) + (which.max(Dmatrix))] + p[6 * (length(Dmatrix)) + (which.max(Dmatrix))] * (log(t))
lamA <- exp(-thetaA) # 1/rho in the paper
lamB <- exp(-thetaB) # 1/mu
pi0 <- exp(eta) / (1 + exp(eta)) # pi
nu0 <- exp(gamma0) # nu
rate0 <- nu0 * lamB
int <- aend - t
val <- ((1 - present) * log(pi0 * lamA * exp(-lamA * int) +
(1 - pi0) * dgamma(int, shape = nu0, rate = rate0)) +
present * log(pi0 * exp(-lamA * int) +
(1 - pi0) * pgamma(int, shape = nu0, rate = rate0, lower.tail = F)))
exp(val)
}
test_that("Weight functions match those in SCCS package", {
p <- c(0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1)
present <- 1
astart <- 1
aend <- 10
start <- 1
end <- 2
Dmatrix <- c(1)
w1 <- SelfControlledCaseSeries:::testEwad(p, present, astart, aend, start, end)
w2 <- integrate(wsmall_ewad2, lower = start, upper = end, p = p, present = present, astart = astart, aend = aend, Dmatrix = Dmatrix)$value
expect_equal(w1, w2, tolerance = 1E-6)
w1 <- SelfControlledCaseSeries:::testEwid(p, present, astart, aend, start, end)
w2 <- integrate(wsmall_ewid2, lower = start, upper = end, p = p, present = present, aend = aend, Dmatrix = Dmatrix)$value
expect_equal(w1, w2, tolerance = 1E-6)
w1 <- SelfControlledCaseSeries:::testEgad(p, present, astart, aend, start, end)
w2 <- integrate(wsmall_egad2, lower = start, upper = end, p = p, present = present, astart = astart, aend = aend, Dmatrix = Dmatrix)$value
expect_equal(w1, w2, tolerance = 1E-6)
w1 <- SelfControlledCaseSeries:::testEgid(p, present, astart, aend, start, end)
w2 <- integrate(wsmall_egid2, lower = start, upper = end, p = p, present = present, aend = aend, Dmatrix = Dmatrix)$value
expect_equal(w1, w2, tolerance = 1E-6)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.