[Rd] Should there be a confint.mlm ?

[Martin Maechler]

>>>>> steven pav 
>>>>>     on Thu, 19 Jul 2018 21:51:07 -0700 writes:

    > It seems that confint.default returns an empty data.frame
    > for objects of class mlm. For example:

> It seems that confint.default returns an empty data.frame for objects of
> class mlm.

Not quite: Note that 'mlm' objects are also 'lm' objects, and so
it is confint.lm() which is called here and fails.

> For example:
> 
> ```
> nobs <- 20
> set.seed(1234)
> # some fake data
> datf <-
> data.frame(x1=rnorm(nobs),x2=runif(nobs),y1=rnorm(nobs),y2=rnorm(nobs))
> fitm <- lm(cbind(y1,y2) ~ x1 + x2,data=datf)
> confint(fitm)
> # returns:
>      2.5 % 97.5 %
> ```
> 
> I have seen proposed workarounds on stackoverflow and elsewhere, but
> suspect this should be fixed in the stats package. 

I agree.
It may be nicer to tweak  confint.lm() instead though.

I'm looking into doing that.

> A proposed implementation would be:
> 
> ```
> # I need this to run the code, but stats does not:
> format.perc <- stats:::format.perc

or better (mainly for esthetical reasons), use

  environment(confint.mlm) <- asNamespace("stats")

after defining  confint.mlm [below]

> # compute confidence intervals for mlm object.
> confint.mlm <- function (object, level = 0.95, ...) {
>   cf <- coef(object)
>   ncfs <- as.numeric(cf)
>   a <- (1 - level)/2
>   a <- c(a, 1 - a)
>   fac <- qt(a, object$df.residual)
>   pct <- format.perc(a, 3)
>   ses <- sqrt(diag(vcov(object)))
           ^^^^^^^^^^^^^^^^^^^^^^^^
BTW --- and this is a diversion ---  This is nice mathematically
    (and used in other places, also in "base R" I think)
but in principle is a waste:  Computing a full 
k x k matrix and then throwing away all but the length-k
diagonal ...
In the past I had contemplated but never RFC'ed or really
implemented a stderr() generic with default method

   stderr.default <- function(object) sqrt(diag(vcov(object)))

but allow non-default methods to be smarter and hence more efficient.

>   ci <- ncfs + ses %o% fac
>   setNames(data.frame(ci),pct)
> }
> 
> # returning to the example above,
> confint(fitm)
> # returns:
>                  2.5 % 97.5 %
> y1:(Intercept) -1.2261 0.7037
> y1:x1          -0.5100 0.2868
> y1:x2          -2.7554 0.8736
> y2:(Intercept) -0.6980 2.2182
> y2:x1          -0.6162 0.5879
> y2:x2          -3.9724 1.5114
> ```

I'm looking into a relatively small patch to confint.lm()
*instead* of the confint.mlm() above

Thank you very much, Steven, for your proposal!
I will let you (and the R-devel audience) know the outcome.

Best regards,
Martin Maechler
ETH Zurich  and  R Core Team

> -- 
> 
> --sep
> 
> 	[[alternative HTML version deleted]]
>