jfox@mcmaster.ca
2000-Mar-07 14:08 UTC
[Rd] update fails after specific sequence of steps (PR#474)
# Your mailer is set to "none" (default on Windows),
# hence we cannot send the bug report directly from R.
# Please copy the bug report (after finishing it) to
# your favorite email program and send it to
#
# r-bugs@biostat.ku.dk
#
######################################################
I stumbled on this error while doing a classroom demonstration. The error is
reproducible,
but seems to require a specific sequence of operations: (1) fitting a linear
model; (2)
calling the summary function; (3) calling a function of mine, called test.terms;
(4) making
a typing error in a call to update -- typing "~.~" instead of
".~." . (5) After this, a
correct call to update fails. If any of the steps 1-4 is omitted, everything
proceeds
normally.
I've attached the following:
o several listings illustrating when the error occurs and when it does not;
o a listing of the dataset;
o a listing of the test.terms function and the functions that it uses.
The dataset, though, seems unremarkable, and test.tests is a straightforward
function
without (intended) side effects. The library fox just contains some functions
and datasets;
I can, of course, supply it if that seems relevant.
I'm not sure that this is a bug in R, but it seemed worth reporting.
Thanks for your trouble,
John
----------------- first listing, showing error ---------------------------
R : Copyright 2000, The R Development Core Team
Version 1.0.0 (February 29, 2000)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type "?license" or "?licence" for distribution details.
R is a collaborative project with many contributors.
Type "?contributors" for a list.
Type "demo()" for some demos, "help()" for on-line help,
or
"help.start()" for a HTML browser interface to help.
Type "q()" to quit R.
> library(fox)
> data(Moore)
> attach(Moore)
> mod.full<-lm(conform~status*fcategory)
> summary(mod.full)
Call:
lm(formula = conform ~ status * fcategory)
Residuals:
Min 1Q Median 3Q Max
-8.6250 -2.9000 -0.2727 2.7273 11.3750
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.8571 1.7307 6.851 3.44e-08 ***
statuslow 0.7679 2.3699 0.324 0.7477
fcategorylow 5.5429 2.6813 2.067 0.0454 *
fcategorymedium 2.4156 2.2140 1.091 0.2819
statuslow.fcategorylow -9.2679 3.4507 -2.686 0.0106 *
statuslow.fcategorymedium -7.7906 3.5728 -2.181 0.0353 *
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1
` ' 1
Residual standard error: 4.579 on 39 degrees of freedom
Multiple R-Squared: 0.3237, Adjusted R-squared: 0.237
F-statistic: 3.734 on 5 and 39 degrees of freedom, p-value: 0.007397
> test.terms(mod.full)
Sum Sq Df F value Pr(>F)
(Intercept) 984.14 1 46.9348 3.436e-08 ***
status 2.20 1 0.1050 0.74767
fcategory 89.67 2 2.1383 0.13147
status:fcategory 175.49 2 4.1846 0.02257 *
Residuals 817.76 39
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1
` ' 1 > mod.m<-update(mod.full, ~.~-status:fcategory)
Error in model.frame.default(formula = ~conform ~ 1, drop.unused.levels = TRUE)
:
object is not a matrix> mod.m<-update(mod.full, .~.-status:fcategory) # fails after summary,
test.terms, & error
Error in parse(file, n, text, prompt) : parse error>
----------------- second listing, omitting error in call to update, everything
works -------
R : Copyright 2000, The R Development Core Team
Version 1.0.0 (February 29, 2000)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type "?license" or "?licence" for distribution details.
R is a collaborative project with many contributors.
Type "?contributors" for a list.
Type "demo()" for some demos, "help()" for on-line help,
or
"help.start()" for a HTML browser interface to help.
Type "q()" to quit R.
> library(fox)
> data(Moore)
> attach(Moore)
> mod.full<-lm(conform~status*fcategory)
> summary(mod.full)
Call:
lm(formula = conform ~ status * fcategory)
Residuals:
Min 1Q Median 3Q Max
-8.6250 -2.9000 -0.2727 2.7273 11.3750
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.8571 1.7307 6.851 3.44e-08 ***
statuslow 0.7679 2.3699 0.324 0.7477
fcategorylow 5.5429 2.6813 2.067 0.0454 *
fcategorymedium 2.4156 2.2140 1.091 0.2819
statuslow.fcategorylow -9.2679 3.4507 -2.686 0.0106 *
statuslow.fcategorymedium -7.7906 3.5728 -2.181 0.0353 *
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1
` ' 1
Residual standard error: 4.579 on 39 degrees of freedom
Multiple R-Squared: 0.3237, Adjusted R-squared: 0.237
F-statistic: 3.734 on 5 and 39 degrees of freedom, p-value: 0.007397
> test.terms(mod.full)
Sum Sq Df F value Pr(>F)
(Intercept) 984.14 1 46.9348 3.436e-08 ***
status 2.20 1 0.1050 0.74767
fcategory 89.67 2 2.1383 0.13147
status:fcategory 175.49 2 4.1846 0.02257 *
Residuals 817.76 39
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1
` ' 1 > mod.m<-update(mod.full, .~.-status:fcategory) # works following summary
& test.terms
> mod.m
Call:
lm(formula = conform ~ status + fcategory)
Coefficients:
(Intercept) statuslow fcategorylow fcategorymedium
14.72356 -4.60667 0.08089 -1.09511
>
----- third listing, omitting call to test.terms, everything works
--------------------
R : Copyright 2000, The R Development Core Team
Version 1.0.0 (February 29, 2000)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type "?license" or "?licence" for distribution details.
R is a collaborative project with many contributors.
Type "?contributors" for a list.
Type "demo()" for some demos, "help()" for on-line help,
or
"help.start()" for a HTML browser interface to help.
Type "q()" to quit R.
> library(fox)
> data(Moore)
> attach(Moore)
> mod.full<-lm(conform~status*fcategory)
> summary(mod.full)
Call:
lm(formula = conform ~ status * fcategory)
Residuals:
Min 1Q Median 3Q Max
-8.6250 -2.9000 -0.2727 2.7273 11.3750
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.8571 1.7307 6.851 3.44e-08 ***
statuslow 0.7679 2.3699 0.324 0.7477
fcategorylow 5.5429 2.6813 2.067 0.0454 *
fcategorymedium 2.4156 2.2140 1.091 0.2819
statuslow.fcategorylow -9.2679 3.4507 -2.686 0.0106 *
statuslow.fcategorymedium -7.7906 3.5728 -2.181 0.0353 *
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1
` ' 1
Residual standard error: 4.579 on 39 degrees of freedom
Multiple R-Squared: 0.3237, Adjusted R-squared: 0.237
F-statistic: 3.734 on 5 and 39 degrees of freedom, p-value: 0.007397
> mod.m<-update(mod.full, ~.~-status:fcategory)
Error in model.frame.default(formula = ~conform ~ 1, drop.unused.levels = TRUE)
:
object is not a matrix> mod.m<-update(mod.full, .~.-status:fcategory) # works following summary
and error
> mod.m
Call:
lm(formula = conform ~ status + fcategory)
Coefficients:
(Intercept) statuslow fcategorylow fcategorymedium
14.72356 -4.60667 0.08089 -1.09511
>
----------- fourth listing, omitting call to summary, everything works
------------------
R : Copyright 2000, The R Development Core Team
Version 1.0.0 (February 29, 2000)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type "?license" or "?licence" for distribution details.
R is a collaborative project with many contributors.
Type "?contributors" for a list.
Type "demo()" for some demos, "help()" for on-line help,
or
"help.start()" for a HTML browser interface to help.
Type "q()" to quit R.
> library(fox)
> data(Moore)
> attach(Moore)
> mod.full<-lm(conform~status*fcategory)
> test.terms(mod.full)
Sum Sq Df F value Pr(>F)
(Intercept) 984.14 1 46.9348 3.436e-08 ***
status 2.20 1 0.1050 0.74767
fcategory 89.67 2 2.1383 0.13147
status:fcategory 175.49 2 4.1846 0.02257 *
Residuals 817.76 39
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1
` ' 1 > mod.m<-update(mod.full, ~.~-status:fcategory)
Error in model.frame.default(formula = ~conform ~ 1, drop.unused.levels = TRUE)
:
object is not a matrix> mod.m<-update(mod.full, .~.-status:fcategory) # works following
test.terms and error
> mod.m
Call:
lm(formula = conform ~ status + fcategory)
Coefficients:
(Intercept) statuslow fcategorylow fcategorymedium
14.72356 -4.60667 0.08089 -1.09511
------------------------ listing of dataset
----------------------------------------
> Moore
subject status conform fcategory fscore
1 1 low 8 low 37
2 2 low 4 high 57
3 3 low 8 high 65
4 4 low 7 low 20
5 5 low 10 low 36
6 6 low 6 low 18
7 7 low 12 medium 51
8 8 low 4 medium 44
9 9 low 13 low 31
10 10 low 12 low 36
11 11 low 4 medium 42
12 12 low 13 high 56
13 13 low 7 low 28
14 14 low 9 medium 43
15 15 low 9 high 65
16 16 low 24 high 57
17 17 low 6 low 28
18 18 low 7 high 61
19 19 low 23 high 57
20 20 low 13 high 55
21 21 low 8 low 17
22 22 low 12 low 35
23 23 high 19 high 68
24 24 high 12 medium 41
25 25 high 21 low 16
26 26 high 9 high 63
27 27 high 23 low 15
28 28 high 7 high 54
29 29 high 17 medium 48
30 30 high 14 medium 42
31 31 high 11 high 59
32 32 high 16 medium 45
33 33 high 15 low 30
34 34 high 20 medium 44
35 35 high 8 medium 42
36 36 high 12 low 22
37 37 high 14 high 52
38 38 high 14 medium 42
39 39 high 17 medium 41
40 40 high 7 medium 50
41 41 high 17 medium 39
42 42 high 13 high 57
43 43 high 16 low 35
44 44 high 10 high 52
45 45 high 15 medium 44>
------------------ function test.terms and functions it calls
------------------------
> test.terms
function (object, ...)
{
UseMethod("test.terms")
}
> test.terms.lm
function (mod, error, ...)
{
if (!missing(error)) {
sumry <- summary(error, corr = FALSE)
s2 <- sumry$sigma^2
error.df <- error$df.residual
error.SS <- s2 * error.df
}
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
Source <- term.names(mod)
n.terms <- length(Source)
SS <- rep(0, n.terms + 1)
df <- rep(0, n.terms + 1)
F <- rep(0, n.terms + 1)
p <- rep(0, n.terms + 1)
assign <- mod$assign
for (term in 1:n.terms) {
subs <- which(assign == term - intercept)
hyp.matrix <- I.p[subs, ]
test <- if (missing(error))
linear.hypothesis(mod, hyp.matrix, ...)
else linear.hypothesis(mod, hyp.matrix, error.SS = error.SS,
error.df = error.df, ...)
SS[term] <- test$SS
df[term] <- test$df[1]
F[term] <- test$F
p[term] <- test$p
}
Source[n.terms + 1] <- "Residuals"
df.res <- if (missing(error))
mod$df.residual
else error.df
sumry <- summary(mod, corr = FALSE)
s2 <- sumry$sigma^2
SS[n.terms + 1] <- if (missing(error))
s2 * df.res
else error.SS
df[n.terms + 1] <- df.res
F[n.terms + 1] <- NA
p[n.terms + 1] <- NA
result <- data.frame(SS, df, F, p)
row.names(result) <- Source
names(result) <- c("Sum Sq", "Df", "F
value", "Pr(>F)")
class(result) <- c("anova", "data.frame")
result
}
> has.intercept
function (object, ...)
{
UseMethod("has.intercept")
}
> term.names
function (object, ...)
{
UseMethod("term.names")
}
> term.names.default
function (model)
{
term.names <- labels(terms(model))
if (has.intercept(model))
c("(Intercept)", term.names)
else term.names
}
> linear.hypothesis
function (object, ...)
{
UseMethod("linear.hypothesis")
}
> linear.hypothesis.lm
function (mod, hypothesis.matrix, rhs = 0, white.adjust = F,
error.SS, error.df)
{
sumry <- summary(mod, corr = FALSE)
s2 <- if (missing(error.SS))
sumry$sigma^2
else error.SS/error.df
V <- if (white.adjust == FALSE)
sumry$cov.unscaled
else hccm(mod, type = white.adjust)/s2
b <- coefficients(mod)
L <- if (is.null(dim(hypothesis.matrix)))
t(hypothesis.matrix)
else hypothesis.matrix
q <- nrow(L)
SSH <- t(L %*% b - rhs) %*% inv(L %*% V %*% t(L)) %*% (L %*%
b - rhs)
F <- SSH/(q * s2)
df <- if (missing(error.df))
mod$df.residual
else error.df
p <- 1 - pf(F, q, df)
list(SS = SSH[1, 1], F = F[1, 1], df = c(q, df), p = p[1,
1])
}
> hccm
function (model, type = "hc3")
{
if (is.null(class(model)) || class(model) != "lm" ||
!is.null(weights(mod))) {
stop("requires unweighted lm")
}
sumry <- summary(model, corr = FALSE)
s2 <- sumry$sigma^2
V <- sumry$cov.unscaled
if (type == FALSE)
return(s2 * V)
e <- residuals(model)
X <- model.matrix(model)
df.res <- df.residual(model)
factor <- switch(type, hc0 = 1, hc1 = n/df.res, hc2 = 1/(1 -
hat(X)), hc3 = 1/(1 - hat(X))^2)
V %*% t(X) %*% apply(X, 2, "*", (e^2)/factor) %*% V
}>
--please do not edit the information below--
Version:
platform = Windows
arch = x86
os = Win32
system = x86, Win32
status =
major = 1
minor = 0.0
year = 2000
month = February
day = 29
language = R
Windows NT 4.0 (build 1381) Service Pack 3
Search Path:
.GlobalEnv, Moore, package:fox, Autoloads, package:base
|----------------------------------------------------|
| John Fox jfox@McMaster.ca |
| Department of Sociology McMaster University |
|----------------------------------------------------|
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel 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-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Peter Dalgaard BSA
2000-Mar-07 14:43 UTC
[Rd] update fails after specific sequence of steps (PR#474)
Strange things happen on Unix too.> source('fox.R') > Moore<-read.table('Moore.txt') > attach(Moore) > mod.full<-lm(conform~status*fcategory) > summary(mod.full)...> mod.m<-update(mod.full, ~.~-status:fcategory)Error in model.frame.default(formula = ~conform ~ 1, drop.unused.levels = TRUE) : object is not a matrix> mod.m<-update(mod.full, .~.-status:fcategory) # works following summary and error > mod.mCall: lm(formula = conform ~ status + fcategory) Coefficients: (Intercept) statuslow fcategorylow fcategorymedium 14.72356 -4.60667 0.08089 -1.09511 #### Well you *thought* it worked, but watch: ####> mod.full<-lm(conform~status*fcategory) > summary(mod.full)Call: lm(formula = conform ~ L׀@ * fcategory) ^^^^^^^ ...and the attached dataframe is corrupted. Must investigate further.> ls("Moore")[1] "conform" "fcategory" "fscore" "L׀@" "subject" ^^^^^^^ -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel 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-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._