STAT226 Multivariate Analysis Assignment

Generate the respective confidence intervals for the Anatomical.csv data using the following methods:
1. Roy-Bose Simultaneous Confidence Intervals
2. Bonferroni Method of Multiple Comparisons

Data Preparation

# Loading the original data
anatomicalraw <- read.csv("anatomical.csv")
summary(anatomicalraw)
##        id          enumerator     chest           waist      
##  Min.   :  1   Atanante1: 1   Min.   :25.25   Min.   :23.00  
##  1st Qu.: 26   Atanante2: 1   1st Qu.:32.00   1st Qu.:27.00  
##  Median : 51   Atanante3: 1   Median :33.00   Median :28.50  
##  Mean   : 51   Atanante4: 1   Mean   :33.71   Mean   :29.28  
##  3rd Qu.: 76   Barles1  : 1   3rd Qu.:36.00   3rd Qu.:31.00  
##  Max.   :101   Barles2  : 1   Max.   :42.00   Max.   :43.00  
##                (Other)  :95                                  
##      wrist            head            hand          forearm     
##  Min.   :4.500   Min.   :19.50   Min.   :6.000   Min.   :13.90  
##  1st Qu.:5.500   1st Qu.:21.70   1st Qu.:6.500   1st Qu.:15.50  
##  Median :6.000   Median :22.00   Median :7.000   Median :16.00  
##  Mean   :6.017   Mean   :22.27   Mean   :7.005   Mean   :16.33  
##  3rd Qu.:6.500   3rd Qu.:23.00   3rd Qu.:7.300   3rd Qu.:17.00  
##  Max.   :8.300   Max.   :25.00   Max.   :8.500   Max.   :20.00  
##                                                                 
##      height          age             sex        
##  Min.   :53.0   Min.   :10.00   Min.   :0.0000  
##  1st Qu.:61.0   1st Qu.:18.00   1st Qu.:0.0000  
##  Median :63.0   Median :19.00   Median :0.0000  
##  Mean   :63.4   Mean   :20.73   Mean   :0.3663  
##  3rd Qu.:65.5   3rd Qu.:20.00   3rd Qu.:1.0000  
##  Max.   :73.0   Max.   :57.00   Max.   :1.0000  
## 
# Subsetting data to only include the necessary columns
anatomical <- anatomicalraw[,colnames(anatomicalraw) %in% c("chest","waist","wrist","head","hand","forearm")]
summary(anatomical)
##      chest           waist           wrist            head      
##  Min.   :25.25   Min.   :23.00   Min.   :4.500   Min.   :19.50  
##  1st Qu.:32.00   1st Qu.:27.00   1st Qu.:5.500   1st Qu.:21.70  
##  Median :33.00   Median :28.50   Median :6.000   Median :22.00  
##  Mean   :33.71   Mean   :29.28   Mean   :6.017   Mean   :22.27  
##  3rd Qu.:36.00   3rd Qu.:31.00   3rd Qu.:6.500   3rd Qu.:23.00  
##  Max.   :42.00   Max.   :43.00   Max.   :8.300   Max.   :25.00  
##       hand          forearm     
##  Min.   :6.000   Min.   :13.90  
##  1st Qu.:6.500   1st Qu.:15.50  
##  Median :7.000   Median :16.00  
##  Mean   :7.005   Mean   :16.33  
##  3rd Qu.:7.300   3rd Qu.:17.00  
##  Max.   :8.500   Max.   :20.00
# Generating the hypothesized mean vector
(anatomical_hypmeans <- c(33,28,6,21,7,16))
## [1] 33 28  6 21  7 16
(anatomical_hypmeans <- as.data.frame(anatomical_hypmeans))
##   anatomical_hypmeans
## 1                  33
## 2                  28
## 3                   6
## 4                  21
## 5                   7
## 6                  16
row.names(anatomical_hypmeans) <- c("chest","waist","wrist","head","hand","forearm")

# Obtaining the values for n (number of observations) and p(number of parameters)
(n <- nrow(anatomical))
## [1] 101
(p <- ncol(anatomical))
## [1] 6
# Generating the sample mean vector
(anatomical_means <- colMeans(anatomical))
##     chest     waist     wrist      head      hand   forearm 
## 33.710396 29.284653  6.016832 22.270792  7.004950 16.327228
(anatomical_means <- as.data.frame(anatomical_means))
##         anatomical_means
## chest          33.710396
## waist          29.284653
## wrist           6.016832
## head           22.270792
## hand            7.004950
## forearm        16.327228
# Generating the sample variance vector
(anatomical_var <- sapply(anatomical,function(x)var(x)))
##      chest      waist      wrist       head       hand    forearm 
##  9.2715156 12.9706871  0.4368139  1.0929133  0.2661253  1.5670762
(anatomical_var <- as.data.frame(anatomical_var))
##         anatomical_var
## chest        9.2715156
## waist       12.9706871
## wrist        0.4368139
## head         1.0929133
## hand         0.2661253
## forearm      1.5670762

Hotelling T^2 Test Implementation

####################################################################
## Hotellings's T^2 Test
####################################################################

# Using manually generated function
T2Test <- function(nxpdata, hypmeans=0){
  nxpdata <- as.matrix(nxpdata)
  n <- nrow(nxpdata)
  p <- ncol(nxpdata)
  df2 <- n-p
  if(df2 < 1L) stop ("Number of rows should be greater than the number of columns!")
  if(length(hypmeans)!= p) hypmeans<-rep(hypmeans[1],p)
  samplemean <- colMeans(nxpdata)
  samplevar <- cov(nxpdata)
  T2 <- n * t(samplemean - hypmeans) %*% solve(samplevar) %*% (samplemean-hypmeans)
  Fstat <- ((n-p)*T2)/(p*(n-1))
  pval <- 1-pf(Fstat,df1=p,df2=n-p)
  return(data.frame(T2=as.numeric(T2), 
                    Fstat=as.numeric(Fstat),
                    df1=p, 
                    df2=n-p, 
                    p.value=formatC((as.numeric(pval)),format="e",digits=2), 
                    row.names=""))
}

(anatomical_T2TestRFunction <- T2Test(anatomical,hypmeans = c(33,28,6,21,7,16)))
##        T2    Fstat df1 df2  p.value
##  233.9452 37.04133   6  95 0.00e+00

Generation of the Roy-Bose Simultaneous Confidence Intervals

####################################################################
## T^2 Post-Hoc Test 1 : Roy-Bose Simultaneous Confidence Intervals
####################################################################

# Obtaining the Hotelling T-Squared value
# Hotelling T-Squared can be expressed as ((n-1)*p/(n-p)*F(alpha,df1,df2)
(hotellingtsquared <- ((n-1)*p/(n-p)*qf(.95, df1=p, df2=n-p)))
## [1] 13.86662
# Generating the function to create the Roy-Bose simultaneous confidence intervals
roybose_confint <- function(n,samplemean,samplevar,hotellingtsquared) {
  roybose_lowerci <- samplemean - sqrt((1/n)*samplevar*hotellingtsquared)
  roybose_upperci <- samplemean + sqrt((1/n)*samplevar*hotellingtsquared)
  roybose_ci <- c(roybose_lowerci,roybose_upperci)
  return(roybose_ci)
}

# Generating the Roy-Bose confidence intervals 
anatomical_roybose_ci <- roybose_confint(n,anatomical_means,anatomical_var,hotellingtsquared)
anatomical_roybose_ci <- as.data.frame(anatomical_roybose_ci)
colnames(anatomical_roybose_ci)<-c("anatomical_roybose_lowerci","anatomical_roybose_upperci")
anatomical_roybose_ci
##   anatomical_roybose_lowerci anatomical_roybose_upperci
## 1                  32.582160                  34.838632
## 2                  27.950192                  30.619115
## 3                   5.771941                   6.261723
## 4                  21.883430                  22.658154
## 5                   6.813803                   7.196098
## 6                  15.863386                  16.791069
# Creating the summary
anatomical_roybose_summary <- cbind(anatomical_hypmeans,anatomical_means,anatomical_var,anatomical_roybose_ci)
colnames(anatomical_roybose_summary) <-c("HypothesizedMean",
                                         "SampleMean",
                                         "SampleVariance",
                                         "RoyBoseLowerCI",
                                         "RoyBoseUpperCI")
anatomical_roybose_summary <- sapply(anatomical_roybose_summary,as.numeric)
anatomical_roybose_summary <- as.data.frame(anatomical_roybose_summary)
anatomical_roybose_summary <- format(round(anatomical_roybose_summary,2),nsmall=2)
anatomical_roybose_summary$Decision <- ifelse(
  anatomical_roybose_summary$HypothesizedMean<anatomical_roybose_summary$RoyBoseLowerCI | 
  anatomical_roybose_summary$HypothesizedMean>anatomical_roybose_summary$RoyBoseUpperCI, 
  "Hypothesized Mean Beyond CI",
  "Hypothesized Mean Within CI")
row.names(anatomical_roybose_summary) <- c("chest","waist","wrist","head","hand","forearm")


# Viewing the summary
anatomical_roybose_summary 
##         HypothesizedMean SampleMean SampleVariance RoyBoseLowerCI
## chest              33.00      33.71           9.27          32.58
## waist              28.00      29.28          12.97          27.95
## wrist               6.00       6.02           0.44           5.77
## head               21.00      22.27           1.09          21.88
## hand                7.00       7.00           0.27           6.81
## forearm            16.00      16.33           1.57          15.86
##         RoyBoseUpperCI                    Decision
## chest            34.84 Hypothesized Mean Within CI
## waist            30.62 Hypothesized Mean Within CI
## wrist             6.26 Hypothesized Mean Within CI
## head             22.66 Hypothesized Mean Beyond CI
## hand              7.20 Hypothesized Mean Within CI
## forearm          16.79 Hypothesized Mean Within CI

Generation of the Bonferroni Method of Multiple Comparisons

####################################################################
## T^2 Post-Hoc Test 2 : Bonferroni Method of Multiple Comparisons
####################################################################

# Setting the alpha
(alpha <- (0.05/6)/2)
## [1] 0.004166667
# Obtaining the t distribution value
(tvalue <- qt(1-alpha, df=n-1))
## [1] 2.691755
# Generating the function to create the Bonferroni confidence intervals
bonferroni_confint <- function(n,samplemean,samplevar,tvalue) {
  bonferroni_lowerci <- samplemean - sqrt((1/n)*samplevar)*tvalue
  bonferroni_upperci <- samplemean + sqrt((1/n)*samplevar)*tvalue
  bonferroni_ci <- c(bonferroni_lowerci,bonferroni_upperci)
  return(bonferroni_ci)
}

# Generating the Bonferroni confidence intervals 
anatomical_bonferroni_ci <- bonferroni_confint(n,anatomical_means,anatomical_var,tvalue)
anatomical_bonferroni_ci <- as.data.frame(anatomical_bonferroni_ci)
colnames(anatomical_bonferroni_ci)<-c("anatomical_bonferroni_lowerci","anatomical_bonferroni_upperci")
anatomical_bonferroni_ci
##   anatomical_bonferroni_lowerci anatomical_bonferroni_upperci
## 1                     32.894847                     34.525945
## 2                     28.320033                     30.249274
## 3                      5.839811                      6.193852
## 4                     21.990786                     22.550798
## 5                      6.866779                      7.143122
## 6                     15.991938                     16.662517
# Creating the summary
anatomical_bonferroni_summary <- cbind(anatomical_hypmeans,anatomical_means,anatomical_var,anatomical_bonferroni_ci)
colnames(anatomical_bonferroni_summary) <-c("HypothesizedMean",
                                            "SampleMean",
                                            "SampleVariance",
                                            "BonferroniLowerCI",
                                            "BonferroniUpperCI")
anatomical_bonferroni_summary <- sapply(anatomical_bonferroni_summary,as.numeric)
anatomical_bonferroni_summary <- as.data.frame(anatomical_bonferroni_summary)
anatomical_bonferroni_summary <- format(round(anatomical_bonferroni_summary,2),nsmall=2)
anatomical_bonferroni_summary$Decision <- ifelse(
  anatomical_bonferroni_summary$HypothesizedMean<anatomical_bonferroni_summary$BonferroniLowerCI | 
  anatomical_bonferroni_summary$HypothesizedMean>anatomical_bonferroni_summary$BonferroniUpperCI, 
  "Hypothesized Mean Beyond CI",
  "Hypothesized Mean Within CI")
row.names(anatomical_bonferroni_summary) <- c("chest","waist","wrist","head","hand","forearm")

# Viewing the summary
anatomical_bonferroni_summary
##         HypothesizedMean SampleMean SampleVariance BonferroniLowerCI
## chest              33.00      33.71           9.27             32.89
## waist              28.00      29.28          12.97             28.32
## wrist               6.00       6.02           0.44              5.84
## head               21.00      22.27           1.09             21.99
## hand                7.00       7.00           0.27              6.87
## forearm            16.00      16.33           1.57             15.99
##         BonferroniUpperCI                    Decision
## chest               34.53 Hypothesized Mean Within CI
## waist               30.25 Hypothesized Mean Beyond CI
## wrist                6.19 Hypothesized Mean Within CI
## head                22.55 Hypothesized Mean Beyond CI
## hand                 7.14 Hypothesized Mean Within CI
## forearm             16.66 Hypothesized Mean Within CI