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]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
>>>>> 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]] > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
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 <[hidden email] > 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]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
In reply to this post by Martin Maechler
> BTW --- and this is a diversion --- This is nice mathematically
A fine idea for sure, but perhaps use a different name in order to avoid
> (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. > conflicts with: > base::stderr function () .Internal(stderr()) <bytecode: 0x32e1e98> <environment: namespace:base> ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel |
Free forum by Nabble | Edit this page |