Nothing
## relevance.R
## two groups and one sample
twosamples <- #f
onesample <- function(x, ...) UseMethod("twosamples")
## ---
twosamples.default <- #f
function(x, y=NULL, paired=FALSE, table=NULL,
hypothesis=0, var.equal=TRUE,
testlevel=getOption("testlevel"),
log=NULL, standardize=NULL,
rlv.threshold=getOption("rlv.threshold"), ...)
{ ## effect (group difference) and relevance
lcheck <-
list(x=list(cnr(na.ok=FALSE), cdf()), y=cnr(na.ok=TRUE), paired=clg(),
table=list(cnr(dim=c(2,2)), cdf()), hypothesis=cnr(), var.equal=clg(na.ok=FALSE),
testlevel=cnr(range=c(0.001,0.5)), log=list(clg(), cch()),
standardize=clg(), rlv.threshold=cnr(range=c(0,Inf))
)
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------
ltlev2 <- testlevel/2
lonegroup <- FALSE
lse <- NA
ltb <- FALSE
if (length(table)) {
table <- as.table(as.matrix(table))
lbin <- ltb <- TRUE
paired <- FALSE
lonegroup <- is.na(ncol(table))
if (lonegroup) {
lnx <- sum(table)
lny <- 0
} else {
lnx <- sum(table[1,])
lny <- sum(table[2,])
}
x <- 0
} else {
if (length(y)==0) {
if (NCOL(x)==2) {
y <- x[,2]
x <- x[,1]
} else {
if (!is.atomic(x))
stop("!twosamples! argument 'x' not suitable")
}
}
if (paired) {
if (length(x)!=length(y))
stop("!twosamples! 'x' and 'y' must have equal length if 'paired' is TRUE")
x <- y-x
lonegroup <- TRUE
}
else if (length(y)==0) lonegroup <- TRUE
## quantitative or proportion?
lbin <- length(table(c(x,y)))<=2 ## all(c(x,y)%in%c(0,1,NA))
## ---
x <- x[is.finite(x)]
lmnx <- mean(x)
ln <- lnx <- length(x)
lny <- 0
lest <- NULL
}
## ---
if (lbin) { ## binary data
lrlvth <-
i.def(rlv.threshold["prop"],
c(rlv.threshold, getOption("rlv.threshold")["prop"])[1])
lx <- if (ltb) table[2] else sum(x)
if (lonegroup) { ## one proportion
if (paired) x <- (x[x!=0]+1)/2 ## McNemar
lt <- binom.test(lx,lnx, conf.level=1-testlevel)
lest <- lt$estimate
## !!! se
leff <- unname(qlogis(lest))
leffci <- qlogis(lt$conf.int) ## c("ci.low","ci.up") ))
lar <-
qlogis( (qintpol(c(ltlev2, 1-ltlev2), lnx, plogis(hypothesis))+0:1) / lnx )
## c("accreg.low","accreg.up") ))
lsg <- (leff-hypothesis)/abs(lar-hypothesis)
lsig <- ifelse(leff<0, lsg[1], lsg[2])
lv <- lest*(1-lest)
effname <- paste("log odds", if(paired) "of changes")
method <- "One proportion, binomial inference"
} else { ## two proportions
if (ltb) {
if (any(dim(table)!=c(2,2)))
stop("!twosamples! argument 'table' not suitable")
} else {
y <- y[is.finite(y)]
lmny <- mean(y)
lny <- length(y)
table <- table(x,y)
}
lt <- fisher.test(table, conf.int=TRUE, conf.level=1-testlevel)
lest <- table[,2]/(table[,1]+table[,2]) ## c(mean(x), mean(y))
## !!! se
leff <- unname(log(lt$estimate))
leffci <- log(lt$conf.int) ## c("ci.low","ci.up") ))
lse <- diff(leffci)/(2*qnorm(1-ltlev2))
lsig <- leff/(leff-ifelse(leff>0, leffci[1], leffci[2]))
effname <- "log odds ratio"
method <- "Two proportions, Fisher's exact inference"
}
lrlvci <- c(leff, leffci)/lrlvth
ltst <- unname(lt$teststatistic)
lpv <- lt$p.value
ltype <- "proportion"
lsigth <- NA
ldf <- NULL
lv <- NA
} else { ## ----------- quantitative data
log <- i.def(log, "", "log", "")
llogrescale <- 1
if(log=="log10") llogrescale <- 1/log(10)
llog <- substr(log,1,3)=="log"
lirl <- ifelse(llog,"rel","stand")
lrlvth <-
i.def(rlv.threshold[lirl],
c(rlv.threshold, getOption("rlv.threshold")[lirl],
relevance.options[["rlv.threshold"]][lirl])[1])
lssx <- sum((x-lmnx)^2)
lpq <- 1-ltlev2
if (lonegroup) { ## one quantitative sample
leff <- lmnx
ldf <- lnx-1
lv <- lssx/ldf
lny <- 0
effname <- paste("mean", if(paired) "of differences")
method <- "One Sample t inference"
lest <- c(mean=lmnx, se=sqrt(lv/lnx))
} else { # two quantitative samples
y <- y[is.finite(y)]
lny <- length(y)
ln <- lnx+lny
lmny <- mean(y)
leff <- lmny-lmnx
lssy <- sum((y-lmny)^2)
lsex2 <- lssx/(lnx-1)/lnx
lsey2 <- lssy/(lny-1)/lny
if (var.equal) {
lv <- ln*(1/lnx+1/lny)*(lssx+lssy)/(ln-2)
ldf <- ln-2
method <- "Two Sample t inference, equal variances assumed"
} else { ## welch
lv <- ln*(lssx/(lnx-1)/lnx+lssy/(lny-1)/lny)
ldf <- (lsex2+lsey2)^2/(lsex2^2/(lnx - 1) + lsey2^2/(lny - 1))
method <- "Two Sample t inference, unequal variances (Welch)"
}
effname <- "difference of means"
lest <- cbind(mean=c(lmnx,lmny), se=sqrt(c(lsex2, lsey2)))
}
lse <- sqrt(lv/ln)
ltst <- leff/lse
lq <- qt(lpq, ldf)
lciwid <- lq*lse
leffci <- leff + c(-1,1)*lciwid
lpv <- 2*pt(-abs(ltst), ldf)
##
lsc <- 1
ltype <- "raw"
if (llog) {
lsc <- llogrescale
ltype <- "relative"
} else if(u.notfalse(standardize)) {
lsc <- sqrt(lv)
ltype <- "standardized"
}
lrlvci <- c(leff,leffci)/lsc/lrlvth
lsig <- leff/lciwid
lsigth <- c((lrlvci[1]-1)*2/diff(lrlvci[2:3]))
}
if (leff<0) lrlvci <- -lrlvci[c(1,3,2)]
rr <- structure(
c(effect=leff, se=lse, teststatistic=ltst, p.value=lpv, Sig0=lsig,
ciLow=leffci[1], ciUp=leffci[2],
Rle=lrlvci[1], Rls=lrlvci[2], Rlp=lrlvci[3],
Sigth=lsigth,
df = ldf, scatter=sqrt(lv)),
class = c("inference", "twosamples"), type="simple",
method = method, effectname = effname,
hypothesis = hypothesis,
n = c(lnx, lny), teststatistic = ltst, ## V = lv, df = ldf,
estimate = lest,
data = if (ltb) table else {if (lonegroup) x else list(x=x, y=y)},
rlv.threshold = lrlvth,
rlv.type = ltype
)
attr(rr, "rlvclass") <- rlvClass(rr)
rr
} # end twosamples.default
## ============================================================================
twosamples.formula <- #f
function (x, data=NULL, subset, na.action, log=NULL, ...)
{ ## adapted from t.test.formula
lcheck <-
list(data=cdf(), na.action=cfn())
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------
if (missing(x) || (length(x) != 3L))
stop("!twosamples! 'formula' must have left and right term")
lfy <- x[[2]]
llog <-
if(length(lfy)>1)
if(as.character(lfy[1])=="log") "log"
else if(as.character(lfy[1])%in%c("log10", "logst")) "log10"
llog <- i.def(log, if(is.null(llog)) "" else llog)
## ---
oneSampleOrPaired <- FALSE
if (length(attr(terms(x[-2L]), "term.labels")) != 1L)
if (x[[3]] == 1L)
oneSampleOrPaired <- TRUE
else stop("!twosamples! 'formula' incorrect")
m <- match.call(expand.dots = FALSE)
if (is.matrix(eval(m$data, parent.frame())))
m$data <- as.data.frame(data)
m[[1L]] <- quote(stats::model.frame)
m$... <- NULL
names(m)[2] <- "formula"
mf <- eval(m, parent.frame())
## DNAME <- paste(paste("`",names(mf),"'",sep=""), collapse = " by ")
names(mf) <- NULL
response <- attr(attr(mf, "terms"), "response")
if (!oneSampleOrPaired) { ## two indep. samples
g <- factor(mf[[-response]])
if (nlevels(g) != 2L)
stop("grouping factor must have exactly 2 levels")
rr <- do.call("twosamples",
c(setNames(split(mf[[response]], g), c("x", "y")),
list(paired=FALSE, log=llog, ...)))
## rr <- twosamples.default(DATA$x, DATA$y, paired=FALSE, log=llog, ...)
names(attr(rr, "n")) <- levels(g)
} else { ## one sample
respVar <- mf[[response]]
if (inherits(respVar, "Pair")) {
rr <- do.call("twosamples",
list(x = respVar[, 1], y = respVar[, 2],
paired = TRUE, log=llog, ...))
} else {
rr <- do.call("twosamples",
list(x = respVar, paired = FALSE, log=llog, ...))
}
}
## attr(rr, "data.name") <- DNAME
ldn <- substitute(data)
structure(rr, formula=x, data.name=if (is.name(ldn)) format.default(ldn))
}
## -------------------------------------------------------------------------
twosamples.table <- #f
function(x, ...) twosamples.default(table=x, ...)
## ========================================================================
correlation <- #f
function(x, y = NULL, method = c("pearson", "spearman"), ## , "kendall"
hypothesis = 0, testlevel=getOption("testlevel"),
rlv.threshold=getOption("rlv.threshold"), ...)
{ ## ---
lcheck <-
list(x=cnr(na.ok=FALSE), y=cnr(na.ok=TRUE), hypothesis=cnr(c(-1,1)),
testlevel=cnr(range=c(0.001,0.5)),
method=cch(c("pearson", "spearman")), ## , "kendall"
rlv.threshold=cnr(range=c(0,Inf))
)
##- lcall <- match.call(expand.dots = FALSE)
##- lcall$... <- NULL
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------
if (length(y)) {
if (length(x)!=length(y))
stop("!correlation! 'x' and 'y' must have the same length")
x <- cbind(x,y)
}
if (NCOL(x)!=2)
stop("!correlation! argument 'x' not suitable")
x <- stats::na.omit(x)
ln <- nrow(x)
if (ln<=3)
stop("!correlation! not enough observations")
## ---------
testlevel = i.def(testlevel, getOption("testlevel"))
lrlvth <-
i.def(rlv.threshold["corr"],
c(rlv.threshold, getOption("rlv.threshold")["corr"],
relevance.options[["rlv.threshold"]]["corr"])[1])
lmethod <- i.def(method, "pearson")
## ---
lt <- cor.test(x[,1], x[,2], method=lmethod[1], conf.level=1-testlevel, ...)
lest <- lt$estimate
lci <- lt$conf.int
lci <- if (length(lci)) atanh(lci)
else atanh(lest) + c(-1,1)*qnorm(1-testlevel/2)/sqrt(ln-3)
## !!! this may not be quite adequate for nonparametric correlation ...
## ---
leffci <- c(atanh(lest),lci)-atanh(hypothesis)
lrlvci <- leffci/lrlvth
lrlvwid <- lrlvci[1]-ifelse(lrlvci[1]>0, lrlvci[2], lrlvci[3])
lsig <- lrlvci[1]/lrlvwid
lsigth <- (lrlvci[1]-1)/lrlvwid
leffnm <- "correlation, z-transformed"
lmeth <- paste("Correlation --",
switch(lmethod[1],
pearson="Pearson product moment c.",
spearman="Spearman's rank c.",
kendall="Kendall's nonparametric c.") )
lneg <- lest<0
# lrl <- unname(lrlvci)
structure(
c(leffci[1], lt$statistic, lt$p.value, lsig, leffci[2:3],
(1-2*lneg)*lrlvci[c(1,2+lneg,3-lneg)], lsigth),
names = c("effect", "teststatistic", "p.value",
"Sig0", "ciLow", "ciUp", "Rle", "Rls", "Rlp", "Sigth"),
class = c("inference"), type="simple",
method = lmeth, effectname = leffnm, hypothesis = hypothesis,
n = length(x)-sumNA(x),
estimate = c(lest, ciLow=lt$conf.int[1], ciUp=lt$conf.int[2]),
data = x, rlv.threshold = lrlvth, rlv.type="correlation"
)
} ## end correlation
## ========================================================================
inference <- #f
function (object=NULL, estimate=NULL, teststatistic=NULL,
se=NA, n=NULL, df=NULL, stcoef=TRUE,
rlv=TRUE, rlv.threshold=getOption("rlv.threshold"),
testlevel = getOption("testlevel"), ...)
{ ## for a coefficients table,
## calculate confidence interval, significance and relevance
lcheck <-
list(object=list(cnr(),cls()), estimate=cnr(),
teststatistic=cnr(), se=cnr(), n=cnr(range=c(2,Inf)),
df=cnr(range=c(1,Inf)), stcoef=list(clg(),cnr()),
rlv=clg(), rlv.trheshold=cnr(range=c(0,Inf)),
testlevel=cnr(range=c(0.001,0.5))
)
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------------------------
if (inherits(object, "inference"))
warning(":inference: The first argument is already an 'inference' object")
if (inherits(object, relevance.modelclasses)) {
lcftb <- getcoeftable(object)
lsum <- summary(object)
lcl <- object$call
ldn <- lcl$data
if(is.symbol(ldn)) ldn <- as.character(ldn)
##- if (length(ldf)==0)
##- ldf <-
##- if (class(object)[1]%in%c("rlm", "coxph"))
##- length(object$residuals)-length(object$coef)
##- else df.residual(object)
## ldfm <- summary(object)$df[1] ## ok for lm, glm, survreg, ...
ltt <- termtable(object, lsum, rlv=rlv, ...)
## special for polr
lic <- NULL
if (inherits(object,"summary.polr")) {
lcf <- coefficients(lsum)
lic <- lcf[(lsum$pc+1):nrow(lcf),1:2]
}
rr <-
structure(
list(summary = lsum, termtable = ltt,
termeffects = termeffects(object, rlv=rlv, ...),
intercepts = lic
## , deviancetable = NULL !!!
),
class=c("inference", "model", "list"), ## type = "model",
method=as.character(lcl[1]), formula=lcl$formula, data.name=ldn,
rlv.threshold=attr(ltt, "rlv.threshold"))
## --- standardized coefficients
if (!u.notfalse(stcoef)) return (rr)
## warning(":ciSigRlv: no standardized coefficients given. No relevances")
if (length(stcoef)==0||u.true(stcoef)) {
if (length(object)==0) {
warning(":inference: argument is not a fitted model. No relevances")
return (rr)
}
if (inherits(object, "nls")) {
warning(":inference: cannot calculate standardized coefficients for 'nls' objects.")
return(rr)
}
lfac <- getcoeffactor(object)
## lcls <- attr(lfac, "fitclass")
stcoef <- lcftb[,1]*lfac[row.names(lcftb)]
}
if (length(stcoef)!=NROW(lcftb)) {
warning(":inference: argument 'stcoef' not suitable. No relevances")
return (rr)
}
lest <- lcftb[,1]
ltq <- qt(1-testlevel/2, attr(lcftb,"df"))
lci <- lest+outer(ltq*lcftb[,2], c(ciLow=-1,ciUp=1))
lfac <- stcoef/lest
lstci <- cbind(effect=stcoef, effLow=lfac*lci[,1], effUp=lfac*lci[,2])
lsum$coefficients <- cbind(lcftb, lstci)
rr$summary <- lsum
## ---
return(structure(rr, summary=lsum))
} ## end model
## --------------------------------------------------
if (length(object) && length(dim(object))==0) {
if(inherits(object, "htest"))
stop("!inference! cannot cope with test results.",
"Call inference on data directly!")
## object <-as.data.frame(rbind(object))
}
ldt <- i.getIfrData(object, estimate=estimate,
teststatistic=teststatistic, se=se, n=n, df=df)
lest <- ldt$estimate
lse <- ldt$se ## NA is ok
if (length(lse)==0) lse <- attr(object, "se")
ltst <- i.def(ldt$teststatistic, lest/lse)
ln <- ldt$n
ldf <- i.def(ldt$df, attr(object, "df")) ## available if estimate is generated by getcoeftable
ldf <- i.def(ldf, ln-1)
ldfm <- 1 ## potentially modified in some cases
ln <- i.def(ln, ldf+ldfm)
if (length(ln)==0 || any(is.na(ln)))
stop("!inference! either 'n' or 'df' must be available")
if (length(lse)==0 || all(is.na(lse))) lse <- lest/ltst
if (length(lse)==0) {
warning(":inference: no standard error available -> NA")
lse <- NA
}
##- df <- i.def(ldf, ln-1, valuefalse=NA)
##- if (i.def(df,0)==0) {
##- warning(":inference: no finite standard errors")
##- return(cbind(object))
##- }
## ---
##- if (!(is.atomic(object)||length(dim(object))==2))
##- stop("!inference! first argument not suitable")
##- if (is.null(se))
##- if (NCOL(object)>1) {
##- se <- object[,2]
##- }
##- if (length(dim(object)))
##- object <- structure(object[,1], names=row.names(object))
ltq <- qt(1-testlevel/2, ldf) ## if (is.finite(ldf)) qt(1-testlevel/2, ldf) else qnorm(1-testlevel/2)
if(length(ldf)==1) ltq <- rep(ltq, length(lest))
lsc <- lse*sqrt(ln)
lest <- i.def(lest, ltst*lse)
lci0 <- outer(ltq, c(estimate=0, ciLow=-1, ciUp=1)) ## matrix
## if(length(ltq)==1) lci0 <- lci0[rep(1,length(lse)),]
lci <- lest+lci0*lse
## ltst <- lest/lse
leff <- ltst/sqrt(ln)
leffci <- structure(leff+lci0/sqrt(ln),
dimnames=list(NULL,c("effect", "effLow", "effUp")))
lsgf <- ltst/ltq
lpv <- 2*pt(-abs(ltst), ldf)
rr <- data.frame(lci, leffci, se=lse,
teststatistic=ltst, p.value=lpv, Sig0=lsgf,
n=ln, scatter=lsc)
## --- relevance
if (!u.notfalse(rlv)) return(rr)
if (length(rlv.threshold)>1) rlv.threshold <- rlv.threshold["coef"]
if (is.na(rlv.threshold)) {
warning(":inference: argument 'rlv.threshold' not suitable. No relevances")
return (rr)
}
lrlv <- leffci/rlv.threshold
li <- which(leffci[,1]<0)
if (length(li)) lrlv[li,] <- - lrlv[li,c(1,3,2)]
colnames(lrlv) <- c("Rle","Rls","Rlp")
lrlvclass <-
rlvClass(lrlv[,"Rle"], lrlv[,c("Rls","Rlp")], relevance=1)
structure(cbind(rr, lrlv, rlvclass=lrlvclass),
row.names=row.names(ldt),
class=c("inference", "data.frame"), type = "simple")
}
## ==========================================================================
termtable <- #f
function (object, summary=summary(object), testtype=NULL, r2x = TRUE,
rlv = TRUE, rlv.threshold = getOption("rlv.threshold"),
testlevel = getOption("testlevel"), ...)
{
## Purpose: generate term table for various models
## --------------------------------------------------------
lcheck <-
list(object=cls(), summary=cls(), testtype=cch(), r2x=clg(na.ok=FALSE),
rlv=clg(), rlv.threshold=cnr(range=c(0,Inf)),
testlevel=cnr(range=c(0.001,0.5))
)
##- lcall <- match.call(expand.dots = FALSE)
##- lcall$... <- NULL
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------------------------
## thresholds
ltlev1 <- 1-testlevel
## lrlthrl <- rlv.threshold[["rel"]]
lrlthcf <- rlv.threshold[["coef"]]
lrlthdr <- rlv.threshold[["drop"]]
lrlthpr <- rlv.threshold[["pred"]]
## --- scatter and threshold for standardized coefficient
## lscatter <- c(object$scatter, summary$scatter, 1)[1]
lf.intercept <- function(object)
length(ln <- names(object$coefficients))&&ln[1]=="(Intercept)" ## !!! polr and other models?
lfamily <- if (is.character(lfm <- object$family)) lfm else lfm$family
ldist <- object$dist
lmethod <- paste(c(setdiff(class(object),"regr")[1],
if (length(lfamily)) c(" ; family = ", lfamily),
if (length(ldist)) c(" ; distribution = ", ldist),
" : Drop-term inference"),
collapse="")
ldname <- attr(object, "data.name")
if (length(ldname)==0) ldname <- object$call$data
ldname <- if (is.language(ldname)) as.character(ldname)
else {
if (!is.character(ldname)) ldname
}
## if ((inherits(object, "glm") &&
## lfamily %in% c("poisson","quasipoisson")) ||
## (inherits(object, "survreg")&&
## ldist %in% c("weibull", "exponential", "lognormal"))
## ) lrlthcf <- lrlthrl
## --- testtype
if (length(testtype)==0) {
testtype <- "LRT"
if (inherits(object, c("lm","lmrob"))) testtype <- "F"
if (inherits(object, "glm")) {
testtype <-
if (lfamily %in% c("quasibinomial","quasipoisson")) "F" else "LRT"
}
if (inherits(object, c("survreg", "polr"))) testtype <- "Chisq"
}
## ---
lterms <- terms(object)
ldfres <- df.residual(object)
if (class(object)[1]%in%c("rlm", "coxph"))
ldfres <- length(object$residuals)-length(object$coef)
if(length(attr(lterms,"term.labels"))==0 | ldfres<1)
return(data.frame(
coef=i.def(object$coefficients,NA), df = NA, se=NA,
teststatistic=NA, p.value=NA, Sig0=NA, ciLow=NA, ciUp=NA,
stcoef=NA, stciLow = NA, stciUp = NA, R2x=NA,
stringsAsFactors=FALSE)
)
## --- coefficients
lcoef <- object$coefficients
lcft <- getcoeftable(object)
colnames(lcft) <-
c("estimate","se","teststatistic","p.value")[1:ncol(lcft)]
lcoeftab <-
inference(lcft, rlv=rlv, rlv.threshold=lrlthcf)
names(lcoeftab)[1] <- "coef"
## --- drop1
ldr1 <-
if (inherits(object, c("lm","lmrob"))&&!inherits(object, "glm")) {
try(drop1Wald(object, test=testtype, scope=lterms), silent=TRUE)
} else {
try(i.drop1(object, scope=lterms, test=testtype), silent=TRUE)
}
if (inherits(ldr1, "try-error")) {
warning(":termtable: drop1 did not work. I return the coefficient table")
lii <- match(c("Rle","Rls","Rlp"), names(lcoeftab), nomatch=0)
names(lcoeftab)[lii] <- c("coefRle","coefRls","coefRlp")[lii!=0]
return(
structure(cbind(coef=lcoeftab[,1,drop=FALSE],df=1,lcoeftab[,-1,drop=FALSE]),
class=c("inference", "data.frame"), type="terms",
formula=formula(object), data.name=ldname
))
}
ldr1 <- ldr1[-1,]
ldr1$RSS <- NULL # same ncol for lm and glm
ldf <- ldr1[,1]
if (inherits(object,"rlm")) ldr1[,4] <- ldr1[,2]/ldf
## -- critical value for test
ltstq <- if (testtype=="F") qf(ltlev1,c(1,pmax(1,ldf)),ldfres) else {
if (testtype%in%c("Chisq", "LRT")) qchisq(ltlev1,c(1,pmax(1,ldf))) else NA }
ltstq[which(ldf==0)+1] <- NA
ltstq1 <- sqrt(ltstq[1]) ## 1 degree of freedom
ltstq <- ltstq[-1]
##- if (inherits(object,"mlm")||inherits(object,"manova"))
##- return(list(termtable=ldr1)) ## !!! needs much more
## ---------------------
## model.matrix
lmmt <- object[["x"]]
if (length(lmmt)==0) lmmt <- model.matrix(object)
lasg <- attr(lmmt,"assign")[!is.na(lcoef)]
## if (class(object)[1] %in% c("polr")) lasg <- lasg[-1] ## ,"coxph"
## terms without factor involvement
lfactors <- attr(lterms,"factors")
lvcont <- !attr(lterms,"dataClasses")[row.names(lfactors)] %in%
c("numeric","logical") ## [...] excludes .weights. and possibly others
## terms only containing continuous variables
lcont <- which( lvcont %*% lfactors ==0 )
## -----------------------------------
## --- r2x
lr2 <- NA
if (r2x)
{
if (inherits(object, c("polr", "coxph")))
warning("!termtable! No R2x for such models")
else {
lvift <- ## lterms: n of levels for each term
## if (u.debug()) vif.regr(object, mmat=lmmt) else
try(vif.regr(object, mmat=lmmt), silent=TRUE)
if (inherits(lvift, "try-error") || length(lvift)==0) {
warning(":termtable: error in the calculation of R2x")
lvif <- NA
} else lvif <- lvift[,3]^2
lr2 <- 1-1/lvif
}
}
## -----------------------------------
## --- prepare table
lpvcol <- pmatch("Pr(",names(ldr1), nomatch=ncol(ldr1))
lpv <- ldr1[,lpvcol]
ldf <- ldr1[,1]
## !!!
lcont <- which(ldf==1)
lnobs1 <- ldfres+sum(ldf)-lf.intercept(object)
## drop effect relevance
ltst <- ldr1[, ifelse(inherits(object, "polr"), 3, 4)]
ldrncci <- rbind(confintF(ltst, ldf, ldfres, testlevel))
ldreff2 <- cbind(ltst*ldf, ldrncci)/lnobs1
ldrrl <- cbind(NA,NA,NA, sqrt(ldreff2)/lrlthdr,
pmax(0.5*log((ldfres+ldreff2*lnobs1)/(ldfres+ldf))/lrlthpr, 0))
dimnames(ldrrl) <-
list(NULL, c(t(outer(c("coef","drop","pred"),c("Rle","Rls","Rlp"),paste, sep=""))))
ltst <- ldr1[,lpvcol-1]
## table, filled partially
ltb <- data.frame(coef=NA, ldr1[,1,drop=FALSE], ## keep name if only 1 coef
se=NA, teststatistic=ltst, p.value=lpv,
Sig0=sqrt(pmax(0,ltst)/ltstq), ciLow=NA, ciUp=NA,
stcoef=NA, stciLow=NA, stciUp=NA, R2x=lr2, ldrrl,
stringsAsFactors=FALSE)
names(ltb)[2] <- "df"
## intercept
ljint <- "(Intercept)"==names(lcoef)[1]
if (ljint) {
ltb <- rbind("(Intercept)"=ltb[1,],ltb)
ltb[1,] <- NA
ltb[1,"df"] <- 1
ltstq <- c(ltstq1, ltstq)
lcont <- c(0, lcont)
} else if (!inherits(object, c("coxph", "polr")))
warning(":termtable: No intercept. Statistics are difficult to interpret.")
lcont1 <- lcont+ljint # row number in dr1
## --- coefficients and statistics for terms with 1 df
if (length(lcont)) { ## lcont refers to assign
## ltlb <- dimnames(ltb)[[1]]
## lclb <- ltlb[lcont1] ## lcont1 is the row in the coef table of summary(object)
ljc <- match(lcont,lasg) # index of coefs for cont variables
## lj <- c("coef","se","teststatistic","ciLow","ciUp","stcoef", "stciLow","stciUp",
lj <- c("coef","se","teststatistic","ciLow","ciUp","effect", "effLow","effUp",
"coefRle","coefRls","coefRlp")
ljj <- sub("coefRl","Rl",lj)
ltb[lcont1,lj] <- lcoeftab[ljc,ljj]
ltb[lcont1,"Sig0"] <- sign(ltb[lcont1,"coef"])*ltb[lcont1,"Sig0"]
}
structure(ltb, class=c("termtable", "data.frame"), type="terms",
testtype=testtype, method=lmethod,
fitclass=class(object), family=lfamily, dist=ldist,
formula=formula(object), data.name=ldname,
rlv.threshold=rlv.threshold[c("coef","drop","pred")])
}
## end termtable
## --------------------------------------------------------------------
i.drop1 <- #F
function (object, scope=drop.scope(object), scatter = 0, test = NULL, k = 2,
sorted = FALSE, ...)
{
## Purpose: drop1/add1 for regr objects
## ----------------------------------------------------------------------
lfam <- i.def(object$distrname, "gaussian")
if (is.null(test))
test <- if (inherits(object, c("lm","roblm"))||
((lfam=="binomial"|lfam=="poisson")&&object$dispersion>1)) {
if (inherits(object,"mlm")) "Wilks" else "F" }
else "Chisq"
if (length(scope)==0) { ## || !is.formula(scope) ## drop.scope is character
warning(":drop1/add1.regr: no valid scope")
return(data.frame(Df = NA, "Sum of Sq" = NA, RSS =NA, AIC = NA,
row.names = "<none>") )
}
class(object) <- setdiff(class(object), "regr")
fcall <- object$funcall
if (!is.null(fcall)) object$call <- fcall
##- if (class(object)[1]%in%c("lmrob")) ## to be expanded
##- drop1Wald(object, test="F", ...) else {
ldata <- object$model
if (is.null(ldata)) ldata <- eval(object$call$data)[,all.vars(formula(object))]
if (is.null(ldata)) stop("!i.drop1! no data found ")
## all predictors must get the same missing observations
object$call$data <- na.omit(ldata)
dr1 <- drop1(object, scope=scope, scatter=scatter, test=test, k=k, ...)
if (sorted)
if (0!=(lsrt <- match(c("AIC","p.value"),colnames(dr1), nomatch=0)))
dr1 <- dr1[order(dr1[, lsrt[1]]), ]
dr1
}
## ===========================================================
termeffects <- #f
function (object, se = 2, df = df.residual(object),
rlv = TRUE, rlv.threshold=getOption("rlv.threshold"), ...)
## --------------------------------------------------------------
{
lcheck <-
list(object=cls(), se=cnr(), df=cnr(),
rlv=clg(), rlv.threshold=cnr(range=c(0,Inf))
)
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------------------------
if (is.atomic(object)||is.null(terms(object)))
stop("!termeffects! inadequate first argument")
## xl <- object$xlevels
Terms <- delete.response(terms(object))
int <- attr(Terms, "intercept")
facs <- attr(Terms, "factors")
if (length(facs)==0) {
message(".termeffects. No terms, no termeffects")
return(NULL)
}
lfamily <- if (is.character(lfm <- object$family)) lfm else lfm$family
lmethod <- paste(c(class(object),
if (length(lfamily)) c(" ; family = ", lfamily),
if (length(ldist <- object$dist)) c(" ; distribution = ", ldist),
" : Term effects"),
collapse="")
tl <- attr(Terms, "term.labels")
dcl <- attr(Terms,"dataClasses")[-1]
## result already available?
allc <- object$termeffects
if ((!is.null(allc))&&length(allc)==length(tl)&&
(is.matrix(allc[[length(allc)]])|!se)) return(allc) ## !!! check!
## ---
## if (rlv) lscatter <- getscalepar(object)
mf <- object$model ##! d.c used all.vars
if (is.null(mf)) {
object$call$model <- object$call$x <- object$call$envir <- NULL
mf <- model.frame(object)
}
xtnm <- dimnames(facs)[[1]] ## names ##! replaces vars
xtlv <- lapply(mf[,xtnm, drop=FALSE],function(x) levels(x)) ## levels
lcontr <- object$contrasts
imat <- which(substr(dcl,1,7)=="nmatrix") ## resulting from bs()
if (length(imat)) {
xtlv[imat] <-
lapply(as.list(dcl[imat]),
function(x) as.character(1:as.numeric(substr(x,9,12))))
## lcontr <- c(lcontr, structure(rep(contr.id,length(tl)), names=tl)[imat])
lctr <- list()
for (li in imat) lctr <- c(lctr, list(diag(length(xtlv[[li]]))))
names(lctr) <- names(dcl)[imat]
lcontr <- c(lcontr, lctr)
}
xtnl <- pmax(sapply(xtlv,length),1) ## number of levels
termnl <- apply(facs, 2L, function(x) prod(xtnl[x > 0])) ##! lterms
nl <- sum(termnl)
## --- df.dummy: data frame of simple terms
args <- setNames(vector("list", length(xtnm)), xtnm)
for (i in xtnm)
args[[i]] <- if (xtnl[[i]] == 1) rep.int(1, nl) else
factor(rep.int(xtlv[[i]][1L], nl), levels = xtlv[[i]])
df.dummy <- as.data.frame(args) # do.call("data.frame", args)
names(df.dummy) <- xtnm
## rnn: names of rows
pos <- 0
rn <- rep.int(tl, termnl)
rnn <- rep.int("", nl)
## fill df.dummy
for (j in tl) {
i <- unlist(xtnm[facs[, j] > 0])
ifac <- i[xtnl[i] > 1]
if (length(ifac) == 0L) {
rnn[pos + 1] <- j
}
else if (length(ifac) == 1L) {
df.dummy[pos + 1L:termnl[j], ifac] <- xtlv[[ifac]]
rnn[pos + 1L:termnl[j]] <- as.character(xtlv[[ifac]])
}
else {
tmp <- expand.grid(xtlv[ifac])
df.dummy[pos + 1L:termnl[j], ifac] <- tmp
rnn[pos + 1L:termnl[j]] <-
apply(as.matrix(tmp), 1L, function(x) paste(x, collapse = ":"))
}
pos <- pos + termnl[j]
}
## attributes
attr(df.dummy,"terms") <- attr(mf,"terms")
lci <- sapply(df.dummy,is.factor)
lcontr <- lcontr[names(lci)[lci]] ## factors with 1 level have disappeared (?)
if (inherits(object, "polr")) {
attr(Terms, "intercept") <- 1
mm <- model.matrix(Terms, df.dummy, contrasts.arg=lcontr, xlev=xtlv)
asgn <- attr(mm, "assign")[-1]
mm <- mm[,-1]
} else {
mm <- model.matrix(Terms, df.dummy, contrasts.arg=lcontr, xlev=xtlv)
asgn <- attr(mm, "assign")
}
if (anyNA(mm)) {
warning("some terms will have NAs due to the limits of the method")
mm[is.na(mm)] <- NA
}
## calculate dummy coefs
coef <- object$coefficients ##!!! cf <-
##- if (!use.na)
##- coef[is.na(coef)] <- 0
names(asgn) <- colnames(mm)
##- lnna <- is.finite(coef)
##- if (any(!lnna)){
##- coef <- coef[lnna]
##- mm <- mm[,lnna]
##- asgn <- asgn[lnna]
##- }
if (se) {
cov <- vcov(object)
if (is.null(cov)) {
warning(":termeffects: no covariance matrix of coefficients found.",
" Returning coefficients only")
se <- FALSE
} else if (inherits(object, "polr"))
cov <- cov[1:length(coef),1:length(coef)]
}
##- licf <- pmatch(colnames(mm), names(coef))
licf <- pmatch(names(coef), colnames(mm))
lasgn <- asgn[licf]
mm <- mm[,licf,drop=FALSE]
## asgn <- asgn[names(coef)] ## !!!
res <- setNames(vector("list", length(tl)), tl)
ljfail <- NULL
for (j in seq_along(tl)) {
mmr <- rn == tl[j] ## rows corresponding to the term
mmc <- lasgn==j ## & !is.na(coef)
##- lcf <- coef[licf[mmc]]
lcf <- coef[mmc]
## mmc <- names(asgn)[asgn == j & !is.na(coef)] ## lcols (logical fails for polr, vcov() too large) !!! was which
##- mmpart <- mm[mmr, mmc, drop=FALSE]
if (all(is.finite(lcf))) {
mmpart <- mm[mmr,mmc, drop=FALSE]
rrj <- setNames(drop(mmpart %*% lcf), rnn[mmr]) ## coef[mmc]
if (se) {
sej <- sqrt(diag(mmpart %*% cov[mmc,mmc] %*% t(mmpart)))
if (any(is.na(rrj))|any(!is.finite(sej))) {
ljfail <- c(ljfail, tl[j])
} else {
lscatter <- getscalepar(object)
rrj <- inference(estimate=rrj, se=sej, df=df,
stcoef=rrj*0.5/lscatter, rlv=rlv, ...)
li <- match(c("Rle","Rlp","Rls"), names(rrj))
names(rrj)[li] <- c("coefRle","coefRlp","coefRls")
}
}
res[[j]] <- rrj
}
}
if (length(ljfail))
warning(":termeffects: error calculating se for terms ",
paste(ljfail, collapse=", "))
if (int > 0) {
res <- c(list(`(Intercept)` = coef[int]), res)
}
##- if (inherits(object, "polr")) {
##- lcfi <- object$intercepts[,1]
##- res <- c(res,
##- list("(Intercepts)"=inference(lcfi, object$intercepts[,2], df=df, stcoef=lcfi) ))
##- }
## class(res) <- "termeffects" ## don't do that:
## want to be able to print the whole table
structure(res, class="termeffects", head=lmethod, rlv.threshold=rlv.threshold)
} ## end termeffects
## -------------------------------------------------------------------------
"[.termeffects" <- #f
function(x, i=NULL) structure(unclass(x)[i], class="termeffects")
## unclass avoids infinite recursion!
## ============================================================================
qintpol <- function(p, par1, par2, dist="binom")
{
## interpolated (theoretical) quantile for discrete distributions
lqfn <- get(paste("q",dist,sep=""))
lpfn <- get(paste("p",dist,sep=""))
lq <- lqfn(p, par1, par2)
lp <- lpfn(lq, par1, par2)
lmaxx <- if (dist=="binom") par1 else Inf
lp0 <- lpfn(pmin(pmax(0,lq-1),lmaxx), par1, par2)
lq - (lp-p)/(lp-lp0)
}
## ------------------------------------------------------------------
confintF <- #f
function(f, df1, df2=Inf, testlevel=0.05) {
## confidence interval for non-centrality of F distribution
p <- testlevel/2
lf.fq <- function(x, fvalue, df1, df2, p) qf(p,df1,df2,x)-fvalue
lf.ciup <- function(fvalue, df, p) { ## upper bound for upper limit
lq <- 1.5*qnorm(p)
lu <- lq^2*2/df
df*(fvalue-1+lu+sqrt(lu*(lu+2*fvalue-1)))
}
ln <- max(length(f), length(df1), length(df2), length(p))
f <- rep(f, length=ln)
df1 <- rep(df1, length=ln)
df2 <- rep(df2, length=ln)
p <- rep(p, length=ln)
p <- pmin(p,1-p)
## ---------------------------
rr <- matrix(NA, ln, 2)
for (li in 1:ln) {
lx <- f[li]
if (!is.finite(lx)) next
if (lx>100)
rr[li,] <- df1[li]*(sqrt(lx)+c(-1,1)*abs(qt(p[li],df2[li]))/sqrt(df1[li]))^2
else {
ldf1i <- df1[li]
if (ldf1i==0) next
rr[li,1] <- { ## lower limit
lf0 <- lf.fq(0, f[li], ldf1i, df2[li], 1-p[li])
if (lf0>=0) 0
else
pmax(0, ## !!! something wrong, should not be needed
uniroot(lf.fq, c(0,df1[li]*f[li]),
fvalue=f[li], df1=ldf1i, df2=df2[li], p=1-p[li])$root
)
}
rr[li,2] <- ## upper limit
if (pf(f[li], ldf1i, df2[li])<=p[li]) 0 ## tiny F value
else
uniroot(lf.fq, interval=c(df1[li]*f[li], lf.ciup(f[li], df1[li], 1-p[li])),
fvalue=f[li], df1=ldf1i, df2=df2[li], p=p[li], extendInt="upX")$root
}
}
if (ln==1) c(rr) else rr
}
## ===========================================================================
print.inference <- #f
function (x, show = getOption("show.inference"), print=TRUE,
digits = getOption("digits.reduced"),
transpose.ok = TRUE, legend = NULL, na.print = getOption("na.print"),
...)
{
lf.format <- function(x, header="")
{ ## format matrix -> character vector including names
x <- format(x)
x <-
if (length(dim(x))==0)
c(header,
paste(apply(format(rbind(names(x),x)), 1, paste, collapse=" "),
collaps="\n"))
else {
lcnm <- colnames(x)
ltb <- rbind(lcnm,x)
lrnm <- c(if(length(lcnm)) "", row.names(x))
c(header, ## , rep("", ncol(x)-1+(length(lrnm)>0))),
paste(apply(format(cbind(lrnm, ltb), justify="right"),
1, paste, collapse=" "), "\n", sep=""))
##- cbind(format(c(header, rep("", nrow(x)))),
##- format(c("",rownames(x))),
##- format(rbind(colnames(x), x)), "\n")
}
unname(x) ## apply(x, 1, paste, collapse=" "))
}
## -------------------------------------
lismodel <- inherits(x, "model")
if (is.list(x)&&!(is.data.frame(x)||lismodel)) return(print(unclass(x)))
## show what?
lsh <- c(show.signif.stars = getOption("show.signif.stars"),
show.symbollegend = getOption("show.symbollegend"),
show.rlv.threshold = getOption("show.rlv.threshold"),
show.rlv.class = getOption("show.Rlv.class"))
legend <- i.def(legend, lsh,
structure(rep(TRUE, length(lsh)), names=names(lsh)))
lItable <- length(dim(x))
type <- i.def(attr(x,"type"), "simple")
lshow <- if (lismodel) {
lsh <- if (is.list(show)) show[1:2] else list(show, show)
list(i.getshow(lsh[[1]], "terms", x=x$termtable),
i.getshow(lsh[[2]], "termeffects", x=x$termeffects))
} else i.getshow(show, type, x=x)
## --- head
ldn <- attr(x, "data.name")
if (!(is.character(ldn)&length(ldn)==1)) ldn <- NULL
lout <- c(
if (length(ldn)) paste("data: ", ldn),
if (length(lfo <- attr(x, "formula")))
paste("target variable: ", format(lfo[[2]]) ) ## was format()
)
lhead <- c(attr(x, "method"),
if (length(ldist <- attr(x, "distribution")))
c(" ; distribution: ", ldist), "\n",
if (length(lout)) paste(paste(lout, collapse=" ; ")))
## lhead <- if (any(nchar(lhead)>2)) c(lhead,"\n")
## ---
lleg <- NULL
lx <- x
## ---------------------------------------------------------------
if (lismodel) {
lout <-
print.coeftable(x$termtable, show=lshow[[1]], digits=digits,
transpose.ok=FALSE, na.print=na.print, print=FALSE)
lte <- x$termeffects
if (!u.true(single)) lte <- lte[sapply(lte, function(x) NROW(x)>1)]
if (length(lte)) {
ltep <- lapply(lte, print.coeftable, show=lshow[[2]], digits=digits,
transpose.ok=transpose.ok, na.print=na.print, print=FALSE)
lout <-
list(termtable=structure(lout, tail="\n"), termeffects=ltep)
}
} else { ## --------------------------------
if (!lItable) {
lout <- c(
if (!is.na(leff <- x[grep("effect", names(x))]))
paste(if(length(leffn <- attr(x, "effectname")))
paste(leffn, ": ",sep="") else "effect: ", format(leff)),
if (getOption("show.confint")) {
lci <- x[c("ciLow","ciUp")]
if (any(!is.na(lci)))
paste(" ; confidence int.: [",
paste(format(lci), collapse=", "),"]\n")
else "\n"
} else "\n"
)
## ---
if (any(c("teststatistic","p.value","Sig0")%in%lshow)) {
lpv <- x["p.value"]
lps <-
if (length(lpv)& ("p.symbol"%in%lshow | getOption("show.signif.stars")) )
symnum(lpv, p.symbols$cutpoint, p.symbols$symbol)
llout <- c(
paste("Test: hypothesis: effect = ", attr(x,"hypothesis"), "\n "),
paste(c(if(length(ltst <- attr(x, "teststatistic")))
paste("teststatistic: ", round(ltst, digits)),
## !!! df
if("p.value"%in%lshow && length(lpv <- x["p.value"]))
paste("p value: ", round(lpv, digits+1), lps),
if("Sig0"%in%lshow&&length(lsig <- x["Sig0"]))
paste("Sig0: ", round(lsig, digits),
if (!"p.value"%in%lshow) lps)), collapse=" ; ")
)
lout <-
c(lout, if (length(llout)) paste(paste(llout, collapse=""), "\n") )
}
## ---
if (any(substring(lshow,1,2)=="Rl")) {
lrs <- if ("Rls.symbol"%in%lshow)
symnum(x["Rls"], rlv.symbols$cutpoint, rlv.symbols$symbol)
llout <- c(
if("Rle"%in%lshow) paste("Rle: ", round(x["Rle"], digits)),
if("Rlp"%in%lshow) paste("Rlp: ", round(x["Rlp"], digits)),
if("Rls"%in%lshow) paste("Rls: ", round(x["Rls"], digits), lrs),
if("Rlv.class"%in%lshow) paste("Rlv.class: ", attr(x, "rlvclass"))
)
if (length(llout)) lout <- c(lout, paste(paste(llout, collapse=" ; "), "\n") )
}
} else { ## ======================== table
lout <- print.coeftable(x, lshow, print=FALSE)
}
}
## --- legend(s)
if(u.notfalse(legend)) {
lsh <- unlist(lshow)
lrlv <- length(intersect(grep("Rl", lsh), grep(".symbol", lsh))) > 0
lleg <- c(
if(u.true(legend["show.signif.stars"]) && "p.symbol"%in%lsh)
paste("Significance codes for p.value: ", i.symleg(getOption("p.symbols")),"\n"),
if(u.true(legend["show.symbollegend"]) && lrlv)
paste("Relevance codes: ", i.symleg(getOption("rlv.symbols")),"\n"),
if(u.true(legend["show.rlv.threshold"]) && lrlv &&
length(ll <- attr(lx, "rlv.threshold")))
paste("Relevance threshold", if(length(ll)>1) "s",": ",
paste(paste(names(ll),as.character(ll), sep=" = "),collapse=", "),"\n", sep="")
)
##- ll, "; type: ",
##- attr(lx,"rlv.type"))
}
ltail <- c(lleg,
if (lismodel)
c("\n", print.modelextras(x$summary, digits=digits, print=FALSE) ))
## -----------
lshe <- getOption("show.estimate")
lshet <- if(is.logical(lshe)) lshe else length(lshe)>0
if (lshet && length(lest <- attr(x, "estimate"))) {
if (is.data.frame(lest)||!is.list(lest)) lest <- list(lest)
ltail <- c(ltail, "\nestimate:\n")
lnm <- names(lest)
for (ll in seq_along(lest)) {
lle <- lest[[ll]]
if (length(lcn <- colnames(rbind(lle))))
if ((!is.logical(lshe))&&is.data.frame(lle))
lle <- lle[, lcn%in%lshe]
if (length(lle))
ltail <-
c(ltail, lf.format(lle,
header= if (length(lnm))
paste(lnm[ll],"\n",sep=":")))
}
}
## ---
rr <- structure(lout, class=c("printInference", class(lout)),
head=lhead, tail=ltail)
if (print) print.printInference(rr)
invisible(rr)
}
## -----------------------------------------------------------------------------
print.coeftable <- #f
function(x, show, print=TRUE,
digits = getOption("digits.reduced"),
transpose.ok = TRUE, na.print = getOption("na.print"), ...)
{
lf.mergesy <- function(lx, x, var, place, concatenate=FALSE)
{
lxx <- x[,var]
ltypep <- substring(var,1,3)=="p.v"
lxx[is.na(lxx)] <- ltypep
lsymb <- if(ltypep) p.symbols else rlv.symbols
ls <- symnum(lxx, lsymb$cutpoint, lsymb$symbol)
lsy <- format(c(" ",ls))[-1]
## trick needed to ensure flush left priniting of the symbols
lcl <- match(c(var, place, "coef", "coef..sy "), colnames(lx), nomatch=0)
lcl <- lcl[lcl>0][1]
if (concatenate) {
rr <- lx
rr[,lcl] <- paste(formatNA(lx[,lcl], na.print=" .", digits=digits),lsy)
names(rr)[lcl] <- paste(names(lx)[lcl], "..sy ", sep="")
} else {
rr <- if(lcl<ncol(lx))
cbind(lx[,1:lcl,drop=FALSE], lsy, lx[,-(1:lcl),drop=FALSE])
else cbind(lx[,1:lcl,drop=FALSE], lsy)
names(rr)[lcl+1] <-
if (ltypep) "p.symbol" else paste(var, "symbol", sep=".")
## paste(substring(var,1,1),if (!ltypep) "R", "sy",sep="")
}
## attr(rr, if(ltypep) "pLegend" else "rlvLegend") <- attr(ls, "legend")
rr
}
## -------------------------------------------------------------
if (length(x)<=1) return(format(x)) ## usually this applies to "(Intercept)"
x <- rbind(x)
show[show=="estimate"] <- "coef"
colnames(x)[colnames(x)=="estimate"] <- "coef"
lshow <- intersect(show, c(colnames(x), relevance.symbnames))
##- if (any(li <- !lshow%in%c(colnames(x),lsymbnm,"test")))
##- warning(":print.inference: ", paste(lshow[li], collapse=", "),
##- " not available")
lcols <- intersect(lshow, colnames(x))
if (length(lcols)==0) {
warning(":print.inference: no existing columns selected")
return(invisible(x))
}
lx <- as.data.frame(x)[,lcols, drop=FALSE]
## --- round some columns to 3 digits
ljrp <- lcols[pmatch(c("R2","p.v"), lcols, nomatch=0)]
if (length(ljrp)) lx[,ljrp] <- round(lx[,ljrp],digits) ## as.matrix()
ljrp <- lcols[c(grep("Rl", lcols),grep("Sig", lcols))]
if (length(ljrp)) lx[,ljrp] <- round(lx[,ljrp],i.last(digits)-1)
## --- paste symbols to numbers
lIsymb <- FALSE
if ("p.symbol"%in%lshow & "p.value"%in%colnames(x)) {
lx <- lf.mergesy(lx, x, "p.value", "Sig0", concatenate=TRUE)
lshow <- setdiff(lshow, "p.symbol")
lIsymb <- TRUE
}
li <- grep("Rls.symbol", lshow)
if (length(li)) {
lrs <- lshow[li]
lrc <- sub(".symbol","",lrs)
for (lii in seq_along(li)) {
lvar <- lrc[lii]
if (!lvar%in%names(x)) {
warning(":print.inference: column ",lvar, " not available")
next
}
lx <- lf.mergesy(lx, x, lvar, c("dropRls", "coefRls", "coef"),
concatenate=TRUE) ## lvar%in%lshow
## should be FALSE for something like Rlp.symbol if Rlp is not in show
lshow <- setdiff(lshow, lvar)
}
lIsymb <- TRUE
}
## result
rr <-
if (transpose.ok && ncol(lx)==1)
setNames(lx[[1]], paste(row.names(lx), if (lIsymb) " "))
else setNames(formatNA(lx, na.print=na.print, digits=digits), names(lx))
## apply(format(lx), 2, function(x) sub("NA", na.print, x))
if (print) print.printInference(rr, quote=FALSE, ...)
invisible(rr) ## structure(rr, show=lshow)
}
## --------------------------------------
print.printInference <- #f
function(x, ...)
{
ltail <- NULL
if (is.list(x)&!is.data.frame(x)) { ## print global head
if (length(lt <- attr(x,"head")))
cat(lt, "\n", sep="")
##---
ltail <- attr(x,"tail")
} else x <- list(x)
lInam <- length(lnam <- names(x))
## -------------------------------------
for (li in seq_along(x)) {
if (lInam) cat("$", lnam[li], "\n")
lx <- x[[li]]
class(lx) <- setdiff(class(lx), c("printInference", "inference"))
## attr(lx, "tail") <- NULL
## print
if (length(lhead <- attr(lx,"head"))) cat(lhead, "\n", sep="")
## attr(lx, "head") <- NULL
if (length(dim(lx))) print(as.data.frame(lx), ...)
## as.data.frame justifies the (character) columns correctly
else {
if(length(names(lx))) print(c(lx), quote=FALSE, ...)
else cat(if(is.character(lx)) lx else format(lx), "\n", sep="")
}
if (length(ltl <- attr(lx,"tail"))) cat(ltl) ## cat("\n", ltl, sep="")
## cat("\n")
}
if (length(ltail)) cat(paste(ltail, collapse="")) ## cat("\n", ltail, sep="")
## cat("\n")
}
## -----------------------------------------------------------
i.getshow <- #f
function(show, type=c("simple", "terms", "termeffects"), x=NULL)
{ ## collect items to be shown by print.inference
if (identical(show, "all")) {
if(0==length(lcn <- colnames(x))) {
warning(":getshow: no columns found. Using default")
show <- NULL
} else return(lcn)
}
lcoll <- intersect(show, c("test", "relevance", "classical"))
if(length(lcoll)) {
Nms <- c("simple", "terms", "termeffects")
ltype <- pmatch(type, Nms, nomatch=0)
## if (length(ltype)==0)
lc <- paste("show", Nms[ltype], lcoll, sep=".")
for (l in lc)
show <- c(show, getOption(l)) # or (as below!) ?? c(getOption(l), show)
}
if (all(c("nocoef","coef")%nin%show)) show <- c("coef", show)
## if (any(c("teststatistic","p.value","Sig0")%in%show)) show <- c(show,"test")
setdiff(show,lcoll)
}
## =============================================================================
getcoeftable <-
function (object)
{ ## get coefficient table from model fit
if (inherits(object, "regr")) ltb <- object$coeftable
else {
if (inherits(object, "survreg")) {
ltb <- summary(object)$table
## ltb <- ltb[-nrow(ltb),]
} else {
if (inherits(object, "coxph")) {
lcoef <- object$coefficients
se <- sqrt(diag(object$var))
ltb <- cbind(lcoef, se, lcoef/se, ## exp(lcoef),
pchisq((lcoef/se)^2, 1, lower.tail = FALSE))
dimnames(ltb) <-
list(names(lcoef), c("coef", "se", "z", "p")) ## "exp(coef)",
} else ltb <- summary(object)$coefficients
}
ldf <-
if (class(object)[1]%in%c("rlm", "coxph"))
length(object$residuals)-length(object$coef)
else df.residual(object)
if (inherits(object, "polr")) ltb <- ltb[names(object$coefficients),]
}
lnm <- names(coefficients(object))
if (any(is.na(match(lnm,row.names(ltb))))) {
rr <- structure(matrix(NA, length(lnm), ncol(ltb)), dimnames=list(lnm,colnames(ltb)))
rr[row.names(ltb),] <- ltb
} else rr <- ltb
structure(rr, df=ldf)
}
## --------------------------------------------------------------------------------
print.termtable <- #f
function (x, show = getOption("show.inference"), ...)
{
show <- i.getshow(show, "terms", x=x)
print.inference(x, show=show, ...)
}
## --------------------------------------------------------------------------------
print.termeffects <- #f
function (x, show = getOption("show.inference"), transpose.ok=TRUE,
single=FALSE, print = TRUE, warn = TRUE, ...)
{
lshowall <- length(show)==1 && show=="all"
show <- i.getshow(show, "termeffects", x=x)
lx <- if (single) x else x[sapply(x, function(x) NROW(x)>1)]
if (length(lx)==0) {
if (warn) warning(":print.termeffects: no termeffects",
if (!single) " with >1 degree of freedom")
return()
}
for (li in seq_along(lx)) {
lxx <- lx[[li]]
lsh <- if(lshowall) c(colnames(lxx), relevance.symbnames) else show
lr <-
if (length(lxx)>1) {
print.inference(lxx, show=lsh, transpose.ok=transpose.ok, print=FALSE, ...)
} else format(lxx) ## e.g., "(Intercept)"
ltail <- attr(lr,"tail")
attr(lr, "head") <- attr(lr,"tail") <- NULL
if (length(lr)==1) lr <- paste(" ", lr)
lx[[li]] <- lr
}
rr <- structure(lx, class="printInference",
head=paste(attr(x, "head"),"\n"), tail=ltail)
if (print) print.printInference(rr, ...)
invisible(rr)
}
## -------------------------------------------------------------
print.modelextras <- #f
function(x, digits=getOption("digits"), print=TRUE, ...)
{
lf.form <- function(x, header="", digits=digits)
paste(header,": ", formatC(x, digits = digits), sep="")
rr <- NULL
## error
rdf <- x$df[2] ## df.residual(x)
if (length(lsig <- x$scatter) && !u.true(attr(lsig,"fixed")))
rr <-
c(rr, scatter = lf.form(lsig, "St.dev.error", digits = digits),
df.residual = paste(" on ", rdf, " degrees of freedom") )
if (length(ldp <- x$dispersion))
rr <-
c(rr, dispersion =
lf.form(ldp, paste("dispersion parameter: ",
if (u.true(attr(ldp,"fixed"))) "fixed at "),
digits=digits) )
if (length(lsc <- x$scatter))
rr <-
c(rr, scatter =
lf.form(lsc, paste("shape parameter ('scatter')",
if (u.true(attr(ldp,"fixed"))) "fixed at "),
digits=digits) )
## if (length(lout)) rr <- c(rr, lout, "\n")
rr <- c(rr, "\n")
lr2 <- x$r.squared
if (length(lr2)&&!is.na(lr2))
rr <-
c(rr,
R2 = paste(
c(
lf.form(lr2, "Multiple R^2", digits=digits),
if (length(lr2a <- x$adj.r.squared))
lf.form(lr2a, "Adjusted R^2", digits = digits) ),
collapse="; "),"\n")
if (length(lAIC <- x$AIC)&&!is.na(lAIC))
rr <- c(rr, AIC =
lf.form(lAIC, "AIC", digits = log10(abs(lAIC))+3) )
if (length(lfst <- x$fstatistic)>0)
rr <-
c(rr, fstatistic =
paste(lf.form(lfst[1], "F-statistic",digits = digits),
" on", x$fstatistic[2], "and", x$fstatistic[3],
"d.f.; ",
lf.form(pf(x$fstatistic[1], x$fstatistic[2], x$fstatistic[3],
lower.tail=FALSE), "p.value", digits = digits),
"\n") )
## cat("\nDistribution: ", x$distrname)
if (print) print(rr, quote=FALSE, ...)
invisible(rr)
}
## =========================================================================
i.prep.plinference <- #f
function(x, overlap)
{
if (is.list(x)&&!is.data.frame(x))
return(lapply(x, i.prep.plinference, overlap=overlap))
if (is.null(dim(x))) x <- rbind(x)
lnm <- colnames(x)
lj <- match(c("Rle","Rls","Rlp"), lnm)
if (any(is.na(lj)))
lj <- match(c("coefRle","coefRls","coefRlp"), lnm)
if (any(is.na(lj)))
stop("!plot.inference! Rle, Rls, Rlp not found")
x <- x[,lj, drop=FALSE]
if (all(is.na(x)))
stop("!plot.inference! no finite values for Rle, Rls, Rlp")
if (nrow(x)>1 & overlap) {
loverlapfactor <-
if (nrow(x)>2)
0.707 else { ## only two intervals to compare
lqse <- abs(x[,3]-x[,1])
sqrt(sum(lqse^2))/sum(lqse)
}
x <- cbind(x, x[,1]+loverlapfactor*(x[,2:3]-x[,1]) )
}
x
}
## -----------------------------------------------------------------
plot.inference <- #f
function(x, pos = NULL, overlap = FALSE, refline = c(0, 1, -1),
xlab="relevance", ...) ## sub=NULL,
{
x <- i.prep.plinference(x, overlap)
plconfint(x, pos=pos, xlab = "relevance", ...)
if (length(refline))
abline(v=refline, lwd=i.def(attr(refline, "lwd"),2),
col=i.def(attr(refline, "col"), "gray70"))
invisible(x)
}
## ----------------------------------------------------------------
plconfint <- #f
function(x, y = NULL, select=NULL, overlap = NULL, pos = NULL,
xlim = NULL, refline = 0, add = FALSE, bty = "L", col = NULL,
plpars=list(lwd=c(2,3,1,4,2), posdiff=0.35,
markheight=c(1,0.6,0.6),
extend=NA, reflinecol = "gray70"),
label = TRUE, label2 = NULL,
xlab="", ...)
{
lcheck <-
list(x=list(cnr(na.ok=FALSE, dim=c(1,3)), cls()), y=cnr(),
select=cnr(range=c(0,Inf)), overlap=clg(), pos=cnr(),
xlim=cnr(), add=clg(),
bty=cch(), col=ccl(), plpars=cls(), xlab=cch()
)
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------
i.extendrange <-
function(range, ext=0.04) range + c(-1,1)*ext*diff(range)
## ---
if (is.list(x) && !is.data.frame(x)) {
if (length(y))
warning(":plconfint: argument 'y' ignored because 'x' is a list")
y <- x[[2]]
x <- x[[1]]
}
if (inherits(x, "inference"))
x <- x[,c("Rle","Rls","Rlp")] ## x[,c("effect", "se")]
lx <- as.matrix(rbind(x))
if (ncol(lx)==2) lx <- lx[,1] + outer(lx[,2], c(0,-1,1))
if (length(select)) lx <- lx[select,]
lnx <- nrow(lx)
## ---
lposlb <- lpos <- if(length(pos)) pos else lnx:1
lcol <- rep(i.def(col,1), length=lnx)
## ---
llb <- label
if (u.true(llb)) llb <- row.names(lx)
ljlb <- length(llb)>0
llb2 <- label2
ljlb2 <- length(llb2)>0
## ---
if (length(y)) {
if (inherits(y, "inference")) y <- y[,c("Rle","Rls","Rlp")]
ly <- as.matrix(rbind(y))
if (length(select)) ly <- ly[select,]
lny <- nrow(ly)
if (lny!=lnx)
stop("!plconfint! arguments 'x' and 'y' do not match")
if (ncol(ly)==2) ly <- ly[,1] + outer(ly[,2], c(0,-1,1))
if (u.notfalse(overlap)) {
lxw <- (lx[,3]-lx[,2])/2
lyw <- (ly[,3]-ly[,2])/2
lovfac <- sqrt(lxw^2+lyw^2)/(lxw+lyw)
lx <- cbind(lx[,1:3,drop=FALSE],
lx[,1] + lovfac*(lx[,2:3,drop=FALSE]-lx[,1]))
ly <- cbind(ly[,1:3,drop=FALSE],
ly[,1] + lovfac*(ly[,2:3,drop=FALSE]-ly[,1]))
}
if (ncol(lx)!=ncol(ly))
stop("!plconfint! number of columns of 'x' and 'y' are different")
lx <- rbind(lx,ly)[c(outer(c(0,lnx), 1:lnx, "+")),]
if (length(lpos)==lnx) {
lmd <- if (lnx>2) min(diff(sort(lpos))) else 1
lpos <- c(outer(lmd*plpars$posdiff/2*c(1,-1), lpos, "+"))
}
lnx <- 2*lnx
lcol <- rep(i.def(col,c("blue","red")), length=lnx)
}
## ---
if (nrow(lx)==0)
stop("!plconfint! no elements left for argument ''x'")
## ---
l2ci <- ncol(lx)>=5
if (!u.true(single)) {
li <- !is.na(lx[,1])
lx <- lx[li,, drop=FALSE]
lpos <- lpos[li]
lcol <- lcol[li]
}
if (length(lx)==0L)
stop("!plconfint! no intervals to plot")
lwd <- rep(i.def(plpars[["lwd"]],2), length=5)
lmh <- 0.1 * rep(c(plpars[["markheight"]],1),length=3) ## * diff(range(lpos))/ln
if (!add) {
## rr <- replication(x,y)
## range
lrg <- i.extendrange(range(c(lx), na.rm=TRUE))
if (length(xlim)) {
if (length(xlim)!=2)
warning(":plconfint: argument 'xlim' not suitable")
else lrg[!is.na(xlim)] <- xlim[!is.na(xlim)]
}
lxt <- i.def(plpars[["extend"]], 1/lnx)
lposlim <- matrix(c(1+lxt, -lxt, -lxt, 1+lxt),2)%*%range(lpos)
lmar <- par("mar")
if (ljlb) lmar[2] <- 0.7*max(nchar(llb))+1
if (ljlb2) lmar[4] <- 0.7*max(nchar(llb2))+1
loldp <- par(mar=lmar)
on.exit(par(loldp))
plot(lrg, lposlim, yaxs="i", type="n", axes=FALSE,
xlab=xlab, ylab="", xaxs="i", yaxs="i", ...) # c(min(0,lrg[1]), max(1,lrg[2]))
box(bty=bty)
axis(1)
}
lrlcol <- i.def(attr(refline, "col"), plpars[["reflinecol"]])
if (length(refline)) abline(v=refline, col=lrlcol)
segments(lx[,2],lpos, lx[,3],lpos, lwd=lwd[1], col=lcol) ## interval line
segments(lx[,1],lpos-lmh[1],lx[,1],lpos+lmh[1], lwd=lwd[2], col=lcol) ## midpoint = estimate
segments(lx[,2], lpos-lmh[2], lx[,2], lpos+lmh[2],
lwd=lwd[3], col=lcol) ## endmarks
segments(lx[,3], lpos-lmh[2], lx[,3], lpos+lmh[2],
lwd=lwd[3], col=lcol) ## endmarks
if (l2ci) {
segments(lx[,4],lpos, lx[,5],lpos, lwd=lwd[4], col=lcol) ## interval line
segments(lx[,4], lpos-lmh[2], lx[,4], lpos+lmh[2],
lwd=lwd[5], col=lcol) ## endmarks
segments(lx[,5], lpos-lmh[2], lx[,5], lpos+lmh[2],
lwd=lwd[5], col=lcol) ## endmarks
}
## ---
if (ljlb) mtext(llb, side=2, at=lposlb, line=1, adj=1, las=1)
if (ljlb2) mtext(llb2, side=4, at=lposlb, line=1, adj=0, las=1)
invisible(
structure(lpos, labelpos=if(length(lpos)!=length(lposlb)) lposlb))
}
## ---------------------------------------------------------------
pltwosamples <- function(x, ...) UseMethod("pltwosamples")
## ---
pltwosamples.default <- #f
function(x, y = NULL, overlap = TRUE, ...) ## , sub=":"
{
lcheck <-
list(x=list(cnr(na.ok=FALSE),cdf()), y=cnr(), overlap=clg()
)
##- lcall <- match.call(expand.dots = FALSE)
##- lcall$... <- NULL
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------
if (is.matrix(x)) x <- as.data.frame(x)
if (inherits(x, "list")) {
## stop("!pltwosamples! First argument not suitable")
lnm <- names(x)
y <- x[[2]]
x <- x[[1]]
}
if (length(y)==0)
stop("!pltwosamples! Argument 'y' not specified")
lx <- onesample(x)
ly <- onesample(y)
plconfint(c(lx[c("effect","ciLow","ciUp","se")],n=attr(lx,"n")[1]),
c(ly[c("effect","ciLow","ciUp","se")],n=attr(ly,"n")[1]),
overlap=overlap, ...)
}
## -------------------------------------------------------------
pltwosamples.formula <- #f
function(formula, data=NULL, ...) ## pos = NULL, col=1,
{
if (length(formula) != 3L)
stop("!twosamples! 'formula' must have left and right term")
oneSampleOrPaired <- FALSE
if (length(attr(terms(formula[-2L]), "term.labels")) != 1L)
if (formula[[3]] == 1L)
oneSampleOrPaired <- TRUE
else stop("!twosamples! 'formula' incorrect")
m <- match.call(expand.dots = FALSE)
if (is.matrix(eval(m$data, parent.frame())))
m$data <- as.data.frame(data)
m[[1L]] <- quote(stats::model.frame)
m$... <- NULL
mf <- eval(m, parent.frame())
ll <- split(mf[,1],mf[,2])
pltwosamples.default(ll, ...)
}
## ---------------------------------------------------------------
plot.termeffects <- #f
function(x, pos = NULL, single=FALSE, overlap = TRUE,
termeffects.gap = 0.2, refline = c(0, 1, -1),
xlim=NULL, ylim=NULL, xlab = "relevance", mar=NA,
labellength = getOption("labellength"), ...)
{
lcheck <-
list(x=cnr(na.ok=FALSE), single=clg(), overlap=clg(),
termeffects.gap=cnr(range=c(-1,10)),
labellength=cnr(range=c(0,50))
)
##- lcall <- match.call(expand.dots = FALSE)
##- lcall$... <- NULL
largs <- check.args(lcheck, envir=parent.frame())
for (lnm in names(largs)) assign(lnm, largs[[lnm]])
## --------------
mardefault <- par("mar")
li <- sapply(x, is.atomic)
x <- x[!li] ## (Intercept)
lx <-
lapply(x, function(x) x <- x[any(is.finite(x[,"effect"])),,drop=FALSE])
llen <- sapply(lx, nrow)
lx <- lx[llen>0]
llen <- llen[llen>0]
if (!single) {
lx <- lx[llen>1]
llen <- llen[llen>1]
}
if (length(llen)==0)
stop("!plot.termeffects! No termeffects", if(!single) " with length >1")
lnx <- length(lx)
lnmx <- names(lx)
lxx <- matrix(,0,3)
llb <- NULL
for (ll in seq_along(lx)) {
llx <- rbind(lx[[ll]])
if (llen[ll] <- nrow(llx)) {
lxs <- llx[,c("coefRle","coefRls","coefRlp")]*sign(llx[,"estimate"])
lxx <- rbind(lxx, lxs)
llb <- c(llb, if(llen[ll]>1) row.names(llx) else lnmx[ll])
}
}
if (is.null(pos)) {
pos <-
rep(termeffects.gap*(1:length(llen)-(lnx==1))+
cumsum(llen>1), llen) + 1:sum(llen)
pos <- max(pos)+1-pos
}
if (length(xlim)==0) xlim <- range(unlist(lxx), na.rm=TRUE)
if (length(ylim)==0) {
lrg <- range(pos) + c(0, llen[1]>1 & lnx>1)
ylim <- i.extendrange(lrg, 0.3*getOption("plottermeffects.extend")/diff(lrg) )
}
if (length(mar)==1) mar <- c(NA, mar, NA,NA)
if (length(mar)&&is.na(mar[2])) {
lmr2 <- max(nchar(llb))
if (lmr2>labellength) {
llb <- sub(")","", sub("log.*\\((.*)\\)","l:\\1", llb))
llb <- shortenstring(llb, labellength)
lmr2 <- labellength
}
mar[2] <- lmr2*0.7 + 1
}
mar <-
if (length(mar))
ifelse(is.finite(mar), mar, mardefault) else mardefault
## ---
lop <- par(mar=mar)
plot(xlim, ylim, xlab=xlab, ylab="", axes=FALSE, type="n")
on.exit(par(lop))
axis(1)
box()
lx0 <- xlim[1]
li0 <- 0
lnmx <- names(x)
for (lt in 1:lnx) {
ltn <- llen[lt]
li <- li0+1:ltn
llx <- lxx[li,]
row.names(llx) <- row.names(lx[[lt]])
if (ltn>1&lnx>1)
text(lx0, pos[li0+1]+1, lnmx[lt], adj=0)
else row.names(llx) <- llb[li]
plconfint(llx, pos=pos[li], add=TRUE)
li0 <- li0+ltn
}
if (length(refline))
abline(v=refline, lwd=i.def(attr(refline, "lwd"),2),
col=i.def(attr(refline, "col"), "gray70"))
}
## =============================================================================
i.logscale <- #f
function(object)
{
substring(as.character(formula(object)[[2]]),1,3)=="log" ||
(inherits(object, "glm") &&
object$family%in%
c("binomial","poisson","quasibinomial","quasipoisson")) ||
inherits(object,"survreg")&&object$dist!="gaussian"
}
## --------------------------------------------------------------------
getscalepar <- #f
function(object)
{ ## get scale parameter of a fit
lsry <- summary(object)
rr <- c(lsry$sigma, lsry$scatter)[1]
if (length(rr)==0) rr <- sqrt(c(lsry$dispersion,1)[1])
rr
}
## -----------------------------------------------------------
getcoeffactor <- #f
function(object, standardize=TRUE)
{
## get factor for converting coef to coef effect
## model matrix
lmmt <- object[["x"]]
if (length(lmmt)==0) object$x <- lmmt <- model.matrix(object)
lfamily <- object$family$family
ldist <- object$dist
lscatter <- if (standardize) getscalepar(object) else 1
lfac <- apply(lmmt, 2, sd)/lscatter
lfac[lfac==0] <- NA
structure(lfac, scatter=lscatter, fitclass=class(object),
family=lfamily, dist=ldist)
}
## -----------------------------------------------------------
rlvClass <- #f
function(effect, ci=NULL, relevance=NA)
{ ## relevance class !!! ciwidth - interval !!!
lrlvth <- i.def(relevance, getOption("rlv.threshold")[1])
if (inherits(effect, "inference")) {
if (length(ci))
warning(":rlvClass: argument 'ci' ignored")
rr <- attr(effect, "rlvclass")
if (length(rr)) return(rr)
rle <- effect["Rle"]
rls <- effect["Rls"]
rlp <- effect["Rlp"]
} else {
if (length(ci)==length(effect))
ci <- effect + outer(c(ci), c(-1,1))
if (length(unlist(ci))!=2*length(effect))
stop("!rlvClass! lengths of 'estimate' and 'ci' not compatible")
rle <- effect/lrlvth
lcis <- rbind(ci/lrlvth)
rls <- lcis[,1]
rlp <- lcis[,2]
names(rls) <- i.def(row.names(effect),
i.def(row.names(lcis), NULL) )
}
rr <- ifelse(rls>=1, "Rlv", "Amb.Sig")
rr <- ifelse(rls<0, "Amb", rr) ## rls may be a vector
rr <- ifelse(rlp<1, "Ngl", rr)
rr <- ifelse(rlp<0, "Ctr", rr)
rr
}
## ====================================================================
relevance.modelclasses <-
c("regr","lm","lmrob","glm","polr","survreg","coxph")
## ,"rlm" : no correlation matrix of coef
## ,"rq"
p.symbols <- list(symbol=c("***", "**", "*", ".", " "),
cutpoint=c(0, 0.001, 0.01, 0.05, 0.1, 1) )
rlv.symbols <- list(symbol=c(" ", ".", "+", "++", "+++"),
cutpoint=c(-Inf,0,1,2,5,Inf) )
i.symleg <-
function(x) paste(c(rbind(as.character(x$cutpoint),c(x$symbol,""))),
collapse=" ")
relevance.symbnames <-
c("p.symbol", "coefRls.symbol", "dropRls.symbol", "predRls.symbol")
## -----------------------------------------------------------
relevance.options <- list(
digits.reduced = 3,
testlevel = 0.05,
rlv.threshold =
c(stand=0.1, rel=0.1, prop=0.1, corr=0.1, coef=0.1, drop=0.1, pred=0.05),
termtable = TRUE,
show.confint = TRUE,
show.estimate =
c("estimate","scatter", "n", "effect","Rle","Rls","Rlp",
"rlvclass","rplclass"),
show.doc = TRUE,
show.inference = "relevance",
show.simple.relevance = c("Rle", "Rlp", "Rls", "Rls.symbol",
"rlvclass", "rplclass"),
show.simple.test = c("Sig0", "p.value", "p.symbol"),
show.simple.classical = c("teststatistic", "p.value", "p.symbol"), ## !!! symbols?
show.terms.relevance = c("df", "R2x", "coefRlp", "coefRls", ## "dropRle",
"dropRls", "dropRls.symbol", "predRle"),
show.terms.test = c("df", "ciLow","ciUp", "R2x", "Sig0", "p.value",
"p.symbol"),
show.terms.classical = c("df", "se", "teststatistic", "p.value", "p.symbol"),
show.termeffects.relevance = c("coef","coefRls.symbol"),
show.termeffects.test = c("coef","p.symbol"),
show.termeffects.classical = c("coef","p.symbol"),
show.symbollegend = TRUE, show.rlv.threshold = TRUE,
show.Rlv.class = TRUE,
plottermeffects.extend = 1,
labellength = 8,
na.print = ". ",
p.symbols = p.symbols,
rlv.symbols = rlv.symbols
)
.onLoad <- function(lib, pkg) options(relevance.options)
## =============================================================================
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.