Code appendix
#' Name, UUN
#' Add your own function definitions on this file.
#' neg_log_lik
#
#' @description Evaluate the negated log-likelihood for model A and B
#' @param beta A vector with the beta parameters
#' @param data A `data.frame` with the same variables as the `filament1` data set.
#' Must have columns `CAD_Weight` and `Actual_Weight`
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional
#' scaling error model
neg_log_lik <- function(beta, data, model){
mu <- beta[1] + beta[2]*data[["CAD_Weight"]]
# distinguish between the two models to find the particular standard deviation for the betas
if(model == "A") {
sigma <- sqrt(exp(beta[3] + beta[4]*data[["CAD_Weight"]]))
}else{
sigma <- sqrt(exp(beta[3])+exp(beta[4]) * (data[["CAD_Weight"]]^2))
}
- sum(dnorm(data[["Actual_Weight"]],
mean = mu,
sd=sigma,
log = TRUE))
}
#' filament_estimate
#
#' @description Estimate filament models with different variance structure
#' @param data A `data.frame` with the same variables as the `filament1` data set.
#' Must have columns `CAD_Weight` and `Actual_Weight`
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional
#' scaling error model
#' @return An estimation object suitable for use with [filament1_predict()]
filament1_estimate <- function(data, model) {
model <- match.arg(model, c("A", "B"))
if (model == "A") {
beta_start <- c(-0.1, 1.07, -2, 0.05)
} else {
beta_start <- c(-0.15, 1.07, -13.5, -6.5)
}
opt <- optim(beta_start,
neg_log_lik,
data = data,
model = model,
hessian = TRUE,
method = "Nelder-Mead",
control = list(maxit = 5000)
)
fit <- list(
model = model,
par = opt$par,
hessian = opt$hessian
)
class(fit) <- c("filament1_estimate", "list")
fit
}
#' filament1_aux_EV
#'
#' @description Evaluate the expectation and variance for model A and B
#' @param beta A vector with the beta parameters
#' @param data A `data.frame` containing the required predictors, including `CAD_Weight`
#' @param model Either "A" for a log-linear variance model, or "B" for a proportional
#' scaling error model
#' @param Sigma_beta : If not NULL, an estimate of the covariance matrix for
# the uncertainty of estimated betas
#' @return A list with four elements:
# E : E(y|beta,x)
# V : Var(y|beta,x)
# VE : Var(E(y|beta,x)|x) or NULL
# EV : E(Var(y|beta,x)|x) or NULL
filament1_aux_EV <- function(beta, data, model = c("A", "B"),
Sigma_beta = NULL) {
model <- match.arg(model)
if (model == "A") {
ZE.0 <- model.matrix( ~ 1 + CAD_Weight, data = data)
ZV.0 <- model.matrix( ~ 1 + CAD_Weight, data = data)
ZE = cbind(ZE.0, ZV.0 * 0)
ZV = cbind(ZE.0 * 0, ZV.0)
VE <- EV <- NULL
if (!is.null(Sigma_beta)) {
# E(Var(y|beta,x)|x)
EV <- exp(ZV %*% beta + rowSums(ZV * (ZV %*% Sigma_beta)) / 2)
# Var(E(y|beta,x)|x)
VE <- rowSums(ZE * (ZE %*% Sigma_beta))
}
out <- list(
E = ZE %*% beta,
V = exp(ZV %*% beta),
VE = VE,
EV = EV
)
} else {
ZE.0 <- model.matrix( ~ 1 + CAD_Weight, data = data)
ZV.0 <- model.matrix( ~ 1 + I(CAD_Weight^2), data = data)
ZE = cbind(ZE.0, ZV.0 * 0)
ZV = cbind(ZE.0 * 0, ZV.0)
VE <- EV <- NULL
if (!is.null(Sigma_beta)) {
# E(Var(y|beta,x)|x)
# (pmin: Ignore large Sigma_beta values)
EV <- ZV %*% exp(beta + pmin(0.5^2, diag(Sigma_beta)) / 2)
# Var(E(y|beta,x)|x)
VE <- rowSums(ZE * (ZE %*% Sigma_beta))
}
out <- list(
E = ZE %*% beta,
V = ZV %*% exp(beta),
VE = VE,
EV = EV
)
}
out
}