R help - Oct 1997 - R-beta: more model.matrix

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. 

 R> fm1 <- lm(y ~ x.1 + x.2 + x.3 + x.4 + x.5, data = e3.7, model = T)
 R> summary(fm1)

 Call:
 lm(formula = y ~ x.1 + x.2 + x.3 + x.4 + x.5, data = e3.7, model = T)

 Residuals:
	Min         1Q     Median         3Q        Max 
 -0.3944667 -0.1184741  0.0005345  0.0831318  0.5623167 

 Coefficients:
	     Estimate Std.Error t Value Pr(>|t|)
 (Intercept)  -2.1561    0.9135 -2.3603   0.0333
 x.1           0.0000    0.0005 -0.0174   0.9864
 x.2           0.0013    0.0013  1.0415   0.3153
 x.3           0.0001    0.0001  1.6618   0.1188
 x.4           0.0079    0.0140  0.5642   0.5815
 x.5           0.0001    0.0001  1.9208   0.0754

 Residual standard error: 0.2618 on 14 degrees of freedom
 Multiple R-Squared: 0.8107,	Adjusted R-squared: 0.743 
 F-statistic: 11.99 on 5 and 14 degrees of freedom,	p-value: 0.0001184 

 R> fm1Sampl <- bsCoefSample(fm1, 500)
 Error: Object "y" not found

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.
-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=