[R] se.contrast ....too hard??? .... Too easy????? .....too trivial???? ...... Too boring.....too????????

[Duncan Mackay]

Hi all,
Regular and avid readers of this column will know that Don Driscoll and
I have recently posted two messages requesting assistance concerning an
apparent failure of "se.contrast" to produce an se for a contrast. So
far, an ominous silence rings in our ears, but read on Gentle Reader,
and see if even the machinations of "debug" doesn't stimulate you to
respond with a revelatory epistle.

Thanks!!!
Duncan



Here is my first message :-

Just to follow up Don Driscoll's earlier post, can anyone please explain
why "se.contrast" fails here??

> shp<-factor(rep(c("reserve","strip"),each=96))
> 
> 
>
site<-factor(rep(c("1g","1p","1t","2g","2p","2t","3g","3p","3t","4g","4p
","4t"),each=16))
> 
> pit<-factor(rep(1:16,12))
> 
>
reptsp<-c(4,5,6,4,6,6,6,7,3,5,2,2,4,8,5,4,2,4,2,2,4,5,2,4,4,4,3,2,3,2,5,
3,5,3,4,4,4,3,4,
+
3,4,4,4,3,4,3,6,3,3,5,4,6,4,4,2,4,2,6,5,5,5,7,4,4,5,1,4,5,6,5,5,2,6,3,5,
6,4,5,
+
4,8,2,4,2,4,2,4,3,3,4,4,3,2,1,3,4,4,2,2,3,2,4,1,2,2,3,4,5,5,3,5,5,4,1,1,
2,1,3,
+
1,4,1,6,1,2,3,2,2,2,1,1,2,2,6,5,3,2,3,5,3,2,3,2,1,3,2,4,4,3,3,3,1,2,4,3,
4,5,6,
+
5,2,3,2,2,5,5,5,2,2,5,2,4,4,3,2,2,3,2,2,2,2,5,4,3,3,5,2,5,4,3,2,2,2,1,2)
> 
> ddata<-data.frame(shp,pit,site,reptsp)
> 
> repmod2<-aov(reptsp~shp/site+ Error(shp/site))
> summary(repmod2)

Error: shp
    Df Sum Sq Mean Sq
shp  1  53.13   53.13

Error: shp:site
         Df Sum Sq Mean Sq
shp:site 10 61.885   6.189

Error: Within
           Df Sum Sq Mean Sq F value Pr(>F)
Residuals 180 318.56    1.77               
> 
> table(ddata$shp)

reserve   strip 
     96      96 
> 
> se.contrast(repmod2, list(shp=="strip", shp=="reserve"),data=ddata)
Error in qr.qty(strata$qr, scontrast) : qr and y must have the same
number of rows
> 
?????????????????????


Thanks,
Duncan


...............and here is what debug says:-

> debug("qr.qty")
> se.contrast(repmod2, list(shp=="strip", shp=="reserve"),data=ddata)
debugging in: qr.qty(e.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"))
    }
    a <- attr(qr, "useLAPACK")
    if (!is.null(a) && is.logical(a) && a) 
        return(.Call("qr_qy_real", qr, as.matrix(y), 1, PACKAGE =
"base"))
    n <- nrow(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]> 
debug: if (!is.qr(qr)) stop("argument is not a QR decomposition")
Browse[1]> 
debug: 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"))
}
Browse[1]> 
debug: a <- attr(qr, "useLAPACK")
Browse[1]> 
debug: if (!is.null(a) && is.logical(a) && a) return(.Call("qr_qy_real",

    qr, as.matrix(y), 1, PACKAGE = "base"))
Browse[1]> 
debug: n <- nrow(qr$qr)
Browse[1]> 
debug: k <- as.integer(qr$rank)
Browse[1]> 
debug: ny <- NCOL(y)
Browse[1]> 
debug: if (NROW(y) != n) stop("qr and y must have the same number of
rows") Browse[1]> ny [1] 1 Browse[1]> NCOL(y) [1] 1 Browse[1]>
nrow(qr$qr) [1] 192 Browse[1]> print(n) [1] 192 Browse[1]> dim(qr$qr)
[1] 192  24 Browse[1]> NROW(y) [1] 192 Browse[1]> NROW(y)==n [1] TRUE
Browse[1]> Q

So, if NROW(y)==n, why do I get the message "qr and y must have the same
number of rows"

Any help gratefully appreciated....
Duncan

**********************************************

R, me bucko

*****************************************
Dr. Duncan Mackay
School of Biological Sciences
Flinders University
GPO Box 2100
Adelaide
S.A.    5001
AUSTRALIA

Ph (08) 8201 2627    FAX (08) 8201 3015

http://www.scieng.flinders.edu.au/biology/people/mackay_d/index.html

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html