pdf("rpd8_5.r") xy <- matrix(c( 1, -1, -1, -1, 53, 1, -1, -1, 1, 54, 1, -1, 1, -1, 40, 1, -1, 1, 1, 37, 1, 1, -1, -1, 84, 1, 1, -1, 1, 76, 1, 1, 1, -1, 40, 1, 1, 1, 1, 50, 1, 0, 0, 0, 50, 1, 1.215, 0, 0, 61, 1, -1.215, 0, 0, 54, 1, 0, 1.215, 0, 39, 1, 0, -1.215, 0, 67, 1, 0, 0, 1.215, 44, 1, 0, 0, -1.215, 61, 2, -1, -1, -1, 50, 2, -1, -1, 1, 42, 2, -1, 1, -1, 31, 2, -1, 1, 1, 28, 2, 1, -1, -1, 57, 2, 1, -1, 1, 78, 2, 1, 1, -1, 49, 2, 1, 1, 1, 54, 2, 0, 0, 0, 50, 2, 1.215, 0, 0, 76, 2, -1.215, 0, 0, 45, 2, 0, 1.215, 0, 33, 2, 0, -1.215, 0, 54, 2, 0, 0, 1.215, 45, 2, 0, 0, -1.215, 38, 3, -1.2500, -1.8867, -0.6350, 46, 3, 0.8600, -2.2200, -0.4250, 66, 3, 1.0000, -2.2400, -0.3100, 68, 3, 2.1165, -2.4167, -0.1450, 75, 3, 2.5825, -2.4900, -0.0800, 75, 3, 3.2475, -2.6667, 0.0800, 68, 3, 1.1760, -1.3333, 0, 78, 3, 1.4700, -1.6667, 0, 93, 3, 1.7640, -2.0000, 0, 96, 3, 2.0580, -2.3333, 0, 66), byrow=T, ncol=5) y <- xy[,5] x01 <- c(rep(1,15),rep(0,25)) x02 <- c(rep(0,15),rep(1,15),rep(0,10)) x03 <- c(rep(0,30),rep(1,10)) x1 <- xy[,2] x2 <- xy[,3] x3 <- xy[,4] x11 <- x1*x1 x22 <- x2*x2 x33 <- x3*x3 x12 <- x1*x2 x13 <- x1*x3 x23 <- x2*x3 xf <- cbind(x01,x02,x03,x1,x2,x3,x11,x22,x33,x12,x13,x23) betaf <- solve(t(xf) %*% xf) %*% t(xf) %*% y yhatf <- xf %*% betaf ssef <- t(y-yhatf) %*% (y-yhatf) dfef <- nrow(xf)-ncol(xf) s2f <- ssef[1,1]/dfef varbetaf <- s2f[1]*solve(t(xf) %*% xf) sebetaf <- sqrt(diag(varbetaf)) tbetaf <- betaf/sebetaf print(cbind(betaf,sebetaf,tbetaf)) print(cbind(s2f,dfef)) xr <- cbind(x01,x02,x03,x1,x2,x22) betar <- solve(t(xr) %*% xr) %*% t(xr) %*% y yhatr <- xr %*% betar sser <- t(y-yhatr) %*% (y-yhatr) dfer <- nrow(xr)-ncol(xr) s2r <- sser[1,1]/dfer varbetar <- s2r[1]*solve(t(xr) %*% xr) sebetar <- sqrt(diag(varbetar)) tbetar <- betar/sebetar print(cbind(betar,sebetar,tbetar)) print(cbind(s2r,dfer)) ffullred <- ((sser-ssef)/(dfer-dfef))/(ssef/dfef) ffrcrit <- qf(.95,dfer-dfef,dfef) probffr <- 1-pf(ffullred,dfer-dfef,dfef) print(cbind(ffullred,ffrcrit,probffr)) dev.off()