~ubuntu-branches/ubuntu/trusty/r-cran-rms/trusty-proposed

« back to all changes in this revision

Viewing changes to R/survplot.survfit.s

• browse files at revision 12
• compare with another revision
• download diff
• download tarball
• view history from revision 12

Show diffs side-by-side

added

removed

R/survplot.survfit.s

1

1

survplot.survfit <-

2

2

  function(fit, xlim, 

3

3

           ylim, xlab, ylab, time.inc,

4

 

           conf=c("bands","bars","none"), add=FALSE, 

 

4

           conf=c("bands","bars","diffbands","none"), add=FALSE, 

5

5

           label.curves=TRUE,

6

6

           abbrev.label=FALSE, levels.only=FALSE,

7

 

           lty,lwd=par('lwd'),

 

7

           lty, lwd=par('lwd'),

8

8

           col=1, col.fill=gray(seq(.95, .75, length=5)),

9

9

           loglog=FALSE, fun, n.risk=FALSE, logt=FALSE,

10

10

           dots=FALSE, dotsize=.003, grid=NULL,

16

16

  conf.int <- fit$conf.int

17

17

  if(!length(conf.int) | conf=="none") conf.int <- 0

18

18

 

 

19

  opar <- par(c('mar', 'xpd'))

 

20

  on.exit(par(opar))

 

21

 

 

22

  fit.orig <- fit

19

23

  units <- fit$units

20

24

  if(!length(units)) units <- "Day"

21

25

  maxtime <- fit$maxtime

96

100

  yd <- ylim[2] - ylim[1]

97

101

  

98

102

  if(n.risk && !add) {

99

 

    mar <- par()$mar

100

 

    if(mar[4]<4) {mar[4] <- mar[4]+2; par(mar=mar)}

 

103

    mar <- opar$mar

 

104

    if(mar[4] < 4) {mar[4] <- mar[4] + 2; par(mar=mar)}

101

105

  }

 

106

  

102

107

  ## One curve for each value of y, excl style used for C.L.

103

108

  lty <- if(missing(lty)) seq(ns+1)[-2] else rep(lty, length=ns)

104

109

  lwd <- rep(lwd, length=ns)

107

112

  if(labelc || conf=='bands') curves <- vector('list',ns)

108

113

  Tim <- Srv <- list()

109

114

  

110

 

  oxpd <- par('xpd')

111

115

  par(xpd=NA)

112

 

  on.exit(par(xpd=oxpd))

113

116

 

114

117

  for(i in 1:ns) {

115

118

    st <- stemp == i

120

123

    if(i==1 & !add) {

121

124

      plot(time, surv, xlab=xlab, xlim=xlim,

122

125

           ylab=ylab, ylim=ylim, type="n", axes=FALSE)

123

 

      mgp.axis(1,at=if(logt)pretty(xlim) else

 

126

      mgp.axis(1, at=if(logt) pretty(xlim) else

124

127

               seq(xlim[1], max(pretty(xlim)), time.inc),

125

128

               labels=TRUE)

126

129

      

138

141

                  ylim[1] + if(n.risk && missing(y.n.risk)) yi else 0,

139

142

                  by = -yi)

140

143

        if(dots)

141

 

          for(tt in xp)symbols(rep(tt, length(yp)), yp,

142

 

                               circles=rep(dotsize, length(yp)),

143

 

                               inches=dotsize, add=TRUE)

 

144

          for(tt in xp)

 

145

            symbols(rep(tt, length(yp)), yp,

 

146

                    circles=rep(dotsize, length(yp)),

 

147

                    inches=dotsize, add=TRUE)

144

148

        else abline(h=yp, v=xp, col=grid, xpd=FALSE)

145

149

      }

146

150

    }

168

172

      }

169

173

    }

170

174

    if(logt) {

171

 

      if(conf != 'bands')

 

175

      if(conf %nin% c('bands', 'diffbands'))

172

176

        lines(tim, srv, type="s", lty=lty[i], col=col[i], lwd=lwd[i])

173

 

      if(labelc || conf=='bands') curves[[i]] <- list(tim,srv)

 

177

      if(labelc || conf %in% c('bands', 'diffbands'))

 

178

        curves[[i]] <- list(tim,srv)

174

179

    }

175

180

          else {

176

181

      xxx <- c(mintime,tim)

177

182

      yyy <- c(fun(1),srv)

178

 

      if(conf != 'bands')

 

183

      if(conf %nin% c('bands', 'diffbands'))

179

184

        lines(xxx, yyy, type="s", lty=lty[i], col=col[i], lwd=lwd[i])

180

 

      if(labelc || conf=='bands') curves[[i]] <- list(xxx, yyy)

 

185

      if(labelc || conf %in% c('bands', 'diffbands'))

 

186

        curves[[i]] <- list(xxx, yyy)

181

187

    }

182

188

    if(pr) {

183

 

      zest <- rbind(time[s],surv[s])

184

 

      dimnames(zest) <- list(c("Time","Survival"),

185

 

                             rep("",sum(s)))

186

 

      if(slevp)cat("\nEstimates for ", slev[i],"\n\n")

 

189

      zest <- rbind(time[s], surv[s])

 

190

      dimnames(zest) <- list(c("Time", "Survival"),

 

191

                             rep("", sum(s)))

 

192

      if(slevp)cat("\nEstimates for ", slev[i], "\n\n")

187

193

      print(zest, digits=3)

188

194

    }

189

195

    if(conf.int > 0) {

190

 

      if(conf=="bands") {

 

196

      if(conf=='bands') {

191

197

        if(logt)

192

198

          polyg(x = c(tim, max(tim), rev(tim)),

193

199

                y = c(blower, rev(bupper), max(bupper)),

197

203

                y = c(fun(1), blower, rev(c(fun(1), bupper))),

198

204

                col = col.fill[i], type = "s")

199

205

      }

 

206

      else if(conf == 'diffbands')

 

207

        survdiffplot(fit.orig, conf=conf, fun=fun)

 

208

 

200

209

      else {

201

210

        j <- if(ns ==1) TRUE else vs == olev[i]

202

211

        tt <- v$time[j]  #may not get predictions at all t

222

231

      text(tt[1], yy, nri[1], cex=cex.n.risk,

223

232

           adj=adj.n.risk, srt=srt.n.risk)

224

233

      text(tt[-1], yy, nri[-1], cex=cex.n.risk, adj=1)

225

 

      if(slevp)text(xlim[2]+xd*.025,

 

234

      if(slevp) text(xlim[2] + xd * .025,

226

235

                    yy, adj=0, sleva[i], cex=cex.n.risk)

227

236

    }

228

237

  }

229

238

  

230

 

  if(conf=='bands') for(i in 1:ns)

231

 

    lines(curves[[i]][[1]], curves[[i]][[2]],

232

 

          lty=lty[i], lwd=lwd[i], col=col[i], type='s')

 

239

  if(conf %in% c('bands', 'diffbands'))

 

240

    for(i in 1:ns)

 

241

      lines(curves[[i]][[1]], curves[[i]][[2]],

 

242

            lty=lty[i], lwd=lwd[i], col=col[i], type='s')

 

243

 

233

244

  if(labelc) labcurve(curves, sleva, type='s', lty=lty, lwd=lwd,

234

245

                      opts=label.curves, col.=col)

235

246

  

238

249

 

239

250

 

240

251

survdiffplot <-

241

 

  function(fit, order=1:2, xlim, ylim, xlab,

 

252

  function(fit, order=1:2, fun=function(y) y, xlim, ylim, xlab,

242

253

           ylab="Difference in Survival Probability", time.inc,

243

 

           conf.int, conf=c("shaded", "bands", "none"),

 

254

           conf.int, conf=c("shaded", "bands", "diffbands", "none"),

244

255

           add=FALSE, lty=1, lwd=par('lwd'), col=1,

245

256

           n.risk=FALSE,  grid=NULL,

246

257

           srt.n.risk=0, adj.n.risk=1,

248

259

 

249

260

  conf <- match.arg(conf)

250

261

  if(missing(conf.int)) conf.int <- fit$conf.int

251

 

  if(!length(conf.int) | conf=="none") conf.int <- 0

 

262

  if(! length(conf.int) | conf == "none") conf.int <- 0

252

263

 

 

264

  opar <- par(c('xpd', 'mar'))

 

265

  on.exit(par(opar))

 

266

  

253

267

  units <- fit$units

254

268

  if(!length(units)) units <- "Day"

255

269

  maxtime <- fit$maxtime

264

278

  }

265

279

  if(n.risk && !length(fit$n.risk)) {

266

280

    n.risk <- FALSE

267

 

    warning("fit does not have number at risk\nIs probably from a parametric model\nn.risk set to F")

 

281

    warning("fit does not have number at risk\nIs probably from a parametric model\nn.risk set to FALSE")

268

282

    }

269

283

 

270

284

  if(missing(xlab)) xlab <- if(units==' ') 'Time' else paste(units, "s", sep="")

277

291

  polyg <- ordGridFun(grid=FALSE)$polygon

278

292

 

279

293

  times <- sort(unique(c(fit$time, seq(mintime, maxtime, by=time.inc))))

280

 

  

 

294

 

 

295

  ## Note: summary.survfit computes standard errors on S(t) scale

281

296

  f <- summary(fit, times=times)

282

297

  

283

298

  slev <- levels(f$strata)

298

313

  b <- h(slev[order[2]], times, f)

299

314

 

300

315

 

301

 

  surv  <- a$surv - b$surv

 

316

  surv  <- if(conf == 'diffbands') (fun(a$surv) + fun(b$surv)) / 2

 

317

   else fun(a$surv) - fun(b$surv)

302

318

  se    <- sqrt(a$se^2 + b$se^2)

303

319

 

304

320

  z  <- qnorm((1 + conf.int) / 2)

 

321

  if(conf == 'diffbands') {

 

322

    lo <- surv - 0.5 * z * se

 

323

    hi <- surv + 0.5 * z * se

 

324

    k <- !is.na(times + lo + hi)

 

325

    polyg(c(times[k], rev(times[k])), c(lo[k], rev(hi[k])),

 

326

           col=gray(.9), type='s')

 

327

    return(invisible(slev))

 

328

  }

305

329

  lo <- surv - z * se

306

330

  hi <- surv + z * se

307

331

 

333

357

           lines(times, lo, col=gray(.7))

334

358

           lines(times, hi, col=gray(.7))

335

359

         },

 

360

         diffbands=NULL,

336

361

         none=NULL)

337

362

  lines(times, surv, type='s', lwd=lwd, col=col)

338

363

  abline(h=0, col=gray(.7))

348

373

    yd         <- ylim[2] - ylim[1]

349

374

    

350

375

    if(!add) {

351

 

      mar <- par()$mar

 

376

      mar <- opar$mar

352

377

      if(mar[4] < 4) {mar[4] <- mar[4] + 2; par(mar=mar)}

353

378

    }

354

 

    oxpd <- par('xpd')

355

379

    par(xpd=NA)

356

 

    on.exit(par(xpd=oxpd))

357

380

    

358

381

    tt <- nrisktimes

359

382

    tt[1] <- xlim[1]

• Older »

Loggerhead is a web-based interface for Breezy
Version: 2.0.1