R/relevance.R

In relevance: Calculate Relevance

## 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)
## =============================================================================