se.contrast: matrix contrast.obj doesn't work as documented (PR#1613)

cmascott@att.net cmascott@att.net
Thu, 30 May 2002 17:55:06 +0200 (MET DST)


The man page for se.contrast, when describing the contrast.obj
parameter, states that "Multiple contrasts should be specified
by a matrix as returned by contrasts."

When doing an unbalanced single factor ANOVA, using a contrast.obj
as returned by contrasts results in the following error from
qr.qty when se.contrast is called:

Error in qr.qty(object$qr, contrast) : qr and y must have the same number of rows

In fact, the matrix contrast.obj parameter that se.contrast
requires is not a contrast matrix as returned by contrasts.
It is a matrix each of whose columns is in the same form as
the vector constructed by se.contrast when contrast.obj is
a list.  This matrix has one row per observation, not one
row per treatment (contrast matrix).  See the session log
below for an example.

I wrote a little function to produce the one-row-per-
observation matrix and se.contrast works correctly when
given this matrix.

-----------------------------------------------------------------
Script started on Thu May 30 11:21:31 2002

R : Copyright 2002, The R Development Core Team
Version 1.5.0  (2002-04-29)

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 more information.

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.

[Previously saved workspace restored]

> debug("qr.qty")
> se.contrast(reading.aov.t, list(group=="1", group=="2"), c(-1,1), reading.df)
debugging in: qr.qty(object$qr, contrast)
debug: {
    if (!is.qr(qr)) 
        stop("argument is not a QR decomposition")
    if (is.complex(qr$qr)) {
        y <- as.matrix(y)
        if (!is.complex(y)) 
            y[] <- as.complex(y)
        return(.Call("qr_qy_cmplx", qr, y, 1, PACKAGE = "base"))
    }
    n <- nrow(qr$qr)
    p <- ncol(qr$qr)
    k <- as.integer(qr$rank)
    ny <- NCOL(y)
    if (NROW(y) != n) 
        stop("qr and y must have the same number of rows")
    storage.mode(y) <- "double"
    .Fortran("dqrqty", as.double(qr$qr), n, k, as.double(qr$qraux), 
        y, ny, qty = y, PACKAGE = "base")$qty
}
Browse[1]> y
       
         [,1]
   [1,] -0.20
   [2,] -0.20
   [3,] -0.20
   [4,] -0.20
   [5,] -0.20
   [6,]  0.25
   [7,]  0.25
   [8,]  0.25
   [9,]  0.25
  [10,]  0.00
  [11,]  0.00
  [12,]  0.00
  [13,]  0.00
  [14,]  0.00
  [15,]  0.00
  [16,]  0.00
  [17,]  0.00
  [18,]  0.00
  [19,]  0.00
  [20,]  0.00
  [21,]  0.00
  [22,]  0.00
  [23,]  0.00
Browse[1]> c
exiting from: qr.qty(object$qr, contrast)
[1] 3.126500
> undebug("qr.qty")
> se.contrast(reading.aov.t, contr.t5, data=reading.df)
Error in qr.qty(object$qr, contrast) : qr and y must have the same number of rows
> contr.t5
   2  3  4  5
1 -1 -1 -1 -1
2  1  0  0  0
3  0  1  0  0
4  0  0  1  0
5  0  0  0  1
> q()
Save workspace image? [y/n/c]: n

Script done on Thu May 30 11:25:49 2002
-----------------------------------------------------------------


--please do not edit the information below--

Version:
 platform = i386--freebsd4.3
 arch = i386
 os = freebsd4.3
 system = i386, freebsd4.3
 status = 
 major = 1
 minor = 5.0
 year = 2002
 month = 04
 day = 29
 language = R

Search Path:
 .GlobalEnv, package:ctest, Autoloads, package:base

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