Nothing
R/generics.R
In hopit: Hierarchical Ordered Probit Models with Application to Reporting Heterogeneity
Defines functions
plot.hopitCP
summary.hopitHI
print.hopitHI
plot.hopitHI
plot.hopitDW
print.summary.hopitDW
summary.hopitDW
print.hopitDW
print.summary.hopit
summary.hopit
print.vcov.hopit
vcov.hopit
print.hopit
coef.hopit
Documented in
coef.hopit
plot.hopitCP
plot.hopitDW
plot.hopitHI
print.hopit
summary.hopit
vcov.hopit
#' Extracting the model coefficients
#'
#' Extract the model coefficients from the \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param aslist a logical indicating whether the model coefficients should be returned
#' as a list of three vectors
#' related to latent variables, threshold lambdas, and threshold gammas.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @usage \method{coef}{hopit}(object, aslist = FALSE, ...)
#' @keywords internal
#' @author Maciej J. Danko
#' @aliases coefficients.hopit
#' @method coef hopit
coef.hopit <- function(object, aslist = FALSE, ...)
if (aslist[1]) object$coef.ls else object$coef
#' Printing basic information about fitted \code{hopit} model
#'
#' Print a \code{hopit} model.
#' @param x a \code{hopit} object.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @usage \method{print}{hopit}(x, ...)
#' @method print hopit
#' @keywords internal
#' @author Maciej J. Danko
print.hopit<-function(x, ...){
cat(hopit_msg(65), deparse(x$latent.formula), fill = TRUE)
cat(hopit_msg(66), deparse(x$thresh.formula), fill = TRUE)
cat(hopit_msg(72), x$link, fill = TRUE)
cat(hopit_msg(73), x$N, fill = TRUE)
cat(hopit_msg(74), toString(levels(x$y_i)), fill = TRUE)
if (x$hasdisp) cat(hopit_msg(77), x$coef[length(x$coef)], fill = TRUE)
cat(hopit_msg(78))
print(x$coef.ls$latent.params)
cat(hopit_msg(79))
print(x$coef.ls$thresh.lambda)
if(length(x$coef.ls$thresh.gamma)){
cat(hopit_msg(80))
print(x$coef.ls$thresh.gamma)
}
#if(length(x$coef.ls$logSigma)){
if (x$hasdisp) cat(hopit_msg(82)) else cat(hopit_msg(81))
cat(exp(x$coef.ls$logSigma),'\n')
#}
invisible(NULL)
}
#' Variance-covariance matrix from the fitted model
#'
#' Returns the variance-covariance matrix of the main parameters of a fitted \code{hopit} model object.
#' @param object a \code{hopit} object.
#' @param robust.vcov a logical indicating whether to use the sandwich estimator to
#' calculate the variance-covariance matrix.
#' If a survey design is detected, then this option is ignored.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @usage \method{vcov}{hopit}(object, robust.vcov, ...)
#' @keywords internal
#' @author Maciej J. Danko
#' @method vcov hopit
vcov.hopit<-function(object, robust.vcov, ...){
z <- object$vcov
mirv <- missing(robust.vcov)
if ("try-error" %in% class(z)) stop(paste(hopit_msg(37),
attr(z,"condition"),sep=''),call.=NULL)
if (!length(z)) stop(hopit_msg(38),call.=NULL)
if (length(object$design)){
if (!mirv) {
if (robust.vcov) {
warning(call. = FALSE, hopit_msg(39))
}
}
robust.vcov <- FALSE
} else {
if (mirv) robust.vcov <- TRUE
if (length(object$weights)) divw <- object$weights else divw <- 1
#check how weights work here, they must be standardized.
if (robust.vcov) z <- (z %*% t(object$estfun) %*%
(object$estfun/divw) %*% (z))
}
attr(z, 'survey.design') <- (length(object$design) > 0L)
attr(z, 'robust.vcov') <- robust.vcov
class(z) <- 'vcov.hopit'
z
}
# Printing the variance-covariance matrix
#
# Print the variance-covariance matrix calculated by the \code{\link{vcov.hopit}}.
# @param x a \code{vcov.hopit} object
# @param digits see \code{\link{print.default}}.
# @param ... further arguments passed to or from other methods.
# @usage \method{print}{vcov.hopit}(x, digits = 3L, ...)
#' @noRd
#' @method print vcov.hopit
#' @keywords internal
#' @export
print.vcov.hopit <- function(x, digits = 3L, ...){
cat(hopit_msg(40))
print.default(x, digits = digits, ...)
if (attr(x, 'survey.design')) cat(hopit_msg(41))
if (!is.na(attr(x, 'robust.vcov')) && attr(x, 'robust.vcov'))
cat(hopit_msg(42))
invisible(NULL)
}
#' Calculate the model summary
#'
#' Summarize a \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param robust.se a logical indicating whether to use robust standard errors based
#' on the sandwich estimator.
#' If a survey design is detected, then this option is ignored.
#' @param control a list with control parameters.
#' See \code{\link{hopit.control}}.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @author Maciej J. Danko
#' @useDynLib hopit
#' @usage \method{summary}{hopit}(object, robust.se, ...)
#' @method summary hopit
#' @keywords internal
#' @importFrom Rcpp evalCpp
summary.hopit <- function(object, robust.se, ...){
if (missing(robust.se)) {
if (length(object$design)) robust.se <- FALSE else robust.se <- TRUE
}
varcov <- vcov.hopit(object, robust.se, ...)
robust.se <- attr(varcov, 'robust.vcov')
dvcov <- diag(varcov)
if (any(dvcov<0)) warning(hopit_msg(43),call.=NA)
SE <- suppressWarnings(sqrt(abs(dvcov)))
if (length(object$design)){
message(hopit_msg(44))
}
if ((!robust.se) && (any(is.na(SE))) && !(length(object$design)))
warning(call. = FALSE, hopit_msg(45))
if (length(object$coef) != length(SE)) stop(hopit_msg(46),call.=NULL)
tstat <- object$coef/SE
pvalue <- pstdnorm(-abs(tstat)) * 2L
table1 <- data.frame(Estimate = object$coef, 'Std. Error' = SE,
'z value' = tstat, 'Pr(>|z|)' = pvalue, check.names = FALSE)
tmp <- list(coef = table1, vcov = varcov, model = object,
robust.se = robust.se)
class(tmp) <- 'summary.hopit'
tmp
}
# Print an object calculated by \code{\link{summary.hopit}}
#
# @param x an object created with \code{\link{summary.hopit}}.
# @param ... further arguments passed to or from other methods.
#' @keywords internal
#' @export
# @usage \method{print}{summary.hopit}(x, ...)
#' @method print summary.hopit
# @author Maciej J. Danko
#' @noRd
print.summary.hopit <- function(x, ...){
model <- x$model
cat(hopit_msg(65), deparse(model$latent.formula), fill = TRUE)
cat(hopit_msg(66), deparse(model$thresh.formula), fill = TRUE)
cat(hopit_msg(72), model$link, fill = TRUE)
cat(hopit_msg(73), model$N, fill = TRUE)
cat(hopit_msg(74), toString(levels(model$y_i)), fill = TRUE)
if(x$robust.se) cat(hopit_msg(71))
cat('\n')
stats::printCoefmat(x = x$coef, P.values = TRUE,
has.Pvalue = TRUE, digits = 4L, dig.tst = 2L)
cat(hopit_msg(70), exp(model$coef.ls$logSigma), fill = TRUE)
cat(hopit_msg(75), model$LL, fill = TRUE)
cat(hopit_msg(76), model$deviance, fill = TRUE)
if (!length(model$design)) cat('AIC:', AIC.hopit(model), fill = TRUE)
cat('\n')
invisible(NULL)
}
#' @keywords internal
#' @export
#' @method print hopitDW
#' @noRd
print.hopitDW <- function(x, digits=4, quote=FALSE, ...){
x1 <- as.matrix(format(round(as.vector(x),digits),nsmall=digits, digits=digits,
scientific = FALSE))
dimnames(x1) <- dimnames(x)
print.default(x1, quote=quote, digits=digits, ...)
}
#' @keywords internal
#' @export
#' @method summary hopitDW
#' @noRd
summary.hopitDW<-function(object, ...){
x <- data.frame(cbind('Std. coef'=as.vector(object), attr(object,'Ptab')),
check.names = FALSE)
class(x) <- c('summary.hopitDW','data.frame')
attr(x,'legend') <- attr(object,'legend')
return(x)
}
#' @keywords internal
#' @export
#' @method print summary.hopitDW
#' @noRd
print.summary.hopitDW <- function(x, digits=4, show.coef.names=TRUE, ...){
if (show.coef.names[1]) {
g <- data.frame('Coefficient name' = x[,2],
'Standardized coeficient'=
format(round(x[,1],digits),nsmall=digits, digits=digits,
scientific = FALSE),
'Pr(>|z|)' = format(round(x[,3],digits),
nsmall=digits, digits=digits, scientific = FALSE),
x[,4],
fix.empty.names = FALSE, check.names = FALSE)
} else {
g <- data.frame('Standardized coeficient'=
format(round(x[,1],digits), nsmall=digits, digits=digits,
scientific = FALSE),
'Pr(>|z|)' = format(round(x[,3],digits),
nsmall=digits, digits=digits, scientific = FALSE),
x[,4],
fix.empty.names = FALSE, check.names = FALSE)
}
rownames(g)<-rownames(x)
print(g)
cat("---\nSignif. codes: ",attr(x,'legend'),'\n')
}
#' Plotting standardized coefficients
#' @keywords internal
#' @export
#' @method plot hopitDW
#' @param x a object generated by \code{\link{standardizeCoef}} function.
#' @param ordered a logical indicating whether to sort the disability weights.
#' @param show.signif show significance codes.
#' @param mar,oma graphic parameters, see \code{\link{par}}.
#' @param ylab,xlab,density,angle,col,las,... arguments passed to \code{\link{barplot}}.
#' @importFrom graphics text
#' @author Maciej J. Danko
plot.hopitDW <- function(x,
ordered = TRUE,
show.signif = TRUE,
mar = c(10, 4, 1.5, 1),
oma = c(0, 0, 0, 0),
ylab = "Disability weight",
xlab = "",
density = 20,
angle = 45,
col = 'orange',
las = 3, ...){
x <- summary.hopitDW(x)
if (ordered[1]) {
oz <- order(x[,1], decreasing = TRUE)
x <- x[oz,]
}
if (length(mar) || length(oma)) {
opar <- graphics::par(c("mar", "oma"))
if (length(mar)) graphics::par(mar = mar)
if (length(oma)) graphics::par(oma = oma)
}
rr <- graphics::barplot(x[,1], las = las, names.arg = rownames(x), ylab = ylab,
angle=angle, col=col, density=density, ...)
if (show.signif[1]) graphics::text(rr,x[,1],paste('',x[,4]),srt=90, xpd=NA, adj=0, col='darkred')
}
#' Plotting Latent Index
#' @keywords internal
#' @export
#' @method plot hopitHI
#' @param x a object generated by \code{\link{latentIndex}} function.
#' @param response X-axis plotting option; choose \code{'data'} for the raw responses and \code{'fitted'} for the responses reclassified by the model.
#' @param xlab a label of the x-axis.
#' @param ylab a label of the y-axis.
#' @param ... further parameters passed to the \code{\link{plot}} function.
#' @author Maciej J. Danko
plot.hopitHI<-function(x,
response = c('data','fitted'),
xlab = '', ylab = 'Latent index', ...){
mi <- attr(x, 'model.info')
response <- tolower(match.arg(response))
if (response=='data') YY <- mi$y_i else
if (response=='fitted') YY <- mi$Ey_i else
stop(hopit_msg(83),call.=NULL)
graphics::plot(YY, x ,las=3, ylab=ylab, xlab=xlab,...)
}
#' @keywords internal
#' @export
#' @method print hopitHI
#' @noRd
print.hopitHI <- function(x, digits = 2, quote=FALSE, ...){
print.default(format(round(as.vector(x), digits),
scientific=FALSE, digits=digits,
nsmall=digits), quote=quote, ...)
}
#' @keywords internal
#' @export
#' @method summary hopitHI
#' @noRd
summary.hopitHI <- function(object, breaks = seq(0,1,length.out = 21), ...){
z<-cut(object, include.lowest = FALSE, breaks=breaks)
z<-table(z)
dimnames(z)<-unname(dimnames(z))
z
}
#' Plotting Cut-Points
#' @keywords internal
#' @export
#' @method plot hopitCP
#' @param decreasing.levels a logical indicating whether self-reported health classes are ordered in decreasing order.
#' @param plotf a logical indicating whether to plot the results.
#' @param XLab,XLab.cex a label of the x axis and it's size.
#' @param YLab,YLab.cex a label of the y axis and it's size.
#' @param border.lwd,border.lty,border.col graphic parameters for vertical lines used to plot cut-points.
#' @param mar,oma graphic parameters, see \code{\link{par}}.
#' @param group.labels.type a position of the legend. One of \code{middle}, \code{border}, or \code{none}.
#' @param ... further plotting arguments.
#' @author Maciej J. Danko
plot.hopitCP<-function(x,
decreasing.levels=x$decreasing.levels,
mar=c(4,4,1,1),
oma=c(0,0,0,0),
XLab='Health index',
XLab.cex=1.1,
YLab='Counts',
YLab.cex=1.1,
border.col=2,
border.lty=2,
border.lwd=1.5,
group.labels.type=c('middle','border','none'),
...){
if (decreasing.levels[1]) dorev <- rev else dorev <- identity
lv <- dorev(as.character(levels(attr(x,'y_i'))))
Nm <- paste(lv[-length(lv)],lv[-1],sep=' | ')
Nm <- sapply(Nm, function(k) bquote(bold(.(k))))
group.labels.type<-tolower(group.labels.type[1])
if (group.labels.type %notin% c('middle','border','none'))
stop(hopit_msg(84),call.=NULL)
opar <- graphics::par(c('mar','oma'))
graphics::par(mar=mar, oma=oma)
z<-graphics::hist(x$h.index, 100,xlab='',ylab='' ,
main='', yaxs='i', col=grDevices::grey(0.4, alpha = 0.5),
border=grDevices::grey(0.4, alpha = 0.5))
if (group.labels.type == 'border') {
for (j in seq_along(Nm)) graphics::text(x=R1[j],y=(1.1*max(z$counts))/2,
labels=Nm[[j]],
srt=90,pos=2,offset=0.67,col=2)
} else if (group.labels.type == 'middle'){
R1 <- x$cutpoints
R11=-diff(c(0,R1,1))/2+c(R1,1)+graphics::strheight('S',units='figure')/2
for (j in seq_along(lv)) graphics::text(x=R11[j],
y=(3*1.1*max(z$counts))/4,
labels=lv[j],
srt=90,
pos=3,
offset=0.67,
col=2)
}
graphics::box()
graphics::abline(v=R1,lwd=border.lwd,col=border.col,lty=border.lty)
graphics::mtext(XLab, 1, cex=XLab.cex, line = 2.5)
graphics::mtext(YLab, 2, cex=YLab.cex, line = 2.5)
suppressWarnings(graphics::par(opar))
}
#' Plotting getLevels object
#' @keywords internal
#' @export
#' @method plot hopitLV
#' @param x aobject generated by \code{\link{getLevels}}.
#' @param mar,oma graphic parameters, see \code{\link{par}}.
#' @param YLab,YLab.cex a label for the y-axis and it's size.
#' @param legbg a legend background color. See \code{bg} parameter in \code{\link{legend}}.
#' @param legbty a legend box type. See \code{bty} parameter in \code{\link{legend}}.
#' @param las the style of axis labels, see \code{\link{par}}.
#' @param ... further plotting arguments.
#' @author Maciej J. Danko
plot.hopitLV<-function(x,
mar = c(7,2,1.5,0.5), oma = c(0,3,0,0),
YLab = 'Fraction [%]',
YLab.cex = 1.1,
legbg = grDevices::adjustcolor('white',alpha.f=0.4),
legbty = 'o',
las = 3,
...){
opar <- graphics::par(c('mar','oma','mfrow'))
graphics::par(mfrow=c(1,2))
graphics::par(mar=mar,oma=oma)
graphics::barplot(t(x$original),las=las,main='Original')
graphics::barplot(t(x$adjusted),las=las,main='Adjusted', legend.text=TRUE,
args.legend = list(x='center', box.col=NA,
bg=legbg, bty=legbty))
graphics::par(mfrow=c(1,1))
graphics::par(mar=mar,oma=rep(0,4))
graphics::mtext(YLab,2,cex=YLab.cex)
suppressWarnings(graphics::par(opar))
}
#' Extracting a log likelihood of the fitted model
#'
#' Extract the log likelihood of a \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param ... additional objects of the class \code{hopit}.
#' @keywords internal
#' @export
#' @usage \method{logLik}{hopit}(object, ...)
#' @author Maciej J. Danko
#' @method logLik hopit
logLik.hopit<-function(object, ...) {
objects <- list(object, ...)
tmp <- deparse(substitute(list(object, ...)))
ob.nam <- gsub(' ', '', strsplit(substring(tmp, 6L,
nchar(tmp) - 1L), ',', fixed = TRUE)[[1L]])
res <- sapply(objects,function(object) object$LL)
names(res) <- ob.nam
res
}
#' Extracting the Akaike Information Criterion from the fitted model
#'
#' Extract the Akaike Information Criterion (AIC) from a fitted \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param k a penalty per parameter to be used; the default k = 2 is the
#' classical AIC.
#' @param ... additional objects of the class \code{hopit}.
#' @keywords internal
#' @export
#' @usage \method{AIC}{hopit}(object, ..., k = 2L)
#' @method AIC hopit
#' @author Maciej J. Danko
AIC.hopit<-function(object, ..., k = 2L) {
objects <- list(object, ...)
tmp <- deparse(substitute(list(object, ...)))
ob.nam <- gsub(' ', '', strsplit(substring(tmp, 6L, nchar(tmp) - 1L), ',',
fixed = TRUE)[[1L]])
res <- sapply(objects,function(object)
if (!length(object$design)) object$AIC else
stop(hopit_msg(47), call=NULL))
names(res) <- ob.nam
res
}
#' Likelihood Ratio Test Tables
#'
#' Perform the likelihood ratio test(s) for two or more \code{hopit} objects.
#' @param object an object containing the results returned by a \code{hopit}.
#' @param ... an additional object(s) of the same type.
#' @param method the method of ordered model comparisons. Choose \code{"sequential"}
#' for 1-2, 2-3, 3-4, ... comparisons or
#' \code{"with.most.complex"} for 1-2, 1-3, 1-4, ... comparisons,
#' where 1 is the most complex model (the least complex for \code{"with.least.complex"}).
#' @param direction determine if the complexity of listed models is
#' \code{"increasing"} or \code{"decreasing"} (default).
# @keywords internal
#' @usage \method{anova}{hopit}(object, ..., method = c("sequential",
#' "with.most.complex", 'with.least.complex'),
#' direction = c("decreasing", "increasing"))
#' @method anova hopit
#' @export
#' @return a vector or a matrix with the results of the test(s).
#' @author Maciej J. Danko
#' @seealso \code{\link{print.lrt.hopit}},
#' \code{\link{lrt.hopit}}, \code{\link{hopit}}.
#' @examples
#' # DATA
#' data(healthsurvey)
#'
#' # the order of response levels decreases from the best health to
#' # the worst health; hence the hopit() parameter decreasing.levels
#' # is set to TRUE
#' levels(healthsurvey$health)
#'
#' # Example 1 ---------------------
#' \donttest{
#' # fitting two nested models
#' model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # a model with an interaction between hypertension and high_cholesterol
#' model2 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # a likelihood ratio test
#' lrt1 <- anova(model1, model2)
#' lrt1
#'
#' # print results in a shorter form
#' print(lrt1, short = TRUE)
#'
#' # or equivalently
#' lrt.hopit(model2, model1)
#' }
#' # Example 2 ---------------------
#' \donttest{
#' # fitting additional nested models
#' model3 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese * diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' model4 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese * diabetes + other_diseases,
#' thresh.formula = ~ sex * ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # sequential likelihood ratio tests
#' # model complexity increases so direction = "increasing"
#' anova(model1, model2, model3, model4,
#' direction = "increasing", method = "sequential")
#'
#' # likelihood ratio tests of the most complex model with the rest of the models
#' anova(model1, model2, model3, model4,
#' direction = "increasing", method = "with.most.complex")
#'
#' # likelihood ratio tests of the least complex model with the rest of the models
#' anova(model1, model2, model3, model4,
#' direction = "increasing", method = "with.least.complex")
#' }
anova.hopit<-function(object, ..., method = c('sequential', 'with.most.complex',
'with.least.complex'),
direction = c('decreasing', 'increasing')){
method <- tolower(match.arg(method))
direction <- tolower(match.arg(direction))
if (length(list(object, ...)) > 1L) {
objects <- list(object, ...)
tmp <- deparse(substitute(list(object, ...)))
ob.nam <- gsub(' ', '', strsplit(substring(tmp, 6L, nchar(tmp) - 1L), ',',
fixed = TRUE)[[1L]])
} else stop(hopit_msg(48))
if (length(objects) == 2L){
if(length(objects[[1L]]$coef)+objects[[1L]]$hasdisp >
length(objects[[2L]]$coef)+objects[[2L]]$hasdisp) {
return(lrt.hopit(objects[[1L]], objects[[2L]]))
} else {
return(lrt.hopit(objects[[2L]], objects[[1L]]))
}
} else {
out <- NULL
rna <- NULL
if (direction == 'increasing') {
objects <- rev(objects)
ob.nam <- rev(ob.nam)
} else if (direction != 'decreasing')
stop(call.=NULL, hopit_msg(50))
for (k in 1L : (length(objects) - 1L)) {
if (method == 'sequential'){
tmp <- lrt.hopit(objects[[k]], objects[[k + 1L]])
rna <- c(rna, paste(ob.nam[k], 'vs.', ob.nam[k + 1L], sep = ' '))
} else if (method == 'with.most.complex') {
tmp <- lrt.hopit(objects[[1L]], objects[[k + 1L]])
rna <- c(rna, paste(ob.nam[1L], 'vs.', ob.nam[k + 1L], sep = ' '))
} else if (method == 'with.least.complex') {
tmp <- lrt.hopit(objects[[k]], objects[[length(objects)]])
rna <- c(rna, paste( ob.nam[k], 'vs.',
ob.nam[length(objects)], sep = ' '))
}
out <- rbind(out, c('Chi^2' = tmp$chisq, df = tmp$df,
'Pr(>Chi^2)' = tmp$pval))
}
rownames(out) <- rna
if (direction == 'increasing') out <- out[dim(out)[1L] : 1L,]
}
out <- list(table = out, objets = objects, names = ob.nam, method = method)
class(out) <- 'anova.hopit'
out
}
# Print an object calculated by \code{\link{anova.hopit}}
#
# @param x an object generated by \code{\link{anova.hopit}}.
# @param ... further arguments passed to or from other methods.
#' @keywords internal
#' @export
# @author Maciej J. Danko
# @usage \method{print}{anova.hopit}(x, ...)
# @seealso \code{\link{anova.hopit}}, \code{\link{hopit}}.
#' @method print anova.hopit
#' @noRd
print.anova.hopit <- function(x, ...){
cat(hopit_msg(49), x$method, '"\n\n', sep = '')
stats::printCoefmat(x$table, signif.stars = TRUE, P.values = TRUE,
has.Pvalue = TRUE, digits = 5L, dig.tst = 3L, tst.ind = 1L)
invisible(NULL)
}
#' Likelihood ratio test for a pair of models
#'
#' @param full,nested models to be compared.
#' @keywords internal
#' @return a vector with the results of the test.
#' @export
#' @author Maciej J. Danko
#' @seealso \code{\link{print.lrt.hopit}}, \code{\link{anova.hopit}},
#' \code{\link{hopit}}.
#' @examples
#' \donttest{
#' # DATA
#' data(healthsurvey)
#'
#' # the order of response levels decreases from the best health to
#' # the worst health; hence the hopit() parameter decreasing.levels
#' # is set to TRUE
#' levels(healthsurvey$health)
#'
#' # Example 1 ---------------------
#'
#' # fitting two nested models
#' model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # model with an interaction between hypertension and high_cholesterol
#' model2 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # Likelihood ratio test
#' lrt1 <- lrt.hopit(full = model2, nested = model1)
#' lrt1
#'
#' # print the results in a shorter form
#' print(lrt1, short = TRUE)
#'
#' # equivalently
#' print(anova(model2, model1), short = TRUE)
#' }
lrt.hopit <- function(full, nested){
if (!identical(full$design, nested$design)) stop(hopit_msg(51),call. = NULL)
if (length(full$coef) + full$hasdisp <= length(nested$coef)+ nested$hasdisp)
stop(hopit_msg(52),call. = NULL)
if (abs(full$LL - nested$LL) < .Machine$double.eps^0.45) {
message(hopit_msg(53))
} else if (full$LL - nested$LL < -.Machine$double.eps^0.45){
warning(call. = FALSE, hopit_msg(54))
}
if (ncol(full$latent.mm) < ncol(nested$latent.mm)) {
cat(hopit_msg(64))
cat("--",hopit_msg(65), deparse(full$latent.formula), fill = TRUE)
cat(hopit_msg(67))
cat("--",hopit_msg(65), deparse(nested$latent.formula), fill = TRUE)
stop(hopit_msg(68))
}
if (ncol(full$thresh.mm) < ncol(nested$thresh.mm)) {
cat(hopit_msg(64))
cat("--",hopit_msg(66), deparse(full$thresh.formula), fill = TRUE)
cat(hopit_msg(67))
cat("--",hopit_msg(66), deparse(nested$thresh.formula), fill = TRUE)
stop(hopit_msg(69))
}
if ((full$hasdisp) < (nested$hasdisp)) stop(hopit_msg(55))
if ((ncol(full$latent.mm)) && (ncol(nested$latent.mm)))
if (!(all(colnames(nested$latent.mm) %in% colnames(full$latent.mm))))
warning(call. = FALSE, hopit_msg(56))
if ((ncol(full$thresh.mm)) && (ncol(nested$thresh.mm)))
if (!(all(colnames(nested$thresh.mm) %in% colnames(full$thresh.mm))))
warning(call. = FALSE, hopit_msg(57))
stat <- 2L*( logLik.hopit(full) - logLik.hopit(nested))
if (!length(full$design)) {
df.diff <- length(full$coef.ls$latent.params) -
length(nested$coef.ls$latent.params) +
length(full$coef.ls$thresh.lambda) -
length(nested$coef.ls$thresh.lambda) +
length(full$coef.ls$thresh.gamma) -
length(nested$coef.ls$thresh.gamma) +
(full$hasdisp) - (nested$hasdisp)
p <- 1L - stats::pchisq(stat, df.diff)
} else {
stop(hopit_msg(58), call=NULL)
}
z <- list(chisq = stat, df = df.diff, pval = p, full = full, nested = nested)
class(z) <- 'lrt.hopit'
z
}
#' Printing an object calculated by \code{\link{lrt.hopit}}
#'
#' @param x an object obtained from \code{\link{lrt.hopit}}.
#' @param short a logical indicating whether to show a shortened description.
#' @param ... further arguments passed to or from other methods.
#' @keywords internal
#' @export
#' @usage \method{print}{lrt.hopit}(x, short = FALSE, ...)
#' @author Maciej J. Danko
#' @method print lrt.hopit
#' @seealso \code{\link{lrt.hopit}}, \code{\link{anova.hopit}},
#' \code{\link{hopit}}.
print.lrt.hopit <- function(x, short = FALSE, ...){
if (!short[1]) {
cat(hopit_msg(64))
cat("--", hopit_msg(65), deparse(x$full$latent.formula), fill = TRUE)
cat("--",hopit_msg(66), deparse(x$full$thresh.formula), fill = TRUE)
cat("--",hopit_msg(70),x$full$hasdisp, fill=TRUE)
cat(hopit_msg(67))
cat("--",hopit_msg(65), deparse(x$nested$latent.formula), fill = TRUE)
cat("--",hopit_msg(66), deparse(x$nested$thresh.formula), fill = TRUE)
cat("--",hopit_msg(70),x$nested$hasdisp, fill=TRUE)
cat('\n')
}
cat(hopit_msg(59))
# if (length(x$df)) {
out <- t(as.matrix(c('Chi^2' = unname(x$chisq), df = unname(x$df),
'Pr(>Chi^2)' = unname(x$pval))))
# out2 <- NULL
# else {
# out <- t(as.matrix(c('Chi^2' = unname(x$chisq),
# 'Pr(>Chi^2)' = unname(x$pval))))
# out2 <- x$scalef
# }
row.names(out) <- ''
stats::printCoefmat(out, signif.stars = TRUE, P.values = TRUE,
has.Pvalue = TRUE, digits = 5L, dig.tst = 3L, tst.ind = 1L)
# if (length(out2)) print(paste(hopit_msg(63),out2))
invisible(NULL)
}
#' Calculate the log likelihood profile for the fitted \code{hopit} model
#'
#' @param fitted a \code{hopit} object (a fitted model).
#' @param scope a value (fraction) defining the plotting range for a coefficient.
#' The range is \code{c(coef \* (1-scope), coef \* (1+scope))}.
#' @param steps at how many equally spaced points the log likelihood
#' function is calculated for each coefficient.
#' @param ... unused now.
#' @export
#' @keywords internal
#' @author Maciej J. Danko
#' @usage \method{profile}{hopit}(fitted, ..., scope = 0.15, steps = 101)
#' @method profile hopit
#' @seealso \code{\link{plot.profile.hopit}}, \code{\link{print.profile.hopit}},
#' \code{\link{hopit}}
#' @examples
#' \donttest{
#' # DATA
#' data(healthsurvey)
#'
#' # the order of response levels decreases from the best health to
#' # the worst health; hence the hopit() parameter decreasing.levels
#' # is set to TRUE
#' levels(healthsurvey$health)
#'
#' # Example 1 ---------------------
#'
#' # fitting the model:
#' model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # check the fit using the profile function (at 51 points)
#' pr <- profile(model1, steps = 51)
#' print(pr, plotf = FALSE)
#'
#' # plot profile
#' plot(pr, relative = FALSE)
#'
#' # alternative plot
#' plot(pr, relative = TRUE)
#' }
profile.hopit<-function(fitted, ..., scope=0.15, steps=101){
steps <- floor(steps/2)*2+1
if (fitted$hasdisp) COEF <- c(fitted$coef, fitted$coef.ls$logSigma) else
COEF <- fitted$coef
sub <- function(x,y) if (x==1) c(y,COEF[2:length(COEF)]) else
if (x==length(COEF)) c(COEF[-length(COEF)],y) else
c(COEF[1:(x-1)],y,COEF[(x+1):length(COEF)])
lo <- COEF*(1-scope)
hi <- COEF*(1+scope)
GG <- function(x) sapply(seq(lo[x],hi[x],length.out=steps),function(y)
hopit_negLL(parameters=sub(x,y),fitted,negative = FALSE))
val <- sapply(seq_along(COEF), function(x) GG(x))
attr(val,'scope') <- scope
attr(val,'steps') <- steps
attr(val,'lo') <- lo
attr(val,'hi') <- hi
colnames(val) <- names(COEF)
class(val) <- c("profile.hopit", "profile")
val
}
#' Plot the log likelihood profile for a profile.hopit object
#'
#' Plot the method for a profile.hopit object.
#' @param x a \code{profile.hopit} object.
#' @param leg.cex a character expansion factor relative to the current
#' \code{par("cex")} (see \code{\link{legend}}).
#' @param leg.col a color used for the legend text.
#' @param ylim see \code{\link{plot}}.
#' @param relative a logical indicating whether \code{ylim} on each panel should be the
#' same (\code{TRUE}) or not (\code{FALSE}).
#' @param ... arguments to be passed to the \code{\link{plot}}() function (see
#' \code{\link{par}}).
#' @export
#' @keywords internal
#' @usage \method{plot}{profile.hopit}(x, ..., ylim = NULL, relative = FALSE,
#' leg.cex = 0.85, leg.col = 'blue4')
#' @author Maciej J. Danko
#' @method plot profile.hopit
#' @seealso \code{\link{profile.hopit}}, \code{\link{print.profile.hopit}},
#' \code{\link{hopit}}.
plot.profile.hopit<-function(x, ..., ylim = NULL, relative = FALSE,
leg.cex = 0.85, leg.col = 'blue4'){
z <- sqrt(ncol(x))
zy <- round(z)
zx <- ceiling(z)
spar <- graphics::par(c('mfrow','mar'))
graphics::par(mfrow=c(zx,zy),mar=c(0,0,0,0))
if (relative[1]) ylim <- range(x)
for (j in seq_len(ncol(x))) {
graphics::plot(x[,j],type='l',axes='F', ylim = ylim, ...)
graphics::abline(v=floor(nrow(x)/2)+1,col=2,lty=2)
graphics::legend('bottom',colnames(x)[j], bty='n',cex=leg.cex,
text.col=leg.col)
graphics::box()
}
suppressWarnings(graphics::par(spar))
}
#' Print method for a profile.hopit object
#'
#' @param x a \code{profile.hopit} object.
#' @param plotf a logical indicating whether to plot the profile.
#' @param ... arguments to be passed to the \code{plot}() function (see
#' \code{\link{plot.profile.hopit}}).
#' @export
#' @keywords internal
#' @usage \method{print}{profile.hopit}(x, ..., plotf = TRUE)
#' @author Maciej J. Danko
#' @method print profile.hopit
#' @seealso \code{\link{profile.hopit}}, \code{\link{plot.profile.hopit}},
#' \code{\link{hopit}}
print.profile.hopit<-function(x, ..., plotf = TRUE){
test <- apply(x,2,which.max)==floor(nrow(x)/2)+1
if(plotf[1]) plot.profile.hopit(x, ...)
if (any(!test)) {
message(hopit_msg(60))
message(paste(hopit_msg(61),paste(names(test)[!test],sep='',
collapse = ', ')))
} else {
cat(hopit_msg(62))
}
}
#' Extract the \code{Sigma} parameter from a \code{hopit} model
#'
#' Extract the \code{Sigma} parameter from a \code{hopit} model
#' @param model a fitted \code{hopit} model.
#' @usage \method{sigma}{hopit}(object, ...)
#' @keywords internal
#' @export
#' @author Maciej J. Danko
#' @method sigma hopit
sigma.hopit <- function(object, ...)
unname(exp(object$coef.ls$logSigma))
Try the hopit package in your browser
Any scripts or data that you put into this service are public.
hopit documentation built on Oct. 1, 2022, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.hopitCP summary.hopitHI print.hopitHI plot.hopitHI plot.hopitDW print.summary.hopitDW summary.hopitDW print.hopitDW print.summary.hopit summary.hopit print.vcov.hopit vcov.hopit print.hopit coef.hopit
Documented in coef.hopit plot.hopitCP plot.hopitDW plot.hopitHI print.hopit summary.hopit vcov.hopit
#' Extracting the model coefficients
#'
#' Extract the model coefficients from the \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param aslist a logical indicating whether the model coefficients should be returned
#' as a list of three vectors
#' related to latent variables, threshold lambdas, and threshold gammas.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @usage \method{coef}{hopit}(object, aslist = FALSE, ...)
#' @keywords internal
#' @author Maciej J. Danko
#' @aliases coefficients.hopit
#' @method coef hopit
coef.hopit <- function(object, aslist = FALSE, ...)
if (aslist[1]) object$coef.ls else object$coef
#' Printing basic information about fitted \code{hopit} model
#'
#' Print a \code{hopit} model.
#' @param x a \code{hopit} object.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @usage \method{print}{hopit}(x, ...)
#' @method print hopit
#' @keywords internal
#' @author Maciej J. Danko
print.hopit<-function(x, ...){
cat(hopit_msg(65), deparse(x$latent.formula), fill = TRUE)
cat(hopit_msg(66), deparse(x$thresh.formula), fill = TRUE)
cat(hopit_msg(72), x$link, fill = TRUE)
cat(hopit_msg(73), x$N, fill = TRUE)
cat(hopit_msg(74), toString(levels(x$y_i)), fill = TRUE)
if (x$hasdisp) cat(hopit_msg(77), x$coef[length(x$coef)], fill = TRUE)
cat(hopit_msg(78))
print(x$coef.ls$latent.params)
cat(hopit_msg(79))
print(x$coef.ls$thresh.lambda)
if(length(x$coef.ls$thresh.gamma)){
cat(hopit_msg(80))
print(x$coef.ls$thresh.gamma)
}
#if(length(x$coef.ls$logSigma)){
if (x$hasdisp) cat(hopit_msg(82)) else cat(hopit_msg(81))
cat(exp(x$coef.ls$logSigma),'\n')
#}
invisible(NULL)
}
#' Variance-covariance matrix from the fitted model
#'
#' Returns the variance-covariance matrix of the main parameters of a fitted \code{hopit} model object.
#' @param object a \code{hopit} object.
#' @param robust.vcov a logical indicating whether to use the sandwich estimator to
#' calculate the variance-covariance matrix.
#' If a survey design is detected, then this option is ignored.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @usage \method{vcov}{hopit}(object, robust.vcov, ...)
#' @keywords internal
#' @author Maciej J. Danko
#' @method vcov hopit
vcov.hopit<-function(object, robust.vcov, ...){
z <- object$vcov
mirv <- missing(robust.vcov)
if ("try-error" %in% class(z)) stop(paste(hopit_msg(37),
attr(z,"condition"),sep=''),call.=NULL)
if (!length(z)) stop(hopit_msg(38),call.=NULL)
if (length(object$design)){
if (!mirv) {
if (robust.vcov) {
warning(call. = FALSE, hopit_msg(39))
}
}
robust.vcov <- FALSE
} else {
if (mirv) robust.vcov <- TRUE
if (length(object$weights)) divw <- object$weights else divw <- 1
#check how weights work here, they must be standardized.
if (robust.vcov) z <- (z %*% t(object$estfun) %*%
(object$estfun/divw) %*% (z))
}
attr(z, 'survey.design') <- (length(object$design) > 0L)
attr(z, 'robust.vcov') <- robust.vcov
class(z) <- 'vcov.hopit'
z
}
# Printing the variance-covariance matrix
#
# Print the variance-covariance matrix calculated by the \code{\link{vcov.hopit}}.
# @param x a \code{vcov.hopit} object
# @param digits see \code{\link{print.default}}.
# @param ... further arguments passed to or from other methods.
# @usage \method{print}{vcov.hopit}(x, digits = 3L, ...)
#' @noRd
#' @method print vcov.hopit
#' @keywords internal
#' @export
print.vcov.hopit <- function(x, digits = 3L, ...){
cat(hopit_msg(40))
print.default(x, digits = digits, ...)
if (attr(x, 'survey.design')) cat(hopit_msg(41))
if (!is.na(attr(x, 'robust.vcov')) && attr(x, 'robust.vcov'))
cat(hopit_msg(42))
invisible(NULL)
}
#' Calculate the model summary
#'
#' Summarize a \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param robust.se a logical indicating whether to use robust standard errors based
#' on the sandwich estimator.
#' If a survey design is detected, then this option is ignored.
#' @param control a list with control parameters.
#' See \code{\link{hopit.control}}.
#' @param ... further arguments passed to or from other methods.
#' @export
#' @author Maciej J. Danko
#' @useDynLib hopit
#' @usage \method{summary}{hopit}(object, robust.se, ...)
#' @method summary hopit
#' @keywords internal
#' @importFrom Rcpp evalCpp
summary.hopit <- function(object, robust.se, ...){
if (missing(robust.se)) {
if (length(object$design)) robust.se <- FALSE else robust.se <- TRUE
}
varcov <- vcov.hopit(object, robust.se, ...)
robust.se <- attr(varcov, 'robust.vcov')
dvcov <- diag(varcov)
if (any(dvcov<0)) warning(hopit_msg(43),call.=NA)
SE <- suppressWarnings(sqrt(abs(dvcov)))
if (length(object$design)){
message(hopit_msg(44))
}
if ((!robust.se) && (any(is.na(SE))) && !(length(object$design)))
warning(call. = FALSE, hopit_msg(45))
if (length(object$coef) != length(SE)) stop(hopit_msg(46),call.=NULL)
tstat <- object$coef/SE
pvalue <- pstdnorm(-abs(tstat)) * 2L
table1 <- data.frame(Estimate = object$coef, 'Std. Error' = SE,
'z value' = tstat, 'Pr(>|z|)' = pvalue, check.names = FALSE)
tmp <- list(coef = table1, vcov = varcov, model = object,
robust.se = robust.se)
class(tmp) <- 'summary.hopit'
tmp
}
# Print an object calculated by \code{\link{summary.hopit}}
#
# @param x an object created with \code{\link{summary.hopit}}.
# @param ... further arguments passed to or from other methods.
#' @keywords internal
#' @export
# @usage \method{print}{summary.hopit}(x, ...)
#' @method print summary.hopit
# @author Maciej J. Danko
#' @noRd
print.summary.hopit <- function(x, ...){
model <- x$model
cat(hopit_msg(65), deparse(model$latent.formula), fill = TRUE)
cat(hopit_msg(66), deparse(model$thresh.formula), fill = TRUE)
cat(hopit_msg(72), model$link, fill = TRUE)
cat(hopit_msg(73), model$N, fill = TRUE)
cat(hopit_msg(74), toString(levels(model$y_i)), fill = TRUE)
if(x$robust.se) cat(hopit_msg(71))
cat('\n')
stats::printCoefmat(x = x$coef, P.values = TRUE,
has.Pvalue = TRUE, digits = 4L, dig.tst = 2L)
cat(hopit_msg(70), exp(model$coef.ls$logSigma), fill = TRUE)
cat(hopit_msg(75), model$LL, fill = TRUE)
cat(hopit_msg(76), model$deviance, fill = TRUE)
if (!length(model$design)) cat('AIC:', AIC.hopit(model), fill = TRUE)
cat('\n')
invisible(NULL)
}
#' @keywords internal
#' @export
#' @method print hopitDW
#' @noRd
print.hopitDW <- function(x, digits=4, quote=FALSE, ...){
x1 <- as.matrix(format(round(as.vector(x),digits),nsmall=digits, digits=digits,
scientific = FALSE))
dimnames(x1) <- dimnames(x)
print.default(x1, quote=quote, digits=digits, ...)
}
#' @keywords internal
#' @export
#' @method summary hopitDW
#' @noRd
summary.hopitDW<-function(object, ...){
x <- data.frame(cbind('Std. coef'=as.vector(object), attr(object,'Ptab')),
check.names = FALSE)
class(x) <- c('summary.hopitDW','data.frame')
attr(x,'legend') <- attr(object,'legend')
return(x)
}
#' @keywords internal
#' @export
#' @method print summary.hopitDW
#' @noRd
print.summary.hopitDW <- function(x, digits=4, show.coef.names=TRUE, ...){
if (show.coef.names[1]) {
g <- data.frame('Coefficient name' = x[,2],
'Standardized coeficient'=
format(round(x[,1],digits),nsmall=digits, digits=digits,
scientific = FALSE),
'Pr(>|z|)' = format(round(x[,3],digits),
nsmall=digits, digits=digits, scientific = FALSE),
x[,4],
fix.empty.names = FALSE, check.names = FALSE)
} else {
g <- data.frame('Standardized coeficient'=
format(round(x[,1],digits), nsmall=digits, digits=digits,
scientific = FALSE),
'Pr(>|z|)' = format(round(x[,3],digits),
nsmall=digits, digits=digits, scientific = FALSE),
x[,4],
fix.empty.names = FALSE, check.names = FALSE)
}
rownames(g)<-rownames(x)
print(g)
cat("---\nSignif. codes: ",attr(x,'legend'),'\n')
}
#' Plotting standardized coefficients
#' @keywords internal
#' @export
#' @method plot hopitDW
#' @param x a object generated by \code{\link{standardizeCoef}} function.
#' @param ordered a logical indicating whether to sort the disability weights.
#' @param show.signif show significance codes.
#' @param mar,oma graphic parameters, see \code{\link{par}}.
#' @param ylab,xlab,density,angle,col,las,... arguments passed to \code{\link{barplot}}.
#' @importFrom graphics text
#' @author Maciej J. Danko
plot.hopitDW <- function(x,
ordered = TRUE,
show.signif = TRUE,
mar = c(10, 4, 1.5, 1),
oma = c(0, 0, 0, 0),
ylab = "Disability weight",
xlab = "",
density = 20,
angle = 45,
col = 'orange',
las = 3, ...){
x <- summary.hopitDW(x)
if (ordered[1]) {
oz <- order(x[,1], decreasing = TRUE)
x <- x[oz,]
}
if (length(mar) || length(oma)) {
opar <- graphics::par(c("mar", "oma"))
if (length(mar)) graphics::par(mar = mar)
if (length(oma)) graphics::par(oma = oma)
}
rr <- graphics::barplot(x[,1], las = las, names.arg = rownames(x), ylab = ylab,
angle=angle, col=col, density=density, ...)
if (show.signif[1]) graphics::text(rr,x[,1],paste('',x[,4]),srt=90, xpd=NA, adj=0, col='darkred')
}
#' Plotting Latent Index
#' @keywords internal
#' @export
#' @method plot hopitHI
#' @param x a object generated by \code{\link{latentIndex}} function.
#' @param response X-axis plotting option; choose \code{'data'} for the raw responses and \code{'fitted'} for the responses reclassified by the model.
#' @param xlab a label of the x-axis.
#' @param ylab a label of the y-axis.
#' @param ... further parameters passed to the \code{\link{plot}} function.
#' @author Maciej J. Danko
plot.hopitHI<-function(x,
response = c('data','fitted'),
xlab = '', ylab = 'Latent index', ...){
mi <- attr(x, 'model.info')
response <- tolower(match.arg(response))
if (response=='data') YY <- mi$y_i else
if (response=='fitted') YY <- mi$Ey_i else
stop(hopit_msg(83),call.=NULL)
graphics::plot(YY, x ,las=3, ylab=ylab, xlab=xlab,...)
}
#' @keywords internal
#' @export
#' @method print hopitHI
#' @noRd
print.hopitHI <- function(x, digits = 2, quote=FALSE, ...){
print.default(format(round(as.vector(x), digits),
scientific=FALSE, digits=digits,
nsmall=digits), quote=quote, ...)
}
#' @keywords internal
#' @export
#' @method summary hopitHI
#' @noRd
summary.hopitHI <- function(object, breaks = seq(0,1,length.out = 21), ...){
z<-cut(object, include.lowest = FALSE, breaks=breaks)
z<-table(z)
dimnames(z)<-unname(dimnames(z))
z
}
#' Plotting Cut-Points
#' @keywords internal
#' @export
#' @method plot hopitCP
#' @param decreasing.levels a logical indicating whether self-reported health classes are ordered in decreasing order.
#' @param plotf a logical indicating whether to plot the results.
#' @param XLab,XLab.cex a label of the x axis and it's size.
#' @param YLab,YLab.cex a label of the y axis and it's size.
#' @param border.lwd,border.lty,border.col graphic parameters for vertical lines used to plot cut-points.
#' @param mar,oma graphic parameters, see \code{\link{par}}.
#' @param group.labels.type a position of the legend. One of \code{middle}, \code{border}, or \code{none}.
#' @param ... further plotting arguments.
#' @author Maciej J. Danko
plot.hopitCP<-function(x,
decreasing.levels=x$decreasing.levels,
mar=c(4,4,1,1),
oma=c(0,0,0,0),
XLab='Health index',
XLab.cex=1.1,
YLab='Counts',
YLab.cex=1.1,
border.col=2,
border.lty=2,
border.lwd=1.5,
group.labels.type=c('middle','border','none'),
...){
if (decreasing.levels[1]) dorev <- rev else dorev <- identity
lv <- dorev(as.character(levels(attr(x,'y_i'))))
Nm <- paste(lv[-length(lv)],lv[-1],sep=' | ')
Nm <- sapply(Nm, function(k) bquote(bold(.(k))))
group.labels.type<-tolower(group.labels.type[1])
if (group.labels.type %notin% c('middle','border','none'))
stop(hopit_msg(84),call.=NULL)
opar <- graphics::par(c('mar','oma'))
graphics::par(mar=mar, oma=oma)
z<-graphics::hist(x$h.index, 100,xlab='',ylab='' ,
main='', yaxs='i', col=grDevices::grey(0.4, alpha = 0.5),
border=grDevices::grey(0.4, alpha = 0.5))
if (group.labels.type == 'border') {
for (j in seq_along(Nm)) graphics::text(x=R1[j],y=(1.1*max(z$counts))/2,
labels=Nm[[j]],
srt=90,pos=2,offset=0.67,col=2)
} else if (group.labels.type == 'middle'){
R1 <- x$cutpoints
R11=-diff(c(0,R1,1))/2+c(R1,1)+graphics::strheight('S',units='figure')/2
for (j in seq_along(lv)) graphics::text(x=R11[j],
y=(3*1.1*max(z$counts))/4,
labels=lv[j],
srt=90,
pos=3,
offset=0.67,
col=2)
}
graphics::box()
graphics::abline(v=R1,lwd=border.lwd,col=border.col,lty=border.lty)
graphics::mtext(XLab, 1, cex=XLab.cex, line = 2.5)
graphics::mtext(YLab, 2, cex=YLab.cex, line = 2.5)
suppressWarnings(graphics::par(opar))
}
#' Plotting getLevels object
#' @keywords internal
#' @export
#' @method plot hopitLV
#' @param x aobject generated by \code{\link{getLevels}}.
#' @param mar,oma graphic parameters, see \code{\link{par}}.
#' @param YLab,YLab.cex a label for the y-axis and it's size.
#' @param legbg a legend background color. See \code{bg} parameter in \code{\link{legend}}.
#' @param legbty a legend box type. See \code{bty} parameter in \code{\link{legend}}.
#' @param las the style of axis labels, see \code{\link{par}}.
#' @param ... further plotting arguments.
#' @author Maciej J. Danko
plot.hopitLV<-function(x,
mar = c(7,2,1.5,0.5), oma = c(0,3,0,0),
YLab = 'Fraction [%]',
YLab.cex = 1.1,
legbg = grDevices::adjustcolor('white',alpha.f=0.4),
legbty = 'o',
las = 3,
...){
opar <- graphics::par(c('mar','oma','mfrow'))
graphics::par(mfrow=c(1,2))
graphics::par(mar=mar,oma=oma)
graphics::barplot(t(x$original),las=las,main='Original')
graphics::barplot(t(x$adjusted),las=las,main='Adjusted', legend.text=TRUE,
args.legend = list(x='center', box.col=NA,
bg=legbg, bty=legbty))
graphics::par(mfrow=c(1,1))
graphics::par(mar=mar,oma=rep(0,4))
graphics::mtext(YLab,2,cex=YLab.cex)
suppressWarnings(graphics::par(opar))
}
#' Extracting a log likelihood of the fitted model
#'
#' Extract the log likelihood of a \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param ... additional objects of the class \code{hopit}.
#' @keywords internal
#' @export
#' @usage \method{logLik}{hopit}(object, ...)
#' @author Maciej J. Danko
#' @method logLik hopit
logLik.hopit<-function(object, ...) {
objects <- list(object, ...)
tmp <- deparse(substitute(list(object, ...)))
ob.nam <- gsub(' ', '', strsplit(substring(tmp, 6L,
nchar(tmp) - 1L), ',', fixed = TRUE)[[1L]])
res <- sapply(objects,function(object) object$LL)
names(res) <- ob.nam
res
}
#' Extracting the Akaike Information Criterion from the fitted model
#'
#' Extract the Akaike Information Criterion (AIC) from a fitted \code{hopit} model.
#' @param object a \code{hopit} object.
#' @param k a penalty per parameter to be used; the default k = 2 is the
#' classical AIC.
#' @param ... additional objects of the class \code{hopit}.
#' @keywords internal
#' @export
#' @usage \method{AIC}{hopit}(object, ..., k = 2L)
#' @method AIC hopit
#' @author Maciej J. Danko
AIC.hopit<-function(object, ..., k = 2L) {
objects <- list(object, ...)
tmp <- deparse(substitute(list(object, ...)))
ob.nam <- gsub(' ', '', strsplit(substring(tmp, 6L, nchar(tmp) - 1L), ',',
fixed = TRUE)[[1L]])
res <- sapply(objects,function(object)
if (!length(object$design)) object$AIC else
stop(hopit_msg(47), call=NULL))
names(res) <- ob.nam
res
}
#' Likelihood Ratio Test Tables
#'
#' Perform the likelihood ratio test(s) for two or more \code{hopit} objects.
#' @param object an object containing the results returned by a \code{hopit}.
#' @param ... an additional object(s) of the same type.
#' @param method the method of ordered model comparisons. Choose \code{"sequential"}
#' for 1-2, 2-3, 3-4, ... comparisons or
#' \code{"with.most.complex"} for 1-2, 1-3, 1-4, ... comparisons,
#' where 1 is the most complex model (the least complex for \code{"with.least.complex"}).
#' @param direction determine if the complexity of listed models is
#' \code{"increasing"} or \code{"decreasing"} (default).
# @keywords internal
#' @usage \method{anova}{hopit}(object, ..., method = c("sequential",
#' "with.most.complex", 'with.least.complex'),
#' direction = c("decreasing", "increasing"))
#' @method anova hopit
#' @export
#' @return a vector or a matrix with the results of the test(s).
#' @author Maciej J. Danko
#' @seealso \code{\link{print.lrt.hopit}},
#' \code{\link{lrt.hopit}}, \code{\link{hopit}}.
#' @examples
#' # DATA
#' data(healthsurvey)
#'
#' # the order of response levels decreases from the best health to
#' # the worst health; hence the hopit() parameter decreasing.levels
#' # is set to TRUE
#' levels(healthsurvey$health)
#'
#' # Example 1 ---------------------
#' \donttest{
#' # fitting two nested models
#' model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # a model with an interaction between hypertension and high_cholesterol
#' model2 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # a likelihood ratio test
#' lrt1 <- anova(model1, model2)
#' lrt1
#'
#' # print results in a shorter form
#' print(lrt1, short = TRUE)
#'
#' # or equivalently
#' lrt.hopit(model2, model1)
#' }
#' # Example 2 ---------------------
#' \donttest{
#' # fitting additional nested models
#' model3 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese * diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' model4 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese * diabetes + other_diseases,
#' thresh.formula = ~ sex * ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # sequential likelihood ratio tests
#' # model complexity increases so direction = "increasing"
#' anova(model1, model2, model3, model4,
#' direction = "increasing", method = "sequential")
#'
#' # likelihood ratio tests of the most complex model with the rest of the models
#' anova(model1, model2, model3, model4,
#' direction = "increasing", method = "with.most.complex")
#'
#' # likelihood ratio tests of the least complex model with the rest of the models
#' anova(model1, model2, model3, model4,
#' direction = "increasing", method = "with.least.complex")
#' }
anova.hopit<-function(object, ..., method = c('sequential', 'with.most.complex',
'with.least.complex'),
direction = c('decreasing', 'increasing')){
method <- tolower(match.arg(method))
direction <- tolower(match.arg(direction))
if (length(list(object, ...)) > 1L) {
objects <- list(object, ...)
tmp <- deparse(substitute(list(object, ...)))
ob.nam <- gsub(' ', '', strsplit(substring(tmp, 6L, nchar(tmp) - 1L), ',',
fixed = TRUE)[[1L]])
} else stop(hopit_msg(48))
if (length(objects) == 2L){
if(length(objects[[1L]]$coef)+objects[[1L]]$hasdisp >
length(objects[[2L]]$coef)+objects[[2L]]$hasdisp) {
return(lrt.hopit(objects[[1L]], objects[[2L]]))
} else {
return(lrt.hopit(objects[[2L]], objects[[1L]]))
}
} else {
out <- NULL
rna <- NULL
if (direction == 'increasing') {
objects <- rev(objects)
ob.nam <- rev(ob.nam)
} else if (direction != 'decreasing')
stop(call.=NULL, hopit_msg(50))
for (k in 1L : (length(objects) - 1L)) {
if (method == 'sequential'){
tmp <- lrt.hopit(objects[[k]], objects[[k + 1L]])
rna <- c(rna, paste(ob.nam[k], 'vs.', ob.nam[k + 1L], sep = ' '))
} else if (method == 'with.most.complex') {
tmp <- lrt.hopit(objects[[1L]], objects[[k + 1L]])
rna <- c(rna, paste(ob.nam[1L], 'vs.', ob.nam[k + 1L], sep = ' '))
} else if (method == 'with.least.complex') {
tmp <- lrt.hopit(objects[[k]], objects[[length(objects)]])
rna <- c(rna, paste( ob.nam[k], 'vs.',
ob.nam[length(objects)], sep = ' '))
}
out <- rbind(out, c('Chi^2' = tmp$chisq, df = tmp$df,
'Pr(>Chi^2)' = tmp$pval))
}
rownames(out) <- rna
if (direction == 'increasing') out <- out[dim(out)[1L] : 1L,]
}
out <- list(table = out, objets = objects, names = ob.nam, method = method)
class(out) <- 'anova.hopit'
out
}
# Print an object calculated by \code{\link{anova.hopit}}
#
# @param x an object generated by \code{\link{anova.hopit}}.
# @param ... further arguments passed to or from other methods.
#' @keywords internal
#' @export
# @author Maciej J. Danko
# @usage \method{print}{anova.hopit}(x, ...)
# @seealso \code{\link{anova.hopit}}, \code{\link{hopit}}.
#' @method print anova.hopit
#' @noRd
print.anova.hopit <- function(x, ...){
cat(hopit_msg(49), x$method, '"\n\n', sep = '')
stats::printCoefmat(x$table, signif.stars = TRUE, P.values = TRUE,
has.Pvalue = TRUE, digits = 5L, dig.tst = 3L, tst.ind = 1L)
invisible(NULL)
}
#' Likelihood ratio test for a pair of models
#'
#' @param full,nested models to be compared.
#' @keywords internal
#' @return a vector with the results of the test.
#' @export
#' @author Maciej J. Danko
#' @seealso \code{\link{print.lrt.hopit}}, \code{\link{anova.hopit}},
#' \code{\link{hopit}}.
#' @examples
#' \donttest{
#' # DATA
#' data(healthsurvey)
#'
#' # the order of response levels decreases from the best health to
#' # the worst health; hence the hopit() parameter decreasing.levels
#' # is set to TRUE
#' levels(healthsurvey$health)
#'
#' # Example 1 ---------------------
#'
#' # fitting two nested models
#' model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # model with an interaction between hypertension and high_cholesterol
#' model2 <- hopit(latent.formula = health ~ hypertension * high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # Likelihood ratio test
#' lrt1 <- lrt.hopit(full = model2, nested = model1)
#' lrt1
#'
#' # print the results in a shorter form
#' print(lrt1, short = TRUE)
#'
#' # equivalently
#' print(anova(model2, model1), short = TRUE)
#' }
lrt.hopit <- function(full, nested){
if (!identical(full$design, nested$design)) stop(hopit_msg(51),call. = NULL)
if (length(full$coef) + full$hasdisp <= length(nested$coef)+ nested$hasdisp)
stop(hopit_msg(52),call. = NULL)
if (abs(full$LL - nested$LL) < .Machine$double.eps^0.45) {
message(hopit_msg(53))
} else if (full$LL - nested$LL < -.Machine$double.eps^0.45){
warning(call. = FALSE, hopit_msg(54))
}
if (ncol(full$latent.mm) < ncol(nested$latent.mm)) {
cat(hopit_msg(64))
cat("--",hopit_msg(65), deparse(full$latent.formula), fill = TRUE)
cat(hopit_msg(67))
cat("--",hopit_msg(65), deparse(nested$latent.formula), fill = TRUE)
stop(hopit_msg(68))
}
if (ncol(full$thresh.mm) < ncol(nested$thresh.mm)) {
cat(hopit_msg(64))
cat("--",hopit_msg(66), deparse(full$thresh.formula), fill = TRUE)
cat(hopit_msg(67))
cat("--",hopit_msg(66), deparse(nested$thresh.formula), fill = TRUE)
stop(hopit_msg(69))
}
if ((full$hasdisp) < (nested$hasdisp)) stop(hopit_msg(55))
if ((ncol(full$latent.mm)) && (ncol(nested$latent.mm)))
if (!(all(colnames(nested$latent.mm) %in% colnames(full$latent.mm))))
warning(call. = FALSE, hopit_msg(56))
if ((ncol(full$thresh.mm)) && (ncol(nested$thresh.mm)))
if (!(all(colnames(nested$thresh.mm) %in% colnames(full$thresh.mm))))
warning(call. = FALSE, hopit_msg(57))
stat <- 2L*( logLik.hopit(full) - logLik.hopit(nested))
if (!length(full$design)) {
df.diff <- length(full$coef.ls$latent.params) -
length(nested$coef.ls$latent.params) +
length(full$coef.ls$thresh.lambda) -
length(nested$coef.ls$thresh.lambda) +
length(full$coef.ls$thresh.gamma) -
length(nested$coef.ls$thresh.gamma) +
(full$hasdisp) - (nested$hasdisp)
p <- 1L - stats::pchisq(stat, df.diff)
} else {
stop(hopit_msg(58), call=NULL)
}
z <- list(chisq = stat, df = df.diff, pval = p, full = full, nested = nested)
class(z) <- 'lrt.hopit'
z
}
#' Printing an object calculated by \code{\link{lrt.hopit}}
#'
#' @param x an object obtained from \code{\link{lrt.hopit}}.
#' @param short a logical indicating whether to show a shortened description.
#' @param ... further arguments passed to or from other methods.
#' @keywords internal
#' @export
#' @usage \method{print}{lrt.hopit}(x, short = FALSE, ...)
#' @author Maciej J. Danko
#' @method print lrt.hopit
#' @seealso \code{\link{lrt.hopit}}, \code{\link{anova.hopit}},
#' \code{\link{hopit}}.
print.lrt.hopit <- function(x, short = FALSE, ...){
if (!short[1]) {
cat(hopit_msg(64))
cat("--", hopit_msg(65), deparse(x$full$latent.formula), fill = TRUE)
cat("--",hopit_msg(66), deparse(x$full$thresh.formula), fill = TRUE)
cat("--",hopit_msg(70),x$full$hasdisp, fill=TRUE)
cat(hopit_msg(67))
cat("--",hopit_msg(65), deparse(x$nested$latent.formula), fill = TRUE)
cat("--",hopit_msg(66), deparse(x$nested$thresh.formula), fill = TRUE)
cat("--",hopit_msg(70),x$nested$hasdisp, fill=TRUE)
cat('\n')
}
cat(hopit_msg(59))
# if (length(x$df)) {
out <- t(as.matrix(c('Chi^2' = unname(x$chisq), df = unname(x$df),
'Pr(>Chi^2)' = unname(x$pval))))
# out2 <- NULL
# else {
# out <- t(as.matrix(c('Chi^2' = unname(x$chisq),
# 'Pr(>Chi^2)' = unname(x$pval))))
# out2 <- x$scalef
# }
row.names(out) <- ''
stats::printCoefmat(out, signif.stars = TRUE, P.values = TRUE,
has.Pvalue = TRUE, digits = 5L, dig.tst = 3L, tst.ind = 1L)
# if (length(out2)) print(paste(hopit_msg(63),out2))
invisible(NULL)
}
#' Calculate the log likelihood profile for the fitted \code{hopit} model
#'
#' @param fitted a \code{hopit} object (a fitted model).
#' @param scope a value (fraction) defining the plotting range for a coefficient.
#' The range is \code{c(coef \* (1-scope), coef \* (1+scope))}.
#' @param steps at how many equally spaced points the log likelihood
#' function is calculated for each coefficient.
#' @param ... unused now.
#' @export
#' @keywords internal
#' @author Maciej J. Danko
#' @usage \method{profile}{hopit}(fitted, ..., scope = 0.15, steps = 101)
#' @method profile hopit
#' @seealso \code{\link{plot.profile.hopit}}, \code{\link{print.profile.hopit}},
#' \code{\link{hopit}}
#' @examples
#' \donttest{
#' # DATA
#' data(healthsurvey)
#'
#' # the order of response levels decreases from the best health to
#' # the worst health; hence the hopit() parameter decreasing.levels
#' # is set to TRUE
#' levels(healthsurvey$health)
#'
#' # Example 1 ---------------------
#'
#' # fitting the model:
#' model1 <- hopit(latent.formula = health ~ hypertension + high_cholesterol +
#' heart_attack_or_stroke + poor_mobility + very_poor_grip +
#' depression + respiratory_problems +
#' IADL_problems + obese + diabetes + other_diseases,
#' thresh.formula = ~ sex + ageclass + country,
#' decreasing.levels = TRUE,
#' control = list(trace = FALSE),
#' data = healthsurvey)
#'
#' # check the fit using the profile function (at 51 points)
#' pr <- profile(model1, steps = 51)
#' print(pr, plotf = FALSE)
#'
#' # plot profile
#' plot(pr, relative = FALSE)
#'
#' # alternative plot
#' plot(pr, relative = TRUE)
#' }
profile.hopit<-function(fitted, ..., scope=0.15, steps=101){
steps <- floor(steps/2)*2+1
if (fitted$hasdisp) COEF <- c(fitted$coef, fitted$coef.ls$logSigma) else
COEF <- fitted$coef
sub <- function(x,y) if (x==1) c(y,COEF[2:length(COEF)]) else
if (x==length(COEF)) c(COEF[-length(COEF)],y) else
c(COEF[1:(x-1)],y,COEF[(x+1):length(COEF)])
lo <- COEF*(1-scope)
hi <- COEF*(1+scope)
GG <- function(x) sapply(seq(lo[x],hi[x],length.out=steps),function(y)
hopit_negLL(parameters=sub(x,y),fitted,negative = FALSE))
val <- sapply(seq_along(COEF), function(x) GG(x))
attr(val,'scope') <- scope
attr(val,'steps') <- steps
attr(val,'lo') <- lo
attr(val,'hi') <- hi
colnames(val) <- names(COEF)
class(val) <- c("profile.hopit", "profile")
val
}
#' Plot the log likelihood profile for a profile.hopit object
#'
#' Plot the method for a profile.hopit object.
#' @param x a \code{profile.hopit} object.
#' @param leg.cex a character expansion factor relative to the current
#' \code{par("cex")} (see \code{\link{legend}}).
#' @param leg.col a color used for the legend text.
#' @param ylim see \code{\link{plot}}.
#' @param relative a logical indicating whether \code{ylim} on each panel should be the
#' same (\code{TRUE}) or not (\code{FALSE}).
#' @param ... arguments to be passed to the \code{\link{plot}}() function (see
#' \code{\link{par}}).
#' @export
#' @keywords internal
#' @usage \method{plot}{profile.hopit}(x, ..., ylim = NULL, relative = FALSE,
#' leg.cex = 0.85, leg.col = 'blue4')
#' @author Maciej J. Danko
#' @method plot profile.hopit
#' @seealso \code{\link{profile.hopit}}, \code{\link{print.profile.hopit}},
#' \code{\link{hopit}}.
plot.profile.hopit<-function(x, ..., ylim = NULL, relative = FALSE,
leg.cex = 0.85, leg.col = 'blue4'){
z <- sqrt(ncol(x))
zy <- round(z)
zx <- ceiling(z)
spar <- graphics::par(c('mfrow','mar'))
graphics::par(mfrow=c(zx,zy),mar=c(0,0,0,0))
if (relative[1]) ylim <- range(x)
for (j in seq_len(ncol(x))) {
graphics::plot(x[,j],type='l',axes='F', ylim = ylim, ...)
graphics::abline(v=floor(nrow(x)/2)+1,col=2,lty=2)
graphics::legend('bottom',colnames(x)[j], bty='n',cex=leg.cex,
text.col=leg.col)
graphics::box()
}
suppressWarnings(graphics::par(spar))
}
#' Print method for a profile.hopit object
#'
#' @param x a \code{profile.hopit} object.
#' @param plotf a logical indicating whether to plot the profile.
#' @param ... arguments to be passed to the \code{plot}() function (see
#' \code{\link{plot.profile.hopit}}).
#' @export
#' @keywords internal
#' @usage \method{print}{profile.hopit}(x, ..., plotf = TRUE)
#' @author Maciej J. Danko
#' @method print profile.hopit
#' @seealso \code{\link{profile.hopit}}, \code{\link{plot.profile.hopit}},
#' \code{\link{hopit}}
print.profile.hopit<-function(x, ..., plotf = TRUE){
test <- apply(x,2,which.max)==floor(nrow(x)/2)+1
if(plotf[1]) plot.profile.hopit(x, ...)
if (any(!test)) {
message(hopit_msg(60))
message(paste(hopit_msg(61),paste(names(test)[!test],sep='',
collapse = ', ')))
} else {
cat(hopit_msg(62))
}
}
#' Extract the \code{Sigma} parameter from a \code{hopit} model
#'
#' Extract the \code{Sigma} parameter from a \code{hopit} model
#' @param model a fitted \code{hopit} model.
#' @usage \method{sigma}{hopit}(object, ...)
#' @keywords internal
#' @export
#' @author Maciej J. Danko
#' @method sigma hopit
sigma.hopit <- function(object, ...)
unname(exp(object$coef.ls$logSigma))
Try the hopit package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.