I would like to write a function that computes Tukey's 1 df for
nonadditivity. Here is a simplified version of the function I'd like to
write: (m is an object created by lm):
tukey.test <- function(m) {
m1 <- update(m, ~.+I(predict(m)^2))
summary(m1)$coef
}
The t-test for the added variable is Tukey's test. This won't work:
data(BOD)
m1 <- lm(demand~Time,BOD)
tukey.test(m1)
Error in predict(m) : object "m" not found
This function doesn't work for two reasons:
1. The statement m1 <- update(m, ~.+I(predict(m)^2)) can't see
'm' in
the call to predict.
2. If in creating m missing values had been present, then predict(m),
even if it could be computed, could be of the wrong length.
Can anyone help?
--
Sanford Weisberg, sandy at stat.umn.edu
Office and mailing address:
University of Minnesota, School of Statistics
312 Ford Hall, 224 Church St. SE, Minneapolis, MN 55455
612-625-8355, FAX 612-624-8868
St. Paul office:
146 Classroom-Office Building, 612-625-8777
On 10/24/2007 3:43 PM, Sandy Weisberg wrote:> I would like to write a function that computes Tukey's 1 df for > nonadditivity. Here is a simplified version of the function I'd like to > write: (m is an object created by lm): > > tukey.test <- function(m) { > m1 <- update(m, ~.+I(predict(m)^2)) > summary(m1)$coef > } > > The t-test for the added variable is Tukey's test. This won't work: > > data(BOD) > m1 <- lm(demand~Time,BOD) > tukey.test(m1) > > Error in predict(m) : object "m" not found > > This function doesn't work for two reasons: > 1. The statement m1 <- update(m, ~.+I(predict(m)^2)) can't see 'm' in > the call to predict.It is looking in the environment of m$terms, and failing to find m there. Here's an ugly way to do what you want; there might be a cleaner way, but I don't know it: tukey.test <- function(m) { # Create a new environment whose parent is where update is currently # looking newenv <- new.env(parent=environment(m$terms)) # Put the new predictor there. You could instead put m there, but # that will create a loop whose evaluation will destroy the universe newenv$predm2 <- predict(m)^2 # Attach the new environment in place of the old one environment(m$terms) <- newenv # Proceed as before m1 <- update(m, ~.+predm2) summary(m1)$coef } Another variation is just to assign predm2 into the existing environment, i.e. environment(m$terms)$predm2 <- predict(m)^2 but that's potentially destructive of an existing variable named predm2 so I didn't do it. Environments are reference objects, so you would kill the real one, and have side effects outside your scope: bad idea.> 2. If in creating m missing values had been present, then predict(m), > even if it could be computed, could be of the wrong length.Not sure about the best approach here, but the prediction vector has names indicating which observations were used, so you can probably do something with those. Duncan Murdoch Duncan
Another way to do this without messing with environments is to update the
data and formula locally within the function and re-run the regression on
the updated model/data:
tukey.test <- function(m) {
ud <- data.frame( m$model, pred2=m$fitted.values^2 )
uf <- update.formula(formula(m$terms), ~ . + pred2)
summary(lm(uf, ud))$coef
}
> data(BOD)
> m1 <- lm(demand~Time,BOD)
> tukey.test(m1)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.5854086 5.0344850 2.4998403 0.0877190
Time 7.5390634 6.1271204 1.2304415 0.3062126
pred2 -0.1096744 0.1148654 -0.9548081 0.4101142
> -----Original Message-----
> From: r-help-bounces at r-project.org
> [mailto:r-help-bounces at r-project.org] On Behalf Of Sandy Weisberg
> Sent: Wednesday, October 24, 2007 3:44 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] scoping problem
>
> I would like to write a function that computes Tukey's 1 df
> for nonadditivity. Here is a simplified version of the
> function I'd like to
> write: (m is an object created by lm):
>
> tukey.test <- function(m) {
> m1 <- update(m, ~.+I(predict(m)^2))
> summary(m1)$coef
> }
>
> The t-test for the added variable is Tukey's test. This won't
work:
>
> data(BOD)
> m1 <- lm(demand~Time,BOD)
> tukey.test(m1)
>
> Error in predict(m) : object "m" not found
>
> This function doesn't work for two reasons:
> 1. The statement m1 <- update(m, ~.+I(predict(m)^2))
> can't see 'm' in the call to predict.
> 2. If in creating m missing values had been present,
> then predict(m), even if it could be computed, could be of
> the wrong length.
>
> Can anyone help?
>
>
>
> --
> Sanford Weisberg, sandy at stat.umn.edu
> Office and mailing address:
> University of Minnesota, School of Statistics
> 312 Ford Hall, 224 Church St. SE, Minneapolis, MN 55455
> 612-625-8355, FAX 612-624-8868
>
> St. Paul office:
> 146 Classroom-Office Building, 612-625-8777
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
See: http://finzi.psych.upenn.edu/R/Rhelp02a/archive/67474.html and also the other posts in that thread. On 10/24/07, Sandy Weisberg <sandy at stat.umn.edu> wrote:> I would like to write a function that computes Tukey's 1 df for > nonadditivity. Here is a simplified version of the function I'd like to > write: (m is an object created by lm): > > tukey.test <- function(m) { > m1 <- update(m, ~.+I(predict(m)^2)) > summary(m1)$coef > } > > The t-test for the added variable is Tukey's test. This won't work: > > data(BOD) > m1 <- lm(demand~Time,BOD) > tukey.test(m1) > > Error in predict(m) : object "m" not found > > This function doesn't work for two reasons: > 1. The statement m1 <- update(m, ~.+I(predict(m)^2)) can't see 'm' in > the call to predict. > 2. If in creating m missing values had been present, then predict(m), > even if it could be computed, could be of the wrong length. > > Can anyone help? > > > > -- > Sanford Weisberg, sandy at stat.umn.edu > Office and mailing address: > University of Minnesota, School of Statistics > 312 Ford Hall, 224 Church St. SE, Minneapolis, MN 55455 > 612-625-8355, FAX 612-624-8868 > > St. Paul office: > 146 Classroom-Office Building, 612-625-8777 > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >