R-beta: more model.matrix

Douglas Bates bates at stat.wisc.edu
Fri Oct 17 19:25:48 CEST 1997


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



More information about the R-help mailing list