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
# 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
####################################################################
## 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
####################################################################
## 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
####################################################################
## 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