hwk741_3

1

z<-c(0.8,0.6,0.4,0.6,0.1,0.4,0.3,0.5)
idtr <- c(1,0,1,1,0,0,1,1)
which(idtr==1)

## [1] 1 3 4 7 8

Ub <- function(b){
  u <- c()
  for(i in which(idtr==1)){
     ui<-z[i] - (sum(z[i:length(z)]*exp(z[i:length(z)]*b)))/(sum(exp(z[i:length(z)]*b)))
     u <-c(u,ui)
  }
  return(sum(u))
}

plot(-2:12, sapply(-2:12,FUN = Ub))

ABSUB<- function(b){
  abs(Ub(b))
}
optimize(f = ABSUB,interval = c(-2,12))

## $minimum
## [1] 5.6908865848836871
## 
## $objective
## [1] 7.0585918393595293e-07

Ub_prime <- function(b){
  u <- c()
  for(i in which(idtr==1)){
    ui<- - (sum(z[i:length(z)]^2*exp(z[i:length(z)]*b)) * sum(exp(z[i:length(z)]*b)) -
              sum(z[i:length(z)]*exp(z[i:length(z)]*b))^2 )/(sum(exp(z[i:length(z)]*b)))^2
    u <-c(u,ui)
  }
  return(sum(u))
}
Ub_prime(5.69)

## [1] -0.060252627044907095

2

library(xlsx)

## Loading required package: rJava

## Loading required package: xlsxjars

d1<-read.xlsx(file = 'Failure_Time/Data Set I.xls',1)
head(d1)

##   treat type time censor status mfd age prior
## 1     1    1   72      1     60   7  69     0
## 2     1    1  411      1     70   5  64    10
## 3     1    1  228      1     60   3  38     0
## 4     1    1  126      1     60   9  63    10
## 5     1    1  118      1     70  11  65    10
## 6     1    1   10      1     20   5  49     0

logrank test

library(survival)

## Warning: package 'survival' was built under R version 3.2.5

survdiff(Surv(time,censor) ~ treat, data = d1 )

## Call:
## survdiff(formula = Surv(time, censor) ~ treat, data = d1)
## 
##          N Observed      Expected         (O-E)^2/E         (O-E)^2/V
## treat=1 69       64 64.5001966636 0.003879006812683 0.008227343202350
## treat=2 68       64 63.4998033364 0.003940117750476 0.008227343202351
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9277272333401

cox model

coxph(Surv(time,censor)~treat + age + prior,data = d1)

## Call:
## coxph(formula = Surv(time, censor) ~ treat + age + prior, data = d1)
## 
##                   coef        exp(coef)         se(coef)        z       p
## treat  0.0036213475721  1.0036279125736  0.1831510837382  0.01977 0.98422
## age    0.0071235764879  1.0071490095144  0.0096749421742  0.73629 0.46155
## prior -0.0135622780308  0.9865292753039  0.0201177459087 -0.67415 0.50022
## 
## Likelihood ratio test=1.09  on 3 df, p=0.7797321823695
## n= 137, number of events= 128