Upload R code

Final_TM_Code.R

+#### VE Model + parameters ####
+D1 <- function(z) {
+  1+z
+}
+
+D2 <- function(z) {
+  0
+}
+
+VE <- function(y,z,eps) {
+  D1(z)*y + D2(z) + eps
+}
+
+r=1e7
+n <- 4
+xbar<- 50
+s <- 5
+y0 <- 50/8.5
+alpha <- 8
+sig0 <- 1/2
+beta<- alpha * sig0^2
+
+#### JCGM 101 method ####
+ziJ <- runif(r,5,10)
+xiJ <- xbar + rt(r,n-1)*s/sqrt(n)
+yiJ <- (xiJ)/(1+ziJ)
+
+#### JCGM 101 + inf. method ####
+ziJi <- runif(r,5,10)
+xiJi <- xbar + rt(r,2*alpha + n-1)*sqrt((2*beta + (n-1)*s^2)/(n*(2*alpha + n-1)))
+yiJi <- (xiJi)/(1+ziJi)
+
+#### MCVE method ####
+yia <- c()
+i=1
+while (i<=r){
+  epsi <- rnorm(4,0,sig0)
+  zi <- runif(1,5,10)
+  xsimi <- c(VE(y0,zi,epsi[1]),VE(y0,zi,epsi[2]),VE(y0,zi,epsi[3]),VE(y0,zi,epsi[4]))
+  xvebar <- mean(xsimi)
+  ci=rchisq(1,df=(n-1))
+  mu <- y0*(zi+1)
+  Tia <- s * (sqrt(n-1))/(sig0*sqrt(ci))
+  yia[i] <- ((1+zi)^(-1)) * (Tia*(mu-xvebar)+(xbar-mu)) + y0
+  i=i+1
+}
+
+#### MCVE inf. method ####
+yib <- c()
+i=1
+while (i<=r){
+  epsi <- rnorm(4,0,sig0)
+  zi <- runif(1,5,10)
+  xsimi <- c(VE(y0,zi,epsi[1]),VE(y0,zi,epsi[2]),VE(y0,zi,epsi[3]),VE(y0,zi,epsi[4]))
+  xvebar <- mean(xsimi)
+  ci=rchisq(1,df=(2*alpha+n-1))
+  mu <- y0*(zi+1)
+  Tib <- sqrt(2*beta + (n-1)*s^2)/(sig0*sqrt(ci))
+  yib[i] <- ((1+zi)^(-1)) * (Tib*(mu-xvebar)+(xbar-mu)) + y0
+  i=i+1
+}
+
+
+#### Comparison ####
+## Mean
+c(mean(yiJ),mean(yiJi),mean(yia),mean(yib))
+
+## Standard Deviation
+c(sd(yiJ),sd(yiJi),sd(yia),sd(yib))
+
+## Shortest 95% Coverage Interval
+# JCGM 101 method
+SortedJ <- sort(yiJ)
+CovIntsJ <- c()
+i=1
+while (i<=(0.05*r)){
+  CovIntsJ[i]<-SortedJ[0.95*r+i]-SortedJ[i]
+  i=i+1
+}
+ShCovIJ <- c(SortedJ[which(CovIntsJ==min(CovIntsJ))],SortedJ[which(CovIntsJ==min(CovIntsJ))+(0.95*r)])
+# JCGM 101 + inf. method
+SortedJi <- sort(yiJi)
+CovIntsJi <- c()
+i=1
+while (i<=(0.05*r)){
+  CovIntsJi[i]<-SortedJi[0.95*r+i]-SortedJi[i]
+  i=i+1
+}
+ShCovIJi <- c(SortedJi[which(CovIntsJi==min(CovIntsJi))],SortedJi[which(CovIntsJi==min(CovIntsJi))+(0.95*r)])
+# MC-VE method
+Sorteda <- sort(yia)
+CovIntsa <- c()
+i=1
+while (i<=(0.05*r)){
+  CovIntsa[i]<-Sorteda[0.95*r+i]-Sorteda[i]
+  i=i+1
+}
+ShCovIa <- c(Sorteda[which(CovIntsa==min(CovIntsa))],Sorteda[which(CovIntsa==min(CovIntsa))+(0.95*r)])
+# MC-VE inf. method
+Sortedb <- sort(yib)
+CovIntsb <- c()
+i=1
+while (i<=(0.05*r)){
+  CovIntsb[i]<-Sortedb[0.95*r+i]-Sortedb[i]
+  i=i+1
+}
+ShCovIb <- c(Sortedb[which(CovIntsb==min(CovIntsb))],Sortedb[which(CovIntsb==min(CovIntsb))+(0.95*r)])
+
+ShCovIJ
+ShCovIJi
+ShCovIa
+ShCovIb
+
+## Coefficient of Variance
+sd(yiJ)/mean(yiJ)
+sd(yiJi)/mean(yiJi)
+sd(yia)/mean(yia)
+sd(yib)/mean(yib)
+
+#### Density Plot ####
+histJ <- hist(ifelse(yiJ>2 & yiJ<10,yiJ,NA),breaks=100,plot=FALSE)
+histJi <- hist(ifelse(yiJi>2 & yiJi<10,yiJi,NA),breaks=100,plot=FALSE)
+hista <- hist(ifelse(yia>2 & yia<10,yia,NA),breaks=100,plot=FALSE)
+histb <- hist(ifelse(yib>2 & yib<10,yib,NA),breaks=100,plot=FALSE)
+
+par(cex=1.5)
+plot(histJi$mids,histJi$density,type="l",lwd=4,xlab="Y",xlim=c(4,9),ylab="Probability Density",col=2)
+lines(histb$mids,histb$density,lwd=4,lty=2,col="yellow")
+lines(histJ$mids,histJ$density,lwd=4,col=1)
+lines(hista$mids,hista$density,lwd=4,lty=2,col=5)
+legend(x="topright", legend=c("JCGM 101","MC-VE", "JCGM 101 + inf.","MC-VE inf."),fill=c(1,5,2,"yellow"))
+
+#### Assessing different values of alpha ####
+## alpha >= 1
+k=1
+meantestJi <- c()
+meantestb <- c()
+sdtestJi <- c()
+sdtestb <- c()
+while (k<=15){
+  alpha <- k
+  sig0 <- 5
+  beta<- alpha * sig0^2
+  #### JCGM 101 + inf. method
+  ziJi <- runif(r,5,10)
+  xiJi <- xbar + rt(r,2*alpha + n-1)*sqrt((2*beta + (n-1)*s^2)/(n*(2*alpha + n-1)))
+  yiJi <- (xiJi)/(1+ziJi)
+  #### MCVE inf. method
+  yib <- c()
+  i=1
+  while (i<=r){
+    epsi <- rnorm(4,0,sig0)
+    zi <- runif(1,5,10)
+    xsimi <- c(VE(y0,zi,epsi[1]),VE(y0,zi,epsi[2]),VE(y0,zi,epsi[3]),VE(y0,zi,epsi[4]))
+    xvebar <- mean(xsimi)
+    ci=rchisq(1,df=(2*alpha+n-1))
+    mu <- y0*(zi+1)
+    Tib <- sqrt(2*beta + (n-1)*s^2)/(sig0*sqrt(ci))
+    yib[i] <- ((1+zi)^(-1)) * (Tib*(mu-xvebar)+(xbar-mu)) + y0
+    i=i+1
+  }
+  meantestJi[k] <- mean(yiJi)
+  meantestb[k] <- mean(yib)
+  sdtestJi[k] <- sd(yiJi)
+  sdtestb[k] <- sd(yib)
+  k=k+1
+}
+meantestJi
+meantestb
+sdtestJi
+sdtestb
+
+## alpha < 1
+k=1
+meansmallJi <- c()
+meansmallb <- c()
+sdsmallJi <- c()
+sdsmallb <- c()
+while (k<=9){
+  alpha <- k/10
+  sig0 <- 5
+  beta<- alpha * sig0^2
+  #### JCGM 101 + inf. method
+  ziJi <- runif(r,5,10)
+  xiJi <- xbar + rt(r,2*alpha + n-1)*sqrt((2*beta + (n-1)*s^2)/(n*(2*alpha + n-1)))
+  yiJi <- (xiJi)/(1+ziJi)
+  #### MCVE inf. method
+  yib <- c()
+  i=1
+  while (i<=r){
+    epsi <- rnorm(4,0,sig0)
+    zi <- runif(1,5,10)
+    xsimi <- c(VE(y0,zi,epsi[1]),VE(y0,zi,epsi[2]),VE(y0,zi,epsi[3]),VE(y0,zi,epsi[4]))
+    xvebar <- mean(xsimi)
+    ci=rchisq(1,df=(2*alpha+n-1))
+    mu <- y0*(zi+1)
+    Tib <- sqrt(2*beta + (n-1)*s^2)/(sig0*sqrt(ci))
+    yib[i] <- ((1+zi)^(-1)) * (Tib*(mu-xvebar)+(xbar-mu)) + y0
+    i=i+1
+  }
+  meansmallJi[k] <- mean(yiJi)
+  meansmallb[k] <- mean(yib)
+  sdsmallJi[k] <- sd(yiJi)
+  sdsmallb[k] <- sd(yib)
+  k=k+1
+}
+meansmallJi
+meansmallb
+sdsmallJi
+sdsmallb
+
+#### JCGM 101 with Perfect PK ####
+xiperf <- rnorm(r,50,5/2)
+ziperf <- runif(r,5,10)
+yiperf <- (xiperf)/(1+ziperf)
+
+#### Line Plot ####
+par(cex=1.5)
+plot(c((1:9)/10,1:15),c(sdsmallJi,sdtestJi),type="l",lwd=4,col=2,xlab="alpha",ylab="Standard Deviation",ylim=c(1.10,1.2))
+lines(c((1:9)/10,1:15),c(sdsmallb,sdtestb),lwd=4,lty=2,col="yellow")
+abline(h=sd(yiJ),col=1)
+abline(h=sd(yiperf),col=4)
+legend(x="right", legend=c("JCGM 101 + inf.","MC-VE inf.","JCGM 101 + No PK", "JCGM 101 + Perfect PK"),fill=c(2,"yellow",1,4))