Monday, September 30, 2013

Time Series II: AR, MA, ARMA, ARIMA

rm(list=ls())
# white noise
set.seed(1)
w <- rnorm(100)
plot(w, type = "l")
x <- seq(-3,3, length = 1000)
hist(rnorm(100), prob = T); points(x, dnorm(x), type = "l")
set.seed(2)
acf(rnorm(100))
# random walk
set.seed(2)
x <- w <- rnorm(1000)
for (t in 2:1000) x[t] <- x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
acf(diff(x))
www <- "http://elena.aut.ac.nz/~pcowpert/ts/HP.txt"
HP.dat <- read.table(www, header = T) ; attach(HP.dat)
plot (as.ts(Price))
DP <- diff(Price) ; plot (as.ts(DP)) ; acf (DP)
mean(DP) + c(-2, 2) * sd(DP)/sqrt(length(DP))
# AR(1)
rho <- function(k, alpha) alpha^k
layout(1:2)
plot(0:10, rho(0:10, 0.7), type = "b")
plot(0:10, rho(0:10, -0.7), type = "b")
layout(1:1)
set.seed(1)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
x.ar <- ar(x, method = "mle")
x.ar$order
x.ar$ar
x.ar$ar + c(-2, 2) * sqrt(x.ar$asy.var)
www = "http://elena.aut.ac.nz/~pcowpert/ts/pounds_nz.dat"
Z = read.table(www, header = T)
Z.ts = ts(Z, st = 1991, fr = 4)
Z.ar <- ar(Z.ts)
mean(Z.ts)
Z.ar$order
Z.ar$ar
Z.ar$ar + c(-2, 2) * sqrt(Z.ar$asy.var)
acf(Z.ar$res[-1])
www = "http://elena.aut.ac.nz/~pcowpert/ts/global.dat"
Global = scan(www)
Global.ts = ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.ar <- ar(aggregate(Global.ts, FUN = mean), method = "mle")
mean(aggregate(Global.ts, FUN = mean))
Global.ar$order
Global.ar$ar
acf(Global.ar$res[-(1:Global.ar$order)], lag = 50)
# MA: correlogram and simulation
# correlogram
rho = function(k, beta){
  q = length(beta) - 1
  if (k > q) ACF = 0
  else{
    s1 = 0
    s2 = 0
    for (i in 1:(q - k + 1)){
      s1 = s1 + beta[i] * beta[i + k]
    }
    for (i in (1:(q + 1))){
      s2 = s2 + beta[i]^2
    }
    ACF = s1/s2
  }
  ACF 
}
 
beta <- c(1, 0.7, 0.5, 0.2) 
rho.k <- rep(1, 10)
for (k in 1:10) rho.k[k] <- rho(k, beta)
plot(0:10, c(1, rho.k), pch = 4, ylab = expression(rho[k])) 
abline(0, 0) 
 
# simulation of MA(3) 
set.seed(1) 
b = c(0.8, 0.6, 0.4) 
x = w = rnorm(1000)
for (t in c(4:1000)){
  for (j in c(1:3)){
    x[t] = x[t] + b[j] * w[t - j]
  }
}
plot(x, type = "l")
acf(x)
x.ma <- arima(x, order = c(0, 0, 3))
x.ma
x.ma <- arima(x, order = c(0, 0, 3), include.mean=FALSE)
x.ma
# fit MA model
www = "http://elena.aut.ac.nz/~pcowpert/ts/pounds_nz.dat"
x <- read.table(www, header = T)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ma
# ARMA simulation
set.seed(1)
x = arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1)))
AIC(arima(x, order = c(1, 0, 1)))
acf(resid(arima(x, order = c(1, 0, 1))))
www = "http://elena.aut.ac.nz/~pcowpert/ts/cbe.dat"
CBE = read.table(www, header = T)
Elec.ts = ts(CBE[,3], start = 1958, freq = 12)
Time = 1:length(Elec.ts)
Imth = cycle(Elec.ts)
Elec.lm = lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth))
summary(Elec.lm)
acf(resid(Elec.lm))
# train ARMA model
best.order = c(0, 0, 0)
best.aic = Inf
for (i in c(0:2)){
  for (j in c(0:2)){
    fit.aic = AIC(arima(resid(Elec.lm), order = c(i, 0, j)))
    if (fit.aic < best.aic){
      best.order = c(i, 0, j)
      best.arma = arima(resid(Elec.lm), order = best.order)
      best.aic = fit.aic           
    }
  }
}
best.order
acf(resid(best.arma))
# prediction
new.time <- seq(length(Elec.ts), length = 36)
new.data <- data.frame(Time = new.time, Imth = rep(1:12, 3))
predict.lm <- predict(Elec.lm, new.data)
predict.arma <- predict(best.arma, n.ahead = 36)
elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991, freq = 12)
ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)
www = "http://elena.aut.ac.nz/~pcowpert/ts/wave.dat"
wave.dat <- read.table(www, header = T)
attach(wave.dat)
layout(1:3)
plot(as.ts(waveht), ylab = 'Wave height')
acf(waveht)
pacf (waveht)
wave.arma <- arima(waveht, order = c(4,0,4))
acf(wave.arma$res[-(1:4)])
pacf(wave.arma$res[-(1:4)])
hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')
# train Non-seasonal ARIMA model
www = "http://elena.aut.ac.nz/~pcowpert/ts/cbe.dat"
CBE = read.table(www, header = T)
Elec.ts = ts(CBE[, 3], start = 1958, freq = 12)
layout(c(1, 1, 2, 3))
plot(Elec.ts)
plot(diff(Elec.ts))
plot(diff(log(Elec.ts)))
layout(c(1, 1))
# ARIMA(1, 1, 1) simulation and inference
set.seed(1)
x = w = rnorm(1000)
for (i in 3:1000){
  x[i] = 0.5 * x[i - 1] + x[i - 1] - 0.5 * x[i - 2] + w[i] + 0.3 * w[i - 1]
}
plot(x, type = "l")
arima(x = x, order  = c(1, 1, 1))
x = arima.sim(model = list(order = c(1, 1, 1), ar = 0.5, ma = 0.3), n = 1000)
arima(x = x, order  = c(1, 1, 1))
# IMA(1, 1)
Beer.ts = ts(CBE[, 2], start = 1958, freq = 12)
Beer.ima = arima(Beer.ts, order = c(0, 1, 1))
Beer.ima
acf(resid(Beer.ima))
# prediction
Beer.1991 = predict(Beer.ima, n.ahead = 12)
plot(Beer.1991$pred)
plot(Beer.1991$se)
Posted by at No comments:

Time Series I: Basics

rm(list=ls())
Fact = function(n) if (n == 1) 1 else n * Fact(n - 1)
Fact(5)
data(AirPassengers)
AP = AirPassengers
AP
class(AP)
start(AP)
end(AP)
frequency(AP)
summary(AP)
plot(AP, ylab = "Passengers (1000's)")
layout(1:2)
plot(aggregate(AP))
boxplot(AP ~ cycle(AP))
www = "http://elena.aut.ac.nz/~pcowpert/ts/Maine.dat"
Maine.month = read.table(www, header = TRUE)
attach(Maine.month)
class(Maine.month)
Maine.month.ts = ts(unemploy, start = c(1996, 1), freq = 12)
Maine.annual.ts = aggregate(Maine.month.ts)/12
layout(1:2)
plot(Maine.month.ts, ylab = "unemployed (%)")
plot(Maine.annual.ts, ylab = "unemployed (%)")
Maine.Feb = window(Maine.month.ts, start = c(1996,2), freq = TRUE)
plot(Maine.Feb)
Maine.Aug = window(Maine.month.ts, start = c(1996,8), freq = TRUE)
plot(Maine.Aug)
Feb.ratio = mean(Maine.Feb) / mean(Maine.month.ts)
Aug.ratio = mean(Maine.Aug) / mean(Maine.month.ts)
Feb.ratio
Aug.ratio
layout(1:1)
www = "http://elena.aut.ac.nz/~pcowpert/ts/USunemp.dat"
US.month = read.table(www, header = T)
attach(US.month)
US.month.ts = ts(USun, start=c(1996,1), end=c(2006,10), freq = 12)
plot(US.month.ts, ylab = "unemployed (%)")
# Multiple time series
www = "http://elena.aut.ac.nz/~pcowpert/ts/cbe.dat"
CBE = read.table(www, header = T)
CBE[1:4, ]
class(CBE)
Elec.ts = ts(CBE[, 3], start = 1958, freq = 12)
Beer.ts = ts(CBE[, 2], start = 1958, freq = 12)
Choc.ts = ts(CBE[, 1], start = 1958, freq = 12)
plot(cbind(Elec.ts, Beer.ts, Choc.ts))
length(AP)
length(Elec.ts)
AP.elec = ts.intersect(AP, Elec.ts)
class(AP.elec)
AP.elec
AP
Elec.ts
AP = AP.elec[,1]
Elec = AP.elec[,2]
layout(1:2)
plot(AP, main = "", ylab = "Air passengers / 1000's")
plot(Elec, main = "", ylab = "Electricity production / MkWh")
layout(1:1)
plot(as.vector(AP), as.vector(Elec),
     xlab = "Air passengers / 1000's",
     ylab = "Electricity production / MWh")
abline(reg = lm(Elec ~ AP))
cor(AP, Elec)
www = "http://elena.aut.ac.nz/~pcowpert/ts/pounds_nz.dat"
Z = read.table(www, header = T)
Z.ts = ts(Z, st = 1991, fr = 4)
plot(Z.ts, xlab = "time / years",
     ylab = "Quarterly exchange rate in $NZ / pound")
Z.92.96 = window(Z.ts, start = c(1992, 1), end = c(1996, 1))
Z.96.98 = window(Z.ts, start = c(1996, 1), end = c(1998, 1))
layout (1:2)
plot(Z.92.96, ylab = "Exchange rate in $NZ/pound",
     xlab = "Time (years)" )
plot(Z.96.98, ylab = "Exchange rate in $NZ/pound",
     xlab = "Time (years)" )

www = "http://elena.aut.ac.nz/~pcowpert/ts/global.dat"
Global = scan(www)
Global.ts = ts(Global, st = c(1856, 1), end = c(2005, 12), fr = 12)
Global.annual = aggregate(Global.ts, FUN = mean)
plot(Global.ts)
plot(Global.annual)
layout(1:1)
New.series = window(Global.ts, start=c(1970, 1), end=c(2005, 12))
New.time = time(New.series)
plot(New.series)
abline(reg=lm(New.series ~ New.time))

# decomposition in R
plot(decompose(Elec.ts))
Elec.decom = decompose(Elec.ts, type = "mult")
plot(Elec.decom)
Trend = Elec.decom$trend
Seasonal = Elec.decom$seasonal
ts.plot(cbind(Trend, Trend * Seasonal), lty = 1:2)

www = "http://elena.aut.ac.nz/~pcowpert/ts/Herald.dat"
Herald.dat = read.table(www, header = T)
attach (Herald.dat)
x = CO; y = Benzoa; n = length(x)
sum((x - mean(x))*(y - mean(y))) / (n - 1)
mean((x - mean(x)) * (y - mean(y)))
cov(x, y)
cov(x,y) / (sd(x)*sd(y))
cor(x,y)
www = "http://elena.aut.ac.nz/~pcowpert/ts/wave.dat"
wave.dat = read.table (www, header=T)
attach(wave.dat)
plot(ts(waveht))
plot(ts(waveht[1:60]))
# acf
acf(waveht)$acf[2]
acf(waveht)$acf[1]
# autocovariance
acf(waveht, type = c("covariance"))$acf[2]
# correlogram
acf(waveht)
acf(AirPassengers)
data(AirPassengers)
AP = AirPassengers
AP.decom = decompose(AP, "multiplicative")
plot(ts(AP.decom$random[7:138]))
acf(AP.decom$random[7:138])
sd(AP[7:138])
sd(AP[7:138] - AP.decom$trend[7:138])
sd(AP.decom$random[7:138])
www = "http://elena.aut.ac.nz/~pcowpert/ts/Fontdsdt.dat"
Fontdsdt.dat = read.table(www, header=T)
attach(Fontdsdt.dat)
plot(ts(adflow), ylab = 'adflow')
acf(adflow, xlab = 'lag (months)', main="")
Build.dat = read.table("http://elena.aut.ac.nz/~pcowpert/ts/ApprovActiv.dat", header=T)
attach(Build.dat)
App.ts = ts(Approvals, start = c(1996,1), freq=4)
Act.ts = ts(Activity, start = c(1996,1), freq=4)
ts.plot(App.ts, Act.ts, lty = c(1,3))

# correlogram and cross-correlogram
acf(ts.union(App.ts, Act.ts))
print(ts.union(App.ts, Act.ts))
print(acf(ts.union(App.ts, Act.ts)))

# Bass model
T79 = 1:10
Tdelt = (1:100) / 10
Sales = c(840,1470,2110,4000, 7590, 10950, 10530, 9470, 7790, 5890)
Cusales = cumsum(Sales)
Bass.nls = nls(Sales ~ M * ( ((P+Q)^2 / P) * exp(-(P+Q) * T79) ) /
                  (1+(Q/P)*exp(-(P+Q)*T79))^2, start = list(M=60630, P=0.03, Q=0.38))
summary(Bass.nls)
Bcoef = coef(Bass.nls)
m = Bcoef[1]
p = Bcoef[2]
q = Bcoef[3]
ngete = exp(-(p+q) * Tdelt)
Bpdf = m * ( (p+q)^2 / p ) * ngete / (1 + (q/p) * ngete)^2
plot(Tdelt, Bpdf, xlab = "Year from 1979",
     ylab = "Sales per year", type='l')
points(T79, Sales)
Bcdf = m * (1 - ngete)/(1 + (q/p)*ngete)
plot(Tdelt, Bcdf, xlab = "Year from 1979",
     ylab = "Cumulative sales", type='l')
points(T79, Cusales)

?HoltWinters
require(graphics)
## Seasonal Holt-Winters
(m <- HoltWinters(co2))
plot(m)
plot(fitted(m))
(m <- HoltWinters(AirPassengers, seasonal = "mult"))
plot(m)
## Non-Seasonal Holt-Winters
x <- uspop + rnorm(uspop, sd = 5)
m <- HoltWinters(x, gamma = FALSE)
plot(m)
## Exponential Smoothing
m2 <- HoltWinters(x, gamma = FALSE, beta = FALSE)
lines(fitted(m2)[,1], col = 3)

Posted by at No comments:

Sunday, September 29, 2013

Some thoughts of subquery efficiency

from http://blog.csdn.net/haiwer/article/details/2826881

1. Use LEFT JOIN instead of NOT IN, NOT EXISTS

Change
SELECT PUB_NAME
FROM PUBLISHERS
WHERE PUB_ID NOT IN
   (SELECT PUB_ID
   FROM TITLES
   WHERE TYPE = 'BUSINESS')

TO
SELECT A.PUB_NAME
FROM PUBLISHERS A LEFT JOIN TITLES B
ON        B.TYPE = 'BUSINESS' AND
          A.PUB_ID=B.PUB_ID
WHERE B.PUB_ID IS NULL

Change
SELECT TITLE
FROM TITLES
WHERE NOT EXISTS
   (SELECT TITLE_ID
   FROM SALES
   WHERE TITLE_ID = TITLES.TITLE_ID)

TO
SELECT TITLE
FROM TITLES LEFT JOIN SALES
ON SALES.TITLE_ID = TITLES.TITLE_ID
WHERE SALES.TITLE_ID IS NULL


2. Use INNER JOIN instead of IN, EXISTS if there are no duplicates in subquery

Change
SELECT PUB_NAME
FROM PUBLISHERS
WHERE PUB_ID IN
   (SELECT PUB_ID
   FROM TITLES
   WHERE TYPE = 'BUSINESS')

To
SELECT DISTINCT A.PUB_NAME
FROM PUBLISHERS A INNER JOIN TITLES B
ON        B.TYPE = 'BUSINESS' AND
          A.PUB_ID=B. PUB_ID


3. Use EXISTS instead of IN in subquery

Change
SELECT PUB_NAME
FROM PUBLISHERS
WHERE PUB_ID IN
   (SELECT PUB_ID
   FROM TITLES
   WHERE TYPE = 'BUSINESS')

To
SELECT PUB_NAME
FROM PUBLISHERS
WHERE EXISTS
   (SELECT 1
   FROM TITLES
   WHERE TYPE = 'BUSINESS' AND
   PUB_ID= PUBLISHERS.PUB_ID)
Posted by at No comments:

Sunday, September 8, 2013

Survival Analysis Part 4: The Proportional Hazards Assumption and Leverage (Based on the book by Hosmer, Lemeshow, May)

rm(list = ls())
setwd("C:\\ASA")
# import data
whas500 = read.table("C:\\ASA\\whas500.DAT", sep = "", header = F)
names(whas500) = c("id",
                   "age",
                   "gender",
                   "hr",
                   "sysbp",
                   "diasbp",
                   "bmi",
                   "cvd",
                   "afb",
                   "sho",
                   "chf",
                   "av3",
                   "miord",
                   "mitype",
                   "year",
                   "admitdate",
                   "disdate",
                   "fdate",
                   "los",
                   "dstat",
                   "lenfol",
                   "fstat")
whas500$bmifp1 = (whas500$bmi/10)^2
whas500$bmifp2 = (whas500$bmi/10)^3
library(survival)
model.coxph = coxph( Surv(lenfol, fstat)~ bmifp1
                      + bmifp2
                      + age
                      + hr
                      + diasbp
                      + gender
                      + chf
                      + age*gender,
                      data = whas500,
                      method = "breslow")
summary(model.coxph)
# Table 6.1
zph.1 = cox.zph(model.coxph, transform = "identity")
zph.1
zph.log = cox.zph(model.coxph, transform = "log")
zph.log
zph.km = cox.zph(model.coxph, transform = "km")
zph.km
zph.rank = cox.zph(model.coxph, transform = "rank")
zph.rank
# Figure 6.1
plot(zph.log[4])
abline(h=0, lty=3)
plot(zph.1[4])
abline(h=0, lty=3)
plot(zph.km[4])
abline(h=0, lty=3)
plot(zph.rank[4])
abline(h=0, lty=3)

# Leverage: Figure 6.4 and Figure 6.5
score = residuals(model.coxph, type="score")
plot(whas500$bmi, score[,1], ylab="Score Residuals", xlab="bmifp1")
plot(whas500$hr, score[,4], ylab="Score Residuals", xlab="heart rate")
plot(whas500$age, score[,3], ylab="Score Residuals", xlab="age")
Posted by at No comments:

Saturday, September 7, 2013

Survival Analysis Part 3: Quartile Design Variable and Fractional Polynomials (Based on the book by Hosmer, Lemeshow, May)


rm(list = ls())
setwd("C:\\ASA")
# import data
whas500 = read.table("C:\\ASA\\whas500.DAT", sep = "", header = F)
names(whas500) = c("id",
                   "age",
                   "gender",
                   "hr",
                   "sysbp",
                   "diasbp",
                   "bmi",
                   "cvd",
                   "afb",
                   "sho",
                   "chf",
                   "av3",
                   "miord",
                   "mitype",
                   "year",
                   "admitdate",
                   "disdate",
                   "fdate",
                   "los",
                   "dstat",
                   "lenfol",
                   "fstat")
library(survival)
model1.coxph = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + bmi + gender + chf,
                      data = whas500, method = "breslow")
# Table 5.6
summary(model1.coxph)
summary(whas500$bmi)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 13.05   23.22   25.95   26.61   29.39   44.84
whas500$bmicat = cut(whas500$bmi, c(13, 23.22, 25.95, 29.39, 44.84))
summary(whas500$bmicat)
# (13,23.2] (23.2,25.9] (25.9,29.4] (29.4,44.8]
# 125         126         124         125
model2.coxph = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + bmicat + gender + chf,
                      data = whas500, method = "breslow")
summary(model2.coxph)
bmi.coeff = data.frame(model2.coxph$coefficients)[4:6,]
(summary(whas500$bmi)[[1]] + summary(whas500$bmi)[[2]])/2
midpt = c((summary(whas500$bmi)[[1]] + summary(whas500$bmi)[[2]])/2,
          (summary(whas500$bmi)[[2]] + summary(whas500$bmi)[[3]])/2,
          (summary(whas500$bmi)[[3]] + summary(whas500$bmi)[[5]])/2,
          (summary(whas500$bmi)[[5]] + summary(whas500$bmi)[[6]])/2
          )
bmi.plot = cbind(midpt, rbind(0,data.frame(bmi.coeff)))
bmi.plot
# Figure 5.1
plot(bmi.plot[,1], bmi.plot[,2], ylab = "Log Hazard", xlab = "BMI", type = "b")
# Table 5.8
library(mfp)
model.null = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf,
                    data = whas500, method = "breslow")
-2*model.null$loglik
model.lin = coxph( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf + bmi,
                   data = whas500, method = "breslow")
-2*model.lin$loglik
# J = 1, df = 2
model.J1 = mfp( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf + fp(bmi, df = 2),
                  data = whas500, method = "breslow", family = cox)
-2*model.J1$loglik
model.J1$powers
# J = 2, df = 4
model.J2 = mfp( Surv(lenfol, fstat)~ age + hr + diasbp + gender + chf + fp(bmi, df = 4),
                data = whas500, method = "breslow", family = cox)
-2*model.J2$loglik
model.J2$powers
# Figure 5.2
whas500$resid = residuals(model.null, type="martingale", data=whas500)
plot(whas500$bmi, whas500$resid, xlab="BMI", ylab="Martingale Residuals")
lines(lowess(whas500$bmi, whas500$resid))
Posted by at No comments:

Thursday, September 5, 2013

Survival Analysis Part 2: Cox Model and Applications (Based on the book by Hosmer, Lemeshow, May)

rm(list = ls())
setwd("C:\\ASA")
# import data
actg320 = read.table("C:\\ASA\\actg320.DAT", sep = "", header = F)
names(actg320) = c("id",
                   "time",
                   "censor",
                   "time_d",
                   "censor_d",
                   "tx",
                   "txgrp",
                   "strat2",
                   "sex",
                   "raceth",
                   "ivdrug",
                   "hemophil",
                   "karnof",
                   "cd4",
                   "priorzdv",
                   "age")
library(survival)
tx.coxph = coxph( Surv(time, censor)~tx, data = actg320, method = "breslow")
# Table 3.1
summary(tx.coxph)
v2.coxph = coxph( Surv(time, censor)~tx + age + sex + cd4 + priorzdv,
                  data = actg320, method = "breslow")
# Table 3.2
summary(v2.coxph)
v3.coxph = coxph( Surv(time, censor)~tx + age + cd4,
                  data = actg320, method = "breslow")
# Table 3.3
summary(v3.coxph)
v4.coxph = coxph( Surv(time, censor)~tx + age + cd4,
                  data = actg320, method = "efron")
summary(v4.coxph)
# baseline cumulative hazard function
H = basehaz(v4.coxph, centered = TRUE)
plot(H[,2], H[,1], xlab = "Time", ylab = "cumulative hazard function", type = "l")
###################################################################
###################################################################
gbcs = read.table("C:\\ASA\\gbcs.DAT", sep = "", header = F)
names(gbcs) = c("id",
                "diagdate",
                "recdate",
                "deathdate",
                "age",
                "menopause",
                "hormone",
                "size",
                "grade",
                "nodes",
                "prog_recp",
                "estrg_recp",
                "rectime",
                "censrec",
                "survtime",
                "censdead")
crude = coxph( Surv(rectime, censrec)~hormone, data = gbcs, method = "breslow")
summary(crude)
adjusted = coxph( Surv(rectime, censrec)~hormone + size, data = gbcs, method = "breslow")
summary(adjusted)
interaction = coxph( Surv(rectime, censrec)~hormone + size + hormone*size,
                     data = gbcs, method = "breslow")
summary(interaction)
hormone.new = data.frame(hormone=c(0,1))
model = coxph( Surv(rectime, censrec)~hormone, data = gbcs, method = "breslow")
plot(survfit(model, newdata=hormone.new), xlab="Time", ylab="Survival Probability")
points(survfit(model, newdata=hormone.new), pch=c(1,2))
legend(200, 0.5, c("Hormone Absent", "Hormone Present"), pch=c(1,2), bg="transparent")
Posted by at No comments:

Wednesday, September 4, 2013

Survival Analysis Part 1: Descriptive Methods for Survival Data (Based on the book by Hosmer, Lemeshow, May)

rm(list = ls())
setwd("C:\\ASA")
# import data
WHAS100 = read.table("C:\\ASA\\whas100.DAT", sep = "", header = F)
names(WHAS100) = c("id",
                   "admitdate",
                   "foldate",
                   "los",
                   "lenfol",
                   "fstat",
                   "age",
                   "gender",
                   "bmi")
head(WHAS100)
# Kaplan-Meier estimate
library(survival)
whas100.surv = survfit(Surv(lenfol, fstat) ~ 1, conf.type = "none", data = WHAS100)
# Table 2.3
summary(whas100.surv)
# Figure 2.2
plot(whas100.surv, xlab="Time", ylab="survival probability")

# Life-Table estimator
t1year = floor(WHAS100$lenfol/365)
fstat = WHAS100$fstat
t1 = data.frame(t1year, fstat)
t1
library(nlme)
die = gsummary(t1, sum, groups = t1year)
total = gsummary(t1, length, groups = t1year)
rm(t1year)
ltab.data = cbind(die[,1:2], total[,2])
names(ltab.data) = c("t1year", "censor", "total")
ltab.data
attach(ltab.data)
lt = length(t1year)
t1year[lt + 1] = NA
die = censor
censored = total[,2] - censor
library(KMsurv)
lifetable = lifetab(t1year, 100, censored, die)
# Table 2.4
lifetable[,1:5]
# nlost = censored
# nevent = die
# Figure 2.3
plot(t1year[1:8], lifetable[,5], type="s",
     xlab="Survival Time (Years)", ylab="Estimated Survival Probability")
# Figure 2.4
plot(t1year[1:8], lifetable[,5], type="b",
     xlab="Survival Time (Years)", ylab="Estimated Survival Probability")
################################
## confidence interval  ########
################################

whas100.surv = survfit(Surv(lenfol, fstat) ~ 1, conf.type = "none", data = WHAS100)
library(km.ci)
ci.log = km.ci(whas100.surv, conf.level = 0.95, tl = NA, tu = NA, method = "log")
ci.loglog = km.ci(whas100.surv, conf.level = 0.95, tl = NA, tu = NA, method = "loglog")
# ci.hallwellner = km.ci(whas100.surv, conf.level = 0.95, tl = NA, tu = NA, method = "hall-wellner")
ci.hallwellner = km.ci(whas100.surv, conf.level = 0.95, tl = NA, tu = NA, method = "loghall")
# look like loghall can recreate the Figure 2.7 instead of hall-wellner
plot(whas100.surv)
plot(ci.log)
plot(ci.loglog)
plot(ci.hallwellner)
log.ci = survfit(Surv(lenfol, fstat) ~ 1, conf.type = "log", data = WHAS100)
loglog.ci = survfit(Surv(lenfol, fstat) ~ 1, conf.type = "log-log", data = WHAS100)
# Figure 2.7
par(cex = 0.8)
plot(ci.hallwellner, lty = 3, lwd = 2)
lines(log.ci, lwd = 2, lty = 1)
lines(log.ci, lwd = 2, lty = 6, conf.int = T)
legend(150, 0.4, c("Kaplan-Meier", "Hall-Wellner", "Pointwise"), lty=c(1, 4, 6))
#######################################
########  Comparison  #################
#######################################
# Table 2.12
diff.logrank = survdiff(Surv(lenfol, fstat) ~ gender, data = WHAS100, rho = 0)
print(diff.logrank)
diff.petopeto = survdiff(Surv(lenfol, fstat) ~ gender, data = WHAS100, rho = 1)
print(diff.petopeto)
# Figure 2.9
diff.drug = survfit(Surv(lenfol, fstat) ~ strata(gender), data = WHAS100, conf.type="log-log")
plot(diff.drug, lty = c(1, 2), xlab = "Time", ylab = "Survival Probability")
diff.drug
legend(200,0.4, c("Male", "Female"), lty = c(1,2))
agecat = cut(WHAS100$age, c(0, 59, 69, 79, 100))
agecat
diff.age = survfit(Surv(lenfol, fstat) ~ strata(agecat), data = WHAS100, conf.type="log-log")
# Table 2.13
print(diff.age)
# Table 2.15
survdiff(Surv(lenfol, fstat) ~ agecat, data = WHAS100, rho = 0)
plot(diff.age, lty=c(1, 2, 3, 4), xlab="Time", ylab="Survival Probability")
legend(200,0.6, lty=c(1, 2, 3, 4),
       c("Age<60", "60<=Age<69", "70<=Age<79", "Age>=80"), bg="transparent")

############################################################
#####  aalen  ################################################
############################################################
aalen = survfit(coxph(Surv(lenfol, fstat) ~ 1, data = WHAS100), type = "aalen")
summary(aalen)
H.aalen = -log(aalen$surv)
aalen.est = cbind(time=aalen$time, d=aalen$n.event, n=aalen$n.risk, H.aalen, s1 = aalen$surv)
KM = survfit(Surv(lenfol, fstat) ~ 1, data = WHAS100)
KM.est = cbind(time=KM$time, s2 = KM$surv)
all = merge(data.frame(aalen.est), data.frame(KM.est), by="time")
all
# Figure 2.12
plot(all$time, all$s1, type="s", xlab="Time", ylab="Survival Probability")
points(all$time, all$s1, pch=1)
lines(all$time, all$s2, type="s")
points(all$time, all$s2, pch=3)
legend(200,0.6, c("Nelson-Aalen", "Kaplan-Meier"), pch=c(1,3))
hazard.rate = all$d/all$n
plot(all$time, hazard.rate, type = "p", pch=20, xlab="Time", ylab="Hazard Ratio")
lines(lowess(all$time, hazard.rate, f=0.75, iter=5))
Posted by at No comments:
Subscribe to: Posts (Atom)