R-beta: more model.matrix

[Thomas Lumley]

On Fri, 17 Oct 1997, Douglas Bates wrote:

> I am trying to show some techniques to my graduate regression class.
> The textbook mentioned using bootstrap samples of regression
> coefficients for assessing variability.  I decided to show them
> reasonably effective ways of doing the resampling.
> 
> The following is a function I wrote to create bootstrap samples of
> coefficients from a fitted linear regression model.
> 
>  bsCoefSample <- 
>    ## Construct a bootstrap sample of coefficients from a
>    ## fitted regression model
>    function(fittedModel, Nsampl, ...)
>  {
>    coef(fittedModel) +
>      qr.coef(qr(model.matrix(fittedModel)),
> 	     matrix(sample(resid(fittedModel),
> 			   length(resid(fittedModel)) * Nsampl,
> 			   repl = T),
> 		    ncol = Nsampl))
>  }
> 
> It works ok with S-PLUS but the R version of model.matrix encounters
> difficulties. 

 
> Can anyone suggest a modified function that would work in both S-PLUS
> and R?
> 
> It appears that model.matrix is a generic in S-PLUS
>  S-PLUS> methods(model.matrix)
>  $SHOME/stat/.Functions $SHOME/stat/.Functions 
>  "model.matrix.default" "model.matrix.lm"     
> 
> Could this be done in R as well?  I am willing to try to do that if
> it seems like a good idea.

We can define 

model.matrix.lm<-function(lmobj){
	model.matrix(delete.response(terms(lmobj)),model.frame(lmobj))
}

so that an explicit call to model.matrix.lm instead of model.matrix would
work in both R and S.  It might well be a good idea to make model.matrix
generic. 

Another option is to use update(,qr=T) to refit the model and return the
qr decomposition.  This is a bit more tricky, since you need to evaluate
things in the right environment.


Thomas Lumley
------------------------------------------------------+------
Biostatistics		: "Never attribute to malice what  :
Uni of Washington	:  can be adequately explained by  :
Box 357232		:  incompetence" - Hanlon's Razor  :
Seattle WA 98195-7232	:				   :
------------------------------------------------------------

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=