From jfox@mcmaster.ca  Tue Mar  7 14:08:03 2000
From: jfox@mcmaster.ca (jfox@mcmaster.ca)
Date: Tue, 7 Mar 2000 15:08:03 +0100 (MET)
Subject: [Rd] update fails after specific sequence of steps (PR#474)
Message-ID: <200003071408.PAA18937@pubhealth.ku.dk>

# 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] update fails after specific sequence of steps (PR#474)
jfox@mcmaster.ca