R/regr.R

Defines functions get_all_vars factor2character colorpale i.multinomfit i.polrfit drop1.default add1.default rrevgumbel qrevgumbelexp qrevgumbel prevgumbel drevgumbel u.nuna IR DB warn getmeth RNAMES extractNames i.main u.merge u.debug u.true nafalse quantNA factor.na is.formula tit doc legendr plcoord asinp logst sumna notna nainf.exclude xNA factorNA last modarg userOptions getUserOption quantilew quinterpol stamp robrange mframe showd subset dropdata plTA.polr plfitpairs plres2x plmboxes plmbox panel.smooth panelDefault plmatrix smoothM smoothMM plresx simresiduals.default simresiduals.glm simresiduals smoothRegr i.plotlws plot.regr plot.xdistResscale xdistResscale xdistResdiff leverage print.modelTable format.modelTable modelTable compareTerms predict.mlm i.findformfac fitcomp linear.predictors residuals.polr nobs.coxph nobs.survreg residuals.regrsurv residuals.regrcoxph residuals.regr condquant predict.polr fitted.polr print.mregr i.add1na add1.mlm drop1.mlm summary.mregr drop1.multinom vif.regr predict.regr fitted.regr terms2order step.regr step drop1Wald add1.regr drop1.regr confint.regr df.residual.regr vcov.regr print.allcoef allcoef summary.regr print.regr contr.wsum ciSignif i.termtable Tobit i.survreg i.polr i.multinomial i.glm i.mlmsum i.lm regr

Documented in add1.mlm add1.regr allcoef asinp colorpale compareTerms condquant contr.wsum DB doc drevgumbel drop1.mlm drop1.multinom drop1.regr drop1Wald dropdata factorNA fitcomp fitted.polr fitted.regr format.modelTable getmeth getUserOption i.glm i.lm i.multinomfit i.multinomial i.polr i.polrfit IR is.formula i.survreg last legendr leverage linear.predictors logst mframe modarg modelTable nainf.exclude notna panelDefault plcoord plfitpairs plmatrix plmbox plmboxes plot.regr plot.xdistResscale plres2x plresx plTA.polr predict.mlm predict.polr predict.regr prevgumbel print.allcoef print.modelTable print.mregr print.regr qrevgumbel quantilew quinterpol regr residuals.polr robrange rrevgumbel showd simresiduals simresiduals.default simresiduals.glm smoothM smoothMM smoothRegr stamp step step.regr subset summary.mregr summary.regr sumna terms2order tit Tobit userOptions warn xdistResdiff xdistResscale xNA xNA

##  regr.R  Functions that are useful for regression, W. Stahel,
## ==========================================================================
regr <- function(formula, data=NULL, tit=NULL, family=NULL, dist=NULL,
                 calcdisp=NULL, suffmean=3,
                 nonlinear = FALSE, start=NULL,
                 robust = FALSE, method=NULL, ## init.reg="f.ltsreg",
                 subset=NULL, weights=NULL, na.action=nainf.exclude,
                 contrasts=getUserOption("regr.contrasts"),
                 model = FALSE, x = TRUE, termtable=TRUE, vif=TRUE,
                 factorNA = TRUE, testlevel = 0.05, hatlim=c(0.99,0.5),...)
{
  ## !!! dispersion: allow to be set.
  ## Purpose:    fit all kinds of regression models
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   formula, data, ...  as with lm
  ##   tit        title (becomes tit attribute of result)
  ##   calcdisp   should dispersion be calculated for
  ##              family=binomial and family=poisson
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  Jan 02
  ## -------------------------------------------------------------------------
  ## --- preparations
  lcall <- match.call()
  ## b. --- convert character formula to formula
  lformula <- as.formula(formula)
  lcall$formula <- lformula
  ## d. --- data
  dataname <- deparse(substitute(data))
  ldata <- as.data.frame(eval(data))
    lnobs <- nrow(ldata)
  if (lnobs==0) stop("!regr! no observations in data")
  ## --- get all variables to be used (as in lm)
##-   lcgetv <- call(quote("get_all_vars"), formula=lformula, data=lcall$data)
  ## nonlinear: drop constants
  lform <- lformula
  if (!(is.null(nonlinear)||as.character(nonlinear)=="FALSE")) {
    lcf <- setdiff(all.vars(formula),names(data))
    if(nonlinear=="lazy") {
      lcfundef <- setdiff(lcf, names(start))
      if (length(lcfundef)) {
        warning(":regr(nonlinear): check if  ",paste(lcfundef, collapse=", "),
                "  are coefficients. Starting values will be 0.")
        lst <- rep(0,length(lcfundef))
        names(lst) <- lcfundef
        start <- c(start,lst)
      }
    }
    lcall$start <- start <- start[names(start)%in%lcf]
    lfo <- setdiff(all.vars(formula),names(start))
    lform <- as.formula(paste("~",paste(lfo,collapse="+")))
    lcall$nonlinear <- nonlinear <- TRUE
  }
  ## get variables
  lcgetv <- call(quote("get_all_vars"), formula=lform, data=lcall$data)
  ## !!!  get_all_vars , when repaired...
  lcgetv$data <- eval(lcall$data, sys.parent()) ## !!! oct 15
##-   lcgetv$formula <- lform[1:2]
  lallvars <- try(eval(lcgetv, parent.frame())) ##
##-   if (length(lform)>2) {
##-     lcgetv$formula <- lform[-2]
##-     lallvar2 <- try(eval(lcgetv, parent.frame())) ##
##-   ## !!! simplify when  get_all_vars  is ok.
  ## variables not found -- which ones?
    if (class(lallvars)=="try-error") {
    lv <- all.vars(formula)
    lvm <- setdiff(lv, names(ldata))
    for (lvv in lvm)  if(exists(lvv)&&!is.function(get(lvv)))
                        lvm <- setdiff(lvm,lvv)
    stop("!regr! undefined variables in formula:  ",
         paste(lvm, collapse=", "))
  }
  ## repair nemes of lallvars
  if (anyNA(names(lallvars))) {
    lnm <- all.vars(lform)
    if (ldot <- match(".",lnm, nomatch=0))
      lnm <- unique(c(lnm[-ldot], names(ldata)))
    if (ncol(lallvars)==length(lnm)) names(lallvars) <- lnm  else
    stop("!regr! bug: variables not correctly identified")
  }
  ## f. --- compose call of fitting function
  lcl <- lcall
  if (!is.null(lcl$weights)) {
    lwgt <- eval(substitute(weights), data)
    lallvars <- cbind(lallvars, .weights.=lwgt)
    lcl$weights <- as.name(".weights.")
  }
  if (!is.null(lcl$offset)) {
    lallvars <- cbind(lallvars, .offset.=eval(substitute(offset), data))
    lcl$offset <- as.name(".offset.")
  }
  if (!is.null(lcl$subset)) {
    li <- eval(substitute(subset), data, environment())
    lallvars <- lallvars[li,, drop=FALSE]
    lcl$subset <- NULL
  }
  ## convert character to factor, drop unused levels, generate .NA. level
  lfacna <- is.character(factorNA) || (is.logical(factorNA)&factorNA)
  lfnalabel <- if(is.character(factorNA)) factorNA else ".NA."
  for (lvn in 1:ncol(lallvars)) {
    lv <- lallvars[[lvn]]
    if (is.character(lv)|is.factor(lv))
      lallvars[[lvn]] <- if (lfacna) factorNA(lv) else factor(lv, lfnalabel)
  }
  ## g. --- check for variables with a single value
  lcgetv[[1]] <- as.name("model.frame")
  lcgetv$data <- lallvars
  lcgetv$formula <-if (length(lform)>2) lform[-2] else lform
  lmodelframe <- eval(lcgetv, parent.frame())
  lv1 <- which( apply(lmodelframe, 2, function(x) all(x==x[1]) ) )
  if (length(lv1)) {  ## adjust formula
    lfac <- attr(terms(lmodelframe),"fac")
    lt1 <- names(which(apply(lfac[lv1,,drop=F],2,any)))
    warning("!regr! formula contains single valued variables: ",
            paste(row.names(lfac)[lv1], collapse=". "),
            "\n     I drop the following terms from the formula:\n    ",
            paste(lt1, collapse=", "))
    lfupd <- as.formula( paste( ".~ .- ",paste(lt1, collapse=" - ") ) )
    lcl$formula <- update(lcl$formula, lfupd)
  }
  ## -------------------------------------------
  ## h. --- response type
  if (length(lformula)==2) { # nonlinear called with formula of type ~...
    ly <- rep(0,NROW(lallvars))
    lytype <- "numeric"
  } else {
    lyf <- model.frame(lformula[1:2], lallvars)
    ## I tried to generate model.frame for x and y together. This failed
    ## because model.frame  needs adequate method (when y is matrix)
    ltrm <- attr(lyf, "terms")
    lytype <- substring(attr(ltrm, "dataClasses"),1,5)
    lysimple <- lytype!="nmatr" ## not a matrix
    ly <- lyf[[1]]
    if (lysimple&&length(unique(ly))==2 &&
        (is.factor(ly[[1]]) || all(ly%in%0:1)))
        lytype <- "binary"
    if (inherits(ly,"Surv"))  {
        lytype <- "survival"
    }
## strange variables
##-   l1v <- sapply(ldta, function(x) all(x==c(x[!is.na(x)],0)[1],na.rm=TRUE) )
##-                                 ## covers case of several or all NAs
##-   if (any(l1v)) {
##-     warning(paste(":regr: variable(s)", paste(lvars[l1v],collapse=", "),
##-                   "has (have) no distinct values")) #  -> dropped.
##-   }
  }
  ## k. --- family and fitting function
  lfam <- if (u.nuna(family)) NULL else as.character(substitute(family))[1]
  if (u.nuna(lfam)) lfam <- dist
  lcl$dist <- NULL
  if (lytype=="survival")
    lfam <- c( lfam, attr(ly,"distribution"))[1]
  if (u.nuna(lfam)) 
    lfam <- switch(substring(lytype,1,5),
                   numer="normal", nmatr="normal", binar="binomial",
                   binco="binomial", order="cumlogit",
                   facto="multinomial", survi="ph", "unknown")
  if (substring(lfam,1,7)=="multinom") lfam <- "multinomial"
  ##
  lfitfun <-
      switch( lfam,
             gaussian="lm", normal="lm", binomial="glm", poisson="glm",
             Gamma="glm",
             cumlogit="polr", multinomial="multinomial",
             weibull="survreg", lognormal="survreg", loggaussian="survreg",
             loglogis="survreg", loglogistic="survreg", extreme="survreg",
             ph="survreg", prop.hazard="survreg",
             "unknown")
  if (lfitfun=="unknown") stop("!regr! Fitting function not identified")
  ## additional checks
  if (lytype=="survival") {
    if (!inherits(ly,"Surv"))
      stop("!regr! bug: convert response to Surv object")
    ## !!! hier machen! lallvars[,1] ersetzen durch Surv davon
    lfitfun <- "survreg"
  }
  else  if (lfitfun=="glm")
    lcl$control <- list(calcdisp=calcdisp, suffmean=suffmean,lcl$control)
  ## 
  lfitname <- paste("i",lfitfun,sep=".")
  if (!exists(lfitname)||!is.function(get(lfitname)))
    stop (paste("!regr! Fitting function",lfitname, "not found"))
  ## m. --- prepare call
  lcl[[1]] <- ## hack --> eval(.) works also when call is source()d ...
      switch(lfitname,
	     "i.lm" = quote(regr0::i.lm),
	     "i.glm" = quote(regr0::i.glm),
	     "i.multinomial" = quote(regr0::i.multinomial),
	     "i.polr" = quote(regr0::i.polr), 
	     "i.smooth" = quote(regr0::i.smooth), ## ??
	     "i.survreg" = quote(regr0::i.survreg),
	     ## default:
	     as.name(lfitname))
##  lcl[[1]] <- as.name(lfitname) ## sonst geht das debuggen nicht.
  lcl$fname <- lfam
  lcl$na.action <- substitute(na.action)
  lcl$data <- ## as.name("lallvars") ## environment(formula)
    ldata <- na.action(lallvars)
  environment(lcl$formula) <- environment() ## !!!
##-  lcl$data <- eval(lcl$data, sys.parent())
  ## problem with environment if different for  data  and  formula
  if (lfitname=="i.survreg") {
    lcl$yy <- ly
    lcl$model <- TRUE  ## model needed, see below
  }
  old.opt <- NULL
  if(is.atomic(contrasts)&&length(contrasts)) {
    if(!is.character(contrasts))
      warning("!regr! invalid contrasts argument") else {
        old.opt <- options(contrasts=c(contrasts,"contr.poly")[1:2])
        lcl$contrasts <- NULL
      }
    if (contrasts[1]=="contr.wsum") lallvars <- contr.wsum(lallvars, ly)
##    lcl$contrasts[1] <- c("contr.sum", lcl$contrasts[2])  ## modify: appply  contr.wfac  to model.frame 
  }
##- ## --------------------------------------------
  ##-   lreg <- eval(lcl, envir=environment(formula))
  lreg <- eval(lcl)
  ## --------------------------------------------
  if (length(old.opt)) options(old.opt)
  lreg$na.action <- attr(ldata, "na.action")
  if (is.null(lreg$distrname)) lreg$distrname <- lfam
##  <<<<<<< .mine
#  lreg$Y <- data.frame(ly) # ly is a model.frame
#  if (ncol(lyy)>1) colnames(lyy) <- colnames(ly[[1]])  
#  lreg$Y <- lyy
## =======  !?!
  lyy <- as.matrix(ly) # ly is a model.frame
  if (ncol(lyy)>1) colnames(lyy) <- colnames(ly[[1]])
  lreg$Y <- lyy
  lreg$response <- ly
  if (nonlinear) lreg$r.squared <- 1-lreg$sigma^2/var(ly)
  ## >>>>>>> .r32
  lreg$allvars <- lallvars ## needed more than $model
  ## since $model contains transformed variables
  ## recover some arguments to effective function call
  lfc <- lreg$call
  largs <- intersect(c("data", "weights", "offset", "subset"), names(lfc))
  ## these arguments should be restored because otherwise,
  ## add1 does not work if they have changed.
  lfc[largs] <- lcall[largs]
  lreg$funcall <- lfc
  lcall$formula <- formula(lreg) # hope this never damages anything
  lreg$call <- lcall
  tit(lreg) <- if (length(tit)==0) attr(data,"tit") else tit
  doc(lreg) <- attr(data,"doc")
  if (model&&length(lreg$model)==0) {
      if (nonlinear) warning(":regr: no $model available for nonlinear regr.")
      else lreg$model <- lm(lformula, data, method="model.frame")
  }
  lterms <- if (nonlinear) NULL else terms(lreg)
  if ((!nonlinear) && is.null(attr(lterms, "predvars")))  ## needed for survreg
    attr(lreg$terms,"predvars") <- attr(attr(lreg$model,"terms"),"predvars")
  lresnm <- colnames(lreg$residuals)
  ## r. --- leverages, standardized res
  lhat <- lreg$leverage
  if (length(lhat)==0 && !nonlinear) {
    lmm <- lreg[["x"]]
    if (length(lmm)==0) lmm <- model.matrix(lterms,lreg$model)
    lreg$leverage <- lhat <- hat(lmm)
    if (length(lhat)==0) warning(":regr: no leverages available")
  }
  if (length(lhat)==NROW(lreg$stres))  lreg$stres <-
      lreg$stres/sqrt(1-pmin(hatlim[1],lhat)) else
      if (length(lreg$stres))
        warning(":regr: bug: leverages and st.res. incompatible")
  ## misc
  if (class(lreg)[1]=="survreg")
    lreg$n.obs <- length(lreg$linear.predictor)
  if (is.null(x) || !x) lreg$x <- NULL
  class(lreg) <- if (class(lreg)[1]=="orig")  ##  nls shall not be regr
    class(lreg)[-1] else c("regr",class(lreg))
  ## result of regr
  lreg
}
## -----------------------------------------------------------------------
i.lm <- function(formula, data, family, fname="gaussian", nonlinear=FALSE,
                 robust=FALSE, method=NULL, control=NULL, 
                 vif=TRUE, termtable=TRUE, testlevel=0.05, ...)
{
  ## Purpose:  internal: fit lm
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 11:18
  ## b. --- method
  lcall <- match.call()
  lmeth <- c(lcall$method,"")[1]
  lfn <- if (nonlinear) {
    lcall$contrasts <- NULL ## lcall[-match("contrasts",names(lcall),nomatch=0)]
    "nls" } else {
    if (lmeth=="lmrob") robust <- TRUE
    if (robust) {
      if (is.null(method)) method <- lcall$method
      if (is.null(method)) method <- c("lmrob","KS")
      method[1]
    } else "lm"}
  if (robust) {
    if (lfn=="lmrob") {
 ##     require(robustbase)  ## !?!
      if (length(method)>1) {
        if (substr(method[2],1,2)=="KS")
          method <- NULL
        lcall$setting <- "KS2011"
      }
      lcall$x <- TRUE
  }
    if (lfn=="rlm") {
##      require(MASS)    ##  !?!
      lcall$method <- c(method,"MM")
      lcall$x.ret <- TRUE
    }
  } else  lcall$x <- TRUE
  if (lmeth=="rq"|lmeth=="quantreg") { # quantile regression
##    require(quantreg) ## !?!
    lfn <- "rq"
    lcall$method <- if(length(lcall$method)>1) lcall$method[-1] else NULL
    lcall$x <- NULL
  }
  ## d. --- call
  lcall$method <- if (length(method)>1) method[-1] else NULL
  ##  method[-1]  produces character(0) which is not NULL!
  mkFn <- function(fn) { ## hack --> eval(.) works also when call is source()d ...## ??? wieso function?
      switch(fn,
	     "lmrob" = quote(robustbase::lmrob),
	     "rlm" = quote(MASS::rlm),
             "rq" = quote(quantreg::rq),
	     ## default:
	     as.name(fn))
  }
  lcall[[1]] <- mkFn(lfn)
  lcall$fname <- lcall$family <- lcall$vif <- lcall$nonlinear <-
    lcall$robust <- NULL
  ## --------------------------
  lreg <- eval(lcall, envir=environment(formula))
  ## --------------------------
  ## f. --- collect results
  lreg$call$formula <- formula
  lreg$fitfun <- lfn
  lreg$distrname <- "gaussian"
  lttype <- switch(lfn,
                   rq="Chisq",
                   rlm="Chisq",
                   "F"
                   )
##-   ## leverage
##-   if (!nonlinear) {
##-     lhat <- pmax(0,hat(lreg$x))
##-     if (length(lhat)==0) warning(":regr/i.lm: no leverages")  ## else {
##- ##-       if (length(lhat)!=NROW(lreg$stres))
##- ##-         if (length(lreg[["w"]])==NROW(lreg$stres))
##- ##-           lhat <- u.merge(lreg$leverage, 0, lreg[["w"]]>0)
##- ##-     }
##-     lreg$leverage <- lhat
##-   }
  ## multivariate
  if (class(lreg)[1]=="mlm")
    return(i.mlmsum(lreg, termtable))
  ## 
  lreg1 <- summary(lreg)
  lsig <- lreg1$sigma
  if (is.null(lsig)) lsig <- lreg$scale ## lmrob
  if (is.null(lsig)) lsig <- sd(resid(lreg)) # !!! used for rq
  lreg$sigma <- lsig
  ## standardized residuals
  if (is.finite(lsig)&&lsig>0) {
    lreg$stres <- lreg$residuals/lsig
    if (length(lreg$weights)) lreg$stres <- lreg$stres*sqrt(lreg$weights)
  } ## standardization by lhat is done in regr
  if (class(lreg)=="lmrob") lreg1$cov.unscaled <- lreg$cov/lsig^2 ## !!!
  ## from summary
  lcomp <- c("r.squared","fstatistic","colregelation","aliased",
             "df","cov.unscaled")
  lreg[lcomp] <- lreg1[lcomp]
  ## degrees of freedom
  if (is.null(lreg$df)) # needed for rq
    lreg$df <- c(length(coef(lreg))-attr(terms(lreg),"intercept"),
                 length(resid(lreg))-length(coef(lreg)))
  lreg$df.residual <- ldfr <- df.residual(lreg)
  if (nonlinear) {
    lcftab <- lreg1$coefficients
    lreg$coefficients <- lcftab[,1]
    ltq <- qt(1-testlevel/2, ldfr)
    lci <- lcftab[,1]*(1+outer(ltq*lcftab[,2], c(ciLow=-1,ciHigh=1)))
    lreg$termtable <- data.frame(coef=lcftab[,1],se=lcftab[,2],lci)
#   lreg$r.squared <- 1-(lsig/lsdy)^2
  }
  lreg$adj.r.squared <- 1-(1-lreg$r.squared)*(length(lreg$residuals)-1)/ldfr
  ## cov of estimates
  lcov <- lreg$cov.unscaled*lsig^2
  lreg$covariance <- lcov
  lse <- sqrt(diag(lcov))
  lreg$correlation <- lcov/outer(lse, lse)
  ## --- table of terms
  if (!nonlinear) {
    if(termtable) {
      ly <- lreg$model[[1]]
      lsdy <- sqrt(var(ly))
      ltt <- i.termtable(lreg, lreg1$coef, data, lcov, lttype, lsdy=lsdy,
                         vif=vif, leverage=TRUE) 
      lcmpn <- c("termtable","allcoef","leverage")
      lreg[lcmpn[lcmpn%in%names(ltt)]] <- ltt
    }  else  class(lreg) <- c("orig",class(lreg))
  }
  ## result of i.lm
  lreg
}
## -----------------------------------------------------------------------
i.mlmsum <- function(object, termtable=TRUE)
{
  ## Purpose:  internal: fit multivariate lm;  called from i.lm() only
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 11:18
##-   lreg <- lm(formula, data=data, weights=.Weights, model=model, ...)
##-   lreg$call$formula <- eval(lreg$call$formula) # patch
  lreg1 <- summary(object)
  lform <- formula(object)
  lts <- ltp <- NULL
  for (ly in 1:length(lreg1)) {
    lrg <- lreg1[[ly]]
    lts <- cbind(lts,c(lrg[["sigma"]],lrg[["r.squared"]],
                       lrg[["fstatistic"]]))
    ltp <- cbind(ltp,lrg[["coefficients"]][,4])
  }
  lmodel <- nrow(lts)>=5  # non-trivial model
  if (lmodel) {
    lts[4,] <- pf(lts[3,],lts[4,],lts[5,], lower.tail=FALSE)
    lts <- lts[1:4,]
  }
  dimnames(object$coefficients)[[2]] <- as.character(lform[[2]])[-1]
  dimnames(ltp) <- dimnames(object$coefficients)
  dimnames(lts) <- list(rep(c("sigma","r.squared","fstatistic","p-value"),
                            length=nrow(lts)), dimnames(ltp)[[2]])
  object$pvalues <- ltp
  object$stats <- lts
  object$sigma <- lsig <- lts["sigma",]
  lres <- residuals(object)
  if (all(lsig>0)) {
    object$stres <- sweep(lres,2,lsig,"/")
    if (length(object$weights))
      object$stres <- object$stres*sqrt(object$weights)
  }
  object$resmd <- mahalanobis(lres,0,var(lres))
  ldfr <- object$df.residual
  object$r.squared <- lr2 <- lts["r.squared",]
  object$adj.r.squared <- 1-(1-lr2)*(nrow(object$resid)-1)/ldfr
  lcomp <- c("aliased","df","cov.unscaled")
  object[lcomp] <- lreg1[[1]][lcomp]
  object$drop1 <- if (lmodel) drop1.mlm(object)
##  class(lreg) <- c("mregr","mlm","lm")
  object
} # i.mlmsum
## -----------------------------------------------------------------------
i.glm <- function(formula, data, family, fname,
                  control=NULL, vif=TRUE, termtable=TRUE, ...)
{
  ## Purpose:  internal: fit glm
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 11:18
  lfamily <- get(fname)
##-   environment(formula) <- environment()
  lcall <- match.call()
  lcall[[1]] <- as.name("glm")
  lcall$family <- lcall$fname
  lcall$calcdisp <- lcall$fname <- lcall$control <- lcall$vif <- NULL
  lcall$x <- TRUE
  ## ---------------
  lreg <- eval(lcall, envir=environment(formula))
  ## ----------------
  lreg$leverage <- pmax(0,hat(lreg$x))
  lreg1 <- summary(lreg)
  lcoeftab <- lreg1$coef
  ly <- as.numeric(lreg$model[,1])
  ldisp <- lreg1$dispersion
  ## ---
  lfcount <- fname=="binomial"|fname=="poisson"
  lcalcdisp <- control$calcdisp
  lsuffmean <- TRUE
  if (lfcount) {
    lsuffmean <- mean(ly)>control$suffmean   # ,na.rm=TRUE
    lcd <- lcalcdisp
    if (length(lcd)==0) lcd <- lsuffmean
    if (lcd) {
      ldisp <- lreg1$deviance/lreg1$df.residual
      if (ldisp>1||length(lcalcdisp)>0) {
        lreg$distrname <- paste("quasi",fname,sep="")
        lcoeftab[,2] <- lcoeftab[,2]*sqrt(ldisp)
        lcoeftab[,3] <- lcoeftab[,3]/sqrt(ldisp)
##-         lcoeftab[,4] <- 2*pnorm(lcoeftab[,3],lower.tail=FALSE)
      }
      else ldisp <- 1
    }
  }  # else calcdisp <- FALSE
  attr(ldisp,"fixed") <- ldisp==1
  lreg$dispersion <- ldisp
  lreg$sigma <- sqrt(ldisp)
  ## ---
  if (ldisp>0) {
      lstr <- residuals(lreg, type="pearson")/sqrt(ldisp)
      lnaa <- lreg$na.action
      if (class(lnaa)=="exclude") lstr <- lstr[-lnaa]
      lreg$stres <- lstr
  }
  ## bug? leverage not taken into account
  lcomp <- c("deviance","aic","df.residual","null.deviance", # "family",
    "df.null","iter","deviance.resid","aliased","df","cov.unscaled")
  lreg[lcomp] <- lreg1[lcomp]
  ## --- deviances
  ltesttype <- ifelse(ldisp==1,"Chisq","F")
  ldev <- unlist(lreg1[c("deviance", "null.deviance")])
  ldf <- lreg1$df[1:2]-c(attr(terms(lreg),"intercept"),0)
  ltbd <- cbind(deviance=c(diff(ldev),ldev), df=c(ldf,sum(ldf)),
                p.value=NA)
  dimnames(ltbd)[[1]] <- c("Model","Residual","Null")
  ltbd[1:2,3] <- pchisq(ltbd[1:2,1], ltbd[1:2,2], lower.tail=FALSE)
  if (!lsuffmean) ltbd[2,3] <- NA
  lreg$devtable <- ltbd
  ## ---
  ## cov of estimates
  ldisp <- lreg$dispersion
  if (is.null(ldisp)) ldisp <- 1
  lreg$covariance <- lcov <- lreg$cov.unscaled*ldisp
  lse <- sqrt(diag(lcov))
  lreg$correlation <- lcov/outer(lse, lse)
  ## ---
  lreg$fitfun <- "glm"
  if (termtable) {
    ltt <- i.termtable(lreg, lcoeftab, data, lcov, ltesttype, lsdy=1, vif=vif)
    lcmpn <- c("termtable","allcoef","leverage")
    lreg[lcmpn[lcmpn%in%names(ltt)]] <- ltt
  }
  ## result of i.glm
  lreg
}
## -----------------------------------------------------------------------
i.multinomial <- function(formula, data, family, fname,
                          model=TRUE, vif=TRUE, termtable=TRUE, ...)
{
  ## Purpose:  internal: fit multinom
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 11:18
##  ltr <- control$trace
##  if (length(ltr)==0) ltr <- trace
##  require(nnet)   ## !?!
  lcall <- match.call()
  lcall[[1]] <- quote(regr0::i.multinomfit)
  lcall$fname <- lcall$family <- lcall$control <- lcall$vif <- NULL
  lcall$trace <- FALSE
  lreg <- eval(lcall, envir=environment(formula))
  ## ---------------
  if (length(lreg$na.action)) {
    lnaact <- attr(lreg$na.action,"class")
    attr(lreg$na.action,"class") <- "omit"
  } else  lnaact <- NULL ##  summary does not work with  exclude
  lreg$call$formula <- formula
  lreg1 <- summary(lreg)
  lreg$dispersion <- lreg$sigma <- 1
  lres <- lreg1$residuals
  lreg$residuals <- lres
  lcf <- lreg1$coefficients
  lreg$coefficients <- lcf
  lreg$aic <- lreg1$AIC
  ldfm <- lreg1$edf-nrow(lcf)
  lreg$df <- c(ldfm,prod(dim(lres)-1)-ldfm,ldfm)
##-   environment(lreg$call$formula) <- environment()
  lreg$fitfun <- "multinom"
  ldr1 <- if (u.debug()) drop1(lreg, test="Chisq", trace=FALSE) else 
              try(drop1(lreg, test="Chisq", trace=FALSE), silent=TRUE)
  if (class(ldr1)[1]=="try-error") {
    warning(paste(":regr/i.multinom: drop1 did not work.",
                  "I return the multinom object"))
    class(lreg) <- c("orig",class(lreg))
    return(lreg)
  } else {
  ldr1 <- ldr1[-1,]}
  names(ldr1) <- c("df", "AIC", "Chisq","p.value") #if(calcdisp) "F" else
  lreg$termtable <- lreg$drop1 <- ldr1
  if (length(lnaact)) attr(lreg$na.action,"class") <- lnaact
## result of i.multinomial
  lreg
}
## -----------------------------------------------------------------------
i.polr <- function(formula, data, family, fname, weights = NULL, 
                   model=TRUE, vif=TRUE, termtable=TRUE, ...)
{
  ## Purpose:  internal: fit ordered y
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 11:18
##  require(MASS)   ## !?!
  lcall <- match.call()
  lcall[[1]] <- quote(MASS::polr) ## quote(i.polrfit)
  lcall$contrasts <- lcall$fname <- lcall$control <- lcall$family <-
    lcall$vif <- NULL
  lcall$Hess <- TRUE
## ---
  lreg <- eval(lcall, envir=environment(formula))
##  lreg$call$formula <- formula
  lreg$w <- data$.weights.
  lreg$leverage <- hat(lreg[["x"]])
  lreg1 <- if (u.debug()) summary(lreg) else
           try(summary(lreg))
  if (class(lreg1)[1]=="try-error") {
    warning(paste(":regr/i.polr: summary did not work.",
                  "I return the polr object"))
##    lreg$call$data <- call$data
    class(lreg) <- c("orig","polr")
    return(lreg)
  }   ## ---
  ## model.matrix
##-   ldata <- eval(data, envir=environment(formula))
##  lreg$x <- model.matrix(formula, ldata)
  lcf <- lreg1$coefficients
  lreg$intercepts <- lcf[(lreg1$pc+1):nrow(lcf),1:2]
  lreg$stres <- NULL
  ## cov of estimates!
  lreg$covariance <- lcov <- vcov(lreg)
  lse <- sqrt(diag(lcov))
  lreg$correlation <- lcov/outer(lse, lse)
  ## --- deviances
  lreg$fitfun <- "polr"
  if (termtable) {
    ltt <- i.termtable(lreg, lreg1$coef, data, lcov, ltesttype="Chisq",
                       lsdy=1, vif=vif, leverage=TRUE)
    lcmpn <- c("termtable","allcoef","leverage")
    lreg[lcmpn[lcmpn%in%names(ltt)]] <- ltt
  }
  lreg$dispersion <- 1
  lreg$residuals <- residuals(lreg)
  lreg$linear.predictors <- fitted.polr(lreg, type="link")
  ## result of i.polr
  lreg
}
## -----------------------------------------------------------------------
i.survreg <-
  function(formula, data, family, yy, fname="ph", method, control,
           vif=TRUE, termtable=TRUE, ...)
{
  ## Purpose:  internal: fit ordered y
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 11:18
##  require(survival)  ## !?!
  lcall <- match.call()
#  lcall <- lcall[-match("contrasts",names(lcall),nomatch=0)]
  ## b. --- method
  if (fname=="ph") {
    lfitfun <- "coxph"
    lcall[[1]] <- quote(survival::coxph)
  } else {
    lfitfun <- "survreg"
    lcall[[1]] <- quote(survival::survreg)
    lcall$dist <- fname
    lcall$method <- lcall$control <- NULL
  }
  lyy <- yy ## lreg$y
  lcall$yy <- lcall$fname <- lcall$family <- lcall$vif <- NULL
  lcall$x <- TRUE
  ## ---
  lreg <- eval(lcall, envir=environment(formula))
  ## ---
  class(lreg$y) <- "Surv"
  attr(lreg$y, "type") <- attr(lyy, "type")
##  lreg$call$formula <- formula
  lreg1 <- if (u.debug()) summary(lreg) else
           try(summary(lreg), silent=TRUE)
  if (class(lreg1)[1]=="try-error") {
    warning(paste(":regr/i.survreg: summary did not work. ",
                  "I return the survreg object"))
##    lreg$call$data <- call$data
    class(lreg) <- c("orig",class(lreg))
    return(lreg)
  }   ## ---
  lreg$resid.orig <- lreg$residuals
  lreg$stres <- NULL
  lcf <- lreg1$coefficients
  ## --- deviances
  ## lreg$scale
  if (lfitfun=="survreg") {
    attr(lreg$scale,"fixed") <- length(lcall$scale)>0
    ldf <- sum(lreg$df) - lreg$idf
    ldfr <- length(lreg$linear.predictors)-sum(lreg$df)
  }
  if (lfitfun=="coxph") {
    lreg1$table <- lreg1$coefficients
    ldf <- length(lreg$coefficients)
    ldfr <- length(lreg$residual)-ldf-1
  }
  lreg$df.residual <- ldfr
  lreg$aic <- extractAIC(lreg)[2]
  lreg$deviance <- -2*lreg$loglik
  lchi <- 2*diff(lreg1$loglik)
  ltbd <- cbind(deviance=c(lchi,-2*lreg1$loglik[2]),
                df=c(ldf, ldf+ldfr),
                p.value=c(pchisq(lchi,ldf,lower.tail=FALSE),NA))
  dimnames(ltbd)[[1]] <- c("Model","Null")
  lreg$devtable <- ltbd
  lreg$covariance <- lcov <- lreg$var
  lsd <- sqrt(diag(lcov))
  lreg$correlation <- lcov/outer(lsd,lsd)
  lreg$fitfun <- lfitfun
  lres <- residuals.regr(lreg)
  ly <- lreg$y
##-   lreg$n.censored <-
##-     if (attr(ly,"type")%in%c("right","left"))
##-       table(ly[,2])[2] else  sum(ly[,2]!=1)    #interval
  ltype <- attr(ly,"type")
  ##-  lreg$n.censored <- sum(lres[,"prob"]>0, na.rm=TRUE)
  ltb <- table(ly[,2])
  lfit <- lres[,"fit"]
  llimit <- attr(ly,"limit")
  lreg$n.censored <- NA
  if (ltype=="left") {
    lreg$n.censored <- structure(ltb[2], names="left")
    if (length(llimit))
      lreg$n.fitout <- structure(sum(lfit<llimit, na.rm=TRUE), names="left")
  } else {
    if (ltype=="right") {
      lreg$n.censored <- structure(ltb[1], names="right")
    if (length(llimit))
    lreg$n.fitout <- structure(sum(lfit>llimit, na.rm=TRUE), names="right")
    } 
  }
  if (termtable) {
    ltt <- i.termtable(lreg, lreg1$table, data, lcov, ltesttype="Chisq",
                       lsdy=1, vif=vif)
    ## log(scale): signif<-NA. no! log(scale)==0 means
    ##    exp.distr for weibull/gumbel
    lcmpn <- c("termtable","allcoef","leverage")
    lreg[lcmpn[lcmpn%in%names(ltt)]] <- ltt
  }
  ## lreg$df <- c(model=ldf, residual=ldfr, original=lreg$df) ## !!! lreg has df = ldf+object$idf !
    ## do not modify before calling i.termtable
  lreg$distrname <- if (lfitfun=="coxph") "prop.hazard" else lreg$dist
  lreg$residuals <- lres
  ## result of i.survreg
  lreg
}

## -----------------------------------------------------------------------
Tobit <- function(data, limit=0, limhigh=NULL, transform=NULL, log=FALSE, ...)
{
  ## Purpose:   create a Surv object for tobit regression
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  1 Jan 2010, 21:49
##  require(survival)  ## !?!
  ltrs <- as.character(substitute(transform))
  data <- pmax(data,limit)
  lright <- !is.null(limhigh)
  if (lright) data <- pmin(data, limhigh)
  if (log[1]) { ## model.frame  evaluates  log  in  data ! Whence [1]
      transform <- logst
      ltrs <- "logst"
  }
  if (!is.null(transform)) {
    if (is.character(transform)) transform <- get(transform)
    if (!is.function(transform))
      stop("!Tobit! argument 'transform' does not yield a function")
    ldt <- transform(c(limit,limhigh,data), ...)
    data <- ldt[-1]
    limit <- ldt[1]
    if (lright) {
      limhigh <- ldt[2]
      data <- data[-1]
    }
  }
  if (sum(data<=limit,na.rm=TRUE)==0)
    warning(":Tobit: no observation <= `limit`")
  if (lright&&sum(data>=limhigh,na.rm=TRUE)<=1)
    warning(":Tobit: no observation >= `limhigh`")
  if (lright) { 
    rr <- survival::Surv(time = data, time2=data,
                         event = (data<=limit) + (data<limhigh),
                         type="interval")
    rr[,2] <- rr[,1]
    } else
      rr <- survival::Surv(data, event = data>limit, type="left")
  structure(rr, distribution="gaussian", transform=ltrs,
            limit=c(limit,limhigh), class=c(class(rr), "matrix"))
}

## -----------------------------------------------------------------------
i.termtable <- function(lreg, lcoeftab, ldata, lcov, ltesttype="F",
                        lsdy, vif=TRUE, leverage=vif)
{
  ## Purpose:  generate term table for various models
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  4 Aug 2004, 15:37
  if(length(attr(terms(lreg),"term.labels"))==0)
    return(list(termtable = data.frame(
      coef=c(lreg$coef,NA)[1], se=NA, ciLow=NA, ciHigh=NA, 
      df=1, testst=NA, signif=NA, p.value=NA, p.symb="", stcoef=NA, R2.x=NA,
      stringsAsFactors=FALSE)
                ))
## degrees of freedom
  ldfr <- df.residual(lreg)
  if (ldfr<1) {
    warning(":regr/i.termtable: no degrees of freedom left.")
    return(list(termtable = data.frame(
      coef=c(lreg$coef,NA)[1], se=NA, ciLow=NA, ciHigh=NA, 
      df=1, testst=NA, signif=NA, p.value=NA, p.symb="", stcoef=NA, R2.x=NA,
      stringsAsFactors=FALSE)
                ))
  }    
  pvCutpoints <- c(0, 0.001, 0.01, 0.05, 0.1, 1)
  pvSymbols <- c("***", "**", "*", ".", " ", "")
  pvLegend <- paste(rbind(pvCutpoints,pvSymbols), collapse="  ")
  ## drop1
  ldr1 <-
      if (class(lreg)[1]%in%c("lm","lmrob")) {
          if (u.debug()) 
              drop1Wald(lreg, test=ltesttype, scope=terms(lreg)) else
              try(drop1Wald(lreg, test=ltesttype, scope=terms(lreg)),
              silent=TRUE) } else {
          if (u.debug()) 
              drop1(lreg, test=ltesttype, scope=terms(lreg)) else
              try(drop1(lreg, test=ltesttype, scope=terms(lreg)),
                  silent=TRUE)
                           }
  if (class(ldr1)[1]=="try-error") {
    warning(paste(":regr: drop1 did not work. I return the table produced by ",
                  lreg$fitfun))
##-     lsum <- summary(lreg)
##-     lcft <- lsum$coef
##-     if (length(lcft)==0) lcft <- lsum$parameters ## nls
##-     return(list(test=lcft)) # !!! noch reparieren
    return(list(termtable=lcoeftab))
  }
  ldr1 <- ldr1[-1,]
  ldr1$RSS <- NULL # same ncol for lm and glm
  if (inherits(lreg,"rlm"))  ldr1[,4] <- ldr1[,2]/ldr1[,1] ## !!!
  if (inherits(lreg,"mlm")||inherits(lreg,"manova"))
    return(list(termtable=ldr1))  ## !!! needs much more
  ltstq <- if (ltesttype=="F") qf(0.95,c(1,ldr1[,1]),ldfr) else {
    if (ltesttype=="Chisq") qchisq(0.95,c(1,ldr1[,1])) else NA }
  ltstq1 <- sqrt(ltstq[1]) ## 1 degree of freedom
  ltstq <- ltstq[-1]
## coefficients
  lcoef <- lreg$coefficients
## model.matrix
  lmmt <- lreg[["x"]]
  if (length(lmmt)==0)
      lmmt <- model.matrix(lreg)
  lasg <- attr(lmmt,"assign")[!is.na(lcoef)]
  if (class(lreg)[1]%in%c("polr","coxph")) lasg <- lasg[-1]
## terms with 1 coef
  lcont <- lasg[licasg <- which(!lasg%in%lasg[duplicated(lasg)])]
  ## vif --> R2.x
  lr2 <- NA
  if (vif) {
    lvift <-     ## lterms: n of levels for each term
        if (u.debug()) vif.regr(lreg, lcov, lmmt) else
        try(vif.regr(lreg, lcov, lmmt), silent=TRUE)
    if (class(lvift)[1]=="try-error" || length(lvift)==0) {
      warning(":regr/i.termtable: error in the calculation of R2.xs")
      lvif <- NA
    } else lvif <- lvift[,3]^2
    lr2 <- 1-1/lvif
  }
## prepare table
  lpvcol <- pmatch("Pr(",names(ldr1), nomatch=ncol(ldr1))
  lpv <- ldr1[,lpvcol]
  ltb <- data.frame(coef=NA, se=NA, ciLow=NA, ciHigh=NA, 
                    df=ldr1[,1], testst=ldr1[,lpvcol-1], signif=NA,
                    p.value=lpv, p.symb="", stcoef=NA, R2.x=lr2,
                    stringsAsFactors=FALSE)
  row.names(ltb) <- row.names(ldr1)
  ## intercept
  ljint <- "(Intercept)"==names(lcoef)[1]
  if (ljint) {
    ltstint <- # if(class(lreg)[1]%in%c("lm","nls","rlm"))
      lcoeftab[1,3]^2 # else lcoeftab[1,3]
    ltb <- rbind(
      "(Intercept)"=
        data.frame(coef=NA, se=NA, ciLow=NA, ciHigh=NA, 
                   df=1, testst=ltstint, signif=NA,
                   p.value=NA, p.symb="", stcoef=NA, R2.x=NA,
                   stringsAsFactors=FALSE),
      ltb)
    ltstq <- c(qf(0.95,1,ldfr), ltstq)
  }
  lcont1 <- lcont+ljint  # row number in dr1
## p.symb and signif
  ltb$signif <- sqrt(pmax(0,ltb$testst)/ltstq)
## 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 lreg1
    ljc <- match(lcont,lasg) # index of coefs for cont variables
    lcf <- lcoef[ljc]
    ## fill in
    ltb$coef[lcont1] <- lcf
    ltb$se[lcont1] <- lse <- lcoeftab[licasg,2]
    lci <- lcf+outer(ltstq1*lse, c(-1,1))
    ## confint(lreg,row.names(ltb)[lcont1]) does not always work...
    ltb[lcont1,c("ciLow","ciHigh")] <- lci
    ltb[lcont1,"signif"] <- sign(lcf)*ltb[lcont1,"signif"]
    ## standardized coefficients
    lstcf <- lcf[lcont>0] *  # exclude intercept term
      sqrt(apply(lmmt[,names(lcf[lcont>0]),drop=FALSE],2,var)) / lsdy
    ltb$stcoef[lcont1[lcont>0]] <- lstcf
  }
  if (row.names(lcoeftab)[nrow(lcoeftab)]=="Log(scale)") { # survreg
    ltsc <- lcoeftab[nrow(lcoeftab),]
    lcont1 <- c(lcont1, nrow(lcoeftab))
    if (!u.true(lreg$dist=="weibull")) ltsc[2:4] <- NA
    ltb <- rbind(ltb,"log(scale)"=
                   c(ltsc[1],ltsc[2],ltsc[1]+c(-1,1)*qnorm(0.975)*ltsc[2],
                     1, ltsc[3], ltsc[3]/qnorm(0.975), ltsc[4], NA, NA, NA))
  }
## p-symbol
  lipv <- as.numeric(cut(ltb$p.value, pvCutpoints))
  ltb[,"p.symb"] <- pvSymbols[lipv]
  attr(ltb, "legend") <- pvLegend
## --- allcoef (dummy coef)
  lallcf <- allcoef(lreg) 
  if (inherits(lreg,"polr")) lreg$coefficients <- c("(Intercept)" = NA, lcoef)
  rr <- list(termtable=ltb, allcoef=lallcf)
  if (leverage) rr <- c(rr, leverage=list(hat(lmmt)))
  rr
}
## --------------------------------------------------------------------------
ciSignif <- function(estimate, se=NULL, df=Inf, testlevel=0.05) {
  if (is.null(se))
    if (NCOL(estimate)>1) {
      se <- estimate[,2]
      estimate <- estimate[,1]
    } else
      stop("!ciSignif! no standard errors found")
  ltq <- qt(1-testlevel/2, df)
  lci <- estimate+outer(ltq*se, c(ciLow=-1,ciHigh=1))
  ltst <- estimate/se
  lsgf <- ltst/ltq
  lpv <- 2*pt(-abs(ltst), df)
  lipv <- as.numeric(cut(lpv, c(0, 0.001, 0.01, 0.05, 0.1, 1)))
  lsst <- c("***", "**", "*", ".", " ")[lipv]
##-     symnum(lpv, corr = FALSE, na = FALSE, 
##-                  cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), 
  ##-                  symbols = c("***", "**", "*", ".", " "))
  data.frame(estimate=estimate, se=se, lci, testst=ltst,
             signif=lsgf, p.value=lpv, p.symb=lsst)
}
  
## ==========================================================================
contr.wsum <- 
  function (n, y=NULL, w = NULL, contrasts = TRUE, sparse = FALSE) 
{  ## provide weighted sum contrasts
  if (is.data.frame(n)) {
    if (!is.null(y)) n <- n[apply(is.finite(cbind(y)), 1, all),, drop=FALSE]
    df <- n
    for (lj in 1:ncol(df)) 
      if(class(df[,lj])[1]=="factor") ## avoid ordered factors (for the time being)
        attr(df[,lj],"contrasts") <- contr.wsum(df[,lj])
    return(df)
  }
  ## argument is a number or ...
  if (is.character(n)) n <- factor(n)
  if (is.factor(n)) { ## ... a factor
    w <- c(table(n))
    n <- levels(n)
  }
  if ((is.numeric(n)&&n==1) || (is.character(n)&&length(n)==1))
             return(matrix(1,1,1)) else
    contr <- contr.sum(n, contrasts = contrasts, sparse = sparse)
  if (is.null(w)) {
    warning(":contr.wsum: no weights provided.",
            " Returning unweighted sum contrasts for\n  ", paste(n,collapse=", "))
    return(contr)
  }
  if (length(w)!=nrow(contr) || anyNA(w) || any(w<0) || all(w==0)) {
    warning(":contr.wsum: weights 'w' not suitable.",
            " Returning unweighted contrasts")
    return(contr)
  }
##-   li0 <- w==0
  ##-   if (any(li0)) contr <- contr[c(!li0,TRUE),!li0]
  n <- nrow(contr)
  if (w[n]==0) 
    warning(":contr.wsum: last weight is zero.",
            " Returning unweighted contrasts") else {
  contr[n,] <- - w[-n]/w[n] }
  structure(contr, w=w)
}
##  ===================================================================
print.regr <- function (x, call=TRUE, correlation = FALSE,
    dummycoef = getUserOption("show.dummycoef"),
    termcolumns = getUserOption("termcolumns"),
    allcoefcolumns = getUserOption("allcoefcolumns"),
    digits = max(3, getUserOption("digits")-2), 
    symbolic.cor = p > 4, signif.stars = getOption("show.signif.stars"),
    residuals=FALSE, niterations=FALSE, ...)
{
##
  ## doc
  ldoc <- getUserOption("doc")
  if (length(ldoc)==0) ldoc <- 1
  if (ldoc>=1) if (length(tit(x)))
    cat("\n ",tit(x),"\n")
  if (ldoc>=2) if (length(doc(x)))
    cat(" ",paste(doc(x),"\n "))
  ## mlm
  if (inherits(x,"mlm")) return(invisible(print.mregr(x, ...)))
  ## preparation
##-   if (length(dummycoef)==0)
##-     dummycoef <- c(getUserOption("show.dummycoef"),TRUE)[1]
  ## call, fitting fn, residuals
  if (call) {
    if(!is.null(x$call)) {
      cat("\nCall:\n")
      cat(paste(deparse(x$call), sep = "\n", collapse = "\n"),"\n", sep = "")
    }
    cat("Fitting function: ",x$fitfun,"\n")
  }
  df <- x$df
    rdf <- c(x$df.resid,df[2])[1]
  if (residuals) {
    resid <- x$residuals
##-     cat(if (!is.null(x$w) && diff(range(x$w)))
##-         "Weighted ", "Residuals:\n", sep = "")
    cat("Residuals:\n")
    if (rdf > 5) {
      nam <- c("Min", "1Q", "Median", "3Q", "Max")
      rq <- if (length(dim(resid)) == 2)
        structure(apply(t(resid), 1, quantile),
                  dimnames = list(nam, dimnames(resid)[[2]]))
        else structure(quantile(resid), names = nam)
      print(rq, digits = digits, ...)
    } else {
      if (rdf > 0) print(resid, digits = digits, ...)
      else  cat("ALL", df[1],
                "residuals are 0: no residual degrees of freedom!\n")
    }
  }
  ## coefficients
    nsingular <- df[3] - df[1]
    if ((!is.na(nsingular))&&nsingular>0)
        cat("\nCoefficients: (", nsingular,
            " not defined because of singularities)\n", sep = "")
    # else {
    if (!is.null(x$sigma))
      if((!is.finite(x$sigma))||x$sigma<=0)
        cat("\n!!! Error variance is 0 !!!")
  ## coef table
  ltc <- x$termtable
  if (length(ltc)>0) {
    lltc <- TRUE
    if (signif.stars) termcolumns <- union(termcolumns, "p.symb")
    if(!is.null(termcolumns)) {
      if (all(termcolumns=="")) lltc <- FALSE else {
        ljp <- match(termcolumns,colnames(ltc), nomatch=0)
        if (sum(ljp)!=0) 
  ##        warning(":print.regr: no valid columns of  termtable  selected") else
          ltc <- ltc[,ljp,drop=FALSE]
      }
    }
    if (lltc) {
      cat("\nTerms:\n")
      ## round R2.x and p.value
      ljrp <- notna(pmatch(c("R2","p.v"),colnames(ltc)))
      if (length(ljrp)) ltc[,ljrp] <- round(as.matrix(ltc[,ljrp]),max(3,digits))
##-     ljps <- match("p.symb",colnames(ltc), nomatch=0)
##-       if (ljps && all(ltc[,ljps]=="")) ltc <- ltc[,-ljps]
      ltcf <- format(ltc)
##-       lsigst <- signif.stars*( signif.stars && !is.na(ljrp[2]) &&
##-                               any((pv <- ltc[,ljrp[2]]) < 0.1, na.rm=TRUE) )
##-       if (lsigst) {
##-         lsignif <- symnum(pv, corr = FALSE, na = FALSE, 
##-                           cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), 
##-                           symbols = c("***", "**", "*", ".", " "))
##-         ltcf <- cbind(ltcf,format(lsignif))
##-         names(ltcf)[ncol(ltcf)] <- " "
##-     }
      print(ltcf, quote=FALSE)
      if (signif.stars>=1) 
        cat("---\nSignif. codes:  ", attr(x$termtable, "legend"),"\n", sep = "")
    } ## end if(lltc)
  ## --- error block
  } else {
    if (length(x$coef)) {
      cat("\nCoefficients:\n")
      print(x$coef)
    }
  }
##-   if (length(x$binlevels)>0) {
##-     cat("\nFactor(s) with two levels converted to 0-1 variable(s):\n")
##-     print(as.matrix(data.frame(x$binlevels,row.names=0:1)))
##-   }
##    cat("\n")
  ## special for polr
  if (length(x$intercepts)) { # polr
    cat("Intercepts:\n")
    print(x$intercepts)
  }
  ## error
  if (length(x$sigma) && !u.true(attr(x$sigma,"fixed")))
    cat("\nSt.dev.error: ", formatC(x$sigma, digits = digits),
        "  on", rdf, "degrees of freedom\n")
  if (length(x$r.squared)&&!is.na(x$r.squared))
    cat("Multiple R^2: ", formatC(x$r.squared, digits = digits),
        "   Adjusted R-squared:",
        formatC(x$adj.r.squared, digits = digits),"\n"
        )
  if (length(x$fstatistic)>0) {
    cat("F-statistic:  ", formatC(x$fstatistic[1],
            digits = digits), "  on", x$fstatistic[2], "and", x$fstatistic[3],
            "d.f.,  p.value:", formatC(pf(x$fstatistic[1],
            x$fstatistic[2], x$fstatistic[3], lower.tail=FALSE),
                                       digits = digits),
         "\n")
    }
  ## deviances
  if (length(x$deviance)>0) {
    if (length(x$devtable)) print(x$devtable)
    if (length(x$n.censored)) {
      lnc <- 100*x$n.censored/x$n.obs
      cat(paste("\ncensored         ",
        paste(paste(names(lnc), "=", round(lnc,1), "%"), collapse=" ;  ")))
    }
    if (length(x$n.fitout)) {
      lnf <- 100*x$n.fitout/x$n.obs
      cat(paste("\nfit outside limit",
        paste(paste(names(lnf), "=", round(lnf,1), "%"), collapse=" ;  ")))
    }
    cat("\nDistribution: ",x$distrname)
    if (length(x$dispersion))
      cat(".  Dispersion parameter: ",
          if ((!is.null(attr(x$dispersion,"fixed")))&&
              attr(x$dispersion,"fixed"))
             "fixed at ", format(x$dispersion))
    else if (length(x$scale))
      cat(".  Shape parameter (`scale`): ",
          if ((!is.null(attr(x$scale,"fixed")))&&
              attr(x$scale,"fixed"))
             "fixed at ", format(x$scale))
    cat("\nAIC: ", format(x$aic, digits = max(4, digits + 1)), "\n", sep = "  ")
    if (niterations&&length(x$iter)>0)
      cat("Number of iterations:", x$iter, "\n")
  }
  ## --- additional coefficients
  if (x$distrname=="multinomial") {
    cat("\nCoefficients:\n")
    print(t(x$coefficients))
  } else {
    if (length(ltc)&u.true(dummycoef)) {        
      lidf <- match("df",colnames(x$termtable))
      if (is.na(lidf)) {
        if (getOption("verbose"))
          warning(":print.regr: df of coef not available")
      } else { ## dummy coefficients
        mterms <-
          unique(c(row.names(x$termtable)[x$termtable[,"df"]>1],
                   names(attr(x$terms,"dataClasses")[-1]%in%
                         c("factor","ordered")) ))
        if (length(mterms)>0 & length(x$allcoef)>0) {
          imt <- mterms%in%names(x$allcoef)
          mt <- x$allcoef[mterms[imt]]
          if (length(mt)>0) {
            cat("\nCoefficients for factors:\n")
            print.allcoef(mt, digits=digits) }
        } ## else  cat("\n")
      }}
  }
  ## ---- correlation
  correl <- x$correlation
  if (length(correl)>0 && correlation) {
    p <- NCOL(correl)
    if (p > 1) {
      cat("\nCorrelation of Coefficients:\n")
      if (symbolic.cor) {
        symbc <- symnum(correl, symbols=c(" ", ".", ",", "+", "*", "H"))
        symbl <- attr(symbc,"legend")
        attr(symbc,"legend") <- NULL
        print(symbc)
        cat("\nSymbols:  ",symbl,"\n")
      } else {
        correl[!lower.tri(correl)] <- NA
        print(correl[-1, -p, drop = FALSE], digits = digits, na = "")
      }
    }
  }
##  cat("\n")
  invisible(x)
}
## ==========================================================================
summary.regr <- function(object, ...)  object ## dispersion=NULL,
## ==========================================================================
allcoef <- function (object, se = 2, # use.na = TRUE, 
                      df = df.residual(object), ...) 
{
  if (is.atomic(object)) stop("!allcoef! inadequate first argument")
  xl <- object$xlevels
  if (!length(xl)) 
    return(as.list(coef(object)))
  Terms <- terms(object)
  tl <- attr(Terms, "term.labels")
  ## result already available?
  allc <- object$allcoef
  if ((!is.null(allc))&&length(allc)==length(tl)&&
      (is.matrix(allc[[length(allc)]])|!se)) return(allc)
  int <- attr(Terms, "intercept")
  facs <- attr(Terms, "factors")[-1, , drop = FALSE]
  Terms <- delete.response(Terms)
  mf <- object$model  ##! d.c used all.vars
  if (is.null(mf)) mf <- model.frame(object)
  xtnm <- dimnames(facs)[[1]]  ## names  ##! replaces vars
  xtlv <- lapply(mf[,xtnm, drop=FALSE],function(x) levels(x)) ## levels
  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")
  lcontr <- object$contrasts
  lci <- sapply(df.dummy,is.factor)
  lcontr <- lcontr[names(lci)[lci]] ## factors with 1 level have disappeared (?) 
  mm <- model.matrix(Terms, df.dummy, contrasts.arg=lcontr, xlev=xl) ## 
  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
  if (se) {
    cov <- vcov(object)
    if (is.null(cov)) {
      warning(":allcoef: no covariance matrix of coefficients found.",
              " Returning coefficients only")
      se <- FALSE
    }
  }
  asgn <- attr(mm, "assign")
  names(asgn) <- colnames(mm)
  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 <- names(asgn)[asgn == j & !is.na(coef)]  ## columns (logical fails for polr, vcov() too large) !!! was  which
    mmpart <- mm[mmr, mmc, drop=FALSE]
    rrj <- setNames(drop(mmpart %*% coef[mmc]), rnn[mmr])
    if (se) {
      if (any(is.na(rrj)))
        warning(":allcoef: missing coef for term '", tl[j],
                "'. no standard errors etc") else {
##-      if (any(is.na(licov))) ljfail <- c(ljfail,j)  else {
          sej <- sqrt(diag(mmpart %*% cov[mmc,mmc] %*% t(mmpart)))
          rrj <- ciSignif(rrj, sej, df)
        }
    }
    res[[j]] <- rrj
  }
  if (length(ljfail))
    warning(":allcoef: error calculating se for terms ",
            paste(ljfail, collapse=", "))
  if (int > 0) {
    res <- c(list(`(Intercept)` = coef[int]), res)
  }
  class(res) <- "allcoef"
  res
}
## --------------------------------------------------------------------
print.allcoef <- function(x, columns=userOptions("allcoefcolumns"),
                          transpose=FALSE, ...) {
  if (is.null(columns)) columns <- "coefsymb"
  columns[columns=="coef"] <- "estimate"
  csymb <- "coefsymb"%in%columns
  if ("all"%in%columns)  columns <-
      if(csymb)
        c("coefsymb", "se", "ciLow", "ciHigh", "testst",
          "signif", "p.value") else
        c("estimate", "se", "ciLow", "ciHigh", "testst",
          "signif", "p.value")
  for (li in seq_along(x)) {
    xi <- x[[li]]
    if (is.null(dim(xi))) next
    if (csymb)
      xi$coefsymb <-
        if ("p.symb"%in%names(xi))
          paste(format(xi[,1],...), format(xi[,"p.symb"])) else  xi[,1]
    xif <- format(xi[,intersect(columns,names(xi)), drop=FALSE],...)
    xif <- if (ncol(xif)==1 || (nrow(xif)>1 & transpose)) t(xif) else xif
    if (nrow(xif)==1) row.names(xif) <- " "  ## drop row name
    if (ncol(xif)==1) colnames(xif) <- " "  ## drop col name
    if (prod(dim(xif))==1) xif <- as.character(xif[1,1])
    x[li] <- list(xif)
  }
  class(x) <- NULL
  print.default(x, quote=FALSE, ...)
}
## -------------------------------------------------------------------------
vcov.regr <- function(object, ...) {
  cov <- object$covariance 
  if (is.null(cov)) {
    class(object) <- setdiff(class(object),"regr")
    vcov(object)
  }
  cov
}
## --------------------------------------------------------------------
df.residual.regr <- function(object, ...) {
  df <- object$df.residual
  if (is.null(df)) df <- object$df
  if (length(df)>=2) df <- df[2]
  if (is.null(df)) {
    sry <- summary(object)
    df <- sry$df
    if (is.null(df))
      df <- NROW(object$residual)-NROW(object$coefficients)
  }
  if (is.null(df)) df <- Inf
  df
}
## ====================================================================
confint.regr <- function(fitted, ...)
{
if (!inherits(fitted, c("glm","nls"))) {
    class(fitted) <- class(fitted)[-1]
    return(confint(fitted, ...))
}
if (inherits(fitted, "glm")) {
      ## confint needs $coefficients from object (a vector) as well as
      ## from its sumary (a matrix containing 'Std. Error"
  summary <- function(fitted) 
    list(coefficients = cbind(fitted$coefficients,
                                 "Std. Error"=sqrt(diag(fitted$covariance))) )
  class(fitted) <- class(fitted)[-1]
} else { ## workaround: call nls again, since  profile.nls  is difficult to adapt...
  call <- fitted$call
  call$start <- fitted$coefficients
  call$nonlinear <- NULL
  call[[1]] <- as.name("nls")
  fitted <- eval(call, parent.frame())
}
  confint(fitted, ...)
}
## ==========================================================================
drop1.regr <-
  function (object, scope=NULL, scale = 0, test = NULL, k = 2,
            sorted = FALSE, add=FALSE, ...)
{
  ## Purpose:    drop1/add1 for regr objects
  ## ----------------------------------------------------------------------
  lfam <- object$distrname
  lres <- object$residuals
  if (is.null(test)) test <- if (is.null(lfam)) "none" else {
    if ((lfam=="gaussian"&&as.character(object$fitfun)%in%c("lm","roblm"))|
        ((lfam=="binomial"|lfam=="poisson")&&object$dispersion>1)) {
          if (inherits(object,"mlm")) "Wilks" else "F" }
    else "Chisq"
  }
  if (length(scope)==0) {
    scope <- if (add) terms2order(object) else drop.scope(object)
  } else {
##-     if (!is.character(scope))
##-             scope <- attr(terms(update.formula(object, scope)),
##-                 "term.labels")
    if (is.character(scope))
      as.formula(paste("~",paste(scope, collapse="+"))) else scope
  }
  if (length(scope)==0) {  ## || !is.formula(scope) ## drop.scope is character
    warning(":drop1/add1.regr: no valid scope")
    ldr1 <- data.frame(Df = NA, "Sum of Sq" = NA, RSS =NA, AIC = NA,
                      row.names = "<none>")
    return(ldr1)
  }
  ##  lform <- update(formula(object), scope !!!)
  ## !!! model.frame for finding the valid rows
  class(object) <- setdiff(class(object), "regr")
  fcall <- object$funcall
  if (!is.null(fcall)) object$call <- fcall
##-   dfm <-  object$df
  ##-   object$df <- dfm[setdiff(names(dfm),"residual")] ## survreg !
##-   if (inherits(object, c("survreg","coxph")))
##-     object$df <- object$df["original"]
##
  dr1 <- if (add) { ## ------------ add
    if (class(object)[1]=="lmrob")
        stop("!add1.regr! 'add1' not (yet) available for 'lmrob' objects")
    ldata <- eval(object$call$data)
    li <- row.names(ldata)%in%RNAMES(object$residuals)
    if (is.null(ldata[li,])) stop("!drop1.regr! no data found ")
    lvars <-unique(c(all.vars(formula(object)),
                     if (is.formula(scope)) all.vars(scope) else scope))
    lvars <- lvars[lvars%in%names(ldata)]
    linotna <- li & !apply(is.na(ldata[,lvars]),1,any)
    lnobs <- sum(linotna)
    lnrd <- sum(li)
    lfc <- object$funcall ## the call to the effective R function
    if (is.null(lfc)) lfc <- object$call
    if (lnobs!= lnrd) {
      warning(":add1.regr: refitting object to ",lnobs," / ",lnrd,
              " observations due to missing values")
      if(!is.null(lsubs <- eval(lfc$subset))) {
        lnsubs <- rep(FALSE,length(linotna))
        lnsubs[lsubs] <- TRUE
        linotna <- linotna &!lnsubs
      }
      lfc$subset <- linotna
      object <- eval(lfc)
##-       object$call[[1]] <-
##-         if (is.null(lfc)) as.name(class(object)[1]) else
##-            lfc[[1]]
##-       object <- update(object, subset=linotna)
#      environment(object$call$formula) <- environment()
##-       class(object) <- setdiff(class(object), "regr")
    }
    if (!all(linotna)) { ## needed if NA's have been generated by transformations
      lfc$subset <- linotna
      object <- eval(lfc)
    }
    add1(object, scope=scope, scale=scale, test=test, k=k, ...)
  } else {  ## ---------------------  drop
    if (class(object)[1]%in%c("lmrob")) ## to be expanded
       drop1Wald(object, test="F", ...) else {
    ldata <- object$allvars # eval(object$call$data)
    if (is.null(ldata)) stop("!drop1.regr! no data found ")
    ## all predictors must get the same missing observations
    lina <- apply(is.na(ldata),1,any)  
    if (any(lina)) ldata[lina,] <- NA
    object$call$data <- ldata
    drop1(object, scope=scope, scale=scale, test=test, k=k, ...)
    }
  }
##-   rnm <- row.names(dr1)
##-   row.names(dr1) <- paste(ifelse(substring(rnm,1,1)=="<","",
##-                                  if (add) "+ " else "- "),rnm,sep="")
  attr(dr1,"drop") <- !add
##-   if(add) attr(dr1,"ndropped") <- lndiff
  if (sorted) {
    lsrt <- notna(match(c("AIC","p.value"),colnames(dr1)))
    if (length(lsrt)) dr1 <- dr1[order(dr1[, lsrt[1]]), ]
  }
  dr1
}
## ==========================================================================
add1.regr <-
  function (object, scope=NULL, scale = 0, test = NULL, k = 2,
            sorted = FALSE, ...)
{
  ## Purpose:    add1 for regr objects
  ## ----------------------------------------------------------------------
  if (!is.null(scope)) {
    if (is.character(scope)) scope <- paste(scope,collapse="+")
    if (is.formula(scope)) scope <- last(as.character(scope))
    scope <- as.formula(paste("~ ",formula(object)[3],"+",scope))
  }
  drop1.regr(object, scope=scope, scale=scale, test=test, k=k,
             sorted=sorted, add=TRUE, ...)
}
## ==========================================================================
drop1Wald <-
  function (object, scope=NULL, scale = 0, test = c("none", "Chisq", "F"),
            k = 2, ...) 
{
    x <- model.matrix(object)
    offset <- model.offset(model.frame(object))
    n <- nrow(x)
    asgn <- attr(x, "assign")
    lterms <- terms(object)
    tl <- attr(lterms, "term.labels")
    attr(lterms, "order") <-  rep(1,length(tl))
    if (is.null(scope)) 
      scope <- tl # drop.scope(lterms)
    else {
        if (!is.character(scope)) 
            scope <- attr(terms(update.formula(object, scope)), 
                "term.labels")
        if (!all(match(scope, tl, 0L) > 0L)) 
            stop("scope is not a subset of term labels")
    }
    ndrop <- match(scope, tl)
    ns <- length(scope)
    rdf <- object$df.residual
    lsig <- c(object$sigma, object$scale)[1]
    chisq <- lsig^2 * rdf
    ## sum(weighted.residuals(object)^2, na.rm = TRUE)
    ## deviance.lm(object)
    dfs <- numeric(ns)
    RSS <- numeric(ns)
    cov <- object$cov.unscaled
    if (is.null(cov)) cov <- object$covariance/lsig^2
    if (length(cov)==0) stop("!drop1Wald! no covariance matrix found")
    cf <- object$coefficients
    ##- jj <- match(names(cf),colnames(cov), nomatch=0)
    ##-     if (!(any(jj==0)&&all(is.na(cf[jj==0]))))
    jj <- match(colnames(cov),names(cf), nomatch=0)
    if (any(jj==0))
      warning(":drop1Wald: coefficient(s) and cov. matrix may not correspond")
    coef <- cf[jj]
    asgn <- asgn[jj]
    if (any(names(coef[!is.na(coef)])%nin%names(coef)))
      stop("!drop1Wald! coefficient(s) not appearing in covariance matrix")
##-     y <- object$residuals + predict(object)
    for (i in 1:ns) {
        ii <- which(asgn==ndrop[i]) ## seq_along(asgn)[asgn == ndrop[i]]
        RSS[i] <- if (length(ii)==1) coef[ii]^2/cov[ii,ii] else
          coef[ii]%*%solve(cov[ii,ii])%*%coef[ii]  ## !!! REPLACE THIS
        dfs[i] <- length(ii)
##-         if (all.cols) 
##-             jj <- setdiff(seq(ncol(x)), ii)
##-         else jj <- setdiff(na.coef, ii)
##-         z <- if (iswt) 
##-             lm.wfit(x[, jj, drop = FALSE], y, wt, offset = offset)
##-         else lm.fit(x[, jj, drop = FALSE], y, offset = offset)
##-         dfs[i] <- z$rank
##-         oldClass(z) <- "lm"
##-         RSS[i] <- deviance(z)
    }
    scope <- c("<none>", scope)
    dfs <- c(c(object$rank,object$df)[1], dfs)
    RSS <- chisq + c(0, RSS)
    if (scale > 0) 
        aic <- RSS/scale - n + k * dfs
    else aic <- n * log(RSS/n) + k * dfs
##-     dfs <- dfs[1] - dfs
##-     dfs[1] <- NA
    aod <- data.frame(Df = dfs, "Sum of Sq" = c(NA, RSS[-1] - 
        RSS[1]), RSS = RSS, AIC = aic, row.names = scope, check.names = FALSE)
    if (scale > 0) 
        names(aod) <- c("Df", "Sum of Sq", "RSS", "Cp")
    test <- match.arg(test)
    if (test == "Chisq") {
        dev <- aod$"Sum of Sq"
        if (scale == 0) {
            dev <- n * log(RSS/n)
            dev <- dev - dev[1]
            dev[1] <- NA
        }
        else dev <- dev/scale
        df <- aod$Df
        nas <- !is.na(df)
        dev[nas] <- pchisq(dev[nas], df[nas], lower.tail = FALSE)
        aod[, "Pr(Chi)"] <- dev
    }
    else if (test == "F") {
        dev <- aod$"Sum of Sq"
        dfs <- aod$Df
        rdf <- object$df.residual
        rms <- aod$RSS[1]/rdf
        Fs <- (dev/dfs)/rms
        Fs[dfs < 1e-04] <- NA
        P <- Fs
        nas <- !is.na(Fs)
        P[nas] <- pf(Fs[nas], dfs[nas], rdf, lower.tail = FALSE)
        aod[, c("F value", "Pr(F)")] <- list(Fs, P)
    }
    head <- c("Single term deletions (Wald test)", "\nModel:",
              deparse(as.vector(formula(object))), 
        if (scale > 0) paste("\nscale: ", format(scale), "\n"))
    class(aod) <- c("anova", "data.frame")
    attr(aod, "heading") <- head
    aod
}

## ==========================================================================
step <- function(object, ...)
  UseMethod("step")

step.default <- stats::step
## step.default <- get("step", pos="package:stats")
#### !!! sollte das anders heissen? step.default <- stats::step  ???

step.regr <- function (object, scope=terms2order(object), scale = 0,
  direction = c("both", "backward","forward"), trace = FALSE, keep = NULL,
  steps = 1000, k = 2, ...)
{
  ## Purpose:
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 15 May 2012, 07:58
    mydeviance <- function(x, ...) {
        dev <- deviance(x)
        if (!is.null(dev))
            dev
        else extractAIC(x, k = 0)[2L]
    }
    cut.string <- function(string) {
        if (length(string) > 1L)
            string[-1L] <- paste0("\n", string[-1L])
        string
    }
    re.arrange <- function(keep) {
        namr <- names(k1 <- keep[[1L]])
        namc <- names(keep)
        nc <- length(keep)
        nr <- length(k1)
        array(unlist(keep, recursive = FALSE), c(nr, nc), list(namr,
            namc))
    }
    step.results <- function(models, fit, object, usingCp = FALSE) {
        change <- sapply(models, "[[", "change")
        rd <- sapply(models, "[[", "deviance")
        if(is.matrix(rd)) rd <- rd[2,]
        dd <- c(NA, abs(diff(rd)))
        rdf <- sapply(models, "[[", "df.resid")
        ddf <- c(NA, diff(rdf))
        AIC <- sapply(models, "[[", "AIC")
        heading <- c("Stepwise Model Path \nAnalysis of Deviance Table",
            "\nInitial Model:", deparse(formula(object)), "\nFinal Model:",
            deparse(formula(fit)), "\n")
        aod <- data.frame(Step = I(change), Df = ddf, Deviance = dd,
            `Resid. Df` = rdf, `Resid. Dev` = rd, AIC = AIC,
            check.names = FALSE)
        if (usingCp) {
            cn <- colnames(aod)
            cn[cn == "AIC"] <- "Cp"
            colnames(aod) <- cn
        }
        attr(aod, "heading") <- heading
        fit$anova <- aod
        fit
      }
    ## end step.results
    Terms <- terms(object)
    object$call$formula <- object$formula <- Terms
##-     ## !!! begin mod by WSt for regr objects
##-     ldata <- eval(object$call$data)
##-     if (is.null(ldata)) stop("!step.regr! no data found ")
##-     lvars <- all.vars(object$formula)
##-     lina <- apply(is.na(ldata[,lvars]),1,sum)>0
##-     lnotna <- sum(lina)
##-     ldt <- ldata[!lina,]
##-     lcdata <- object$call$data
##-     object$call$data <- ldt
##-     ## !!! end
    md <- missing(direction)
    direction <- match.arg(direction)
    backward <- direction == "both" | direction == "backward"
    forward <- direction == "both" | direction == "forward"
    if (missing(scope)) {
        fdrop <- numeric()
        fadd <- attr(Terms, "factors")
        if (md)
            forward <- FALSE
    }
    else {
        if (is.list(scope)) {
            fdrop <- if (!is.null(fdrop <- scope$lower))
                attr(terms(update.formula(object, fdrop)), "factors")
            else numeric()
            fadd <- if (!is.null(fadd <- scope$upper))
                attr(terms(update.formula(object, fadd)), "factors")
        }
        else {
            fadd <- if (!is.null(fadd <- scope))
                attr(terms(update.formula(object, scope)), "factors")
            fdrop <- numeric()
        }
    }
    models <- vector("list", steps)
    if (!is.null(keep))
        keep.list <- vector("list", steps)
    n <- nobs(object, use.fallback = TRUE)
    fit <- object
    bAIC <- extractAIC(fit, scale, k = k, ...)
    edf <- bAIC[1L]
    bAIC <- bAIC[2L]
    if (is.na(bAIC))
        stop("AIC is not defined for this model, so `step` cannot proceed")
    nm <- 1
    if (trace) {
      cat("Start:  AIC=", format(round(bAIC, 2)), "\n",
          cut.string(deparse(formula(fit))), "\n\n", sep = "")
      utils::flush.console()
    }
  models[[nm]] <-
    list(deviance = mydeviance(fit), df.resid = n - edf, change = "",
         AIC = bAIC)
  if (!is.null(keep))
    keep.list[[nm]] <- keep(fit, bAIC)
  usingCp <- FALSE
  ## ------------------------
  while (steps > 0) {
    steps <- steps - 1
    AIC <- bAIC
    ffac <- attr(Terms, "factors")
    scope <- factor.scope(ffac, list(add = fadd, drop = fdrop))
    aod <- NULL
    change <- NULL
    ## backward
    if (backward && length(scope$drop)) {
      aod <- drop1(fit, scope$drop, scale = scale, trace = trace,
                   k = k, sorted=FALSE, ...)
      rn <- row.names(aod)
      row.names(aod) <- c(rn[1L], paste("-", rn[-1L], sep = " "))
      if (any(aod$Df == 0, na.rm = TRUE)) {
        zdf <- aod$Df == 0 & !is.na(aod$Df)
        change <- rev(rownames(aod)[zdf])[1L]
      }
    }
    ## forward
    if (is.null(change)) {
      if (forward && length(scope$add)) {
        aodf <- add1(fit, scope$add, scale = scale, trace = trace,
                     k = k, ...)
        rn <- row.names(aodf)
        row.names(aodf) <- c(rn[1L], paste("+", rn[-1L], sep = " "))
        aod <- if (is.null(aod)) aodf  else {
          names(aodf) <- names(aod)
          rbind(aod, aodf[-1, , drop = FALSE])
        }
      }
    }
    ## backward or forward
    attr(aod, "heading") <- NULL
    nzdf <- if (!is.null(aod$Df))
              aod$Df != 0 | is.na(aod$Df)
    aod <- aod[nzdf, ]
    if (is.null(aod) || ncol(aod) == 0)
      break
    nc <- match(c("Cp", "AIC"), names(aod))
    nc <- nc[!is.na(nc)][1L]
    o <- order(aod[, nc])
    if (trace)
      print(aod[o, ])
    if (o[1L] == 1)  break
    change <- rownames(aod)[o[1L]]
    ## update
    usingCp <- match("Cp", names(aod), 0L) > 0L
    ##    if (is.null(change))  break  else {
    fit <- update(fit, paste("~ .", change), evaluate = FALSE)
    fit <- eval.parent(fit)
    nnew <- nobs(fit, use.fallback = TRUE)
    if (all(is.finite(c(n, nnew))) && nnew != n) {
      warning(":step.regr: number of rows in use has changed: \n  ",
              nnew," observations instead of ", n)
      n <- nnew
    }
    Terms <- terms(fit)
    bAIC <- extractAIC(fit, scale, k = k, ...)
    edf <- bAIC[1L]
    bAIC <- bAIC[2L]
    ## output
    if (trace) {
      cat("\nStep:  AIC=", format(round(bAIC, 2)), "\n",
          cut.string(deparse(formula(fit))), "\n\n", sep = "")
      utils::flush.console()
    }
    ##        if (bAIC >= AIC + 1e-07)
    ##  else   break
    nm <- nm + 1
    models[[nm]] <- list(deviance = mydeviance(fit), df.resid = n - edf,
                         change = change, AIC = bAIC)
    if (!is.null(keep))
      keep.list[[nm]] <- keep(fit, bAIC)
  }
  if (!is.null(keep))
    fit$keep <- re.arrange(keep.list[seq(nm)])
  step.results(models = models[seq(nm)], fit, object, usingCp)
}
## ==========================================================================
terms2order <- function(object, squared = TRUE, interactions = TRUE)
{
  ## Purpose:
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 15 Oct 2009, 16:01
  ltadd <- NULL
  if (squared) {
#    lvrs <- all.vars(formula(object)[-2])
    if (!is.list(object)||length(lterms <- object$terms)==0)
      stop("!terms2order! first argument has no 'terms' component")
    ltl <- attr(lterms, "term.labels")
    lts <- c(grep("\\^",ltl),grep(":",ltl))
    if (length(lts)) ltl <- ltl[-lts]
    ltn <- ltl[attr(object$terms,"dataClasses")[ltl]=="numeric"]
    if (length(ltn)) ltadd <- paste(paste("I(", ltn,"^2)",sep=""), collapse="+")
  }
  if (interactions) {
    ltint <- "(.)^2"
    if (is.null(ltadd)) ltadd <- ltint else ltadd <- paste(ltadd, ltint, sep=" + ")
  }
  if (is.null(ltadd)) {
    warning(":terms2order: nothing to add")
    return(formula(object))
  }
  ltadd <- paste("~.+", ltadd)
##-     attr(terms(update.formula(object, ltadd)), "term.labels")
  update.formula(object, ltadd)
}
## ==========================================================================
fitted.regr <-
  function (object, type=NULL, ...)
{
  if (is.null(type)&&pmatch("fitted",names(object),nomatch=0))
    return( naresid(object$na.action, object$fitted) )
  lres <- object$residuals
  if (inherits(lres, "condquant"))
    structure(lres[,"fit"], names=row.names(lres))  else  {
      class(object) <- setdiff(class(object), "regr")
      predict(object, type=type, ...)
    }
}
## ==========================================================================
predict.regr <-
function (object, newdata = NULL, scale = object$sigma, type = NULL, ...)
  ## bug: if used with NULL newdata, predictions will be produced
  ## for obs in model.matrix, which excludes those with missing values
{
##-   lglm <- inherits(object,"glm")
##-   lmeth <- object$call$method
##-   lnls <- length(lmeth)>0 && lmeth=="nls"
  if (length(type)==0)
    type <- if (inherits(object,"glm")) "link" else
  if (inherits(object, "polr")) "link" else "response"
  ## !!!
  if (object$fitfun=="rlm")
    if (!is.matrix(object[["x"]]))
      object$x <- model.matrix(formula(object), data=object$allvars) ## !!! was $model
##-   if (length(scale)==0) scale <- c(object$sigma,1)[1]
  class(object) <- class(object)[-1]
##-   if (missing(newdata) || length(newdata)==0) {
##-     lpred <- if (lglm)
##-       predict.glm(object, type=type, se.fit=se.fit,
##-                   dispersion=lscale^2, terms=terms, na.action = na.action )
##-     else {
##-       if (lnls) predict.nls(object, type=type, se.fit=se.fit,
##-                             na.action = na.action )  else
##-       predict.lm(object, type=type, se.fit=se.fit,
##-                     terms=terms, na.action = na.action )
##-     }
##-     lpred <- predict(object, type=type, se.fit=se.fit,
##-                   dispersion=lscale^2, terms=terms, na.action = na.action )
##-   } else {
    ldt <- newdata
    ## analyze variables
  if (!is.null(ldt)) {
    for (lvn in names(ldt)) {
      lv <- ldt[[lvn]]
      ## binary factors
      if (is.factor(lv)&&match(lvn,names(object$binlevels),nomatch=0)>0)
        ldt[[lvn]] <- match(lv,object$binlevels[[lvn]])-1
      else  if (is.factor(lv))   ldt[[lvn]] <- factor(lv)
    }
  }
  if (is.null(ldt))
    predict(object, type=type, scale=object$sigma,
            dispersion=object$dispersion^2, ... )  else
    predict(object, newdata=ldt, type=type, scale=object$sigma,
            dispersion=object$dispersion^2, ... )
}
## ==========================================================================
##- extractAIC.regr <- function (fit, scale = 0, k = 2, ...)
##- {  ##- #  AIC, divided by n to allow for comparing models with different n
##- ##-   lres <- fit$residuals
##- ##-   if (is.null(lres)) {
##- ##-     lfit <- fit$fitted
##- ##-     fit$fitted.values <- notna(lfit)
##- ##-   } else  fit$residuals <- notna(lres)
##-   class(fit) <- setdiff(class(fit),"regr")
##-   extractAIC(fit, scale = scale, k = k, ...)
##- }
## ==========================================================================
vif.regr <- function(mod, cov, mmat)
{
  ## Purpose:   vif.lm  of library  car
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: modified by Werner Stahel, Date: 11 Mar 2005, 09:18
##-     v <- mod$sigma^2* mod$cov.unscaled
##-     cls <- dimnames(model.matrix(mod))[[2]]%in%dimnames(v)[[2]]
##-                                         # needed for singular cases
##-     assign <- attributes(model.matrix(mod))$assign[cls]
    cls <- dimnames(mmat)[[2]]%in%dimnames(cov)[[2]]
##-                                         # needed for singular cases
    assign <- attr(mmat,"assign")[cls]
    terms <- labels(terms(mod))
    n.terms <- length(terms)
    if (n.terms < 2) {
##-         stop("model contains fewer than 2 terms")
      return(matrix(1,1,3))
    }
    if (length(cov)==0) { # ||n.terms!=nrow(cov)|nrow(cov)!=ncol(cov)
      warning(":vif.regr: mod$cov.unscaled  is inappropriate. no vifs")
      return(matrix(NA,n.terms,3))
    }
    if (names(coefficients(mod)[1]) == "(Intercept)") {
        cov <- cov[-1, -1]
        assign <- assign[-1]
    }
    else if (mod$fitfun%nin%c("polr","coxph","survreg"))
      warning("No intercept: vifs may not be sensible.")
    sd <- 1/sqrt(diag(cov))
    if (any(!is.finite(sd))) {
      warning(":vif.regr: zero variances of estimates. no R2x")
      return(NULL)
    }
    R <- cov/outer(sd,sd)
    result <- matrix(0, n.terms, 3)
    rownames(result) <- terms
    colnames(result) <- c("GVIF", "Df", "GVIF^(1/2Df)")
    for (term in 1:n.terms) {
      subs <- which(assign == term)
      result[term, 1] <- det(as.matrix(R[subs, subs])) *
        det(as.matrix(R[-subs,-subs]))/det(R)
      result[term, 2] <- length(subs)
    }
    result[, 3] <- result[, 1]^(1/(2 * result[, 2]))
    result
}
## =================================================================
drop1.multinom <-
  function (object, scope, test = c("Chisq","none"), ...)
{
    if (!inherits(object, "multinom"))
        stop("Not a multinom fit")
    if (missing(scope))
        scope <- drop.scope(object)
    else {
        if (!is.character(scope))
            scope <- attr(terms(update.formula(object, scope)),
                "term.labels")
        if (!all(match(scope, attr(object$terms, "term.labels"),
            nomatch = FALSE)))
            stop("scope is not a subset of term labels")
    }
    ns <- length(scope)
    ans <- matrix(nrow = ns + 1, ncol = 2, dimnames = list(c("<none>",
        scope), c("Df", "AIC")))
    ans[1, ] <- c(object$edf, object$AIC)
    if (test[1]=="Chisq") ans <- cbind(ans,Chisq=NA, p.value=NA)
    env <- environment(formula(object))
    for(i in seq(ns)) {
      tt <- scope[i]
##        cat("trying -", tt, "\n")
      nfit <- update(object, as.formula(paste("~ . -", tt)),
                     evaluate = FALSE)
      nfit <- eval(nfit, envir=env) # was  eval.parent(nfit)
##-        nobject <- update(object, paste("~ . -", tt))
      if (nfit$edf == object$edf)
        nfit$AIC <- NA
      ans[i+1, ] <- c(nfit$edf, nfit$AIC,
                    if (test[1]=="Chisq") unlist(anova(object,nfit)[2,6:7]))
    }
    as.data.frame(ans)
}
## =================================================================
## FIXME: this is nowhere called, nor documented ...
summary.mregr <- function(object, ...)
{
  ## Purpose:   collect results for mlm object
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 15 Feb 2005, 09:11
  lsum <- summary(object)
  lts <- ltp <- NULL
  for (ly in 1:length(lsum)) {
    lts <- cbind(lts,c(lsum[[ly]][["sigma"]],lsum[[ly]][["r.squared"]],
                       lsum[[ly]][["fstatistic"]]))
    ltp <- cbind(ltp,lsum[[ly]][["coefficients"]][,4])
  }
  lts[4,] <- pf(lts[3,],lts[4,],lts[5,], lower.tail=FALSE)
  lts <- lts[1:4,]
  dimnames(ltp) <- dimnames(object$coefficients)
  dimnames(lts) <- list(c("sigma","r.squared","fstatistic","p-value"),
                        dimnames(ltp)[[2]])
  list(coefficients=object$coefficients, pvalues=ltp, stats=lts)
}
## =================================================================

##- @title drop1() method for multivariate lm()`s
##- @param object
##- @param scope
##- @param test
##- @param total
##- @param add
##- @param ...
##- @note R`s  stats:::drop1.mlm() just has stop("no `drop1` method for \"mlm\" models")
drop1.mlm <- function (object, scope = NULL,
                       test = c("Wilks", "Pillai", "Hotelling-Lawley", "Roy"),
                       total=TRUE, add=FALSE, ...)
{
  ## add=TRUE   do add instead of drop1
    Pillai <- function(eig, q, df.res) {
        test <- sum(eig/(1 + eig))
        p <- length(eig)
        s <- min(p, q)
        n <- 0.5 * (df.res - p - 1)
        m <- 0.5 * (abs(p - q) - 1)
        tmp1 <- 2 * m + s + 1
        tmp2 <- 2 * n + s + 1
        c(test, (tmp2/tmp1 * test)/(s - test), s * tmp1, s *
            tmp2)
    }
    Wilks <- function(eig, q, df.res) {
        test <- prod(1/(1 + eig))
        p <- length(eig)
        tmp1 <- df.res - 0.5 * (p - q + 1)
        tmp2 <- (p * q - 2)/4
        tmp3 <- p^2 + q^2 - 5
        tmp3 <- if (tmp3 > 0)
            sqrt(((p * q)^2 - 4)/tmp3)
        else 1
        c(test, ((test^(-1/tmp3) - 1) * (tmp1 * tmp3 - 2 * tmp2))/p/q,
            p * q, tmp1 * tmp3 - 2 * tmp2)
    }
    HL <- function(eig, q, df.res) {
        test <- sum(eig)
        p <- length(eig)
        m <- 0.5 * (abs(p - q) - 1)
        n <- 0.5 * (df.res - p - 1)
        s <- min(p, q)
        tmp1 <- 2 * m + s + 1
        tmp2 <- 2 * (s * n + 1)
        c(test, (tmp2 * test)/s/s/tmp1, s * tmp1, tmp2)
    }
    Roy <- function(eig, q, df.res) {
        p <- length(eig)
        test <- max(eig)
        tmp1 <- max(p, q)
        tmp2 <- df.res - tmp1 + q
        c(test, (tmp2 * test)/tmp1, tmp1, tmp2)
    }
##-     if (!is.null(object$drop1)) return (object$drop1)
    if (!(inherits(object, "maov") || inherits(object, "mlm")))
        stop("object must be of class \"maov\" or \"mlm\"")
    test <- match.arg(test)
    asgn <- object$assign[object$qr$pivot[1:object$rank]]
    tl <- attr(object$terms, "term.labels")
## scope
    if (is.null(scope))
        scope <- if (add) attr(terms(update.formula(object, ~(.)^2)),
                "term.labels") else
          drop.scope(object)
    else {
        if (!is.character(scope))
            scope <- attr(terms(update.formula(object, scope)),
                "term.labels")
##-         if (!all(match(scope, tl, FALSE)))
        if (!(add||all(match(scope, tl, FALSE))))
            stop("!drop1.mlm! scope is not a subset of term labels")
    }
    ns <- length(scope)
    rdf <- object$df.residual
    res <- resid(object)
## ::: needed for finding the data later
    ldata <- eval(object$call$data,
                 envir=environment(formula(object)))
    ladd <- 1
    if (add) {
      ladd <- -1
      lna <- i.add1na(object, scope)
      if (!is.null(lna)) res[lna,1] <- NA
    }
    res <- nainf.exclude(res)
    lna <- attr(res,"na.action")
    if (!is.null(lna)) ldata <- ldata[-lna,]
## full model
    rss <- crossprod(as.matrix(res))
    rss.qr <- qr(rss)
    if (rss.qr$rank < NCOL(res))
      stop(paste("!drop1.mlm! residuals have rank", rss.qr$rank, "<",
                 ncol(res)))
    stats <- matrix(NA,length(scope),4)
    dimnames(stats) <- list(scope,c(test,"F.stat","dfnum","dfden"))
    tstfn <- switch(test, Pillai = Pillai, Wilks = Wilks,
                    "Hotelling-Lawley" = HL, HL = HL, Roy = Roy)
    object$call[[1]] <- as.name("lm")
## loop through scope
    for (lsc in scope) {
      lfo <- as.formula(paste(if (add) "~.+" else "~.-",lsc))
      lrg <- update(object, lfo, data=ldata, model=FALSE) # ,data=data
      dfj <- ladd * (lrg$df.residual - rdf)
      bss <- ladd * (crossprod(resid(lrg))-rss)
      eigs <- Re(eigen(qr.coef(rss.qr, bss),symmetric = FALSE)$values)
      stats[lsc,] <- tstfn(eigs, dfj, rdf)
    }
    ldf <- stats[1,3:4]
    names(ldf) <- c("numerator","denominator")
    if (total) {
      lpr <- predict(object)
      if (length(lna)) lpr <- lpr[-lna,] # drop rows with NA
      yy <- scale(res + lpr, scale=FALSE)
      bss <- crossprod(yy)-rss
      eigs <- Re(eigen(qr.coef(rss.qr, bss), symmetric = FALSE)$values)
      stats <- rbind(stats, "<total>"= tstfn(eigs, object$df[1], rdf))
    }
    data.frame(stats,
               p.value = pf(stats[,2],stats[,3],stats[,4], lower.tail = FALSE))
##    attr(stats,"df") <- ldf
##    stats
} ## {drop1.mlm}
## ==========================================================================
add1.mlm <-
  function (object, scope=NULL, test = c("Wilks", "Pillai", "Hotelling-Lawley", "Roy"), ...)
{
  ## Purpose:    add1 for regr objects
  ## ----------------------------------------------------------------------
  drop1.mlm(object, scope=scope, test=test, total=FALSE, add=TRUE, ...)
}
## ===========================================================================
i.add1na <- function (object, scope)
{
  ##  determine rows with NA`s in model.frame for expanded model
    Terms <-
      terms(update.formula(object, paste("~.+",paste(scope, collapse = "+"))))
    fc <- object$call
    fc$formula <- Terms
    fob <- list(call = fc, terms = Terms)
    class(fob) <- oldClass(object)
    m <- model.frame(fob, xlev = object$xlevels)
    r <- cbind(resid(object))
    if (nrow(r)!=nrow(m)) {
      warning(gettextf("!add1! using the %d/%d rows from a combined fit",
                nrow(m), nrow(r)), domain = NA)
      lna <- !row.names(r)%in%row.names(m)
    }
    else lan <- NULL
  }
## ==========================================================================
## currently only called from print.regr():
print.mregr <- function(x, ...)
{
  ## Purpose:   collect results for mregr object
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: Feb 2008
  f.prv <- function(x) paste(paste(names(x),x,sep=" = "),collapse=", ")
  cat("\nCall:\n")
  cat(paste(deparse(x$call), sep = "\n", collapse = "\n"),"\n", sep = "")
  cat("\nCoefficients:\n")
  print(x$coefficients)
  cat("\nP-values:\n")
  print(round(x$pvalues,4))
  cat("\nStatistics for individual response models:\n")
  print(x$stats)
  cat("\nResidual degrees of freedom: ",x$df,"\n")
  ldr <- x$drop1
  if (!is.null(ldr)) {
    cat("\nMultivariate tests for all responses\n  Degrees of freedom: ",
        f.prv(attr(ldr,"df")),"\n")
    print(ldr[,])
  }
  invisible(x)
}

## ===================================================
fitted.polr <- function(object, type="link", na.action=object, ...) {
  if (pmatch(type,"link",nomatch=0)) {
    lfit <- object$linear.predictor
    if (length(lfit)==0) { # if called by original polr
      Terms <- delete.response(object$terms)
      environment(Terms) <- environment() ## ! WSt
      ## from predict.polr
      m <- object$model
      if (length(m)==0)
        m <- model.frame(Terms, object$data, na.action = function(x) x,
                         xlev = object$xlevels)
      if (!is.null(cl <- attr(Terms, "dataClasses")))
        .checkMFClasses(cl, m)
      X <- model.matrix(Terms, m, contrasts = object$contrasts)
      xint <- match("(Intercept)", colnames(X), nomatch = 0)
      if (xint > 0)  X <- X[, -xint, drop = FALSE]
      lfit <- drop(X %*% object$coefficients)
    }
  } else 
    lfit <- object$fitted
  if (type=="class")
    lfit <- factor(max.col(lfit), levels = seq_along(object$lev),
            labels = object$lev)
  naresid(na.action,lfit)
}
## ===================================================
predict.polr <-
  function (object, newdata=NULL,
            type = c("class", "probs", "link"), ...)
  ## type link added by WSt, newdata=NULL
{
    if (!inherits(object, "polr"))
        stop("not a \"polr\" object")
    type <- match.arg(type)
    if (length(newdata)==0) {
      if (type=="link")
        eta <- fitted(object, type="link", na.action=NULL)
        Y <- object$fitted
      na.action <- object$na.action
    }
    else {
      na.action <- NULL
        newdata <- as.data.frame(newdata)
        Terms <- delete.response(object$terms)
        environment(Terms) <- environment() ## ! WSt
        m <- model.frame(Terms, newdata, na.action = function(x) x,
            xlev = object$xlevels)
        if (!is.null(cl <- attr(Terms, "dataClasses")))
            .checkMFClasses(cl, m)
        X <- model.matrix(Terms, m, contrasts = object$contrasts)
        xint <- match("(Intercept)", colnames(X), nomatch = 0)
    ## changed by WSt
##-         if (xint > 0)
##-             X <- X[, -xint, drop = FALSE]
        eta <- drop(X[,names(object$coefficients)] %*% object$coefficients)
        n <- nrow(X)
        q <- length(object$zeta)
##      pgumbel <- function(q) exp(pweibull(log(q))) # ???
        pfun <- switch(object$method, logistic = plogis, probit = pnorm,
            cloglog = prevgumbel, cauchit = pcauchy)
        cumpr <- matrix(pfun(matrix(object$zeta, n, q, byrow = TRUE) -
            eta), , q)
        Y <- t(apply(cumpr, 1, function(x) diff(c(0, x, 1))))
        dimnames(Y) <- list(rownames(X), object$lev)
    }
##-     if (newdata) && !is.null(object$na.action))
##-         Y <- napredict(object$na.action, Y)
    switch(type, class = {
        Y <- factor(max.col(Y), levels = seq_along(object$lev),
            labels = object$lev)
    }, probs = {
    }, link = { Y <- eta })
    Y <- napredict(na.action,Y)
    drop(Y)
}
## ===================================================================
condquant <- function(x, dist="normal", sig=1, randomrange=0.9)
{
  ## Purpose:   conditional quantiles and random numbers
  ## works only for centered scale families
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 24 Oct 2008, 09:30
##-   fp <- switch(dist, normal="pnorm(x,0,sig)", gaussian="pnorm(x,0,sig)",
##-                logis="plogis(x,0,sig)", revgumbel="prevgumbel(x,0,sig)",
##-                # lognormal="pnorm(log(x),0,sig)"
##-                ## x is log regsidual ?
##-                )
  if (length(x)==0) stop("!condquant! bug: no data")
  fp <- switch(dist, normal=pnorm, gaussian=pnorm, unif=function (x) x,
               logis=plogis, logistic=plogis, revgumbel=prevgumbel)
  if (is.null(fp)) stop(paste("!condquant! distribution ", dist, " not known"))
  fq <- switch(dist, normal=qnorm, gaussian=qnorm, unif=function (x) x,
               logis=qlogis, logistic=qlogis, revgumbel=qrevgumbel)
##  if (NCOL(x)>=2) stop("!condquant! x must have 2 columns")
  lx <- t(apply(rbind(x),1,sort))
  lp <- fp(lx/sig)
  lpp <- rbind(rbind(lp)%*%rbind(c(0.5,0.75,0.25),c(0.5,0.25,0.75)))
  lprand <- lp[,1]+(lp[,2]-lp[,1])*
      runif(nrow(lp),(1-randomrange)/2,(1+randomrange)/2)
  ##  <<<<<<< .mine
  lr <- cbind(cbind(median=fq(lpp[,1]),lowq=fq(lpp[,2]),
              uppq=fq(lpp[,3]), random=fq(lprand))*sig,
              prob=lp[,2]-lp[,1])
  if (any(lp0 <- lp[,2]<=0)) lr[lp0,1:4] <- matrix(lx[lp0,2],sum(lp0),4)
  if (any(lp1 <- lp[,1]>=1)) lr[lp1,1:4] <- matrix(lx[lp1,1],sum(lp1),4)
##  dimnames(lr)[[1]] <- row.names(x)
  class(lr) <- c("condquant", "matrix")
  lr
##- =======  !!!???!!!
##-   structure(cbind(sig*cbind(median= fq(lpp[,1]), lowq  = fq(lpp[,2]),
##- 			    uppq  = fq(lpp[,3]), random= fq(lprand)),
##- 		  prob=abs(lp[,2]-lp[,1])),
##- 	    class = c("condquant", "matrix")) # "matrix", e.g., for head()
##  >>>>>>> .r32
}
## ===================================================================
residuals.regr <- function(object, type=NULL, ...)
{
##-   if (!is.na(pmatch(type,"condquant"))) {
##-     ## this seems to apply only if residuals.regr is called explicitly
  lcall <- match.call()
  lcall$type <- if (is.null(type)||is.na(type)) NULL else type
  lff <- object$fitfun
  if (lff=="glm" && !is.null(type) && substr(type,1,4)=="cond") lff <- "polr"
  lcall[[1]] <-
    switch(as.character(lff),
           "polr" = quote(regr0:::residuals.polr),
           "survreg" = quote(regr0:::residuals.regrsurv),
           "coxph" = quote(regr0:::residuals.regrcoxph),
           quote(residuals))
##-   lres <-
##-       if (object$fitfun%in%c("polr","glm"))
##-         residuals.polr(object, type=type, na.action=na.action, ...)  else {
##-           if (object$fitfun=="coxph")
##-             residuals.regrcoxph(object, type=type, na.action=na.action, ...) else
##-           residuals.regrsurv(object, type=type, na.action=na.action, ...)
##-         }
##- ##-     return(lres)
##- ##-   }
##-   if (is.null(type)) return( naresid(na.action$na.action, object$residuals) )
  class(object) <- setdiff(class(object), "regr")
  lcall$object <- object
  structure( eval(lcall, envir=parent.frame()), type=type)
}
## ===================================================================
residuals.regrcoxph <- function(object, type=NULL, na.action=object, ...)
{
  ## Purpose:    conditional quantiles and random numbers for censored obs
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: Aug 2010
  if(u.nuna(type)) type <- "CoxSnellMod"
  if (type!="CoxSnellMod") {
    object$residuals <- object$resid.orig
    return(structure( survival:::residuals.coxph(object, type=type, ...),
                     type=type) )
  }
  lres <- object$residuals
  if (inherits(lres, "condquant")) return(lres)
  ly <- object$y
  lst <- ly[,2]  # status
  li <- lst!=1
  ##  martingale --> coxsnell
  lres <- lst - lres
  lrs <- qnorm(exp(-lres)) # --> normalscores
  ## fill matrix with values adequate for non-censored obs
  lrr <- matrix(lrs,length(lrs),5)
  dimnames(lrr) <- list(row.names(ly),c("median","lowq","uppq","random","prob"))
  lrr[,"prob"] <- 0
  ## censoring
  if (any(li)) {
    llim <- cbind(lrs[li],Inf)
    lr <- condquant(llim, "normal")
    lrr[li,] <- lr
  }
  lrr <- cbind(lrr, fit=object$linear.predictor)
##  class(lrr) <- "condquant"
  structure( naresid(na.action$na.action, lrr),
            class=c("condquant", "matrix"), type=type)
}
## ===================================================================
residuals.regrsurv <- function(object, type=NULL, na.action=object, ...)
{
  ## Purpose:    conditional quantiles and random numbers for censored obs
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 24 Oct 2008, 10:16
  if(u.nuna(type)) type <- "condquant"
  if (type!="condquant") 
    return(structure( survival:::residuals.survreg(object, type, ...),
                     type=type) )
  ly <- object$y  ## log for weibull
  lsig <- object$sigma
  if (length(lsig)==0) lsig <- summary(object)$scale
  if (length(lsig)==0) {
    warning("!residuals.regrsurv! no sigma found. Setting =1")
    lsig <- 1
  }
  lfit <- object$linear.predictors
##  lres <- ly[,1]-lfit
  lres <- ly[,1]-lfit
  ldist <- if (length(object$dist)>0) object$dist  else  "normal"
  li <- match(ldist, c("weibull","lognormal","loglogistic"))
  if (!is.na(li)) ldist <- c("revgumbel","normal","logistic")[li]
  ## for user-defined survreg.distributions with transformation,
  ##   this is not enough.
  ## fill matrix with values adequate for non-censored obs
  lrr <- matrix(lres,length(lres),5)
  dimnames(lrr) <- list(row.names(ly),c("median","lowq","uppq","random","prob"))
  lrr[,"prob"] <- 0
  ## censoring
  lst <- ly[,2]  # status
##   ltt <- attr(object$response[[1]], "type")
  ltt <- attr(object$response, "type")
##-   if (length(ltt)>0&&ltt=="left") lst[lst==0] <- 2
  ltl <- length(ltt)>0&&ltt=="left"
  li <- lst!=1
  if (any(li)) {
##-     llim <- cbind(lres[li],c(Inf,NA,-Inf)[lst[li]+1]) #
      llim <- if(ltl) cbind(-Inf,lres[li]) else cbind(lres[li],Inf)
      lr <- condquant(llim, ldist, lsig)
      lrr[li,] <- lr
  }
  lrr <- cbind(lrr, fit=lfit)
  ##  class(lrr) <- "condquant"
  structure( naresid(na.action$na.action, lrr),
            class=c("condquant", "matrix"), type=type)
}
## ==============================================================
nobs.survreg <- function(object, use.fallback = TRUE) {
  lnobs <- length(object$linear.predictors)
  if (lnobs==0) lnobs <- NROW(residuals(object))
  lnobs
}
nobs.coxph <- function(object, use.fallback = TRUE) {
  object$n
}
## ===================================================================
residuals.polr <- function(object, na.action=object, ...)
{
  ## Purpose:   residuals for cumulative logit regression
  ## ----------------------------------------------------------------------
  ## Arguments:
  ##   object   result of polr
  ## Value:     list with components
  ##   median     "conditional median" residual = median of conditional distr.
  ##            of latent response variable, given observed response
  ##            will probably be replaced by mean in the near future
  ##            in order to allow for adequate smoothing
  ##   lowq, uppq  lower ans upper quartiles of this cond. distribution
  ## Remark:    experimental function !!!
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  2 Oct 2007, 11:31
  if (!((lpolr <- inherits(object, "polr"))|
        (lbin <- inherits(object, "glm")&&object$family=="binomial")))
    stop ("!residuals.polr! unsuitable first argument")
  lyr <- as.character(formula(object)[2])
  ldi <- max(match(c("model","allvars"), names(object), nomatch=0))
  if (ldi) ldt <- object[[ldi]]  else {
      ldt <- if (u.debug())  model.frame(object) else
      try(model.frame(object),silent=TRUE)
      if (class(ldt)=="try-error") 
          stop ("!residuals.polr! no data found")
  }          
  ly <- ldt[,lyr]
  if (length(ly)==0) stop ("!residuals.polr! bug: no response values found")
  if (length(dim(ly))) {
    warning(":residuals.polr: returning simple deviance residuals for non-binary (grouped) data")
    return(residuals(object, type="deviance"))
  }
  #if (lpolr)
  ly <- as.numeric(ordered(ly))
  lfit <- fitted(object, type="link", na.action=list(na.action=NULL))
  lthres <- c(-100, if (lpolr) object$zeta else 0, 100)
  llim <- cbind(lthres[ly],lthres[ly+1])-lfit
  lr <- cbind(condquant(llim,"logis"),fit=lfit,y=ly)
##-   lp <- cbind(plogis(),plogis(lthres[ly+1]-lfit))
##-   lpp <- lp%*%rbind(c(0.5,0.25,0.75),c(0.5,0.75,0.25))
##-   lprand <- lp[,1]+(lp[,2]-lp[,1])*runif(length(lfit))
##-   lr <- cbind(median=qlogis(lpp[,1]),lowq=qlogis(lpp[,2]),
##-               uppq=qlogis(lpp[,3]), random=qlogis(lprand), fit=lfit,y=ly)
  dimnames(lr)[[1]] <- row.names(ldt)
##-   class(lr) <- c("condquant", "matrix") 
  structure( naresid(na.action$na.action, lr), class=c("condquant", "matrix"))
}
## ===========================================================================
linear.predictors <- function(object) {
  llp <- object$linear.predictors
  if (is.null(llp)) llp <- object$fitted.values
  if (is.null(llp))
    stop("linear.predictors! no component linear predictor")
  naresid(object$na.action, llp)
}
linpred <- linear.predictors
## ===========================================================================
fitcomp <- function(object, data=NULL, vars=NULL, se=FALSE, 
                    xm=NULL, xfromdata=FALSE, noexpand=NULL, nxcomp=51)
{
  ## Purpose:    components of a fit
  ## ----------------------------------------------------------------------
  ## !!! make nxcomp >= maximal number of factor levels !!!
  ## !!! why vars??? can possibly be much simpler!!!
  ## ----------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  6 Jan 2004, 22:24
##-   lf.predict <-
##-     function(object, newdata, type="linear.predictor", se.fit = FALSE) {
##-     if (type=="linear.predictor") {
##-       lmod <- model.matrix(formula(object)[-2], newdata)
##-       if (inherits(object,"polr")) lmod <- lmod[,colnames(lmod)!="(Intercept)"]
##-       rr <- lmod%*%object$coefficients[colnames(lmod)]
##-     } else rr <- predict(object, newdata, type=type, se.fit=se.fit)
##-     rr
##-   }
  if (length(data)==0) {
    data <- object$allvars
    if (length(data)==0)
      data <- eval(object$call$data)
  }
  lnls <- c(eval(object$call$nonlinear),FALSE)[1]
  lform <- formula(object)[-2]
  lvars <- all.vars(lform)
  if (lnls)
    lvars <- lvars[match(lvars,names(coefficients(object)),nomatch=0)==0]
  lvmiss <- setdiff(lvars,names(data))
  if (length(lvmiss)) 
    stop(paste("!fitcomp! variable(s)", paste(lvmiss,collapse=", "),
               "not in data"))
  if (length(lvars)==0)
    stop("!fitcomp! no variables selected")
  ldata <- data[,lvars,drop=FALSE]
  for (lj in 1:length(ldata))
    if (is.character(ldata[[lj]])|is.factor(ldata[[lj]]))
      ldata[,lj] <- factor(ldata[,lj])
  lformfac <- NULL
  if (!lnls) {
    if (length(c(grep("factor *\\(", format(lform)),
        grep("ordered *\\(", format(lform))))) {
      warning(
        ":fitcomp: Using 'factor(...)' or 'ordered(...)' in the formula ",
        "is hazardous for fitcomp.\n",
        "  : I try to continue.")
      lformfac <- i.findformfac(lform)
##          return(structure(list(comp=NULL), class="try-error") )
    }
  }
  ## generate means  xm  if needed
  if (length(xm)>0) {
    if ((!is.data.frame(xm))||any(names(xm)!=names(ldata))) {
      warning(":fitcomp: arg. xm  not suitable -> not used")
      xm <- NULL } else xm <- xm[1,]
  }
## median point and prediction for it
  if (length(xm)==0) {
    xm <- ldata[1,,drop=FALSE]
    for (lj in 1:length(ldata)) {
      lv <- ldata[,lj]
      if (is.character(lv)) lv <- factor(lv)
      lnhalf <- ceiling(sum(!is.na(lv))/2)
      xm[1,lj] <-
        if (is.factor(lv)) {
          levels(lv)[
          if (is.ordered(lv)) sort(as.numeric(lv))[lnhalf] else
                     which.max(table(as.numeric(lv)))
                   ]
        } else   ## median(as.numeric(lv),na.rm=TRUE)
        sort(lv)[lnhalf]
        ## median should be attained in >=1 cases
    }
  }
  if (is.null(attr(terms(object), "predvars"))) { # from  model.frame
    lterms <- attr(lm(formula(object), data=data, method="model.frame"),"terms")
    attr(object$terms,"predvars") <- attr(lterms,"predvars")
  }
  ltype <- if (object$fitfun=="coxph") "lp"  else NULL
  lprm <- c(predict(object, newdata=xm, type=ltype)) # lf.
  lny <- length(lprm)
##  expand to matrix
  if (xfromdata) {
    lx <- ldata
  } else {
    lnxj <- sapply(ldata,
                   function(x) if (is.factor(x)) length(levels(x)) else 0)
    if(!is.null(noexpand) && is.numeric(noexpand))
      noexpand <- names(noexpand)[noexpand>0]
    noexpand <- c(noexpand, lformfac) ## 
    lvconv <- names(ldata) %in% noexpand
    names(lvconv) <- names(ldata)
    if (any(lvconv)) lnxj[lvconv] <-
      sapply(ldata[lvconv], function(x) length(unique(x)) )
    lnxc <- max(nxcomp, lnxj)
    lx <- ldata[1,,drop=FALSE][1:lnxc,,drop=FALSE]
    row.names(lx) <- 1:lnxc
  }
  ##  lxm: data.frame of suitable dimension filled with "median"
  lxm <- lx
  for (lv in names(lxm)) lxm[,lv] <- xm[,lv]
  ##
  lvcomp <- names(ldata)
  if (!is.null(vars)) lvcomp <- intersect(lvcomp, vars)
  if (is.null(lvcomp)) {
    warning(":fitcomp: no variables found. Result is NULL")
    return(NULL)
  }
##  components
  lcomp <- array(dim=c(nrow(lx), length(lvcomp), lny)) 
  dimnames(lcomp) <- list(dimnames(lx)[[1]], lvcomp, names(lprm))
  lcse <- if (se) lcomp  else NULL
  for (lv in lvcomp) {
    if (xfromdata) {
      ld <- lxm
      ld[,lv] <- ldata[,lv]
      lfc <- sapply(ld,is.factor) # eliminate extra levels of factors
      if (any(lfc)) ld[lfc] <- lapply(ld[lfc], factor)
    } else { # +++
      ldv <- ldata[,lv]
      if (lnxj[lv]) { # factor levels
        ldx <- if (lvconv[lv]) sort(unique(ldv)) else factor(levels(ldv))
        lnl <- length(ldx)
        ld <- lxm[1:lnl,,drop=FALSE]
        ld[,lv] <- ldx
##-         lx[,lv] <- factor(c(1:lnl,rep(NA,lnxc-lnl)),labels=levels(ldv))
        lx[,lv] <- c(ldx,rep(NA,lnxc-lnl))
        ##
        lpr <- try( predict(object, newdata=ld, se.fit = se),
                   silent=TRUE)
        if (class(lpr)=="try-error") {
          warning(":fitcomp: no fitcomp for variable  ", lv)
          ## predict finds new levels of formfac variables
          next
        }
        if (se) {
          lc <- lpr$fit
          lcse[1:lnl,lv,] <- lpr$se.fit
        } else lc <- lpr
        lcomp[1:lnl,lv,] <- lc
        next # end for loop
      } else { # continuous var
        ld <- lxm
        lx[,lv] <- ld[,lv] <-
          seq(min(ldv,na.rm=TRUE),max(ldv,na.rm=TRUE),length=lnxc)
      } # ---
    } # +++
    ## continuous variable or xfromdata
    lpr <- predict(object, newdata=ld, se = se) # lf.
    if (se) {
      lcomp[,lv,] <- lpr$fit
      lcse[,lv,] <- lpr$se.fit
    } else lcomp[,lv,] <- lpr
  }
  if (lny==1) {
    dim(lcomp) <- dim(lcomp)[1:2]
    dimnames(lcomp) <- dimnames(lx)
    lcomp <- lcomp-lprm
    if (se) {
      dim(lcse) <- dim(lcse)[1:2]
      dimnames(lcse) <- dimnames(lx)
    }
  } else   lcomp <- sweep(lcomp,3,lprm)
  list(comp=lcomp, x=lx[,lvars,drop=FALSE], xm=xm[,lvars,drop=FALSE], se=lcse)
}
## ==========================================================================
i.findformfac <- function(formula) {
  ## find variable involved in explicit factor terms in formula
  lfo <- format(formula)
  lmf <- c(gregexpr("(factor *\\([^)]*\\))", lfo),
           gregexpr("(ordered *\\([^)]*\\))", lfo) )
  lf <- function(x) 
    if(x[1]!=-1) substring(lfo, x, x+attr(x,"match.length"))
  all.vars(as.formula(
        paste("~",paste(unlist(lapply(lmf, lf)), collapse="+"))))
}
## ==========================================================================
predict.mlm <-
  function (object, newdata=NULL, se.fit = FALSE, scale = NULL, df = Inf,
    interval = c("none", "confidence", "prediction"), level = 0.95,
    type = c("response", "terms"), terms = NULL, na.action = na.pass,
##-     pred.var = res.var/wgts, wgts = 1, ...)
    pred.var = NULL, weights = 1, ...) ## ... to absorb unused args
  ## predict.lm, extended for mlm
{
    tt <- terms(object)
    if (missing(newdata) || is.null(newdata)) {
        mm <- X <- model.matrix(object)
        mmDone <- TRUE
        offset <- object$offset
    }
    else {
        Terms <- delete.response(tt)
        m <- model.frame(Terms, newdata, na.action = na.action,
            xlev = object$xlevels)
        if (!is.null(cl <- attr(Terms, "dataClasses")))
            .checkMFClasses(cl, m)
        X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
        offset <- if (!is.null(off.num <- attr(tt, "offset")))
            eval(attr(tt, "variables")[[off.num + 1]], newdata)
        else if (!is.null(object$offset))
            eval(object$call$offset, newdata)
        mmDone <- FALSE
    }
    r <- cbind(object$residuals)
    n <- nrow(r)
    m <- ncol(r)
    ynm <- colnames(r)
    p <- object$rank
    p1 <- seq_len(p)
    piv <- object$qr$pivot[p1]
    if (p < ncol(X) && !(missing(newdata) || is.null(newdata)))
        warning("prediction from a rank-deficient fit may be misleading")
##-     beta <- object$coefficients
    beta <- cbind(object$coefficients)
##-     predictor <- drop(X[, piv, drop = FALSE] %*% beta[piv])
    predictor <- drop(X[, piv, drop = FALSE] %*% beta[piv,,drop=FALSE])
    if (!is.null(offset))
        predictor <- predictor + offset
    interval <- match.arg(interval)
    if (interval == "prediction") {
        if (missing(newdata))
            warning("Predictions on current data refer to _future_ responses\n")
        if (missing(newdata) && missing(weights)) {
            w <- weights(object) ## .default
            if (!is.null(w)) {
                weights <- w
                warning("Assuming prediction variance inversely proportional to weights used for fitting\n")
            }
        }
        if (!missing(newdata) && missing(weights) && !is.null(object$weights) &&
##-             missing(pred.var)
            is.null(pred.var)
            )
            warning("Assuming constant prediction variance even though model fit is weighted\n")
        if (inherits(weights, "formula")) {
            if (length(weights) != 2L)
                stop("`weights` as formula should be one-sided")
            d <- if (missing(newdata) || is.null(newdata))
                model.frame(object)
            else newdata
            weights <- eval(weights[[2L]], d, environment(weights))
        }
    }
    type <- match.arg(type)
    if (se.fit || interval != "none") {
        res.var <- if (is.null(scale)) {
##            r <- object$residuals
            w <- object$weights
##-             rss <- sum(if (is.null(w)) r^2 else r^2 * w)
            rss <- apply( if (is.null(w)) r^2 else r^2 * w ,2,sum)
            df <- n - p
            rss/df
        }  else {
##-         scale^2
        if (length(scale)==m) scale^2 else
           stop("!predict.lm! argument scale has wrong length")
        }
        res.var.mx <- matrix(res.var, p, m, byrow=TRUE)
        if (type != "terms") {
            if (p > 0) {
                XRinv <- if (missing(newdata) && is.null(w))
                  qr.Q(object$qr)[, p1, drop = FALSE]
                else X[, piv] %*% qr.solve(qr.R(object$qr)[p1, p1])
##-                 ip <- drop(XRinv^2 %*% rep(res.var, p))
                ip <- drop(XRinv^2 %*% res.var.mx)
            }
            else ip <- rep(0, n)
        }
    }
    if (type == "terms") {
        if (!mmDone) {
            mm <- model.matrix(object)
            mmDone <- TRUE
        }
        aa <- attr(mm, "assign")
        ll <- attr(tt, "term.labels")
        hasintercept <- attr(tt, "intercept") > 0L
        if (hasintercept)
            ll <- c("(Intercept)", ll)
        aaa <- factor(aa, labels = ll)
        asgn <- split(order(aa), aaa)
        if (hasintercept) {
            asgn$"(Intercept)" <- NULL
            if (!mmDone) {
                mm <- model.matrix(object)
                mmDone <- TRUE
            }
            avx <- colMeans(mm)
            termsconst <- sum(avx[piv] * beta[piv])
        }
        nterms <- length(asgn)
        if (nterms > 0) {
##-             predictor <- matrix(ncol = nterms, nrow = NROW(X))
            predictor <- array(dim=c(NROW(X),nterms,m))
##-             dimnames(predictor) <- list(rownames(X), names(asgn))
            dimnames(predictor) <- list(rownames(X), names(asgn), ynm)
            if (se.fit || interval != "none") {
##-                 ip <- matrix(ncol = nterms, nrow = NROW(X))
##-                 dimnames(ip) <- list(rownames(X), names(asgn))
                ip <- predictor
                Rinv <- qr.solve(qr.R(object$qr)[p1, p1])
            }
            if (hasintercept)
                X <- sweep(X, 2L, avx, check.margin = FALSE)
            unpiv <- rep.int(0L, NCOL(X))
            unpiv[piv] <- p1
            for (i in seq.int(1L, nterms, length.out = nterms)) {
                iipiv <- asgn[[i]]
                ii <- unpiv[iipiv]
                iipiv[ii == 0L] <- 0L
##-                 predictor[, i] <- if (any(iipiv > 0L))
##-                   X[, iipiv, drop = FALSE] %*% beta[iipiv]
                predictor[, i,] <- if (any(iipiv > 0L))
                  X[, iipiv, drop = FALSE] %*% beta[iipiv,]
                else 0
                if (se.fit || interval != "none")
##-                   ip[, i] <- if (any(iipiv > 0L))
##-                     as.matrix(X[, iipiv, drop = FALSE] %*% Rinv[ii,, drop = FALSE])^2 %*%
##-                       rep.int(res.var,p)
                  ip[, i,] <- if (any(iipiv > 0L))
                    as.matrix(X[, iipiv, drop = FALSE] %*% Rinv[ii,, drop = FALSE])^2 %*%
                      res.var.mx
                  else 0
            }
            if (!is.null(terms)) {
                predictor <- predictor[, terms, drop = FALSE]
                if (se.fit)
                  ip <- ip[, terms, drop = FALSE]
            }
        }
        else {
            predictor <- ip <- matrix(0, n, 0)
        }
        attr(predictor, "constant") <- if (hasintercept)
            termsconst
        else 0
    }
    if (interval != "none") {
        tfrac <- qt((1 - level)/2, df)
        if (is.null(pred.var)) pred.var <- res.var.mx/weights ## !!!
        hwid <- tfrac * switch(interval, confidence = sqrt(ip),
            prediction = sqrt(ip + pred.var))
        if (type != "terms") {
            if (m==1) {  ## changed
              predictor <- cbind(predictor, predictor + hwid %o% c(1, -1))
              colnames(predictor) <- c("fit", "lwr", "upr")
            } else {  ## changed
              predictor <- array(c(predictor, predictor - hwid, predictor + hwid),
                               dim=c(n,m,3))
              dimnames(predictor)[[3]] <- c("fit", "lwr", "upr")
            }
        }
        else {
            lwr <- predictor + hwid
            upr <- predictor - hwid
        }
    }
    if (se.fit || interval != "none")
        se <- sqrt(ip)
##-     if (missing(newdata) && !is.null(na.act <- object$na.action)) { ## !!! not yet extended
##-         predictor <- napredict(na.act, predictor)
##-         if (se.fit)
##-             se <- napredict(na.act, se)
##-     }
    if (m==1) { ## !!!
      if (length(dim(predictor))==3) predictor <- predictor[,,1]
      if (se.fit) if (length(dim(se))==3) se <- se[,,1]
    }
    if (type == "terms" && interval != "none") {
##-         if (missing(newdata) && !is.null(na.act)) { # !!! not yet extended
##-             lwr <- napredict(na.act, lwr)
##-             upr <- napredict(na.act, upr)
##-         }
        list(fit = predictor, se.fit = se, lwr = lwr, upr = upr,
            df = df, residual.scale = sqrt(res.var))
    }
    else if (se.fit)
        list(fit = predictor, se.fit = se, df = df, residual.scale = sqrt(res.var))
    else predictor
}
## ==========================================================================
compareTerms <-
  function(..., list=NULL, seq=NULL)
{
  ## Purpose:   compare terms of several models
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   models   character vector of names of model objects
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  2010
  lcl <- match.call()[-1]
  if (!is.null(names(lcl))) lcl <- lcl[names(lcl)!="seq"]
  ## !!! bug: list does not work !!!
##-   if (length(list)) lcl <- c(as.list(lcl[names(lcl)!="list"]),
##-                              c(lcl["list"]))
  lnmod <- length(lcl)
  lmnm <- names(lcl)
  if (is.null(lmnm)) lmnm <- as.character(lcl) else
    lmnm[lmnm==""] <- as.character(lcl)[lmnm==""]
  lterms <- list()
  for (lmd in 1:lnmod) {
    lmd <- eval(lcl[[lmd]])
    if (is.character(lmd)) lmd <- get(lmd)
    ltr <- if(is.list(lmd)) attr(terms(lmd),"term.labels") else NULL
    if (is.null(ltr)) stop("!compareTerms! inadequate argument", lmnm[lmd])
    lterms <- c(lterms, list(ltr))
  }
  ltrm <- unique(unlist(lterms))
  rr <- sapply(lterms, function(x) ltrm%in%x )
  dimnames(rr) <- list(ltrm,lmnm)
  if (!is.null(seq)) rr <- rr[order(match(ltrm,seq,nomatch=length(seq)+1)),]
  rr
}
## ==========================================================================
modelTable <-function(models, seq=NULL) 
{
  ## Purpose:   collect several models into a table
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   models   character vector of names of model objects
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  9 Aug 2000, 11:50
  t.islist <- is.list(models)
  if ((!t.islist)&&!(is.atomic(models)&is.character(models)))
    stop("!modelTable! Inadequate argument 'model'")
  t.nm <- length(models)
  t.mnm <- if (t.islist) names(models) else models
  if (length(t.mnm)==0)  names(models) <- t.mnm <- paste("model",1:t.nm,sep=".")
  t.ls <- t.trm <- t.cf <- setNames(vector("list",t.nm), t.mnm)
  t.nobs <- t.df <- t.sig <- t.fitfun <- setNames(rep(NA,t.nm), t.mnm)
  ## -----------------------------------------------------------
  for (li in t.mnm) {
    lr <- if (t.islist) models[[li]] else get(li,envir=parent.frame())
##-     lfitfun <- NULL
##-     for (lc in c("lm","lmrob","glm","multinom","polr","survreg"))
##-       if (inherits(lr,lc)) lfitfun <- lc
##-     if (is.null(lfitfun))
##-       stop(paste("!modelTable! Model ",li," is not an adequate model"))
    ##-     t.fitfun[li] <- lfitfun
    if (!inherits(lr, "regr"))
      stop("!modelTable! ... only programmed for 'regr' objects")
    t.fitfun[li] <- t.ff <- lr$fitfun
    t.nobs[li] <- lnr <-
      NROW(if(t.ff=="survreg") lr$linear.predictors else lr$fitted.values)
    t.df[li] <- ldf <- lnr-df.residual(lr)
    lt <- terms(lr)
    ltnm <- c( if(attr(lt,"intercept")) "(Intercept)", attr(lt, "term.labels"))
    t.cf[[li]] <-
      lr$termtable[match(ltnm,row.names(lr$termtable),nomatch=0),
                  c("coef","p.value")]
    t.trm[[li]] <- ltnm
##    t.trmc <- c(t.trmc, attr(terms(lr),"dataClasses"))
    lsig <- if (t.ff=="survreg") lr$scale else summary(lr)$sigma
    t.sig[li] <- c(lsig,NA)[1]
  }
  if (length(unique(t.nobs))>1)
    warning(":modelTable: models have different numbers of observations")
## --- collect
  t.tr <- unique(unlist(t.trm))
  ## --- coefs and p values
  t.nt <- length(t.tr)
  t.pr <- t.coef <- matrix(NA,t.nt,t.nm, dimnames=list(t.tr,t.mnm))
  ## ---
  for (li in t.mnm) {
    t.t <- t.trm[[li]]
    if (length(t.t)) {
      lcf <- t.cf[[li]]
      t.coef[t.t,li] <- lcf[,1]
      t.pr[t.t,li] <- lcf[,2]
    }
  }
  ## reorder
  if (length(seq)>0) {
    li <- match(seq, t.tr)
    t.t <- c(t.tr[notna(li)],t.tr[!t.tr%in%seq])
    t.coef <- t.coef[t.t,]
    t.pr <- t.pr[t.t,]
    t.sd <- t.sd[t.t]
  }
#  attr(t.coef,"standardized") <- t.trn
  if (all(is.na(t.sig))) t.sig <- NULL
  t.r <- list(coef=t.coef, p=t.pr, sigma=t.sig, nobs=t.nobs,
              df=t.df, fitfun=t.fitfun) # , sd.terms=t.sd
  class(t.r) <- "modelTable"
  t.r
}
## ==========================================================================
"[.modelTable" <- function(object,rows=NULL,cols=NULL, reduce=TRUE) {
  if (is.null(rows)) rows <- 1:nrow(object$coef)
  if (is.null(cols)) cols <- 1:ncol(object$coef)
  lp <- object$p[rows,cols,drop=FALSE]
  li <- if(reduce) !apply(is.na(lp),1,all)  else 1:nrow(lp)
  if (length(li)==0) stop("![.modelTable! no terms left")
  lsd <- if (length(object$sd.terms)) object$sd.terms[rows][li]
  lsig <- if (length(object$sigma)) object$sigma[rows][li]
  rr <- list(coef=object$