load(file="colestiramina.RDa")
## In this situation, we clearly have paired observation. We study 164
## subjects before and after a treatment on the cholesterol level.
# Question 1 --------------------------------------------------------------
## We create D as the difference between z and y
D<-cholost$z-cholost$y
cholost<-data.frame(cholost,D)
## If D>0, it means a worsening of the cholesterol
# COMMENT FROM THE TEACHER ---> Inadequate problem statement, this will be the approach
# to perform a means comparison test in a paired data setting, not to answer this
# exercise
set.seed(080615)
B<-9999
# perform a one-sample randomization test on D
# for the null hypothesis H0: cor > 0 vs H1 cor <= 0
# here the test statistic is the corelation between before and after the treatment
# COMMENT FROM THE TEACHER ---> Again the same objection, this is a mixture of the
# permutation significance test for the means in paired data (with errors) with a correlation test.
listcor <- replicate(B,cor(cholost$z,(rbinom(length(D),1,.5)*2-1)*D,
method="pearson"))
p.value<-sum(listcor>0)/length(listcor)
p.value #p-value=0.4917492
## Here cor>0, the hypotesis is confirmed
# COMMENT FROM THE TEACHER ---> basic statistical error: this (wrong) p-value
# would not confirm the hypothesis cor > 0
# Question 2 --------------------------------------------------------------
B <- 9999
n<-length(D)
alpha <- 0.05
confLevel = 1 - alpha
## B values of r* the correlation coefficient for each permutation
r.boot <- replicate(B, cor(cholost[sample(1:n, replace=TRUE),])[3,1])
# COMMENT FROM THE TEACHER ---> Again, incorrect initial setting
## Confidence interval percentil boostrap 95%
icBoot.perc = quantile(r.boot, probs = c(alpha/2, 1 - alpha/2))
names(icBoot.perc) = ""
attr(icBoot.perc, "conf.level") = confLevel
icBoot.perc # [0.4511519 ; 0.6393776]
## The confidence interval boostrap percentil for the correlation coefficient
## is [0.4511519 ; 0.6393776]