rm(list =ls())
#Part 1, verify the result in textbook page 210.
dat1 <- data.frame(list(
O = c(-1,1,-1,1,-1,1,-1,1),
H = c(-1,-1,1,1,-1,-1,1,1),
C = c(-1,-1,-1,-1,1,1,1,1),
y = c(5.9,4.0,3.9,1.2,5.3,4.8,6.3,0.8)
))
lm1 <- lm(y~O*H*C,data = dat1)
summary(lm1)
##
## Call:
## lm(formula = y ~ O * H * C, data = dat1)
##
## Residuals:
## ALL 8 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.025 NA NA NA
## O -1.325 NA NA NA
## H -0.975 NA NA NA
## C 0.275 NA NA NA
## O:H -0.725 NA NA NA
## O:C -0.175 NA NA NA
## H:C 0.225 NA NA NA
## O:H:C -0.525 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 7 and 0 DF, p-value: NA
## orthogonal desgin
X <- model.matrix(lm1)
t(X)%*%X
## (Intercept) O H C O:H O:C H:C O:H:C
## (Intercept) 8 0 0 0 0 0 0 0
## O 0 8 0 0 0 0 0 0
## H 0 0 8 0 0 0 0 0
## C 0 0 0 8 0 0 0 0
## O:H 0 0 0 0 8 0 0 0
## O:C 0 0 0 0 0 8 0 0
## H:C 0 0 0 0 0 0 8 0
## O:H:C 0 0 0 0 0 0 0 8
solve(t(X)%*%X)%*%t(X)%*%dat1$y
## [,1]
## (Intercept) 4.025
## O -1.325
## H -0.975
## C 0.275
## O:H -0.725
## O:C -0.175
## H:C 0.225
## O:H:C -0.525
###################################################
#relationship with respect to effects in textbook
####################################################
#1 , main effect, interaction effects estimates in textbook are twice the
# number returned by the fitting
# mean = 4.025 , O = 2.65, H = -1.95 ...
#2, standard deviation from effects in text book doubled, too.
## Comments : better run through hand calculation
lm2 <- lm(y~O+H+C,data = dat1)
summary(lm2)
##
## Call:
## lm(formula = y ~ O + H + C, data = dat1)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -0.15 0.60 -0.20 -0.25 -1.30 0.85 1.65 -1.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0250 0.4697 8.569 0.00102 **
## O -1.3250 0.4697 -2.821 0.04778 *
## H -0.9750 0.4697 -2.076 0.10653
## C 0.2750 0.4697 0.585 0.58967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.329 on 4 degrees of freedom
## Multiple R-squared: 0.7592, Adjusted R-squared: 0.5785
## F-statistic: 4.203 on 3 and 4 DF, p-value: 0.09958
X <- model.matrix(lm2)
t(X)%*%X
## (Intercept) O H C
## (Intercept) 8 0 0 0
## O 0 8 0 0
## H 0 0 8 0
## C 0 0 0 8
solve(t(X)%*%X)%*%t(X)%*%dat1$y
## [,1]
## (Intercept) 4.025
## O -1.325
## H -0.975
## C 0.275
## intesting observation, all std are the same.
## because t(X)%*%X = 8 * identity
#Part 2,
library (BHH2) ## the library with our text book.
library(xtable) ## can help create table from R for latex to use.
xtable(dat1)
## % latex table generated in R 3.2.2 by xtable 1.7-4 package
## % Wed Feb 24 18:30:09 2016
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & O & H & C & y \\
## \hline
## 1 & -1.00 & -1.00 & -1.00 & 5.90 \\
## 2 & 1.00 & -1.00 & -1.00 & 4.00 \\
## 3 & -1.00 & 1.00 & -1.00 & 3.90 \\
## 4 & 1.00 & 1.00 & -1.00 & 1.20 \\
## 5 & -1.00 & -1.00 & 1.00 & 5.30 \\
## 6 & 1.00 & -1.00 & 1.00 & 4.80 \\
## 7 & -1.00 & 1.00 & 1.00 & 6.30 \\
## 8 & 1.00 & 1.00 & 1.00 & 0.80 \\
## \hline
## \end{tabular}
## \end{table}