# t.freg <- f.reg(log10(ersch)~log10(dist)+log10(ladung)+stelle,data=d.spreng14)
## ==========================================================================
f.reg <- function(formula, data, tit=NULL, method="lm",
                 init.reg="f.ltsreg", ...)
{
  ## Purpose:    lm with adequate handling of NAs and more results
  ## -------------------------------------------------------------------------
  ## Arguments:  formula, data, ...  as with lm;
  ##             omit: observations to be omitted from the analysis
  ##             tit:  title (becomes tit attribute of result)
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date:  7 May 93, 13:28
  ## -------------------------------------------------------------------------
##rfr  on.exit(browser())		#
  ## calculations
  dataname <- deparse(substitute(data))
  lcall <- match.call()
  lform <- formula
    ldata <- data
    lnobs <- nrow(ldata)
    lrl <- row.names(ldata)
    row.names(ldata) <- as.character(1:lnobs)
#  synchronize() # in S
  lreg <- lm(lform, qr=T, na.action=na.omit, data=ldata) # , ...
  lnna <- as.character(1:lnobs) %in% names(lreg$residuals)
  if (method=="robust") { lqr <- lreg$qr
       lreg <- mmreg(lform, na.action=na.omit, data=ldata,
                   init.reg=get(t.init.reg), ...)
       lreg$qr <- lqr }
  if (method=="rlm") { lqr <- lreg$qr
       lreg <- rlm(lform, na.action=na.omit, data=ldata, ...)
       lreg$qr <- lqr }
##-   lreg$call$formula <- lform
##-     lreg$call$data <- as.name(dataname)
    lreg$call <- lcall
    ltit <- if(is.null(tit)) attr(data,"tit") else tit
    ldoc <- attr(data,"doc") 
## collecting results
  tit(lreg) <- ltit
    doc(lreg) <- ldoc
  lreg1 <- summary(lreg)
    lhat <- u.merge(hat(lreg$qr),NA,lnna, names=lrl)
  lreg$h <- lhat
  lreg$residuals <- u.merge(residuals(lreg),NA,lnna, names=lrl)
  if (method=="robust") lreg$rweights <- u.merge(lreg$w,NA,lnna, names=lrl)
  lreg$fitted.values <- u.merge(lreg$fitted.values,NA,lnna, names=lrl)
  if (!is.null(lreg$weights))
    lreg$weights <- u.merge(lreg$weights,NA,lnna, names=lrl)
  lreg$sigma <- lreg1$sigma
  lreg$r.squared <- lreg1$r.squared
  lreg$fstatistic <- lreg1$fstatistic
##-   lreg$df <- lreg1$df
  lcov <- lreg1$cov.unscaled * lreg1$sigma
  lreg$covariance <- lcov
  lse <- sqrt(diag(lcov))
  lreg$correlation <- lcov/outer(lse, lse)
  lreg$colregelation <- lreg1$colregelation
  lreg$df <- lreg1$df
##- terms table
  ltq <- qt(0.975,lreg$df.residual)
  ltb <- drop1(lreg, test="F", scope=terms(lreg))[-1,]
  lcont <- ltb$Df==1
  ltlb <- dimnames(ltb)[[1]]
  lclb <- ltlb[lcont]
  lnt <- nrow(ltb)
  ltb <- c(list(coef=rep(NA,lnt), stcoef=rep(NA,lnt),
                testty=rep("F test",lnt)), ltb)
  lcf <- coefficients(lreg)[lclb]
  lsd <- sqrt(apply(lreg$model[,c(T,lcont)],2,var))
  lstcf <- lcf*lsd[1]/lsd[-1]
  ltb$coef[lcont] <- lcf
  ltb$stcoef[lcont] <- lstcf
  ltb$testty[lcont] <- "t-ratio"
  ltb[["F value"]][lcont] <- sign(lcf)*sqrt(ltb[["F value"]][lcont])/ltq
  ltb <- data.frame(ltb[c(1:4,8:9)])
  names(ltb) <- c("coef","stcoef","testty","df","testst","p-value")
  lint <- coefficients(lreg)["(Intercept)"]
  row.names(ltb) <- ltlb
  if (length(lint)) ltb <- rbind("(Intercept)"=c(lint,NA,NA,1,NA,NA),ltb)
  lreg$testcoef <- ltb
  lreg$allcoef <- dummy.coef(lreg)
  lreg$stres <- lreg$res/(lreg1$sigma*sqrt(1- ifelse(lhat<1,lhat,0)))
  class(lreg) <- c("freg",class(lreg1),class(lreg))
##rfr  if (!Browse) on.exit()
## result
  lreg
}

## ==========================================================================
print.freg <- function (x, correlation = FALSE,
    digits = max(3, getOption("digits") - 3), 
    symbolic.cor = p > 4, signif.stars = getOption("show.signif.stars"),
    residuals=FALSE, ...) 
{
##
##rfr    on.exit(browser())
    cat("\nCall:\n")
    cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), 
        "\n", sep = "")
    df <- x$df
    rdf <- df[2]
    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")
      }
    }
    if (nsingular <- df[3] - df[1]) 
        cat("\nCoefficients: (", nsingular,
            " not defined because of singularities)\n", sep = "")
    else {
      cat("\nTerms:\n")
      ltc <- x$testcoef
      ltc[,"p-value"] <- round(ltc[,"p-value"],max(3,digits))
      print(ltc)
      mterms <- x$testcoef[,"df"]>1
      if (any(mterms)) {
      cat("\nCoefficients of blocks:\n")
      print(x$allcoef[mterms]) }
    }
    cat("\nStandard deviation of error", format(signif(x$sigma, 
        digits)), "on", rdf, "degrees of freedom\n")
    if (!is.null(x$fstatistic)) {
        cat("Multiple R-Squared:", formatC(x$r.squared, digits = digits))
        cat(",\tAdjusted R-squared:", formatC(x$adj.r.squared, 
            d = digits), "\nF-statistic:", formatC(x$fstatistic[1], 
            digits = digits), "on", x$fstatistic[2], "and", x$fstatistic[3], 
            "df,  p-value:", formatC(1 - pf(x$fstatistic[1], 
                x$fstatistic[2], x$fstatistic[3]), dig = digits), 
            "\n")
    }
    correl <- x$correlation
    if (correlation && !is.null(correl)) {
        p <- NCOL(correl)
        if (p > 1) {
            cat("\nCorrelation of Coefficients:\n")
            if (symbolic.cor) 
                print(symnum(correl)[-1, -p])
            else {
                correl[!lower.tri(correl)] <- NA
                print(correl[-1, -p, drop = FALSE], digits = digits, 
                  na = "")
            }
        }
    }
    cat("\n")
##rfr    if(!Browse) on.exit()
    invisible(x)
}
## ==========================================================================
summary.freg <- function(object)  object
## ==========================================================================
plot.freg <-
function(object, nmark=NULL, lab=NULL, labcex=0.7, ask = NULL,
         mfrow = c(2,2), oma = 2*(prod(mfrow)>1), lty = c(1,3,2,4,5),
         colors = c(1,5,2,3,4),
         main =NULL, cex.title=1, lowess.it=NA, 
	 qq = TRUE, yfit = F, yfit.lowess=T, 
	 ta = !yfit, ta.lowess = TRUE, show.2sigma=TRUE,
         tascale = TRUE, tascale.lowess = TRUE, 
	 hat = TRUE, xplot = TRUE, x.se=TRUE, x.lowess = TRUE, x.ylim = NULL,
         ...)
{
## Purpose:  more plots for residual analysis
## -------------------------------------------------------------------------
## Arguments:
##  object     result of lm or glm (or nls?)
##  nmark      how many extreme positive and negative residuals should be
##             labeled?
##  lab        labels for points
##  labcex     character expansion for labels
##  ask        graphics pararmeter: should R ask before erasing the screen?
##  mfrow, oma arguments passed to par
##  lty, colors  to be used for different elements of the plots:
##             [1] for points,
##             [2] for reference lines
##             [3] for smooths
##             [4] for terms in g.termsplot
##             [5] for confidence bands of terms
##  main       main title, defaults to the formula of  object
##  cex.title  character expansion for the main title
##  qq         if TRUE, show normal plot of standardized residuals 
##  yfit       if TRUE, show response against fitted
##  ta         if TRUE, show Tukey-Anscombe plot (residuals against fit)
##  ta.lowess  if TRUE, show lowess smooth in TA plot
##  show.2sigma  if TRUE, show 2 sigma limit in TA plot
##  tascale    if TRUE, show absolute residuals against fit
##  tascale.lowess  if TRUE, show lowess smooth in tascale plot
##  hat        if TRUE, show residuals against diagonal of hat matrix
##  xplot      if TRUE, show residuals against terms in the model
##             if a formula or a character vector, it contains the terms
##             against which the residuals are plotted
##             may contain variables that are not in the model (in future!!!)
##  x.lowess   if TRUE, show lowess smooth in term plots
##  x.ylim     y range(s) for terms plots
##  x.lowess   if TRUE, show lowess smooth in x and xadd plots
##  !!! the case with unequal weights is not yet properly treated.
## -------------------------------------------------------------------------
## Author: Werner Stahel, Date:  7 May 93 / 2002
  is.formula <- function(object)
    length(class(object))>0 && class(object)=="formula"
##rfr  on.exit(browser())	#
## prepare graphical elements
  if (is.null(ask)) ask <- !last(c("",names(dev.list)))=="postscript"
  loma <- if (length(oma)==4) oma else c(0,0,oma[1],0)
  op <- par(ask = ask, mfrow = mfrow, oma = loma)
  lform <- formula(object)
  lftext <- deparse(lform)
  if (is.null(main))  main <- if(exists("c.env")) c.env$project else ""
  if (!(is.logical(main)&&!main)) main <- paste(main,lftext,sep="  ")
  if (is.null(cex.title)) cex.title <- max(0.5, min(1,
      par("mfg")[4]*par("pin")[1]/(par("cin")[1]*nchar(main))))
  lf.title <- function(main, cex = cex.title, out = loma[3]>0) 
    mtext(main, 3,1-out, cex = cex, outer = out)
#       	object$fitted.values <- object$fit
  lf <- fitted(object) # predict(object)  does not work for nls
  lr <- residuals(object)
  ly <- lf + lr
  lyname <- deparse(lform[[2]])
  lnna <- !is.na(lr)
  ln <- length(lr)
  lsigma <- sum(lr^2)/object$df.residual
  lh <- object$h
  if(is.null(lh)) lh <- hat(model.matrix(object))
  lstres <- object$stres
  if(is.null(lstres)) lstres <- lr/(lsigma*sqrt(1-lh))
  lrname <- paste("res(", lyname, ")")
  lstrname <- paste("st", lrname, sep = ".")
  if (is.na(lowess.it)) lowess.it <- 3*(c(class(object),"")[1]!="glm")
## - labels
  lrown <- names(lr)
  if (is.null(lrown)) lrown <- as.character(1:ln)
  llab <- if (is.null(lab)) lrown else rep(as.character(lab), length = ln)
  if (is.null(nmark)||is.na(nmark))
    nmark <- if (length(lab)>1) ln else ceiling(sqrt(ln)/2)
  if (nmark==0) llab <- par("pch")
  else  if (nmark<ln/2) {
      lor <- order(lstres)
      li <- lor[(nmark+1):(ln-nmark)]
      lpch <- par("pch")
      if (is.numeric(llab)) {
        llab[li] <- if(is.numeric(lpch)) lpch else 1
      } else  llab[li] <- if (is.numeric(lpch)) "+" else lpch
  }
  ltxt <- is.character(llab)&&any(nchar(lab)>1)
  lpty <- ifelse(ltxt,"n","p")
## --------------------------------------------------------------------------
## start plots
  if(qq) { # normal plot
    qqnorm(lstres, ylab = lstrname, main="", col=colors[1])
    lquart <- quantile(lstres,c(0.25,0.75))
    abline(0, diff(lquart)/(2*qnorm(0.75)), lty = lty[2], col=colors[2])
    lf.title(main)
    g.stamp(sure=FALSE) }
## --- 
  if(yfit) { # y against fit
    plot(lf, ly, xlab = "fitted", ylab = lyname, type=lpty, pch=llab,
         col=colors[1], ...)
    if (ltxt) text(lf, ly, llab, cex=labcex, lty = lty[1], col=colors[1])
    abline(0, 1, lty = lty[2], col=colors[2])
    if(yfit.lowess)
      lines(lowess(lf[lnna], ly[lnna]), lty = lty[3], col=colors[3])
    lf.title(main)
    g.stamp(sure=FALSE)
  }
## --- 
  if(ta) { 
##-     TA.plot(object, fit=lf[lnna], res=lr[lnna], labels=lab, 
##- 	    main="", draw.smooth=ta.lowess, # ylab=rname, 
##- 	    show.call=F, show.2sigma=show.2sigma, lo.iter=lowess.it, ...) 
    plot(lf, lr, xlab = "fitted", ylab = lrname, type=lpty, pch=llab, ...)
    if (ltxt) text(lf, lr, llab, cex=labcex)
    abline(0, 0, lty = lty[2], col=colors[2])
    if(ta.lowess)
      lines(lowess(lf[lnna], lr[lnna]), lty = lty[3], col=colors[3])
    lf.title(main)
    g.stamp(sure=FALSE)
  }
## --- 
  lra <- abs(lstres)
  if(tascale) {
    plot(lf, lra, xlab = "fitted", ylab = "|standardized res|", 
         type=lpty, pch=llab, ylim=c(0,1.05*max(lra)), yaxs="i", ...)
    if (ltxt) text(lf, lra, llab, cex=labcex)
    if(tascale.lowess) {
      lras <- sqrt(lra)
      llo <- lowess(lf[lnna], lras[lnna], iter=lowess.it)
      lines(llo$x, llo$y^2, lty = lty[3], col=colors[3]) }
    lf.title(main)
    g.stamp(sure=FALSE)
  }
## --- 
  if(hat) {
    plot(lh, lr, xlab = "hat diagonal", ylab = lrname, type=lpty, pch=llab, ...)
    if (ltxt) text(lh, lr, llab, cex=labcex)
    abline(0, 0, lty = lty[2], col=colors[2])
    lf.title(main)
    g.stamp(sure=FALSE)
  }
## --- 
  if (is.logical(xplot)) xplot <- if(xplot) lform else NULL
  if(length(xplot) > 0) {
    if (is.character(xplot))
      xplot <- formula(paste("y ~",paste(xplot,collapse="+")))
    xplot <- update(lform, xplot)
    lterms <- as.character(attr(terms(xplot),"variables")[-(1:2)])
    ldataname <- as.character(object$call[3])
    ldata <- get(ldataname)
    if (nrow(ldata)!=length(lr)) {
      li <- match(names(lr),row.names(ldata))
      if (any(is.na(li))) 
        stop("cannot draw res on x's since original y contains NAs")
      else ldata <- ldata[li,]
    }
    g.termplot(object, data=ldata, terms=lterms, se=x.se, partial=TRUE,
      ylim = x.ylim, main=main, cex.title=cex.title, lty=lty, lwd.term=1.5,
      colors = colors)
  }
##rfr  if (!Browse) on.exit()
  "plot.freg done"
}
## --- plot against explanatory variables and variables in  xadd
## ==========================================================================
g.termplot <-
  function (model, data = model.frame(model), partial.resid = TRUE, 
  rug = FALSE, terms = NULL, se = FALSE, addterm = FALSE,
  smooth = lowess, colors = c(1,5,2,3,4), lty = c(1,3,2,4,5),
  xlabs = NULL, ylabs = NULL, main = NULL, cex.title = 1, ylim = NULL,
  lwd.term = 2, cex = 1, pch = par("pch"),
  ask = interactive() && nb.fig < n.tms && .Device != "postscript", 
  ...) 
{
## Purpose:   termplot, adapted
## -------------------------------------------------------------------------
## Arguments:
##   addterm
##   smooth
##   ylim
## -------------------------------------------------------------------------
## Author: Werner Stahel, Date: 12 Jan 2002, 06:44
  ## lty.term = if (addterm) 1 else 3, 
##rfr  on.exit(browser())
  terms <- if (is.null(terms)) 
      predict(model, type = "terms", se = se)
  else predict(model, type = "terms", se = se, terms = terms)
  if (is.logical(smooth)&&smooth) smooth <- lowess
  addterm <- as.logical(addterm)
  if (is.na(addterm)) addterm <- FALSE
  tmssign <- if (addterm) 1 else -1
  tms <- tmssign * as.matrix(if (se) terms$fit else terms)
  n.tms <- ncol(tms)
  mf <- model.frame(model)
  nmt <- colnames(tms)
  cn <- parse(text = nmt)
  if (is.null(ylabs)) 
      ylabs <- if (addterm) paste("Partial for", nmt) else "Residuals"
  ylabs <- rep(ylabs, length=n.tms)
  if (!is.null(ylim)) {
    if (!NROW(ylim)==2) {
      ylim <- NULL
      warning("inadequate argument  ylim") }
    else if (!NCOL(ylim)==n.tms) {
      if (NCOL(ylim)==1) ylim <- matrix(ylim,2,n.tms)
      else warning("inadequate number of columns in  ylim")
    }
  }
  df <- model$Df.residuals
  qnt <- if (is.null(df)) 2 else qt(0.975,df)
## -
  lty <- rep(c(lty,1:5),length=5)
  colors <- rep(c(colors,1:5),length=5)
## - main 
  lmain <- ""
  if (is.logical(main)&&main) lmain <- deparse(model$call)
  if (is.character(main)) lmain <- main
  lonetitle <- length(lmain)==1
  loma <- par("oma")
  lf.title <- function(main, i=1, cex=cex.title, out=lonetitle&&loma[3]>0) {
    lmain <- if (lonetitle) main else main[i]
    if (!(any(par("mfg")[1:2]>1) && lonetitle)) 
    mtext(lmain, 3, 1-out, cex = cex, outer = out) }
## -
  pf <- parent.frame()
  carrier <- function(term) {
    if (length(term) > 1) carrier(term[[2]])
    else eval(term, data, enclos = pf)
  }
  carrier.name <- function(term) {
    if (length(term) > 1) carrier.name(term[[2]])
    else as.character(term)
  }
## ---
  if (is.null(xlabs)) 
      xlabs <- unlist(lapply(cn, carrier.name))
  if (partial.resid) {
    res <- residuals(model, "response")
    if (addterm) pres <- residuals(model, if (addterm) "partial") }
  is.fac <- sapply(nmt, function(i) is.factor(mf[, i]))
## ---
  nb.fig <- prod(par("mfcol"))
  if (ask) {
      op <- par(ask = TRUE)
      on.exit(par(op))
  }
## --- loop ---
  for (i in 1:n.tms) {
## - plot range
    if (partial.resid)
      rr <- if (addterm) pres[,i] else res
    if (is.null(ylim)) {
      if (partial.resid) ylims <- range(rr, na.rm = TRUE)
      else  ylims <- range(tms[, i], na.rm = TRUE)
      if (rug)  ylims[1] <- ylims[1] - 0.07 * diff(ylims)
    }  else  ylims <- ylim[,i]
## -
    if (is.fac[i]) {
      ff <- mf[, nmt[i]]
      ll <- levels(ff)
      xlims <- range(seq(along = ll)) + c(-0.5, 0.5)
      xx <- codes(ff)
##-       if (rug) {
##-           xlims[1] <- xlims[1] - 0.07 * diff(xlims)
##-           xlims[2] <- xlims[2] + 0.03 * diff(xlims)
##-       }
      plot(1, 0, type = "n", xlab = xlabs[i], ylab = ylabs[i], 
           xlim = xlims, ylim = ylims, main = "", ...)
      lf.title (lmain, i, cex=cex.title)
      for (j in seq(along = ll)) {
        ww <- which(ff == ll[j])[c(1, 1)]
        xj <- j + c(-0.4, 0.4)
        lines(xj, tms[ww, i], col = colors[4], lwd = lwd.term, 
          lty=lty.term,  ...)
        if (se) {
          wid <- qnt * terms$se.fit[ww, i]
          lines(xj, addterm * tms[ww, i] + wid, lty = lty[5], 
          col = colors[5])
          lines(xj, addterm * tms[ww, i] - wid, lty = lty[5], 
          col = colors[5])
        }
      }
    }
    else {
      xx <- carrier(cn[[i]])
      xlims <- range(xx, na.rm = TRUE)
      if (rug) 
          xlims[1] <- xlims[1] - 0.07 * diff(xlims)
      io <- order(xx)
      plot(xx[io], tms[io, i], type = "l", xlab = xlabs[i], 
          ylab = ylabs[i], xlim = xlims, ylim = ylims, 
          main = "", col = colors[4], lwd = lwd.term, lty = lty[4],
          ...)
      lf.title(lmain)
      if (se) {
        wid <- qnt * terms$se.fit[io, i]
        lines(xx[io], addterm * tms[io, i] + wid, lty = lty[5], 
              col = colors[4])
        lines(xx[io], addterm * tms[io, i] - wid, lty = lty[5], 
              col = colors[4])
      }
      if (!is.null(smooth)) {
        ys <- smooth(xx[io], res[io])
        if (is.list(ys)) ys <- ys$y
        if (addterm) ys <- ys + tms[io, i] 
        lines(xx[io], ys, lty=lty[3], col=colors[3])
      }
    }
    if (partial.resid) 
        points(xx, rr, cex = cex, pch = pch, col = colors[1])
    if (rug) {
        n <- length(xx)
        lines(rep(jitter(xx), rep(3, n)), rep(ylims[1] + 
            c(0, 0.05, NA) * diff(ylims), n))
##-         if (partial.resid) 
##-             lines(rep(xlims[1] + c(0, 0.05, NA) * diff(xlims), 
##-               n), rep(rr, rep(3, n)))
    }
    if (!addterm) abline(h=0, lty=lty[2], col=colors[2])
    g.stamp(sure=FALSE)
  }
##rfr  if (!Browse) on.exit()
  invisible(n.tms)
}
## ===========================================================================
g.res2x <- function(formula=NULL, data, restricted=NULL, size = 1,
  slwd = 1, scol = 2, xlab = NULL, ylab= NULL, xlim=NULL, ylim=NULL, ...)
{
## Purpose:  plot residuals vs. two x's
## Author:   ARu , Date:  11/Jun/91
## Aenderungen: MMae, 30/Jan/92, Dez.94
## --------------------------------------------------------------------------
## Arguments:
##   formula    z~x+y, where
##              x, y  coordinates of points given by two vector arguments.
##              z     gives orientation (by sign)
##                   and size (by absolute value) of symbol.
##   data       data.frame or regression results.
##   restricted absolute value which truncates the size.
##              The corresponding symbols are marked by stars.
##   size       the symbols are scaled so that 'size' is the size of
##              the largest symbol in cm.
##   ...        additional arguments for the S-function 'plot'
## the function currently only plots  z  for the first two terms of the
## right hand side of  formula
## --------------------------------------------------------------------------
##rfr  on.exit(browser())
  t.form <- as.formula(formula)
  if (inherits(data,"lm")) {
    t.data <- get(as.character(data$call[3]))
    if (is.null(formula)) t.form <- formula(data)[c(1,3)]
    if (length(t.form)<3) {
      t.form <- update.formula(t.form,residuals~.)
      t.data <- f.merge1(t.data,resid(data),namefrom="residuals") } 
  } else  t.data <- data
  if (!is.data.frame(t.data)) {
    if(is.matrix(data)) t.data <- as.data.frame(data) else
      stop("data is not a data.frame") }
  t.d <- model.frame(t.form,t.data)
  t.d <- na.omit(t.d)
  z <- t.d[,1]
  x <- as.numeric(t.d[,2])
  y <- as.numeric(t.d[,3])
  if (is.null(xlim)) xlim <- range(x)
  if (is.null(ylim)) ylim <- range(y)
##--- restrict z values: ---
  if(length(restricted)==0)   restr <- F else {
    restr <- abs(z) > restricted
    z <- pmin( pmax( z, -restricted), restricted) }
##--- fix plot region: ---
  pcm <- par("pin") * 2.54              #damit in cm
##--- damit im Plot das Symbol wirklich die Groesse size hat:
  size <- size/(2 * sqrt(2))
  fx <- (size * diff(xlim))/(pcm[1] - 2 * size)/2
  fy <- (size * diff(ylim))/(pcm[2] - 2 * size)/2
##--
  if(is.null(xlab)) xlab <- names(t.d)[2]
  if(is.null(ylab)) ylab <- names(t.d)[3]
  plot(x, y, xlim = xlim + c(-1,1)* fx, ylim = ylim + c(-1,1)* fy, pch = ".",
       xlab=xlab, ylab=ylab, ...)
##---------------
##--- draw symbols: ---
  z <- z/max(abs(z), na.rm = T)
  usr <- par("usr")
  sxz <-     diff(usr[1:2])/pcm[1] * size * z
  syz <- abs(diff(usr[3:4])/pcm[2] * size * z)
  segments(x - sxz, y - syz,  x + sxz, y + syz, lwd = slwd, col=scol)
##--- mark restricted observations: ---
  if(any(restr)) {
    points((x - sxz)[restr], (y - syz)[restr], pch= 8, mkh = 1/40)
    points((x + sxz)[restr], (y + syz)[restr], pch= 8, mkh = 1/40)
  }
  g.stamp(sure=F)
##rfr  if (!Browse) on.exit()
"g.res2x done"
}
## ============================================================
u.merge <- function(dd1, dd2 = NA, which=NULL, after=NULL,
                    length=NULL, names=NULL)
{
## Purpose:   merge two vectors or expand a vector by NA s
## -------------------------------------------------------------------------
## Arguments:
##   dd1      first vector or matrix or data.frame (?),
##   dd2      second vector, ...
##   which    is T for indices for which first vector is used
##   after    elements of  dd2  will be inserted after "after" in dd1
##   length   length of the result (will be expanded if necessary)
##   names    names of the result (if length is adequate)
## -------------------------------------------------------------------------
## Author: Werner Stahel, Date: 11 Mar 93, 13:50, and later
##rfr  on.exit(browser())
  llen <- length
  n1 <- length(dd1)
  nc1 <- ncol(dd1)
  nc2 <- ncol(dd2)
  if (!is.null(nc1)) {
    n1 <- nrow(dd1)
    if (!(length(dd2)==1 | (nc2==nc1)))
      stop("unsuitable second argument")
    }
## --- generate  which  vector for all cases
  if (is.null(which)) {
## - after  specified
      if (is.null(after)) stop("specify either  which  or  after")
      if (is.logical(after))  after <- which(after)
      wh <- rep(TRUE,n1+length(after))
      wh[after+1:length(after)] <- FALSE }
  else {
## - which  specified
    if(is.logical(which)) wh <- which  
    else {
      if(is.null(llen))  llen <- n1+length(which)
        wh <- rep(TRUE, llen)
        wh[which] <- FALSE }
  }
## --- merge
  nn <- length(wh)
  n2 <- nn-n1
  if (!(is.null(names)|length(names)==nn))
    warning("argument  names  not used (unsuitable length)")
  if (!is.null(nc1)) {
    if (!(length(dd2)==1 | ((!is.null(nc2))&nrow(dd2)==n2)))
      stop("unsuitable number of rows") 
    rr <- matrix(NA,nn,nc1)
    rr[wh,] <- dd1
    rr[!wh,] <- dd2
    if (!is.null(names)) row.names(rr) <- names
  }
  else {
    rr <- rep(NA,nn)
    rr[wh] <- dd1
    rr[!wh] <- dd2
    if (!is.null(names)) names(rr) <- names
  }
  on.exit()
  rr
}
## ============================================================
f.merge1 <- function(data, from, var=NULL, by=NULL, byfrom=by, namefrom=NULL)
{
  ## Purpose:   add data from another data set
  ## -------------------------------------------------------------------------
  ## Arguments:
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 17 Jul 2000, 11:08
##rfr  on.exit(browser())
  t.v <- rep(NA,nrow(data))
  t.fromby <- if (is.null(byfrom)) {
    if (is.data.frame(from)) row.names(from) else
    if (is.matrix(from)) dimnames(from)[[1]] else names(from) } else
    from[,byfrom]
  t.databy <- if (is.null(by)) row.names(data) else data[,by]
  if (is.null(t.fromby)|is.null(t.databy))
    stop("by  not contained in both frames")
  t.i <- match(t.fromby,t.databy)
  t.fr <- as.data.frame(from)
  if (length(var)>0) t.fr <- t.fr[,var,drop=F]
  if (length(namefrom)>0) names(t.fr) <- namefrom
  if (all(length(var)!=c(0,length(t.fr)))|any(is.na(names(t.fr))))
    stop("wrong arg.  var : not all contained in data.frame  from") 
  if (any(is.na(t.i))) {
    warning(paste(sum(is.na(t.i)),"rows of  from  not in  data"))
    t.fr <- t.fr[!is.na(t.i),]
    t.i <- na.omit(t.i) }
  t.var <- names(t.fr)
  for (l.var in t.var) {
    t.v[t.i] <- t.fr[,l.var]
    data[,l.var] <- t.v }
  on.exit()
  data
}
## ==========================================================================
g.stamp <- function(sure=TRUE, outer = NULL, project=c.env$project, 
           step=c.env$step, stamp=c.env$stamp, ...)
{
## Purpose:   plot date and project information
## -------------------------------------------------------------------------
## Arguments:
##   sure     if F, the function only plots its thing if  c.env$stamp>0
##   outer    if T, the date is written in the outer margin
##   project  project title
##   step     title of step of data analysis
##   ...      arguments to  mtext , e.g., line=3
## -------------------------------------------------------------------------
## Author: Werner Stahel, Date: 13 Aug 96, 09:00
##  on.exit(browser())
  if (is.null(outer)) outer <- par('oma')[4]>0 
  t.txt <- date()
  t.txt <- paste(substring(t.txt,5,10),",",substring(t.txt,22,23),"/",
                 substring(t.txt,13,16),sep="")
    if (!is.null(project)) t.txt <- paste(t.txt,project,sep=" | ")
    if (!is.null(step)) t.txt <- paste(t.txt,step,sep=" | ")
  if( sure | stamp==2 | ( stamp==1 & (
##     last figure on page
     { t.tfg <- par("mfg") ; all(t.tfg[1:2]==t.tfg[3:4]) }
       || (is.logical(outer)&&outer) ))  ) 
       mtext(t.txt, 4, cex = 0.6, adj = 0, outer = outer, ...) 
  invisible(t.txt)
}
if (!exists("c.env")) c.env <- 
  list(project="",step="",stamp=1,out=1,mar=c(3,3,1,1)+0.1,mgp=c(2,0.8,0))
## ===========================================================================
last <- function(data,length.out = 1)
{
  ## Purpose: last element(s) of a vector
  ## Author: Werner Stahel, Date:  Tue Jan 21 17:29:42 1992
  ## ----------------------------------------------------------------
  ## Arguments:
  ## data: vector
  ## length: if positive, the result is the  length  last elements of data
  ##         if negative, the last  length  elements are dropped
  ## ----------------------------------------------------------------
  n <- length(data)
  data[sign(length.out)*(n-abs(length.out)+1):n]
}
## ===========================================================================
"doc<-" <- function(x, value)
{
  ##-- Create doc attribute or  PREpend  new doc to existing one.
  attr(x, "doc") <- c(value, attr(x, "doc"))
  x
}
"tit<-" <- function(x, value)
{
  attr(x, "tit") <- value
  x
}
## ===========================================================================