Example 6.7

df = read.table("http://users.stat.umn.edu/~sandy/courses/8053/Data/Wichern_data/T6-4.dat")
Factor1 = df$V1
Factor1[df$V1 == 0] <-"Low"
Factor1[df$V1 == 1] <- "High"
Factor2 = df$V2
Factor2[df$V2 == 0] <-"Low"
Factor2[df$V2 == 1] <- "High"
film = data.frame(df$V3,df$V4,df$V5,Factor1,Factor2)
names(film)[names(film) == "df.V3"] <- "X1"
names(film)[names(film) == "df.V4"] <- "X2"
names(film)[names(film) == "df.V5"] <- "X3"
attach(film)
film
model1 = lm(formula = X1~Factor1+Factor2+Factor1*Factor2,data=film,
    contrasts = list(Factor1 = contr.sum,Factor2 = contr.sum))
anova(model1)
model2 = lm(formula = X2~Factor1+Factor2+Factor1*Factor2,data=film,
    contrasts = list(Factor1 = contr.sum,Factor2 = contr.sum))
anova(model2)
model3 = lm(formula = X3~Factor1+Factor2+Factor1*Factor2,data=film,
    contrasts = list(Factor1 = contr.sum,Factor2 = contr.sum))
anova(model3)
DV = cbind(X1,X2,X3)
output = lm(DV~ Factor1*Factor2,data = film,
            contrasts = list(Factor1 = contr.sum,Factor2 = contr.sum))
library(car)
manova_out = Manova(output,type = "III")
summary(manova_out,multivariate = TRUE)
## 
## Type III MANOVA Tests:
## 
## Sum of squares and products for error:
##        X1     X2     X3
## X1  1.764  0.020 -3.070
## X2  0.020  2.628 -0.552
## X3 -3.070 -0.552 64.924
## 
## ------------------------------------------
##  
## Term: (Intercept) 
## 
## Sum of squares and products for the hypothesis:
##           X1        X2       X3
## X1  920.7245 1264.0455 533.9795
## X2 1264.0455 1735.3845 733.0905
## X3  533.9795  733.0905 309.6845
## 
## Multivariate Tests: (Intercept)
##                  Df test stat approx F num Df den Df     Pr(>F)    
## Pillai            1    0.9992 5950.906      3     14 < 2.22e-16 ***
## Wilks             1    0.0008 5950.906      3     14 < 2.22e-16 ***
## Hotelling-Lawley  1 1275.1941 5950.906      3     14 < 2.22e-16 ***
## Roy               1 1275.1941 5950.906      3     14 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: Factor1 
## 
## Sum of squares and products for the hypothesis:
##         X1      X2      X3
## X1  1.7405 -1.5045  0.8555
## X2 -1.5045  1.3005 -0.7395
## X3  0.8555 -0.7395  0.4205
## 
## Multivariate Tests: Factor1
##                  Df test stat approx F num Df den Df   Pr(>F)   
## Pillai            1 0.6181416 7.554269      3     14 0.003034 **
## Wilks             1 0.3818584 7.554269      3     14 0.003034 **
## Hotelling-Lawley  1 1.6187719 7.554269      3     14 0.003034 **
## Roy               1 1.6187719 7.554269      3     14 0.003034 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: Factor2 
## 
## Sum of squares and products for the hypothesis:
##        X1     X2     X3
## X1 0.7605 0.6825 1.9305
## X2 0.6825 0.6125 1.7325
## X3 1.9305 1.7325 4.9005
## 
## Multivariate Tests: Factor2
##                  Df test stat approx F num Df den Df   Pr(>F)  
## Pillai            1 0.4769651 4.255619      3     14 0.024745 *
## Wilks             1 0.5230349 4.255619      3     14 0.024745 *
## Hotelling-Lawley  1 0.9119183 4.255619      3     14 0.024745 *
## Roy               1 0.9119183 4.255619      3     14 0.024745 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: Factor1:Factor2 
## 
## Sum of squares and products for the hypothesis:
##        X1     X2     X3
## X1 0.0005 0.0165 0.0445
## X2 0.0165 0.5445 1.4685
## X3 0.0445 1.4685 3.9605
## 
## Multivariate Tests: Factor1:Factor2
##                  Df test stat approx F num Df den Df  Pr(>F)
## Pillai            1 0.2228942 1.338522      3     14 0.30178
## Wilks             1 0.7771058 1.338522      3     14 0.30178
## Hotelling-Lawley  1 0.2868261 1.338522      3     14 0.30178
## Roy               1 0.2868261 1.338522      3     14 0.30178

Example 6.15

\(\hat{\beta_l} = (B'S^{-1}_{pooled}B)^{-1}B'S^{-1}_{pooled}\overline{X_l}\), for l = 1,2,…,g

where \(S_{pooled} = \dfrac{1}{(N-g)}((n_1-1)S_1 + ... +(n_g-1)S_g) = \dfrac{1}{N-g}W\)

\(N = \sum\limits_{l=1}^{g}n_l\).

control = read.table("http://users.stat.umn.edu/~sandy/courses/8053/Data/Wichern_data/T6-5.dat")
head(control)
treatment =read.table("http://users.stat.umn.edu/~sandy/courses/8053/Data/Wichern_data/T6-6.dat")
head(treatment)
B = matrix(c(1,0,0,1,1,1,1,2,4,1,3,9),nrow = 4,byrow = T)
B
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    1    1    1
## [3,]    1    2    4
## [4,]    1    3    9
S_1 = cov(control)
S_1
##          V1       V2       V3       V4
## V1 92.11886 86.11057 73.36229 74.58900
## V2 86.11057 89.07638 72.95552 71.77276
## V3 73.36229 72.95552 71.89067 63.59176
## V4 74.58900 71.77276 63.59176 75.44410
S_2 = cov(treatment)
S_2
##          V1        V2       V3       V4
## V1 98.17450  97.01342 89.48242 86.11108
## V2 97.01342 100.59596 88.14246 88.20946
## V3 89.48242  88.14246 86.34963 80.55062
## V4 86.11108  88.20946 80.55062 81.41563
n1 = dim(control)[1] 
n2 = dim(treatment)[1]
N = n1+n2
g = 2
S_pooled = 1/(N-g)*((n1-1)*S_1 + (n2-1)*S_2)
S_pooled
##          V1       V2       V3       V4
## V1 95.25109 91.74997 81.70028 80.54870
## V2 91.74997 95.03478 80.81083 80.27450
## V3 81.70028 80.81083 79.36944 72.36359
## V4 80.54870 80.27450 72.36359 78.53282
S_pooled_inv = solve(S_pooled)
S_pooled_inv
##             V1          V2          V3          V4
## V1  0.19574037 -0.11339613 -0.05422005 -0.03489314
## V2 -0.11339613  0.17332679 -0.02662306 -0.03633220
## V3 -0.05422005 -0.02662306  0.12510876 -0.03245537
## V4 -0.03489314 -0.03633220 -0.03245537  0.11556612
A = round(solve(t(B)%*%S_pooled_inv%*%B),digits = 4)
A
##         [,1]    [,2]    [,3]
## [1,] 93.1744 -5.8368  0.2184
## [2,] -5.8368  9.5699 -3.0240
## [3,]  0.2184 -3.0240  1.1051
mean_control = matrix(c(72.38,73.29,72.47,64.79),nrow = 4)
mean_treatment = matrix(c(69.29,70.66,71.18,64.53),nrow = 4)
beta_1_hat = round(A%*%t(B)%*%S_pooled_inv%*%mean_control,digits = 4)
beta_1_hat
##         [,1]
## [1,] 73.0699
## [2,]  3.6381
## [3,] -2.0248
beta_2_hat = A%*%t(B)%*%S_pooled_inv%*%mean_treatment
beta_2_hat
##           [,1]
## [1,] 70.140048
## [2,]  4.090316
## [3,] -1.853808
t = c(0:10)
control_group= beta_1_hat[1] + beta_1_hat[2]*t + beta_1_hat[2]*t^2
treatment_group=  beta_2_hat[1] + beta_2_hat[2]*t + beta_2_hat[2]*t^2
library(reshape2)
data = data.frame(t,control_group,treatment_group)
df = melt(data,id.vars = 't')
library(ggplot2)
ggplot(data = df,aes(x=t,y=value,colour = variable))+
  geom_line()+theme_bw()

Example 6.16

x1 = c(5.0,4.5,6.0,6.0,6.2,6.9,6.8,5.3,6.6,7.3,4.6,4.9,4.0,3.8,6.2,5.0,5.3,7.1,5.8,6.8)
x2 = c(3.0,3.2,3.5,4.6,5.6,5.2,6.0,5.5,7.3,6.5,4.9,5.9,4.1,5.4,6.1,7.0,4.7,6.6,7.8,8.0)
Group = c(rep(1,10),rep(2,10))
 
data = data.frame(x1,x2,Group)
data
av1 = aov(data$x1 ~ data$Group)
summary(av1)
##             Df Sum Sq Mean Sq F value Pr(>F)
## data$Group   1  2.521   2.521   2.459  0.134
## Residuals   18 18.449   1.025
av2 =aov(data$x2 ~ data$Group)
summary(av2)
##             Df Sum Sq Mean Sq F value Pr(>F)
## data$Group   1   5.10   5.100   2.678  0.119
## Residuals   18  34.29   1.905

\(T^2 = [\overline{x_1}- \overline{x_2} - (\mu_1-\mu_2)]'\left[\dfrac{1}{n_1}S_1 + \dfrac{1}{n_2}S_2\right]^{-1}[\overline{x_1}- \overline{x_2} - (\mu_1-\mu_2)] > \dfrac{(n_1+n_2-2)p}{(n_1+n_2-p-1)}F_{p,n_1+n_2-p-1}(\alpha)\)

Group1 = subset(data,data$Group == 1)
S1 = cov(cbind(Group1$x1,Group1$x2))
Group2 = subset(data,data$Group == 2)
S2 = cov(cbind(Group2$x1,Group2$x2))
p = 2
n1 = dim(Group1)[1]
n2 = dim(Group2)[1]
S_pooled = ((n1-1)/(n1+n2-2))*S1 + ((n2-1)/(n1+n2-2))*S2 
x1_bar = cbind(mean(Group1$x1),mean(Group1$x2))
x2_bar = cbind(mean(Group2$x1),mean(Group2$x2))
T_square = (x1_bar-x2_bar)%*%solve((1/n1 + 1/n2)*S_pooled)%*%t(x1_bar-x2_bar)
T_square
##          [,1]
## [1,] 17.28755
((n1+n2-2)*p)/(n1+n2-p-1)*qf(1-0.01,df1 =p,df2=n1+n2-p-1)
## [1] 12.9433

\(T^2\) = 17.29 > \(c^2\) = \(\dfrac{(18)(2)}{17}F_{2,17}(0.01) = 12.94\).

Example 6.17

data = read.table("http://users.stat.umn.edu/~sandy/courses/8053/Data/Wichern_data/T6-7.dat")
head(data)
log_Mass = log(data$V1)
log_SVL = log(data$V2)
Group = c(rep("C",20),rep("S",40))
log_data = data.frame(log_Mass,log_SVL,Group)
head(log_data)
library(ggplot2)
ggplot(data = log_data,aes(x=log_SVL,y=log_Mass))+
  geom_point(aes(col = Group))+theme_bw()

t.test(log_data$log_Mass~log_data$Group)
## 
##  Welch Two Sample t-test
## 
## data:  log_data$log_Mass by log_data$Group
## t = -0.73633, df = 44.817, p-value = 0.4654
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4789905  0.2225469
## sample estimates:
## mean in group C mean in group S 
##        2.239919        2.368140
t.test(log_data$log_SVL~log_data$Group)
## 
##  Welch Two Sample t-test
## 
## data:  log_data$log_SVL by log_data$Group
## t = 1.7766, df = 47.402, p-value = 0.08205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01140387  0.18407453
## sample estimates:
## mean in group C mean in group S 
##        4.394427        4.308091

\(T^2 = [\overline{x_1}- \overline{x_2} - (\mu_1-\mu_2)]'\left[\dfrac{1}{n_1}S_1 + \dfrac{1}{n_2}S_2\right]^{-1}[\overline{x_1}- \overline{x_2} - (\mu_1-\mu_2)] > \chi^2_{p}(\alpha)\)

Group_C = subset(log_data,log_data$Group == "C")
Group_S = subset(log_data,log_data$Group == "S")
S1 = round(cov(cbind(Group_C$log_Mass,Group_C$log_SVL)),digits = 5)
S2 = round(cov(cbind(Group_S$log_Mass,Group_S$log_SVL)),digits = 5)
n1 = dim(Group_C)[1]
n2 = dim(Group_S)[1]
mu1 = cbind(mean(Group_C$log_Mass),mean(Group_C$log_SVL))
mu2 = cbind(mean(Group_S$log_Mass),mean(Group_S$log_SVL))
T_square = (mu1-mu2)%*%solve((1/n1)*S1 + (1/n2)*S2)%*%t(mu1-mu2)
T_square
##          [,1]
## [1,] 225.7378
p=2
qchisq(1-0.05,df = p)
## [1] 5.991465

\(T^2\) = 225.7378 > \(c^2\) = \(\chi^{2}_{2}(0.05) = 5.99\)