Nothing
#' Analysis of Penalized Deviance for \code{logistf} Models
#'
#'This method compares hierarchical and non-hierarchical logistf models using
#'penalized likelhood ratio tests. It replaces the function logistftest of former
#'versions of logistf.
#'
#'Comparing models fitted by penalized methods, one must consider that the penalized likelihoods
#'are not directly comparable, since a penalty is involved. Or in other words, inserting zero for
#'some regression coefficients will not lead to the same penalized likelihood as if the corresponding
#'variables are simply "unknown" to a model. The anova method takes care that the same penalty is
#'used for two hierarchically nested models, and if the models are not hierarchically nested, it will
#'first relate each penalized likelihood to its null penalized likelihood, and only compare the resulting
#'penalized likelihod ratio statistics. The chi-squared approximation for this latter method (PLR) is
#'considered less accurate than that of the nested method. Nevertheless, it is the only way to go for
#'comparison of non-nested models.
#'
#' @param object A fitted \code{logistf} model object
#' @param fit2 Another fitted \code{logistf} model object, to be compared with \code{object}
#' @param formula Alternatively to \code{fit2}, a formula which specifies terms to omit from the object model fit.
#' @param method One of c("nested","PLR"). nested is the default for hierarchically nested
#' models, and will compare the penalized likelihood ratio statistics (minus twice
#' the difference between maximized penalized log likelihood and null penalized
#' log likelihood), where the null penalized log likelihood is computed from the
#' same, hierarchically superior model. Note that unlike in maximum likelihood
#' analysis, the null penalized likelihood depends on the penalty (Jeffreys prior)
#' which itself depends on the scope of variables of the hierarchically superior
#' model. PLR compares the difference in penalized likelihood ratio between the
#' two models, where for each model the null penalized likelihood is computed
#' within the scope of variables in that model. For PLR, the models need not be
#' hierarchically nested.
#' @param ... Further arguments passed to the method.
#'
#' @return An object of class \code{anova.logistf} with items
#' \item{chisq}{the chisquared statistic for the model comparison}
#' \item{df}{The degrees of freedom}
#' \item{pval}{The p-value}
#' \item{call}{The function call}
#' \item{method}{The method of comparison (input)}
#' \item{model1}{The first model}
#' \item{model2}{The second model which was compared to the first model}
#' \item{PLR1}{The PLR statistic of the first model}
#' \item{PLR2}{the PLR statistic of the second model; for the nested method, this will be the drop in chi-squared due to setting the coefficients to zero}
#' @export
#'
#' @examples
#' data(sex2)
#' fit<-logistf(data=sex2, case~age+oc+dia+vic+vicl+vis)
#'
#' #simultaneous test of variables vic, vicl, vis:
#' anova(fit, formula=~vic+vicl+vis)
#'
#' #test versus a simpler model
#' fit2<-logistf(data=sex2, case~age+oc+dia)
#' # or: fit2<-update(fit, case~age+oc+dia)
#' anova(fit,fit2)
#'
#' # comparison of non-nested models (with different df):
#' fit3<-logistf(data=sex2, case~age+vic+vicl+vis)
#' anova(fit2,fit3, method="PLR")
#'
#'
#' @rdname anova
#' @method anova logistf
#' @exportS3Method anova logistf
anova.logistf<-function(object, fit2, formula, method="nested", ...){
# methods: "PLR": take difference in PLR, "nested": proper method for nested models
# needed in logistf class: $firth, $data
mf <- match.call(expand.dots =FALSE)
m <- match(c("object","fit2","formula","method"), names(mf), 0L)
mf <- mf[c(1, m)]
object <- eval(mf$object, parent.frame())
fit2 <- eval(mf$fit2, parent.frame())
fit1<-object
if(missing(formula)){
out2 <- fit2$formula
if(fit1$df<fit2$df) {
ff0<-fit1 #swap
fit1<-fit2
fit2<-ff0
out2 <- fit2$formula
}
} else {
out2 <- formula
}
if(method=="PLR"){
if(fit1$df==fit2$df) stop("Models not comparable (equal df).\n")
df<-abs(fit1$df-fit2$df)
PLR1<- -2 * (fit1$loglik['null']-fit1$loglik['full'])
PLR2<- -2 * (fit2$loglik['null']-fit2$loglik['full'])
chisq<-PLR1-PLR2
if (chisq<0) chisq<-0
pval<-1-pchisq(chisq,df)
model2<-as.character(fit2$formula)
} else if(method=="nested"){
f1<-fit1$formula
a<-attr(terms(fit1),"term.labels")
#check if just intercept fitted:
if(missing(formula)){
f2<-fit2$formula
b<-attr(terms(fit2),"term.labels")
# find out about which model is nested in the other
upper<-f1
lower<-f2
if(!any(is.na(match(a,b)))) {
lower<-f1
upper<-f2
ab<-a
a<-b
b<-ab
} else if(any(is.na(match(b,a)))) {
stop("Models are not nested. Try method=PLR on non-nested models.\n")
}
aug<-a[is.na(match(a,b))]
f3<-paste("~",aug[1])
if(length(aug)>1) {
for(j in 2:length(aug)) {
f3<-paste(f3, aug[j], sep="+")
}
}
f3<-paste(f3, "-1")
f3<-as.formula(f3, env = environment(object$formula))
} else {
f3<-as.formula(paste(paste(as.character(formula), collapse=""),"-1",collapse=""), env = environment(object$formula))
}
test<-logistftest(object=fit1, test = f3, firth=fit1$firth, weights=model.weights(model.frame(fit1)),...)
chisq<- -2 * (test$loglik['null']-test$loglik['full'])
PLR1<- -2 * (fit1$loglik['null']-fit1$loglik['full'])
PLR2<-PLR1-chisq
df<-test$df
pval<-test$prob
model2<-as.character(f3)
}
res<-list(chisq=chisq, df=df, pval=pval, call=match.call(), method=method, model1=as.character(fit1$formula), model2=out2, PLR1=PLR1, PLR2=PLR2)
attr(res,"class")<-"anova.logistf"
return(res)
}
#' @method anova flic
#' @exportS3Method anova flic
anova.flic<-function(object, fit2, formula, method="nested", ...){
object$firth <- TRUE
if(missing(formula)){
if(inherits(fit2, "flic")){
fit2$firth <- TRUE
anova.logistf(object, fit2 = fit2, ...)
}
} else {
anova.logistf(object, formula = formula, ...)
}
}
#' @exportS3Method print anova.logistf
print.anova.logistf <-function(x,...){
cat("Comparison of logistf models:\n")
obj<-x
pdat<-data.frame(Formula=c(as.character(obj$model1),as.character(obj$model2)), ChiSquared=c(obj$PLR1,obj$PLR2))
pdatrows<-capture.output(print(pdat))
for(i in 1:3) cat(pdatrows[i],"\n")
cat("\nMethod: ", obj$method, "\n")
cat("Chi-Squared: ", obj$chisq, " df=",obj$df," P=", obj$pval,"\n")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.