Skip to content

Commit

Permalink
v0.4.5 Fixed Kenward-Roger! Couple other bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@samuel-watson
samuel-watson committed Jul 25, 2023
1 parent 3d1532d commit cc4ddb9
Show file tree
Hide file tree
Showing 6 changed files with 170 additions and 78 deletions.
109 changes: 67 additions & 42 deletions R/R6Model.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -718,12 +718,16 @@ Model <- R6::R6Class("Model",
#'between iterations at which to stop the algorithm.
#'@param max.iter Integer. The maximum number of iterations of the MCML algorithm.
#'@param se String. Type of standard error to return. Options are "gls" for GLS standard errors (the default),
#' "robust" for Huber robust sandwich estimator, and "kr" for Kenward-Roger bias corrected standard errors.
#' "robust" for Huber robust sandwich estimator, "kr" for Kenward-Roger bias corrected standard errors, "bw" to use
#' GLS standard errors with a between-within correction to the degrees of freedom, "bwrobust" to use robust
#' standard errors with between-within correction to the degrees of freedom.
#'@param sparse Logical indicating whether to use sparse matrix methods
#'@param usestan Logical whether to use Stan (through the package `cmdstanr`) for the MCMC sampling. If FALSE then
#'the internal Hamiltonian Monte Carlo sampler will be used instead. We recommend Stan over the internal sampler as
#'it generally produces a larger number of effective samplers per unit time, especially for more complex
#'covariance functions.
#'@param se.theta Logical. Whether to calculate the standard errors for the covariance parameters. This step is a slow part
#' of the calculation, so can be disabled if required in larger models. Has no effect for Kenward-Roger standard errors.
#'@return A `mcml` object
#'@seealso \link[glmmrBase]{Model}, \link[glmmrBase]{Covariance}, \link[glmmrBase]{MeanFunction}
#'@examples
Expand Down Expand Up @@ -763,7 +767,8 @@ Model <- R6::R6Class("Model",
max.iter = 30,
se = "gls",
sparse = FALSE,
usestan = TRUE){
usestan = TRUE,
se.theta = TRUE){
private$verify_data(y)
private$set_y(y)
Model__use_attenuation(private$ptr,private$attenuate_parameters)
Expand Down Expand Up @@ -859,7 +864,6 @@ Model <- R6::R6Class("Model",
if(trace==2)t3 <- Sys.time()
if(trace==2)cat("\nModel fitting took: ",t3-t2,"s")
if(verbose){
#cat("\ntheta:",theta[all_pars])
cat("\nBeta: ", beta_new)
cat("\nTheta: ", theta_new)
if(var_par_family)cat("\nSigma: ",var_par_new)
Expand All @@ -882,12 +886,20 @@ Model <- R6::R6Class("Model",
if(verbose)cat("\n\nCalculating standard errors...\n")
self$var_par <- var_par_new
u <- Model__u(private$ptr, TRUE)
if(se == "gls"){
if(se == "gls" || se == "bw"){
M <- Matrix::solve(Model__obs_information_matrix(private$ptr))[1:length(beta),1:length(beta)]
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
} else if(se == "robust"){
if(se.theta){
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
} else {
SE_theta <- rep(NA, ncovpar)
}
} else if(se == "robust" || se == "bwrobust"){
M <- Model__sandwich(private$ptr)
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
if(se.theta){
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
} else {
SE_theta <- rep(NA, ncovpar)
}
} else if(se == "kr"){
Mout <- Model__kenward_roger(private$ptr)
M <- Mout[[1]]
Expand All @@ -912,21 +924,24 @@ Model <- R6::R6Class("Model",
p = NA,
lower = NA,
upper = NA)
dof <- rep(self$n(),length(beta))
if(se == "kr"){
for(i in 1:length(all_pars_new)){
for(i in 1:length(beta)){
if(!is.na(res$SE[i])){
bvar <- 1/(res$SE[i]^2)
ef <- 1 + bvar^2*Mout$b
vf <- 1 + bvar^2*Mout$a+6*bvar^2*Mout$b
rho <- vf/(2*ef^2)
m <- 4 + 3/(rho - 1)
lambda <- m/(ef*(m-2))
res$t[i] <- (res$est[i]/res$SE[i])*sqrt(lambda)
res$p[i] <- 2*(1-stats::pt(abs(res$t[i]),m))
res$lower[i] <- res$est - qt(1-0.05/2,m)*res$SE[i]
res$upper[i] <- res$est + qt(1-0.05/2,m)*res$SE[i]
res$t[i] <- (res$est[i]/res$SE[i])#*sqrt(lambda)
res$p[i] <- 2*(1-stats::pt(abs(res$t[i]),Mout$dof[i],lower.tail=FALSE))
res$lower[i] <- res$est[i] - stats::qt(0.975,Mout$dof[i],lower.tail=FALSE)*res$SE[i]
res$upper[i] <- res$est[i] + stats::qt(0.975,Mout$dof[i],lower.tail=FALSE)*res$SE[i]
}
dof[i] <- Mout$dof[i]
}
} else if(se=="bw" || se == "bwrobust"){
res$t <- res$est/res$SE
bwdof <- sum(repar_table$count) - length(beta)
res$p <- 2*(1-stats::pt(abs(res$t),bwdof,lower.tail=FALSE))
res$lower <- res$est - qt(1-0.05/2,bwdof,lower.tail=FALSE)*res$SE
res$upper <- res$est + qt(1-0.05/2,bwdof,lower.tail=FALSE)*res$SE
dof <- rep(bwdof,length(beta))
} else {
res$t <- res$est/res$SE
res$p <- 2*(1-stats::pnorm(abs(res$t)))
Expand Down Expand Up @@ -958,6 +973,7 @@ Model <- R6::R6Class("Model",
link = self$family[[2]],
re.samps = u,
iter = iter,
dof = dof,
P = length(self$mean$parameters),
Q = length(self$covariance$parameters),
var_par_family = var_par_family)
Expand All @@ -981,9 +997,13 @@ Model <- R6::R6Class("Model",
#'@param method String. Either "nloptim" for non-linear optimisation, or "nr" for Newton-Raphson (default) algorithm
#'@param verbose logical indicating whether to provide detailed algorithm feedback (default is TRUE).
#'@param se String. Type of standard error to return. Options are "gls" for GLS standard errors (the default),
#' "robust" for Huber robust sandwich estimator, and "kr" for Kenward-Roger bias corrected standard errors.
#' "robust" for Huber robust sandwich estimator, "kr" for Kenward-Roger bias corrected standard errors, "bw" to use
#' GLS standard errors with a between-within correction to the degrees of freedom, "bwrobust" to use robust
#' standard errors with between-within correction to the degrees of freedom.
#'@param max.iter Maximum number of algorithm iterations, default 20.
#'@param tol Maximum difference between successive iterations at which to terminate the algorithm
#'@param se.theta Logical. Whether to calculate the standard errors for the covariance parameters. This step is a slow part
#' of the calculation, so can be disabled if required in larger models. Has no effect for Kenward-Roger standard errors.
#'@return A `mcml` object
#' @seealso \link[glmmrBase]{Model}, \link[glmmrBase]{Covariance}, \link[glmmrBase]{MeanFunction}
#'@examples
Expand Down Expand Up @@ -1012,7 +1032,8 @@ Model <- R6::R6Class("Model",
verbose = FALSE,
se = "gls",
max.iter = 40,
tol = 1e-2){
tol = 1e-4,
se.theta = TRUE){
private$verify_data(y)
private$set_y(y)
Model__use_attenuation(private$ptr,private$attenuate_parameters)
Expand Down Expand Up @@ -1065,12 +1086,20 @@ Model <- R6::R6Class("Model",
self$var_par <- var_par_new
u <- Model__u(private$ptr,TRUE)
if(verbose)cat("\n\nCalculating standard errors...\n")
if(se == "gls"){
if(se == "gls" || se =="bw"){
M <- Matrix::solve(Model__obs_information_matrix(private$ptr))[1:length(beta),1:length(beta)]
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
} else if(se == "robust"){
if(se.theta){
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
} else {
SE_theta <- rep(NA, ncovpar)
}
} else if(se == "robust" || se == "bwrobust"){
M <- Model__sandwich(private$ptr)
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
if(se.theta){
SE_theta <- tryCatch(sqrt(diag(solve(Model__infomat_theta(private$ptr)))), error = rep(NA, ncovpar))
} else {
SE_theta <- rep(NA, ncovpar)
}
} else if(se == "kr"){
Mout <- Model__kenward_roger(private$ptr)
M <- Mout[[1]]
Expand All @@ -1096,29 +1125,24 @@ Model <- R6::R6Class("Model",
lower = NA,
upper = NA)

dof <- rep(self$n(),length(beta))
if(se == "kr"){
for(i in 1:length(beta)){
if(!is.na(res$SE[i])){
bvar <- 1/(res$SE[i]^4)
a1 <- bvar*Mout$a
a2 <- bvar*Mout$b
ef <- 1/(1 - a2)
g <- (2*a1 - 5*a2)/(3*a2)
c1 <- g/(3+2*(1-g))
c2 <- (1-g)/(3+2*(1-g))
c3 <- (3-g)/(3+2*(1-g))
b <- 0.5*(a1 + 6*a2)
vf <- 2*((1+c1*b)/((1-c2*b)^2*(1-c3*b)))
rho <- vf/(2*ef^2)
m <- 4 + 3/(rho - 1)
lambda <- m/(ef*(m-2))
print(m)
res$t[i] <- (res$est[i]/res$SE[i])*sqrt(lambda)
res$p[i] <- 2*(1-stats::pt(abs(res$t[i]),m))
res$lower[i] <- res$est[i] - stats::qt(0.975,m)*res$SE[i]
res$upper[i] <- res$est[i] + stats::qt(0.975,m)*res$SE[i]
res$t[i] <- (res$est[i]/res$SE[i])#*sqrt(lambda)
res$p[i] <- 2*(1-stats::pt(abs(res$t[i]),Mout$dof[i],lower.tail=FALSE))
res$lower[i] <- res$est[i] - stats::qt(0.975,Mout$dof[i],lower.tail=FALSE)*res$SE[i]
res$upper[i] <- res$est[i] + stats::qt(0.975,Mout$dof[i],lower.tail=FALSE)*res$SE[i]
}
dof[i] <- Mout$dof[i]
}
} else if(se=="bw" || se == "bwrobust"){
res$t <- res$est/res$SE
bwdof <- sum(repar_table$count) - length(beta)
res$p <- 2*(1-stats::pt(abs(res$t),bwdof,lower.tail=FALSE))
res$lower <- res$est - qt(1-0.05/2,bwdof,lower.tail=FALSE)*res$SE
res$upper <- res$est + qt(1-0.05/2,bwdof,lower.tail=FALSE)*res$SE
dof <- rep(bwdof,length(beta))
} else {
res$t <- res$est/res$SE
res$p <- 2*(1-stats::pnorm(abs(res$t)))
Expand Down Expand Up @@ -1150,6 +1174,7 @@ Model <- R6::R6Class("Model",
link = self$family[[2]],
re.samps = u,
iter = iter,
dof = dof,
P = length(self$mean$parameters),
Q = length(self$covariance$parameters),
var_par_family = var_par_family)
Expand Down
26 changes: 16 additions & 10 deletions R/printfunctions.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,14 +29,23 @@ print.mcml <- function(x, ...){

cat("\nFixed effects formula :",x$mean_form)
cat("\nCovariance function formula: ",x$cov_form)
cat("\nFamily: ",x$family,", Link function:",x$link,"\n")
cat("\nFamily: ",x$family,", Link function:",x$link)
setype <- switch(
x$se,
"gls" = "GLS",
"robust" = "Robust",
"kr" = "Kenward Roger",
"bw" = "GLS with between-within correction",
"bwrobust" = "Robust with between-within correction"
)
cat("\nStandard error: ",setype,"\n")
if(x$method%in%c("mcem","mcnr"))cat("\nNumber of Monte Carlo simulations per iteration: ",x$m," with tolerance ",x$tol,"\n\n")


dim1 <- dim(x$re.samps)[1]
pars <- x$coefficients[1:(length(x$coefficients$par)-dim1),2:7]
colnames(pars) <- c("Estimate","Std. Err.","z value","p value","2.5% CI","97.5% CI")
if(x$se == "robust") colnames(pars)[2] <- "Robust Std. Err."
if(x$se == "kr") colnames(pars)[2:3] <- c("Bias corr. Std. Err.", "t value")
if(x$se == "bw" || x$se == "bwrobust" || x$se == "kr")colnames(pars)[3] <- "t value"
rnames <- x$coefficients$par[1:(length(x$coefficients$par)-dim1)]
if(any(duplicated(rnames))){
did <- unique(rnames[duplicated(rnames)])
Expand All @@ -45,26 +54,23 @@ print.mcml <- function(x, ...){
}
}
rownames(pars) <- rnames
pars <- apply(pars,2,round,digits = digits)

total_vars <- x$P+x$Q
if(x$var_par_family)total_vars <- total_vars + 1
if(x$se %in% c("kr","bw","bwrobust"))pars$DoF <- c(x$dof, rep(NA,total_vars - x$P))
pars <- apply(pars,2,round,digits = digits)

cat("\nRandom effects: \n")
print(pars[(x$P+1):(total_vars),c(1,2)])

cat("\nFixed effects: \n")
print(pars[1:x$P,])


cat("\ncAIC: ",round(x$aic,digits))
cat("\nApproximate R-squared: Conditional: ",round(x$Rsq[1],digits)," Marginal: ",round(x$Rsq[2],digits))
cat("\nLog-likelihood: ",round(x$logl,digits))
if(!x$converged)cat("\nMCML ALGORITHM DID NOT CONVERGE!")
# #messages
# if(x$permutation)message("Permutation test used for one parameter, other SEs are not reported. SEs and Z values
# are approximate based on the p-value, and assume normality.")
#if(!x$hessian&!x$permutation)warning("Hessian was not positive definite, standard errors are approximate")
#if(!x$converged)warning("Algorithm did not converge")

return(invisible(pars))
}

Expand Down
17 changes: 17 additions & 0 deletions inst/include/glmmr/general.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,4 +173,21 @@ struct matrix_matrix{
};
};

struct kenward_data{
public:
MatrixXd vcov_beta;
MatrixXd vcov_theta;
VectorXd dof;
VectorXd lambda;
kenward_data(int n1, int m1, int n2, int m2): vcov_beta(n1,m1), vcov_theta(n2,m2), dof(n1), lambda(n1) {};
kenward_data(const kenward_data& x) : vcov_beta(x.vcov_beta), vcov_theta(x.vcov_theta), dof(x.dof), lambda(x.lambda) {};
kenward_data& operator=(kenward_data x){
vcov_beta = x.vcov_beta;
vcov_theta = x.vcov_theta;
dof = x.dof;
lambda = x.lambda;
return *this;
};
};

#endif
Loading

0 comments on commit cc4ddb9

Please sign in to comment.