# Installing useful packages
library(survival)
library(numDeriv)  # differential function
library(ggplot2)
library(magrittr)  # Using pipe operator
library(kableExtra)

Problem 1.

According to a general result, if T is the random variable associated with survival time of an individual, and \(S(t)\) is the corresponding survival function.
Show that the random variable \(U = -log\left(S(t)\right)\) has an exponential distribution with unit mean.
< Proof >
We know \[ \begin{align} S(t) = P(T>t) = 1- F(t) \\ \end{align} \] Thus \[ \begin{align} F_{U}(u) &= P(U \leq u) = P\bigg(-log\left(S(t)\right) \leq u \bigg) \\ &= P\left(S(t) \geq e^{-u} \right) \\ &= P\left(1-F(t) \geq e^{-u} \right) \\ &= P\left(F(t) \leq 1-e^{-u} \right) \\ &= P\left(T \leq F^{-1}(1-e^{-u}) \right) \\ &= F_{T}\left({F_{T}}^{-1}(1-e^{-u}) \right) \\ &= 1-e^{-u} \\ \Rightarrow f_{U}(u) &= \frac{dF_{U}(u)}{du} = e^{-u},\quad u>0 \\ \Rightarrow U &\sim exp(1) \quad Q.E.D \end{align} \]

Problem 2.

Let T be a random variable with conditional survival function S:
\[ S(t|x) = exp\left\{-\Lambda_{0}(t)e^{x\beta} \right\} \cdots \cdots\cdots (1) \] where \(\Lambda_{0}(t)\) corresponds to its cumulative baseline hazard function.

  1. Show that the survival time T of the Cox model (1) can be expressed as \[ T = \Lambda_{0}^{-1}\left\{-(logU)*exp(-x\beta) \right\},\quad where\ U\sim U[0,1] \] < Proof >
    We know \[ \begin{align} S(t|x) &= exp\left\{-\Lambda_{0}(t)e^{x\beta} \right\} \\ \Rightarrow -log\left(S(t|x) \right) &= \Lambda_{0}(t)e^{x\beta} \\ \Rightarrow T &= \Lambda_{0}^{-1}\{-log\left(U \right)*e^{-x\beta}\},\quad let\ U = S(t|x) \sim U[0,1] \end{align} \] where, we know \[ S(t|x)=P(T>t|x) = 1-F(t|x) \] Let \(U = S(t|x)\), thus \[ \begin{align} F_{U}(u) &= P(U \leq u) = P\left(S(t|x) \leq u \right) \\ &= P\left(1-F(t|x) \leq u \right) \\ &= P\left(F(t|x) \geq 1-u \right) \\ &= 1 - P\left(F(t|x) \leq 1-u \right) \\ &= 1 - F_{T}\left({F_{T}}^{-1}(1-u|x) \right) \\ &= 1 - ( 1 - u ) \\ &= u \\ \Rightarrow f_{U}(u) &= \frac{dF_{U}(u)}{du} = 1,\quad u \in [0,1] \\ \Rightarrow U &\sim U[0,1] \quad Q.E.D \end{align} \]

  2. Considering baseline density function \(f_{0}(t)=\alpha vt^{v-1}exp(-\alpha t^{v})\), which is a Weibull density function.
    Show that the survival time T has explicit form \[ T = \left(-\frac{logU}{\alpha*e^{x\beta}}\right)^{\frac{1}{v}}, \] and \[ \lambda(t|x) = \alpha e^{x\beta}vt^{v-1} \] < Proof >
    We know T follows a Weibull distribution, namely \[ \begin{align} f_{0}(t) &= \alpha vt^{v-1}exp(-\alpha t^{v}),\quad t>0,and\ \alpha,v\ >0 \\ \Rightarrow F_{0}(t) &= \int_{t}^{\infty} \alpha vu^{v-1}exp(-\alpha u^{v}) du = \alpha v \int_{t}^{\infty} u^{v-1}exp(-\alpha u^{v}) du \\ &\overset{Let\ w = \alpha u ^{v}}{=} \int_{\alpha t^{v}}^{\infty} exp(-w) dw \\ &= 1 - exp(-\alpha t ^{v}) \\ \Rightarrow S(t) &= P(T>t) = 1 - F_{0}(t) = exp(-\alpha t ^{v}) \\ \Rightarrow \lambda(t) &= \frac{f(t)}{S(t)}= \frac{\alpha vt^{v-1}exp(-\alpha t^{v})}{exp(-\alpha t ^{v})} = \alpha v t^{v-1} \\ \Rightarrow \Lambda(t) &= \int_{0}^{t} \alpha v u^{v-1} du = \alpha t^{v} \end{align} \] So, when T follows a Weibull distribution, then Cox model is \[ \begin{align} \lambda(t|x) &= \alpha vt^{v-1}*e^{x\beta}=\lambda_{0}(t)*e^{x\beta} \\ \Rightarrow S(t|x) &= exp\bigg\{ -\int_{0}^{t} \lambda(u|x) du \bigg\} \\ &= exp\bigg\{ -\alpha v e^{x\beta}\int_{0}^{t} u^{v-1} du \bigg\} \\ &= exp\bigg\{ -\alpha t^{v} e^{x\beta} \bigg\} \\ \Rightarrow -log\left(S(t|x)\right) &= \alpha t^{v} e^{x\beta} \\ \Rightarrow T &= \left(-\frac{logU}{\alpha*e^{x\beta}}\right)^{\frac{1}{v}},\quad where\ U=S(t|x) \sim~U[0,1]\quad Q.E.D \end{align} \]

  3. From (b), we can find that T follows a Weibull distribution. Write down its corresponding shape and scale parameter. We know
    由b小題證明可得 \[ T = \left(-\frac{logU}{\alpha*e^{x\beta}}\right)^{\frac{1}{v}},\quad U \sim U[0,1] \] 可由上述等式觀察到 \[ \left\{ \begin{aligned} \alpha &, Scale\ parameter \\ v &, Shape\ parameter \end{aligned} \right. \]

  4. Using R, assume sample size n=1000, \[ \alpha=0.1,\ \beta=1,\ v=2,\ X_{i}=\left\{ \begin{aligned} 1 &, i=1,3,...999 \\ 0 &, i=2,4,...1000 \\ \end{aligned} \right. \] Let censoring variable \(C_{i} \sim U[0,10]\). Base on the setting, we can generate observed data\((Y_{i},\delta_{i},X_{i})\),
    where \(Y_{i}=min(T_{i},C_{i}),\delta_{i}=I(T_{i}\leq C_{i}),i=1,2...1000.\)

    # Generate observed data
set.seed(1000)
n <- 1000
x <- rep(c(1,0), n/2)
U <- runif(n)
t <- ( -log(U)/ (0.1 * exp(x*1) ) )^(0.5)
C <- runif(n, min = 0, max =10)
data.hw2 <- data.frame(t, C)
data.hw2$Y <- apply(data.hw2, 1, min)
data.hw2$d <- ifelse(data.hw2$t <= data.hw2$C, 1, 0 )
data.hw2$x <- x   # length() = 1000
head(data.hw2)   # dim() = 1000 * 5 and sum(data.hw2$d) = 753

    d-1. Since we know T follows a Weibull distribution, based on your simulated data using MLE approach to estimate \(\beta\) (Denote by \(\hat{\beta}_{M}\))
    Given,T follows a Weibull distribution, thus,
    \[ \begin{align} L = L(\alpha,v,\beta) &= \Pi_{i=1}^{1000} \left[f(y_{i};\alpha,v,\beta) \right]^{\delta_{i}} * \left[S(y_{i};\alpha,v,\beta) \right]^{1-\delta_{i}} \\ &= \Pi_{i=1}^{1000} \left[\lambda(y_{i};\alpha,v,\beta) \right]^{\delta_{i}} * S(y_{i};\alpha,v,\beta) \\ &= \Pi_{i=1}^{1000} \left[\alpha*v*y_{i}^{v-1}*e^{\beta x_{i}} \right]^{\delta_{i}} * e^{-\alpha*y_{i}^{v}*e^{\beta x_{i}}} \\ \Rightarrow lnL &= \sum_{i=1}^{1000} \bigg\{ \delta_{i}*\left[ln\alpha +lnv+(v-1)*lny_{i}+\beta x_{i}\right] - \alpha*y_{i}^{v}* e^{\beta*x_{i}} \bigg\} \\ \Rightarrow &\left\{ \begin{aligned} \frac{\partial lnL}{\partial \alpha} &= \sum_{i=1}^{1000}\bigg\{ \frac{\delta {i}}{\alpha}-y_{i}^{v}e^{\beta*x_{i}}\bigg\} \\ \frac{\partial lnL}{\partial v} &= \sum_{i=1}^{1000}\bigg\{ \delta {i}\left[ \frac{1}{v}+lny_{i}\right] - \alpha y_{i}^{v}*lny_{i}*e^{\beta*x_{i}}\bigg\} \\ \frac{\partial lnL}{\partial \beta} &= \sum_{i=1}^{1000}\bigg\{\delta_{i}*x_{i} - \alpha*y_{i}^{v}* e^{\beta*x_{i}} *x_{i} \bigg\} \\ \end{aligned} \right. \\ \Rightarrow S(\alpha,v,\beta) &=\begin{pmatrix} \frac{\partial lnL}{\partial \alpha} \\ \frac{\partial lnL}{\partial v} \\ \frac{\partial lnL}{\partial \beta} \\ \end{pmatrix}_{3 * 1} \\ &= \begin{pmatrix} \sum_{i=1}^{1000}\bigg\{ \frac{\delta {i}}{\alpha}-y_{i}^{v}e^{\beta*x_{i}}\bigg\} \\ \sum_{i=1}^{1000}\bigg\{ \delta {i}\left[ \frac{1}{v}+lny_{i}\right] - \alpha y_{i}^{v}*lny_{i}*e^{\beta*x_{i}}\bigg\} \\ \sum_{i=1}^{1000}\bigg\{\delta_{i}*x_{i} - \alpha*y_{i}^{v}* e^{\beta*x_{i}} *x_{i} \bigg\} \\ \end{pmatrix}_{3 * 1} \in \mathbb{R}^{3*1}\quad (Score\ function) \\ \Rightarrow &\left\{ \begin{aligned} \frac{\partial^2 lnL}{\partial \alpha^2} &= \sum_{i=1}^{1000} \frac{-\delta {i}}{\alpha^2} \\ \frac{\partial^2 lnL}{\partial \alpha \partial v} &= \sum_{i=1}^{1000} -y_{i}^{v}*lny_{i}*e^{\beta x_{i}} \\ \frac{\partial^2 lnL}{\partial \alpha \partial \beta} &= \sum_{i=1}^{1000} -y_{i}^{v}*e^{\beta x_{i}}*x_{i} \\ \frac{\partial^2 lnL}{\partial v^2} &= \sum_{i=1}^{1000}\bigg\{ \delta {i}\left[ \frac{-1}{v^2}\right] - \alpha y_{i}^{v}*(lny_{i})^2*e^{\beta*x_{i}}\bigg\} \\ \frac{\partial^2 lnL}{\partial v \partial \beta} &= \sum_{i=1}^{1000} -\alpha*y_{i}^{v}*lny_{i}*e^{\beta x_{i}}*x_{i} \\ \frac{\partial^2 lnL}{\partial \beta^2} &= \sum_{i=1}^{1000} - \alpha*y_{i}^{v}* e^{\beta*x_{i}} *x_{i}^2 \\ \end{aligned} \right. \\ \Rightarrow H(\alpha,v,\beta) &= \begin{pmatrix} \frac{\partial^2 lnL}{\partial \alpha^2} & \frac{\partial^2 lnL}{\partial \alpha \partial v} & \frac{\partial^2 lnL}{\partial \alpha \partial \beta} \\ \frac{\partial^2 lnL}{\partial v \partial \alpha} & \frac{\partial^2 lnL}{\partial v^2} & \frac{\partial^2 lnL}{\partial v \partial \beta} \\ \frac{\partial^2 lnL}{\partial \beta \partial \alpha} & \frac{\partial^2 lnL}{\partial \beta \partial v} & \frac{\partial^2 lnL}{\partial \beta^2} \\ \end{pmatrix}_{3 * 3},\quad(Symmetry) \\ &= \begin{pmatrix} \sum_{i=1}^{1000} \frac{-\delta {i}}{\alpha^2} & \sum_{i=1}^{1000} -y_{i}^{v}*lny_{i}*e^{\beta x_{i}} & \sum_{i=1}^{1000} -y_{i}^{v}*e^{\beta x_{i}}*x_{i} \\ & \sum_{i=1}^{1000}\bigg\{ \delta {i}\left[ \frac{-1}{v^2}\right] - \alpha y_{i}^{v}*(lny_{i})^2*e^{\beta*x_{i}}\bigg\} & \sum_{i=1}^{1000} -\alpha*y_{i}^{v}*lny_{i}*e^{\beta x_{i}}*x_{i} \\ & & \sum_{i=1}^{1000} - \alpha*y_{i}^{v}* e^{\beta*x_{i}} *x_{i}^2 \\ \end{pmatrix}_{3 * 3} \\ &\in \mathbb{R}^{3*3}\quad (Hessian\ matrix) \end{align} \] So,the Newton-Raphson algorithm is \[ \begin{align} \hat{\theta}^{(s+1)} &= \hat{\theta}^{(s)} - H^{-1}(\hat{\theta}^{(s)})*S(\hat{\theta}^{(s)}) \in \mathbb{R}^{3*1}, where,\ \hat{\theta}^{(s)} =\begin{pmatrix} \hat{\alpha}^{(s)} \\ \hat{v}^{(s)} \\ \hat{\beta}^{(s)} \\ \end{pmatrix}_{3 * 1} \\ \end{align} \] 其中, \(\hat{\theta}^{(s)}\)表第s次迭代
    將Newton-Raphson algorithm寫成code如下

    # Estimate coefficient by MLE(full likelihood)
MLE.NR <- function(theta = c(a, v, b), tol1 = 1e-8, tol2 = 1e-6){
  d <- data.hw2$d
  y <- data.hw2$Y
  x <- data.hw2$x
  B <- t(t(theta))    # dim() = 3 * 1  Initial value
  for(i in 1:1000){   # Fixed iteration 10000 times
    a <- B[1]
    v <- B[2]
    b <- B[3]
    # Calculate Score function
    S1 <- sum(d/a - y^{v}*exp(b*x))
    S2 <- sum(d*((1/v)+log(y)) - a*y^{v}*log(y)*exp(b*x))
    S3 <- sum(d*x - a*y^{v}*exp(b*x)*x)
    S <- rbind(S1, S2, S3)     # Score function
    # Calculate Hessian matrix
    H11 <- sum(-d/(a^2))
    H12 <- H21 <- sum(-y^{v}*log(y)*exp(b*x))
    H13 <- H31 <- sum(-y^{v}*exp(b*x)*x)
    H22 <- sum((-d/(v^2)) - a*y^{v}*(log(y))^{2}*exp(b*x))
    H23 <- H32 <- sum(-a*y^{v}*log(y)*exp(b*x)*x)
    H33 <- sum(-a*y^{v}*exp(b*x)*x^{2})
    H <- matrix(c(H11, H12, H13,
                  H21, H22, H23,
                  H31, H32, H33), nrow = 3, byrow = TRUE) # Hessian matrix
    Bn <- B - solve(H) %*% S 
    if( all( abs(Bn-B)  <  tol1 * ( abs(B) + tol2 )   )) {
      break
      }
    B <- Bn
    }
  MLE.NR.list <- list(Iteration = i, Coefficients = B)
  rownames(MLE.NR.list[[2]]) <- c("a", "v", "b")
  colnames(MLE.NR.list[[2]]) <- c("")
  return(MLE.NR.list)
  }
# Assume initial value is  a = 0.1, v = 1, b = 1.4
MLE.NR(theta = c(a = 0.1, v = 1, b = 1.4))
    ## $Iteration
## [1] 8
## 
## $Coefficients
##             
## a 0.09077915
## v 2.06731562
## b 1.10416102

    故在起始值(Initial value)取 \[ \begin{align} \theta = \begin{pmatrix} \alpha = 0.1 \\ v = 1 \\ \beta = 1.4 \\ \end{pmatrix}_{3 * 1} \end{align} \]
    用Newton-Raphson algorithm可得 \[ \begin{align} MLE\ of\ \theta= \hat{\theta} = \begin{pmatrix} \hat{\alpha} = 0.0908 \\ \hat{v} = 2.0673 \\ \hat{\beta} = 1.1042 \\ \end{pmatrix}_{3 * 1} \Rightarrow \hat{\beta}_{M} = 1.1042 \end{align} \]

    Note:Use numDeriv Package also can find the score and Hession matrix

    # Writing log likelihood 
log.likelihood <- function(theta){
  y <- data.hw2$Y
  d <- data.hw2$d
  x <- data.hw2$x
  a <- theta[1]
  v <- theta[2]
  b <- theta[3]
  value <- sum( d * ( log(a) + log(v) + (v-1)*log(y) + b*x ) 
            - ( a * y^{v} * exp(b*x) ))
  return(value) 
}
# Newton-Raphson algorithm
MLE.NR.package <- function(theta, tol1 = 1e-8, tol2 = 1e-6){
  for(i in 1: 1000){
    a <- theta[1]
    v <- theta[2]
    b <- theta[3]
    theta <- t(t(theta))
    S <- t(t(grad(log.likelihood, c(a, v, b)))) # Score function
    H <- hessian(log.likelihood, c(a, v, b)) # Hessian matrix 
    theta.new <- theta - solve(H) %*% S
    if(all(abs(theta - theta.new ) < tol1*(abs(theta) + tol2))){
      break
      }else{
        theta <- theta.new
      }
    }
  MLE.NR.list <- list(Iteration = i, Coefficients = t(t(theta)))
  rownames(MLE.NR.list[[2]]) <- c("a", "v", "b")
  colnames(MLE.NR.list[[2]]) <- c("")
  return(MLE.NR.list)
  }
MLE.NR.package(theta = c(a = 0.1, v = 1.2, b = 1.4))
    ## $Iteration
## [1] 9
## 
## $Coefficients
##             
## a 0.09077915
## v 2.06731562
## b 1.10416102

    d-2. Pretend that we only know T follows a proportional hazard model(PHM), based on your siumlated data using partial likelihood approach to estimate \(\beta\) (Denote by \(\hat{\beta}_{P}\)) (In R:coxph function)

    # Estimate coefficient by partial likelihood
# 1. 將資料轉換成survival object
data.hw2$SurvObj <- with(data.hw2, Surv(data.hw2$Y, data.hw2$d)) # use Surv() to build the standard survival object
head(data.hw2, 7)

    由上述表格最後一行(column)可看出,資料已成功轉換成survival object(right censoring)

    # 2.Use coxph function
model.coxph <- coxph(SurvObj ~ x, data =  data.hw2, method = 'breslow')
model.coxph %>% summary %>% coef %>% round(., 4) %>% 
  kable(., "html", 
        caption = "Cox proportional hazards model") %>% 
  kable_styling(bootstrap_options = "striped", 
                full_width = F) %>% 
  footnote(general = paste("n =", model.coxph$n))
    Cox proportional hazards model
    coef exp(coef) se(coef) z Pr(>|z|)
    x 1.0805 2.9462 0.0811 13.3158 0
    Note:
    n = 1000

    由 coxph function可知, ( partial likelihood ) \[ The\ MLE\ of\ \beta=\hat{\beta}_{p}=1.0805 \]

  5. Pretend that we only know T follows a proportional hazard model(PHM), plot Cox-snell residuals(x-axis) versus its cumulative hazard function (y-axis) to see if slope is 1.
    Define \[ \begin{align} &r_{Mi}\ be\ Martingale\ Residuals\ (R\ defalt) \\ &r_{Ci}\ be\ Cox-Snell\ Residuals\ \end{align} \] 且上述兩者具有如下關係式 \[ r_{Mi} = \delta_{i}-r_{Ci},\quad where\ \delta_{i}=\left\{ \begin{aligned} 1 &, uncensored \\ 0 &, censored \\ \end{aligned} \right. \\ \Rightarrow r_{Ci} = \delta_{i}-r_{Mi} \] 將Cox-snell residuals計算寫成code如下所示

    # Cox-snell residuals
cox.snell <- data.hw2$d - resid(model.coxph) 
sv.cox.snell <- survfit(Surv(cox.snell, data.hw2$d) ~ 1)
ggplot() +
  geom_point(aes(x = sv.cox.snell$time, 
                 y = -log(sv.cox.snell$surv))) +
  labs(title = "Cumulative Hazards of Cox-Snell Residuals",
       x = expression(r[Ci] ),
       y = expression(-log(S[r[C]](r[Ci])) ), 
       subtitle = "A way of checking model by Cox-Snell residuals") +
  geom_abline(intercept = 0, slope = 1)  # add 45 degree line

    由上述圖型可看出,因為大部分資料皆靠著45度線,故此模型配適程度是可以接受的.

  6. Regenerate 100 simulating data sets like (d), calculate the mean and 95% C.I of \(\hat{\beta}_{M}\) and \(\hat{\beta}_{P}\) .

    # 1. 寫下生成樣本的函數?
hw2.sample <- function(n, i ){
  set.seed(i)
  x <- rep(c(1,0), n/2)
  U <- runif(n)
  t <- sqrt( -log(U)/ (0.1 * exp(x) ) )
  C <- runif(n, min = 0, max =10)
  data.hw2 <- data.frame(t, C)
  data.hw2$Y <- apply(data.hw2, 1, min)
  data.hw2$d <- ifelse(data.hw2$t <= data.hw2$C, 1, 0 )
  data.hw2$x <- x 
  data.hw2$SurvObj <- with(data.hw2, Surv(data.hw2$Y, data.hw2$d))
  return(data.hw2)
}
    # 2. Create table by 100 simulating data sets
table.coefficient <- matrix(NA, ncol = 2, nrow = 120)
colnames(table.coefficient) <- c("Full", "Partial")
# z <- numeric()
for(i in 1: 120){
  data.hw2 <- hw2.sample(1000, i = 1000 + i )
  # Newton-Raphson
  temp <- try(MLE.NR(theta = c(a = 0.1, v = 1, b = 1.4))[[2]][3], silent = TRUE)
  if (class(temp) != "try-error"){
    table.coefficient[i,1] <- temp
  }
  # coxph
  model.coxph <- coxph(SurvObj ~ x, data =  data.hw2)
  table.coefficient[i,2] <- model.coxph$coefficients
}

# Remove all NA
table.coefficient <- na.omit(table.coefficient)
table.coefficient <- table.coefficient[1:100,1:2]
head(data.frame(table.coefficient), 3)  # dim()  = 100 * 2

    上述表格表每次生成資料時各自的?\(\hat{\beta}\)

    # Mean and Variance
rbind(apply(table.coefficient, 2, mean), 
      apply(table.coefficient, 2, var)) %>% 
  `rownames<-`(., c("Mean", "Variance")) %>% 
  round(., 4) %>% 
  kable(., "html") %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
    Full Partial
    Mean 1.0006 1.0038
    Variance 0.0058 0.0062

    由上述表格可看出,用Full likelihood 求出來的\(\beta\)?其變異數較用Partial likelihood低,表 \(\hat{\beta}_{M}較\hat{\beta}_{P}有效\).
    故其95%C.I為
    \[ \begin{align} \hat{\beta} \pm1.96*S.E(\hat{\beta}) \end{align} \] \[ \begin{align} \Rightarrow \left\{ \begin{aligned} \hat{\beta}_{M} \pm1.96*S.E(\hat{\beta}_{M}) &= 1.000604 \pm 1.96*\sqrt{0.005808805} \\ &= ( 0.8512219,\ 1.149987) \\ &\approx (0.8512,\ 1.15) \\ \hat{\beta}_{P} \pm1.96*S.E(\hat{\beta}_{P}) &= 1.003840 \pm 1.96*\sqrt{0.006189142 } \\ &=(0.8496449,\ 1.158036) \\ &\approx (0.8496,\ 1.158) \end{aligned} \right. \end{align} \]

  7. Describe what you find in (f).
    由上述比較可看出，\(\hat{\beta}_{M}\)之信賴區間較窄,其原因為估計\(\hat{\beta}_{M}\)是使用到Full likelihood去進行估計,有使用到較多資訊,故其variance會較使用Partial likelihood去估計的\(\hat{\beta}_{P}\)較低,即 \[ \begin{align} Var(\hat{\beta}_{M}) \leq Var(\hat{\beta}_{P}) \end{align} \] 故\(\hat{\beta}_{M}\)之信賴區間較窄合理.

  8. From (f), recoding how many simlating data sets reject the null hypothesis \(\beta\) =1 and describe what you find (Just considering partial likelihood approach). \[ \left\{ \begin{aligned} H_0 : \beta = 1 \\ H_1 : \beta\neq 1 \\ \end{aligned} \right. \] 因為 \[ \hat{\beta}_{P}之95\%C.I\approx (0.8496,\ 1.158) \\ \Rightarrow 1 \in (0.8496,\ 1.158) \] 表在\(\alpha=0.05\)下我們沒有充分證據eject\(H_0:\beta=1\).

    # 95%C.I for partial likelihood 
LCI <- mean(table.coefficient[,2]) - 1.96*sqrt(var(table.coefficient[,2])) # Lower C.I
UCI <- mean(table.coefficient[,2]) + 1.96*sqrt(var(table.coefficient[,2])) # Upper C.I
z <- ifelse(table.coefficient[,2] < LCI | UCI < table.coefficient[,2], 1, 0)
z
    ##   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [36] 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    sum(z)
    ## [1] 2

    由上述程式計算知道,在\(\hat{\beta}_{P}之95\%C.I\approx(0.8496,1.158)\)下,(f)的data set中,有兩個data set 會reject null hypothesis\(\beta=1\)