R/estimate_dispersion.R

In kschmaus/gammmbest: Method Of Moments Based Estimation for Generalized Additive Mixed Models

#' Estimate model dispersion.
#'
#' @param xz_svd_list A list of compact singular vector decomposition objects.
#' @param y_list A list of numeric vectors representing the response.
#' @param eta_list A list of numeric vectors encoding fixed and random effects.
#' @param family A description of the error distribution and link function to be
#' used in the model. See \code{\link{glm}}.
estimate_dispersion <- function(xz_svd_list, y_list, eta_list, family) {
	if (family$family == "gaussian") {
		return(estimate_dispersion_gaussian(xz_svd_list, y_list, eta_list))
	} else if (family$family == "binomial") {
		return(1)
	}

	stop("family not supported")
}

estimate_dispersion_gaussian <- function(xz_svd_list, y_list, eta_list) {
	U_list <- Map(getElement, xz_svd_list, "U")
	D_list <- Map(getElement, xz_svd_list, "D")

	residuals <- Map(gaussian_residual, U_list, D_list, y_list, eta_list)
	residuals <- unlist(residuals)
	ssr <- crossprod(residuals)

	total_rank <- map_reduce_sum(getElement, xz_svd_list, "rank")

	sum(residuals**2) / (length(residuals) - total_rank)
}

gaussian_residual <- function(U, D, y, eta) {
	y - U %*% D %*% eta
}