H:\public_html\sta6208\rpd>C:\R-2.9.0\bin\Rterm --vanilla R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > 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)) sebetaf x01 61.8439867 3.440772 17.97386986 x02 56.5106534 3.440772 16.42383008 x03 70.4079421 8.227833 8.55728856 x1 9.1271660 1.771973 5.15084855 x2 -9.8516000 1.854669 -5.31178429 x3 0.2632813 1.861929 0.14140246 x11 -1.2601941 1.464480 -0.86050616 x22 -6.4976734 2.014588 -3.22531071 x33 -2.9849162 2.952339 -1.01103449 x12 -0.9339081 1.510510 -0.61827327 x13 2.2419511 2.150452 1.04254858 x23 -0.1392224 2.137808 -0.06512388 > print(cbind(s2f,dfef)) s2f dfef [1,] 76.53377 28 > > > > 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)) sebetar x01 59.034845 2.429201 24.302170 x02 53.701512 2.429201 22.106660 x03 70.560136 7.603604 9.279828 x1 9.212551 1.474434 6.248197 x2 -9.820102 1.759983 -5.579658 x22 -6.895506 1.558674 -4.423957 > print(cbind(s2r,dfer)) s2r dfer [1,] 69.08662 34 > > > 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)) ffrcrit [1,] 0.4486023 2.445259 0.8397384 > > > dev.off() null device 1 > > > > > >