Re: [Rd] Should there be a confint.mlm ?
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. A fine idea for sure, but perhaps use a different name in order to avoid conflicts with: > base::stderr function () .Internal(stderr()) __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Should there be a confint.mlm ?
I would welcome the fixes suggested by Martin. I did not think my
implementation was pretty, but didn't want to suggest a bugfix without
submitting code.
Regarding the idea of a `stderr` function, I am all for potential
efficiency gains, but I suspect that in many situations `vcov` is the
result of a matrix inversion. This might mean an 'efficient computation of
only the diagonal of the inverse of a matrix' function would be required to
achieve widespread use. (I am not aware of the algorithm for that, but I am
ignorant.)
Again, thanks for patching `confint.lm`.
On Fri, Jul 20, 2018 at 9:06 AM, Martin Maechler wrote:
> > 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]]
> >
>
--
--sep
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Should there be a confint.mlm ?
> 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]]
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Should there be a confint.mlm ?
It seems that confint.default returns an empty data.frame for objects of
class mlm. 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. A proposed
implementation would be:
```
# I need this to run the code, but stats does not:
format.perc <- stats:::format.perc
# 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)))
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
```
--
--sep
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
