##### Multivariate analysis of Repeated measurements
rm(list=(ls()))
trt<- 4
aa<- matrix(c(30.9,31.9,31.3,32.1,30.9,31.3,31.3,32.1,30.3,32.2,30.7,31.6,31.1,31,31.2,31.7,31.8,33,30.9,32.1,30.9,31.6,31,31.7,30.5,31.4,31.8,31.7,30.8,32.2, 30.9,31.7,31.3,31.3,30.8,31.2,31.7,31.5,30.6,32.4),10,4,byrow = F);aa
## [,1] [,2] [,3] [,4]
## [1,] 30.9 30.7 30.9 30.9
## [2,] 31.9 31.6 31.6 31.7
## [3,] 31.3 31.1 31.0 31.3
## [4,] 32.1 31.0 31.7 31.3
## [5,] 30.9 31.2 30.5 30.8
## [6,] 31.3 31.7 31.4 31.2
## [7,] 31.3 31.8 31.8 31.7
## [8,] 32.1 33.0 31.7 31.5
## [9,] 30.3 30.9 30.8 30.6
## [10,] 32.2 32.1 32.2 32.4
bb <- matrix(c(31.5,31.2,31.3,30.4,30.7,29.8,31.4,30.9,31.1,31.3,30.6,31.2,31.3,30.8,30.9,30.8,32,32.4,31.3,31.5,30.8, 31.1,31.5,30.4,30.9,30.9,31.7,31.8,31.2,31.6,31,31.3,31.4,30.2,30.9,30.8,31.6,31.9,31.2,31.7),10,4,byrow=F);bb
## [,1] [,2] [,3] [,4]
## [1,] 31.5 30.6 30.8 31.0
## [2,] 31.2 31.2 31.1 31.3
## [3,] 31.3 31.3 31.5 31.4
## [4,] 30.4 30.8 30.4 30.2
## [5,] 30.7 30.9 30.9 30.9
## [6,] 29.8 30.8 30.9 30.8
## [7,] 31.4 32.0 31.7 31.6
## [8,] 30.9 32.4 31.8 31.9
## [9,] 31.1 31.3 31.2 31.2
## [10,] 31.3 31.5 31.6 31.7
naa <- nrow(aa);naa
## [1] 10
nbb <- nrow(bb);nbb
## [1] 10
xbar1s <- colMeans(aa);xbar1s
## [1] 31.43 31.51 31.36 31.34
xbar2s <- colMeans(bb);xbar1s
## [1] 31.43 31.51 31.36 31.34
c <- matrix(c(-3,1,-1,-1,-1,3,1,-1,-3,3,1,1),3,4);c
## [,1] [,2] [,3] [,4]
## [1,] -3 -1 1 3
## [2,] 1 -1 -1 1
## [3,] -1 3 -3 1
Ybar1s <- c%*%xbar1s;Ybar1s
## [,1]
## [1,] -0.42
## [2,] -0.10
## [3,] 0.36
Ybar2s <- c%*%xbar2s;Ybar2s
## [,1]
## [1,] 0.63
## [2,] -0.31
## [3,] 0.51
poold <- ((naa-1)*var(aa)+(nbb-1)*var(bb))/(naa+nbb-2);poold
## [,1] [,2] [,3] [,4]
## [1,] 0.3402778 0.1832778 0.1993333 0.2104444
## [2,] 0.1832778 0.3969444 0.2312222 0.2270000
## [3,] 0.1993333 0.2312222 0.2462778 0.2347778
## [4,] 0.2104444 0.2270000 0.2347778 0.2613333
round(poold,4)
## [,1] [,2] [,3] [,4]
## [1,] 0.3403 0.1833 0.1993 0.2104
## [2,] 0.1833 0.3969 0.2312 0.2270
## [3,] 0.1993 0.2312 0.2463 0.2348
## [4,] 0.2104 0.2270 0.2348 0.2613
T2<- ((naa*nbb)/(naa+nbb))*(t(xbar1s-xbar2s)%*%solve(poold)%*%(xbar1s-xbar2s));T2
## [,1]
## [1,] 4.410313
one <- round(c%*%poold%*%t(c),4);one
## [,1] [,2] [,3]
## [1,] 1.7576 -0.2998 0.0661
## [2,] -0.2998 0.4394 -0.6816
## [3,] 0.0661 -0.6816 1.8574
two <- round(t(one),4);two
## [,1] [,2] [,3]
## [1,] 1.7576 -0.2998 0.0661
## [2,] -0.2998 0.4394 -0.6816
## [3,] 0.0661 -0.6816 1.8574
Fcal <-T2*((naa+nbb-trt)/((trt-1)*(naa+nbb-2)));Fcal
## [,1]
## [1,] 1.306759
qf(0.05,(trt-1),(naa+nbb-trt),lower.tail = F)
## [1] 3.238872