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