# This helper funciton returns just the coefficient from the WeMix results object.
#' @export
coef.WeMixResults <- function(object, ...) {
object$coef
}
# This helper function prints the coefficient from the WeMix results object.
#' @export
print.WeMixResults <- function(x, ...) {
print(x$coef)
}
# This helper funciton creates a coefficient matrix from WeMix results.
#' @importFrom stats cov2cor
#' @export
summary.WeMixResults <- function(object, ...) {
x0 <- object
#for adaptive quad models
if(object$is_adaptive){
VC <- makeSandwich(object) # calculates the sandwhich estimator of variance
object$vc <- VC$VC #sets variances returned by model to variances calculated by sandwich estimator
se <- VC$se[1:length(x0$coef)]
re_se <- VC$se[-(1:length(x0$coef))] #not used currently
x0$varVC <- list(NULL,VC$VC[names(x0$vars),names(x0$vars)]) #This is hard coded for 2 levels, with a NULL because the are no refs at level 1
# build the var mat output, first adding SEvcov
varDF <- x0$varDF
varDF$SEvcov <- NA
names(re_se) <- gsub(":",".",names(re_se),fixed=TRUE) # change : to .
for(i in 1:nrow(varDF)) {
if(varDF$fullGroup[i] %in% names(re_se)) {
varDF$SEvcov[i] <- re_se[varDF$fullGroup[i]]
}
}
} else{
se <- object$SE
re_se <- rep(0,length(x0$vars))
# build the var mat output from this
varDF <- x0$varDF
}
object$coef <- as.matrix(data.frame(Estimate=x0$coef, "Std. Error"=se, "t value"=x0$coef/se))
# make.names breaks/fixes colnames with spaces, un-break/-fix them.
colnames(object$coef)[2:3] <- c("Std. Error", "t value")
rownames(object$coef) <- names(x0$coef)
# build the var mat output
varsmat <- varDF[is.na(x0$varDF$var2),c("level","grp", "var1","vcov","SEvcov")]
varsmat$st <- sqrt(varsmat$vcov)
# add variance estimates
colnames(varsmat) <- c("Level", "Group", "Name", "Variance", "Std. Error", "Std.Dev.")
for(li in 2:x0$levels) {
vc <- as.matrix(x0$varVC[[li]])
cr <- cov2cor(vc)
if(ncol(vc)>1) {
for(i in 2:ncol(cr)) {
for(j in 1:(i-1)){
varsmat[varsmat$Level==li & varsmat$Name==rownames(cr)[i],paste0("Corr",j)] <- cr[i,j]
}
}
}
}
object$varsmat <- varsmat
object$vars <- cbind(x0$vars, re_se, sqrt(x0$vars))
class(object) <- "summaryWeMixResults"
return(object)
}
# This helper funciton creates a coefficient matrix from WeMix results.
#' @export
#' @importFrom stats printCoefmat
print.summaryWeMixResults <- function(x, ...) {
cat("Call:\n")
print(x$call)
cat("\nVariance terms:\n")
varsmat <- x$varsmat
i <- 1
while(paste0("Corr",i) %in% colnames(varsmat)) {
# round correlations to two digits
varsmat[[paste0("Corr",i)]] <- as.character(round(varsmat[[paste0("Corr",i)]],2))
i <- i + 1
}
print(varsmat, na.print="", row.names=FALSE, digits=3)
cat("Groups:\n")
print(x$ngroups)
cat("\nFixed Effects:\n")
printCoefmat(x$coef, ...)
cat("\nlnl=",format(x$lnl,nsmall=2),"\n")
if(length(x$ICC) > 0) {
cat("Intraclass Correlation=",format(x$ICC, nsmall=3, digits=3),"\n")
}
}
#' This function calculates robust standard errors using the sandwich estimator of the standard errors
#' following Rabe-Hesketh & Skrondal 2006. For linear models, robust standard errors are returned from the main estimation function.
#' This function is only used for post-hoc estimation for non-linear models.
#' @importFrom numDeriv jacobian
#' @keywords internal
makeSandwich <- function(fittedModel) {
#make same estimator based on
par <- c(fittedModel$coef, fittedModel$vars)
FisherInfInverse <- solve(fittedModel$invHessian)
if (fittedModel$robustSE == TRUE){
L <- jacobian(fittedModel$lnlf, par, top=FALSE)
nGroups <- nrow(L)
nParams <- length(par)
J <- matrix(0,ncol=nParams, nrow=nParams)
for (i in 1:nGroups){
J <- J + L[i,] %*% t(L[i,])
}
J <- (nGroups/(nGroups-1))*J
SE <- FisherInfInverse%*%J%*%FisherInfInverse
colnames(SE) <- rownames(SE)<-names(par)
se <- sqrt(diag(SE))
#return NA for obs with variance below threshold
se[which(log(fittedModel$vars) < -3.6) + length(fittedModel$coef) ] <- NA
diag(SE) <- se^2
names(se) <- colnames(SE)
} else {
SE <- -1* FisherInfInverse
colnames(SE) <- rownames(SE)<-names(par)
se <- sqrt(diag(SE))
names(se) <- names(par)
}
return(list(VC=SE,se=se))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.