################################################### ### chunk number 1: preliminaries ################################################### options(width=85, show.signif.stars = FALSE, lattice.theme = function() canonical.theme("pdf", color = FALSE), str = strOptions(strict.width = "cut")) ################################################### ### chunk number 2: strFormaldehyde ################################################### str(Formaldehyde) ################################################### ### chunk number 3: lm1 ################################################### (X <- model.matrix(lm1 <- lm(optden ~ 1 + carb, Formaldehyde))) ################################################### ### chunk number 4: classqr ################################################### class(qrlm1 <- lm1$qr) ################################################### ### chunk number 5: R ################################################### (R <- qr.R(qrlm1)) ################################################### ### chunk number 6: Q1 ################################################### (Q1 <- qr.Q(qrlm1)) ################################################### ### chunk number 7: Q1 ################################################### (Q1R <- Q1 %*% R) ################################################### ### chunk number 8: badcomp ################################################### all(X == Q1R) ################################################### ### chunk number 9: goodcomp ################################################### all.equal(X, Q1R, check.attr = FALSE) ################################################### ### chunk number 10: nmsLm1 ################################################### names(lm1) ################################################### ### chunk number 11: yvec ################################################### (y <- model.response(model.frame(lm1))) ################################################### ### chunk number 12: Qmat ################################################### (Q <- qr.Q(qrlm1, complete=TRUE)) ################################################### ### chunk number 13: checkeffects ################################################### str(cbind(lm1$effects)) str(crossprod(Q, y)) all.equal(cbind(lm1$effects), crossprod(Q, y), check.attr=FALSE) ################################################### ### chunk number 14: XtX ################################################### crossprod(X) ################################################### ### chunk number 15: checkeffects ################################################### all.equal(lm1$effects, as.vector(crossprod(Q, y)), check.attr=FALSE) ################################################### ### chunk number 16: qrqty ################################################### all.equal(lm1$effects, qr.qty(qrlm1, y), check.attr=FALSE) ################################################### ### chunk number 17: Q1tQ1 ################################################### all.equal(crossprod(Q1), diag(nrow=ncol(Q1))) all.equal(crossprod(Q), diag(nrow=nrow(Q))) all.equal(tcrossprod(Q), diag(nrow=nrow(Q))) ################################################### ### chunk number 18: Q1tQ1print ################################################### crossprod(Q1) ################################################### ### chunk number 19: Q1tQ1print ################################################### zapsmall(crossprod(Q1)) zapsmall(crossprod(Q)) ################################################### ### chunk number 20: coef ################################################### coef(lm1) backsolve(R, crossprod(Q1, y)) all.equal(coef(lm1), as.vector(backsolve(R, crossprod(Q1, y))), check.attr=FALSE) ################################################### ### chunk number 21: coef1 ################################################### qr.coef(qrlm1, y) ################################################### ### chunk number 22: rankDeficient ################################################### set.seed(1234) # allow for reproducible "random" numbers badDat <- within(data.frame(x1=1:20, x2=rnorm(20,mean=6,sd=0.2), x4=rexp(20,rate=0.02), y=runif(20,min=18,max=24)), x3 <- x1 + 2*x2) # create linear combination (summary(lm2 <- lm(y ~ x1 + x2 + x3 + x4, badDat))) (lm2qr <- lm2$qr)$rank qr.R(lm2qr) # the columns are rearranged lm2qr$pivot # the permutation vector ################################################### ### chunk number 23: InsectSprays ################################################### str(InsectSprays) unique(mm <- model.matrix(lm3 <- lm(count ~ spray, InsectSprays))) attr(mm, "assign") qr.R(lm3[["qr"]]) ################################################### ### chunk number 24: getOptionscontrasts ################################################### getOption("contrasts") ################################################### ### chunk number 25: ################################################### all.equal(det(R), prod(diag(R))) det(crossprod(Q1)) det(crossprod(Q)) ################################################### ### chunk number 26: betahat ################################################### backsolve(R, crossprod(Q1, y)) ################################################### ### chunk number 27: betahat ################################################### qr.coef(qrlm1, y) ################################################### ### chunk number 28: fittedlm1 ################################################### qr.fitted(qrlm1, y) all.equal(qr.fitted(qrlm1, y), fitted(lm1)) ################################################### ### chunk number 29: residlm1 ################################################### qr.resid(qrlm1, y) all.equal(qr.resid(qrlm1, y), residuals(lm1)) ################################################### ### chunk number 30: P1 ################################################### (P1 <- tcrossprod(Q1)) ################################################### ### chunk number 31: P1 ################################################### (P2 <- tcrossprod(Q[, -(1:2)])) ################################################### ### chunk number 32: idempotent ################################################### all.equal(P1 %*% P1, P1) all.equal(P2 %*% P2, P2) ################################################### ### chunk number 33: P1X ################################################### all.equal(P1 %*% X, X, check.attr=FALSE) ################################################### ### chunk number 34: P1X ################################################### all.equal(P2 %*% X, 0 * X, check.attr=FALSE) ################################################### ### chunk number 35: diaghat ################################################### diag(P1) ################################################### ### chunk number 36: diaghat1 ################################################### rowSums(Q1^2) ################################################### ### chunk number 37: diaghat1 ################################################### hatvalues(lm1) ################################################### ### chunk number 38: rankX ################################################### lm1$rank qrlm1$rank ################################################### ### chunk number 39: cholcrossprod ################################################### chol(crossprod(X)) ################################################### ### chunk number 40: svdx ################################################### str(Xsv <- svd(X)) ################################################### ### chunk number 41: svdreconstruct ################################################### Xsv$u %*% (Xsv$d * t(Xsv$v)) all.equal(Xsv$u %*% (Xsv$d * t(Xsv$v)), X, check.attr=FALSE) ################################################### ### chunk number 42: UVorth ################################################### zapsmall(crossprod(Xsv$u)) zapsmall(crossprod(Xsv$v)) ################################################### ### chunk number 43: eigenXtX ################################################### str(ev <- eigen(crossprod(X))) Xsv$d^2 all.equal(ev$values, Xsv$d^2) ev$vectors Xsv$v all.equal(-Xsv$v, ev$vectors) ################################################### ### chunk number 44: kappa ################################################### svd(Q, nu=0, nv=0)$d kappa(Q) svd(Q1, nu=0, nv=0)$d kappa(Q1) kappa(Xsv$u) ################################################### ### chunk number 45: kappaX ################################################### Xsv$d (kappaX <- Xsv$d[1]/Xsv$d[length(Xsv$d)]) ################################################### ### chunk number 46: kappaX ################################################### kappa(X) kappa(X, exact=TRUE) ################################################### ### chunk number 47: eps ################################################### .Machine$double.eps ################################################### ### chunk number 48: circ ################################################### str(rad <- seq(0, 2*pi, len=201)) str(circ <- rbind(cos(rad), sin(rad))) ################################################### ### chunk number 49: circX ################################################### fits <- X %*% circ ################################################### ### chunk number 50: ellipses ################################################### library(lattice) twoD <- rbind(data.frame(t(crossprod(Xsv$u, fits))), data.frame(t(crossprod(Q1, fits)))) twoD$mat <- factor(rep(c("U","Q"), each = ncol(fits)), levels=c("U","Q")) print(xyplot(X2 ~ X1|mat, twoD,, type=c("g","l"), xlab="", ylab="", aspect="iso", layout=c(2,1)))