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=TRUE,
rlv.threshold=getOption("rlv.threshold"), ...)
{ ## effect (group difference) and relevance
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 <- NULL
} 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 <- NULL
lest <- NULL
}
## ---
if (lbin) { ## binary data
lrlvth <-
i.def(rlv.threshold["prop"],
c(rlv.threshold, relevance.options[["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
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])
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))
leff <- unname(log(lt$estimate))
## lse <- # !!!
leffci <- log(lt$conf.int) ## c("ci.low","ci.up") ))
lsig <- NA
effname <- "log odds ratio"
method <- "Two proportions, Fisher's exact inference"
}
lrlvci <- c(leff, leffci)/lrlvth
ltst <- lt$statistic
lpv <- lt$p.value
ltype <- "proportion"
lsigth <- NA
ldf <- lv <- NULL
} else { ## quantitative data
log <- i.def(log, FALSE, TRUE, FALSE)
llogrescale <- 1
if (is.character(log)) {
if(log=="log10") llogrescale <- 1/log(10)
log <- substr(log,1,3)=="log"
}
lirl <- ifelse(log,"rel","stand")
lrlvth <-
i.def(rlv.threshold[lirl],
c(rlv.threshold, relevance.options[["rlv.threshold"]]["stand"])[1])
lssx <- sum((x-lmnx)^2)
lpq <- 1-ltlev2
if (lonegroup) { ## one quantitative sample
leff <- lmnx
ldf <- lnx-1
lv <- lssx/ldf
lny <- NULL
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 (log) {
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)]
structure(
c(effect=leff, se=lse, statstic=ltst, p.value=lpv, Sig0=lsig,
ciLow=leffci[1], ciUp=leffci[2],
Rle=lrlvci[1], Rls=lrlvci[2], Rlp=lrlvci[3],
Sigth=lsigth),
class = c("inference", "twosamples"), type="simple",
method = method, effectname = effname, hypothesis = hypothesis,
n = c(lnx, lny), estimate=lest, statistic = ltst, V = lv, df = ldf,
data = list(x=x, y=y),
rlv.threshold = lrlvth,
rlv.type = ltype
)
}
## ============================================================================
twosamples.formula <-
function (formula, data, subset, na.action, ...)
{ ## adapted from t.test.formula
if (missing(formula) || (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())
## 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")
DATA <- setNames(split(mf[[response]], g), c("x", "y"))
rr <- do.call("twosamples", c(DATA, paired=FALSE, list(...)))
names(attr(rr, "n")) <- levels(g)
} else { ## one sample
respVar <- mf[[response]]
if (inherits(respVar, "Pair")) {
DATA <- list(x = respVar[, 1], y = respVar[, 2])
rr <- do.call("twosamples", c(DATA, paired = TRUE, list(...)))
} else {
DATA <- list(x = respVar)
rr <- do.call("twosamples", c(DATA, paired = FALSE, list(...)))
}
}
## attr(rr, "data.name") <- DNAME
ldn <- substitute(data)
structure(rr, formula=formula, data.name=if (is.name(ldn)) format.default(ldn))
}
## -------------------------------------------------------------------------
twosamples.table <- #f
function(x, ...) twosamples.default(table=x, ...)
## ========================================================================
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)
}
## ===========================================================================
##
inference <-
function (estimate=NULL, se=NULL, n=NULL, df=Inf, testlevel=0.05, stcoef=TRUE,
rlv=TRUE, rlv.threshold=0.1, object=NULL)
{ ## for a coefficients table,
## calculate confidence interval, significance and relevance
if (length(estimate)==0) estimate <- object
ln <- i.def(n, NA)
ldf <- df
ldfm <- 1
if (inherits(estimate, relevance.modelclasses)) {
object <- estimate
estimate <- getcoeftable(object)
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, ...
}
df <- i.def(ldf, ln-1, valuefalse=NA)
ln <- i.def(ln, ldf+ldfm)
if (!(is.atomic(estimate)||length(dim(estimate))==2))
stop("!inference! first argument not suitable")
if (is.null(se))
if (NCOL(estimate)>1) {
se <- estimate[,2]
estimate <- structure(estimate[,1], names=row.names(estimate))
}
if (is.null(se)) se <- attr(estimate, "se")
if (is.null(se))
stop("!inference! no standard errors found")
if (df==0||!any(is.finite(se))) {
warning("!inference! no finite standard errors")
return(cbind(estimate))
}
ltq <- if (is.finite(ldf)) qt(1-testlevel/2, ldf) else qnorm(1-testlevel/2)
lci <- estimate+outer(ltq*se, c(ciLow=-1,ciUp=1))
ltst <- estimate/se
lestst <- ltst*sqrt(ln)
lsgf <- ltst/ltq
lpv <- 2*pt(-abs(ltst), df)
rr <- data.frame(estimate=estimate, se=se, est.st=lestst,
statistic=ltst, p.value=lpv, Sig0=lsgf, lci)
## --- 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 <- estimate*lfac[names(estimate)]
}
if (length(stcoef)!=length(estimate)) {
warning(":inference: argument 'stcoef' not suitable. No relevances")
return (rr)
}
lfac <- stcoef/estimate
lstci <- cbind(stcoef=stcoef, stciLow=lfac*lci[,1], stciUp=lfac*lci[,2])
rr <- cbind(rr, lstci)
## --- 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 <- lstci/rlv.threshold
li <- which(lstci[,1]<0)
if (length(li)) lrlv[li,] <- - lrlv[li,c(1,3,2)]
colnames(lrlv) <- c("Rle","Rls","Rlp")
cbind(rr, lrlv)
}
## ===========================================================================
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
## --------------------------------------------------------
## thresholds
ltlev1 <- 1-testlevel
## lrlthrl <- rlv.threshold[["rel"]]
## lrlthcf <- rlv.threshold[["coef"]]
lrlthdr <- rlv.threshold[["drop"]]
lrlthpr <- rlv.threshold[["pred"]]
## --- sigma and threshold for standardized coefficient
## lsigma <- c(object$sigma, summary$sigma, 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,
statistic=NA, p.value=NA, Sig0=NA, ciLow=NA, ciUp=NA,
stcoef=NA, stciLow = NA, stciUp = NA, R2x=NA,
stringsAsFactors=FALSE)
)
## --- coefficients
lcoef <- object$coefficients
lcoeftab <- inference(object=object)
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, statistic=ltst, p.value=lpv,
Sig0=sqrt(pmax(0,ltst)/ltstq), ciLow=NA, ciUp=NA,
stcoef=NA, stciLow=NA, stciUp=NA,
testst=ltst, 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","statistic","ciLow","ciUp","stcoef", "stciLow","stciUp",
"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)
}
## end termtable
## --------------------------------------------------------------------
i.drop1 <- #F
function (object, scope=drop.scope(object), scale = 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, scale=scale, 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
}
## ===========================================================
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
}
## ===========================================================================
termeffects <- #f
function (object, se = 2, df = df.residual(object), rlv = TRUE)
## --------------------------------------------------------------
{
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) lsigma <- 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))) {
##- warning(":termeffects: missing coef or non-finite standard error for term '",
##- tl[j], "'. no standard errors etc")
ljfail <- c(ljfail, tl[j])
} else {
lsigma <- getscalepar(object)
rrj <- inference(rrj, sej, df=df, stcoef=rrj*0.5/lsigma, 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)
} ## end termeffects
## -------------------------------------------------------------------------
"[.termeffects" <- #f
function(x, i=NULL) structure(unclass(x)[i], class="termeffects")
## unclass avoids infinite recursion!
## -------------------------------------------------------------------------
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.mergesy <- function(lx, x, var, place)
{
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"), colnames(lx), nomatch=0)
lcl <- lcl[lcl>0][1]
rr <- if(lcl<ncol(lx)) cbind(lx[,1:lcl,drop=FALSE], lsy, lx[,-(1:lcl),drop=FALSE])
else cbind(lx[,1:lcl,drop=FALSE], rsy=lsy)
names(rr)[lcl+1] <- paste(substring(var,1,1),if (!ltypep) "R", "sy",sep="")
attr(rr, if(ltypep) "pLegend" else "rlvLegend") <- attr(ls, "legend")
rr
}
lf.format <- function(x, header="")
{ ## format matrix -> character vector including names
x <- format(x)
x <-
if (length(dim(x))==0)
cbind(format(c(header, "")),
rbind(format(names(x)),x), "\n")
else
cbind(format(c(header, rep("", nrow(x)))),
format(c("",rownames(x))),
format(rbind(colnames(x), x)), "\n")
unname(apply(x, 1, paste, collapse=" "))
}
## -------------------------------------
## show what?
lsh <- c(show.signif.stars = getOption("show.signif.stars"),
show.symbollegend = getOption("show.symbollegend"),
show.rlv.threshold = getOption("show.rlv.threshold"))
legend <- i.def(legend, lsh, structure(rep(TRUE, length(lsh)), names=names(lsh)))
lItable <- length(dim(x))
type <- i.def(attr(x,"type"), "simple")
lsymbnm <-
c("p.symbol",
if (lItable) c("coefRls.symbol", "dropRls.symbol", "predRls.symbol")
else "Rls.symbol")
lshow <- if (identical(show, "all")) c(colnames(x), lsymbnm)
else i.getshow(show, type)
if (all(c("nocoef","coef")%nin%lshow)) lshow <- c("coef", lshow)
if (any(c("statistic","p.value","Sig0")%in%lshow)) lshow <- c(lshow,"test")
## ---
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: ", as.character(lfo[[2]]) ) ## was format()
)
lhead <- c(
paste(attr(x, "method"),"\n",collapse=" "),
if (length(lout)) paste(paste(lout, collapse=" ; "), "\n") ## " ; "
)
if (getOption("show.estimate") && length(lest <- attr(x, "estimate")))
lhead <- c(lhead, lf.format(lest, header="estimates:"))
## ---
lleg <- NULL
if (!lItable) {
lx <- x
lout <- c(
if (!is.na(leff <- x["effect"]))
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 ("test"%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, "statistic")))
paste("statistic: ", 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 (length(llout)) lout <- c(lout, paste(paste(llout, collapse=" ; "), "\n") )
}
} else { ## ======================== table
lshow[lshow=="estimate"] <- "coef"
colnames(x)[colnames(x)=="estimate"] <- "coef"
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(as.matrix(lx[,ljrp]),digits)
ljrp <- lcols[c(grep("Rl", lcols),grep("Sig", lcols))]
if (length(ljrp)) lx[,ljrp] <- round(as.matrix(lx[,ljrp]),i.last(digits)-1)
## --- paste symbols to numbers
if ("p.symbol"%in%lshow & "p.value"%in%colnames(x))
lx <- lf.mergesy(lx, x, "p.value", "Sig0")
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"))
}
}
lout <-
if (transpose.ok &&
((lnc1 <- ncol(lx)==1)|| ncol(lx)==2 && length(grep(".symbol", lshow))) ) {
format(
if (lnc1) lx[[1]]
else setNames(paste(format(lx[[1]]),lx[[2]]),
paste(row.names(lx)," ")))
} else
apply(format(lx), 2, function(x) sub("NA", na.print, x))
}
## --- legend(s)
if(u.notfalse(legend)) {
lleg <- c(
if(u.true(legend["show.signif.stars"]) &&
length(ll <- attr(lx, "pLegend")))
paste("\nSignificance codes for p.value: ", ll),
if(u.true(legend["show.symbollegend"]) &&
length(ll <- attr(lx, "rlvLegend")))
paste("\nRelevance codes: ", ll),
if(u.true(legend["show.rlv.threshold"]) &&
length(ll <- attr(lx, "rlv.threshold")))
paste("\nRelevance threshold: ", ll, "; type: ",attr(lx,"rlv.type"))
)
}
rr <- structure(lout, class=c("printInference", class(lout)), head=lhead, tail=lleg)
if (print) print.printInference(rr)
invisible(rr)
}
## --------------------------------------
print.printInference <- #f
function(x, ...)
{
ltail <- NULL
if (is.list(x)&!is.data.frame(x)) {
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"))
lhead <- attr(lx,"head")
attr(lx, "head") <- NULL
ltl <- attr(lx,"tail")
attr(lx, "tail") <- NULL
if (length(lhead)) cat(lhead, "\n", sep="")
if (length(dim(lx))) print(unclass(lx), quote=FALSE)
## unclass needed because ohterwise, the class attribute is printed!
else {
if(length(names(lx))) print(c(lx), quote=FALSE) else cat(lx, sep="")
}
if (length(ltl)) cat(ltl, sep="")
## cat("\n")
}
if (length(ltail)) cat(ltail, sep="")
cat("\n")
}
## -----------------------------------------------------------
i.getshow <- #f
function(show, type=c("simple", "terms", "termeffects"))
{ ## collect items to be shown by print.inference
lcoll <- intersect(show, c("test", "relevance", "classical"))
if(length(lcoll)) {
ltype <- pmatch(type, c("simple", "terms", "termeffects"), nomatch=0)
## if (length(ltype)==0)
lc <- paste("show", c("simple","terms","termeffects")[ltype], lcoll, sep=".")
for (l in lc)
show <- c(show, getOption(l))
}
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
}
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
rr
} else ltb
}
## --------------------------------------------------------------------------------
print.termtable <- #f
function (x, show = getOption("show.inference"), ...)
{
show <- if (length(show)==1 && show=="all") colnames(x)
else unique(c("coef", i.getshow(show, "terms")))
print.inference(x, ...)
}
## --------------------------------------------------------------------------------
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 <- unique(c("coef", i.getshow(show, "termeffects")))
lx <- x[sapply(x, function(x) NROW(x)>1-single)]
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) colnames(lxx) 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)
}
## =========================================================================
plot.inference <- #f
function(x, pos = NULL, overlap = FALSE, reflines = c(0, 1, -1),
xlab="relevance", ...) ## sub=NULL,
{
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 (nrow(x)>1 & overlap) {
loverlapfactor <-
if (nrow(x)>2)
0.707 else { lqse <- abs(x[,3]-x[,1])
sqrt(sum(lqse^2))/sum(lqse)
}
x <- cbind(x, x[,1]+loverlapfactor*(x[,2:3]-x[,1]) )
}
plconfint(x, xlab = "relevance", ...)
if (length(reflines))
abline(v=reflines, lwd=i.def(attr(reflines, "lwd"),2),
col=i.def(attr(reflines, "col"), "gray70"))
invisible(x)
}
## -------------------------------------------------------------------------
plconfint <- #f
function(x, pos = NULL, xlim = NULL, add = FALSE, bty = "L", col = 1,
plpars=list(lwd=c(2,3,1,2,2), markheight=c(1,0.7,0.85), extend=NA,
reflinecol = "gray70"),
xlab="", ...)
{
i.extendrange <- function(range, ext=0.04) range + c(-1,1)*ext*diff(range)
lx <- as.matrix(rbind(x))
ln <- nrow(lx)
loverlap <- ncol(x)>=5
ly <- seq(ln,1)
if (length(pos)) {
if (length(pos)!=ln | any(is.na(pos)) | any(duplicated(pos)))
warning(":plconfint: unsuitable argument 'pos'")
else ly <- pos
}
lcol <- rep(i.def(col,1), length=ln)
li <- !is.na(lx[,1])
lx <- lx[li,, drop=FALSE]
ly <- ly[li]
lcol <- lcol[li]
lnmeff <- i.def(if (ln>1) row.names(lx), "")
lwd <- rep(i.def(plpars[["lwd"]],2), length=5)
lmh <- 0.1 * rep(c(plpars[["markheight"]],1),length=3) ## * diff(range(ly))/ln
if (!add) {
## 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/ln)
lylim <- matrix(c(1+lxt, -lxt, -lxt, 1+lxt),2)%*%range(ly)
lmar <- par("mar")
lnch <- max(nchar(lnmeff))
lmar[2] <- 0.7*lnch+1
loldp <- par(mar=lmar)
on.exit(par(loldp))
plot(c(min(0,lrg[1]), max(1,lrg[2])), lylim, yaxs="i", type="n", axes=FALSE,
xlab=xlab, ylab="", xaxs="i", yaxs="i", ...)
lrlcol <- plpars[["reflinecol"]]
box(bty=bty, col=lrlcol)
axis(1, col=lrlcol)
}
segments(lx[,2],ly, lx[,3],ly, lwd=lwd[1], col=lcol) ## interval line
segments(lx[,1],ly-lmh[1],lx[,1],ly+lmh[1], lwd=lwd[2], col=lcol) ## midpoint = estimate
segments(lx[,2:3],rep(ly,2)-lmh[2],lx[,2:3],rep(ly,2)+lmh[2],
lwd=lwd[3], col=lcol) ## endmarks
if (loverlap)
segments(lx[,4:5],rep(ly,2)-lmh[3],lx[,4:5],rep(ly,2)+lmh[3],
lwd=lwd[4], col=lcol) ##
## ---
mtext(lnmeff, side=2, at=ly, line=1, adj=1, las=1)
}
## ---------------------------------------------------------------
pltwosamples <- function(x, ...) UseMethod("pltwosamples")
## ---
pltwosamples.default <- #f
function(x, y, overlap = TRUE, ...) ## , sub=":"
{
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]]
}
lx <- onesample(x)
ly <- onesample(y)
lci <- rbind(lx[c("effect","ciLow","ciUp")],
ly[c("effect","ciLow","ciUp")])
lse <- c(lx["se"],ly["se"])
if (overlap) {
loverlapfactor <- sqrt(sum(lse^2))/sum(lse)
lci <- cbind(lci, lci[,1]+loverlapfactor*(lci[,2:3]-lci[,1]) )
}
plconfint(lci, ...)
}
## -------------------------------------------------------------
pltwosamples.formula <- #f
function(formula, data=NULL, ...) ## pos = NULL, col=1,
{
if (missing(formula) || (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, ...) ## , sub=":"
{
li <- sapply(x, is.atomic)
x <- x[!li] ## (Intercept)
llen <- sapply(x, nrow)
if (!single) {
x <- x[llen>1]
llen <- llen[llen>1]
}
if (length(llen)==0)
stop("!plot.termeffects! No termeffects", if(single) "with length >1")
##- if (length(llen)==1) {
##- plot.inference(x[[1]], pos=pos, ...) ## sub=sub,
##- return()
##- }
## lx <- t(sapply(x, function(x) x[,c("coefRle","coefRls","coefRlp")]))
lx <- matrix(,0,3)
for (ll in seq_along(x)) {
lxx <- x[[ll]]
lxs <- lxx[,c("coefRle","coefRls","coefRlp")]*sign(lxx[,"estimate"])
lx <- rbind(lx, lxs)
}
row.names(lx) <- unlist(lapply(x, row.names))
if (is.null(pos)) {
pos <- rep(termeffects.gap*(1:length(llen))+cumsum(llen>1), llen) + 1:sum(llen)
pos <- max(pos)+1-pos
}
plot.inference(lx, pos=pos, overlap=overlap, ...) ## sub=sub,
li <- llen>1
if (any(li)) {
lii <- c(0,i.last(cumsum(llen),-1))[llen>1]+1
mtext(names(x)[li], side=2, at=pos[lii]+1, las=1)
}
}
## ----------------------------------------------------------------
##- shortenstring <- function (x, n=50, endstring="..", endchars=NULL)
##- { ## from plgraphics
##- if (length(endchars)==0) endchars <- pmin(3,ceiling(n/10))
##- if (any(li <- 1 >= (ncut <- n-nchar(endstring)-endchars))) {
##- warning(":shortenstring: argument 'n' too small for given 'endstring' and 'endchar'")
##- endstring <- ifelse(li, ".", endstring)
##- endchars <- ifelse(li, 0, endchars)
##- ncut <- n-nchar(endstring)-endchars
##- }
##- if (length(x) && any(n < (lnc <- nchar(x))))
##- ifelse(n<lnc & ncut>1, paste(substring(x, 1, ncut), endstring,
##- substring(x, lnc-endchars+1, lnc), sep=""), x)
##- else x
##- }
##- ## ======================================================================
getscalepar <- #f
function(object)
{ ## get scale parameter of a fit
lsry <- summary(object)
sigma <- c(lsry$sigma, lsry$scale)[1]
if (length(sigma)==0) sigma <- sqrt(c(lsry$dispersion,1)[1])
sigma
}
## -----------------------------------------------------------
getcoeffactor <- #f
function(object)
{
## 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
lsigma <- getscalepar(object)
##- lsigma <-
##- if (inherits(object,"glm")&&lfamily%in%c("binomial", "quasibinomial"))
##- 1.6683*lsigma ## qlogis(pnorm(1)) ???
##- else {
##- if (inherits(object, c("lm","lmrob","rlm"))||
##- (inherits(object,"survreg")&&ldist=="gaussian"))
##- lsigma else 1
##- }
lfac <- apply(lmmt, 2, sd)/lsigma
lfac[lfac==0] <- NA
## lnm <- names(object$coefficients)
##- if (any(lna <- is.na(match(lnm,names(lfac))))) {
##- warning(":getcoeffactor: error, possibly singular case", paste(lnm[lna], collapse=", "))
##- lfac <- 1
##- } else lfac <- lfac[lnm] ## needed for singular designs
structure(lfac, sigma=lsigma, fitclass=class(object),
family=lfamily, dist=ldist)
}
## ====================================================================
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) )
## -----------------------------------------------------------
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 = TRUE, show.doc = TRUE,
show.inference = "relevance",
show.simple.relevance = c("Rle", "Rlp", "Rls", "Rls.symbol"),
show.simple.test = c("Sig0", "p.value", "p.symbol"),
show.simple.classical = c("statistic", "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", "statistic", "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,
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.