1
2svykm<-function(formula, design, se=FALSE, ...) UseMethod("svykm",design)
3
4svykm.survey.design<-function(formula, design,se=FALSE, ...){
5  if (!inherits(formula,"formula")) stop("need a formula")
6  if (length(formula)!=3) stop("need a two-sided formula")
7  mf<-model.frame(formula, model.frame(design), na.action=na.pass)
8  mf<-na.omit(mf)
9  drop<-attr(mf,"na.action")
10  if (!is.null(drop))
11    design<-design[-drop,]
12  y<-model.response(mf)
13  if (!is.Surv(y) || attr(y,"type")!="right")
14    stop("response must be a right-censored Surv object")
15
16  if (ncol(mf)==1) {
17    if (se)
18      s<-km.stderr(y,design)
19    else
20      s<-svykm.fit(y,weights(design))
21  } else {
22    x<-mf[,-1]
23    if (NCOL(x)>1)
24      groups<-do.call(interaction,x)
25    else
26      groups<-as.factor(x)
27
28    if (se){
29      lhs<-formula
30      lhs[[3]]<-1
31      s<-lapply(levels(groups), function(g) svykm(lhs, subset(design,groups==g),se=TRUE))
32    }else{
33      s<-lapply(levels(groups), function(g) svykm.fit(y[groups==g],weights(design)[groups==g]))
34    }
35    names(s)<-levels(groups)
36    class(s)<-"svykmlist"
37  }
38  call<-match.call()
39  call[[1]]<-as.name(.Generic)
40  attr(s,"call")<-call
41  attr(s, "formula")<-formula
42  attr(s, "na.action")<-drop
43  return(s)
44}
45
46
47svykm.svyrep.design<-function(formula, design,se=FALSE, ...){
48  if (!inherits(formula,"formula")) stop("need a formula")
49  if (length(formula)!=3) stop("need a two-sided formula")
50  mf<-model.frame(formula, model.frame(design), na.action=na.pass)
51  mf<-na.omit(mf)
52  drop<-attr(mf,"na.action")
53  if (!is.null(drop))
54    design<-design[-drop,]
55  y<-model.response(mf)
56  if (!is.Surv(y) || attr(y,"type")!="right")
57    stop("response must be a right-censored Surv object")
58
59  if (ncol(mf)==1) {
60    if (se)
61      stop("SE not yet available")
62    else
63      s<-svykm.fit(y,weights(design,"sampling"))
64  } else {
65    x<-mf[,-1]
66    if (NCOL(x)>1)
67      groups<-do.call(interaction,x)
68    else
69      groups<-as.factor(x)
70
71    if (se){
72      lhs<-formula
73      lhs[[3]]<-1
74      s<-lapply(levels(groups), function(g) svykm(lhs, subset(design,groups==g),se=TRUE))
75    }else{
76      s<-lapply(levels(groups), function(g) svykm.fit(y[groups==g],weights(design)[groups==g]))
77    }
78    names(s)<-levels(groups)
79    class(s)<-"svykmlist"
80  }
81  call<-match.call()
82  call[[1]]<-as.name(.Generic)
83  attr(s,"call")<-call
84  attr(s, "formula")<-formula
85  attr(s, "na.action")<-drop
86  return(s)
87}
88
89
90svykm.fit<-function(y,w){
91  t<-y[,"time"]
92  s<-y[,"status"]
93  nn<-rowsum(cbind(s,1)*w,t)
94  tt<-sort(unique(t))
95  N<-c(sum(w),sum(w),sum(w)-cumsum(nn[-nrow(nn),2]))
96  d<-c(0,nn[,1])
97  surv<-pmax(0,cumprod(1-d/N))
98  rval<-list(time=c(0,tt), surv=surv)
99  class(rval)<-"svykm"
100  rval
101}
102
103km.stderr<-function(survobj,design){
104  time<-survobj[,'time']
105  status<-survobj[,'status']
106  ## Brute force and ignorance: compute Y and dN as totals, use delta-method
107  keep<-which((status==1) & (weights(design)!=0))
108  y<-outer(time,time[keep],">=")
109  dN<-diag(status)[,keep,drop=FALSE]
110  oo<-order(time[keep], -status[keep])
111  okeep<-keep[oo]
112  ntimes<-length(oo)
113  ttime<-time[okeep]
114  sstatus<-status[okeep]
115
116  totals<-svytotal(cbind(dN[,oo,drop=FALSE],y[,oo,drop=FALSE]), design)
117  rm(dN)
118  y<-coef(totals)[-(1:ntimes)]
119  dNbar<-coef(totals)[1:ntimes]
120
121  h<-cumsum(dNbar/y)
122
123  dVn<- vcov(totals)[(1:ntimes),(1:ntimes)]/outer(y,y)
124  dVy <- vcov(totals)[-(1:ntimes),-(1:ntimes)]*outer(dNbar/y^2,dNbar/y^2)
125  dCVny<- -vcov(totals)[(1:ntimes),-(1:ntimes)]*outer(1/y,dNbar/y^2)
126  dV<-dVn+dVy+dCVny+t(dCVny)
127
128  V<-numeric(ntimes)
129  if (ntimes>0)
130      V[1]<-dV[1,1]
131  if (ntimes>1)
132      for(i in 2:ntimes) V[i]<-V[i-1]+sum(dV[1:(i-1),i])+sum(dV[i,1:i])
133
134  rval<-list(time=ttime,surv=exp(-h),varlog=V)
135  class(rval)<-"svykm"
136  rval
137}
138
139
140plot.svykm<-function(x,xlab="time",ylab="Proportion surviving",ylim=c(0,1),ci=NULL,lty=1,...){
141  if (is.null(ci))
142    ci<-!is.null(x$varlog)
143
144  plot(x$time,x$surv,xlab=xlab,ylab=ylab, type="s",ylim=ylim,lty=lty,...)
145
146  if (ci){
147    if (is.null(x$varlog))
148      warning("No standard errors available in object")
149    else{
150      lines(x$time,exp(log(x$surv)-1.96*sqrt(x$varlog)),lty=2,type="s",...)
151      lines(x$time,pmin(1,exp(log(x$surv)+1.96*sqrt(x$varlog))),lty=2,type="s",...)
152    }
153  }
154  invisible(x)
155}
156
157
158lines.svykm<-function(x,xlab="time",type="s",ci=FALSE,lty=1,...){
159  lines(x$time,x$surv, type="s",lty=lty,...)
160  if (ci){
161    if (is.null(x$varlog))
162      warning("no standard errors available in object")
163    else {
164      lines(x$time,exp(log(x$surv)-1.96*sqrt(x$varlog)),lty=2,type="s",...)
165      lines(x$time,pmin(1,exp(log(x$surv)+1.96*sqrt(x$varlog))),lty=2,type="s",...)
166    }
167  }
168  invisible(x)
169}
170
171plot.svykmlist<-function(x, pars=NULL, ci=FALSE,...){
172  if (!is.null(pars)) pars<-as.data.frame(pars,stringsAsFactors=FALSE)
173
174  if(is.null(pars))
175    plot(x[[1]],ci=ci,...)
176  else
177    do.call(plot,c(list(x[[1]]),pars[1,,drop=FALSE],ci=ci,...))
178
179  m<-length(x)
180  if(m==1) return()
181  for(i in 2:m){
182    if(is.null(pars))
183      lines(x[[i]],ci=ci,...)
184    else
185      do.call(lines,c(list(x[[i]]),pars[i,,drop=FALSE],ci=ci,...))
186  }
187  invisible(x)
188}
189
190print.svykm<-function(x, digits=3,...,header=TRUE){
191  if (header) {cat("Weighted survival curve: ")
192               print(attr(x,"call"))}
193  suppressWarnings({iq1<-min(which(x$surv<=0.75))
194                    iq2<-min(which(x$surv<=0.5))
195                    iq3<-min(which(x$surv<=0.25))})
196  if (is.finite(iq1)) q1<-x$time[iq1] else q1<-Inf
197  if (is.finite(iq2)) q2<-x$time[iq2] else q2<-Inf
198  if (is.finite(iq3)) q3<-x$time[iq3] else q3<-Inf
199 cat("Q1 =",round(q1,digits)," median =",round(q2,digits)," Q3 =",round(q3,digits),"\n")
200 invisible(x)
201}
202
203print.svykmlist<-function(x, digits=3,...){
204  cat("Weighted survival curves:\n")
205  print(attr(x,"call"))
206  for(i in 1:length(x)){
207    cat(names(x)[i],": ")
208    print(x[[i]],digits=digits,header=FALSE)
209  }
210  invisible(x)
211}
212
213quantile.svykm<-function(x, probs=c(0.75,0.5,0.25),ci=FALSE,level=0.95,...){
214
215  iq<-sapply(probs, function(p) suppressWarnings(min(which(x$surv<=p))))
216  qq<-sapply(iq, function(i) if (is.finite(i)) x$time[i] else Inf)
217  names(qq)<-probs
218  if (ci){
219    if(is.null(x$varlog)){
220      warning("no confidence interval available.")
221    } else {
222      halfalpha<-(1-level)/2
223      z<-qnorm(halfalpha, lower.tail=FALSE)
224      su<-exp(log(x$surv)+z*sqrt(x$varlog))
225      iu<-sapply(probs, function(p) suppressWarnings(min(which(su<=p))))
226      qu<-sapply(iu, function(i) if (is.finite(i)) x$time[i] else Inf)
227      sl<-exp(log(x$surv)-z*sqrt(x$varlog))
228      il<-sapply(probs, function(p) suppressWarnings(min(which(sl<=p))))
229      ql<-sapply(il, function(i) if (is.finite(i)) x$time[i] else Inf)
230      ci<-cbind(ql,qu)
231      rownames(ci)<-probs
232      colnames(ci)<-format(c(halfalpha,1-halfalpha),3)
233      attr(qq,"ci")<-ci
234    }
235  }
236  qq
237}
238
239
240confint.svykm<-function(object, parm, level=0.95,...){
241  if (is.null(object$varlog)) stop("no standard errors in object")
242
243  parm<-as.numeric(parm)
244  idx<-sapply(parm, function(t) max(which(object$time<=t)))
245  z<-qnorm((1-level)/2)
246  ci<-exp(log(object$surv[idx])+outer(sqrt(object$varlog[idx]),c(z,-z)))
247  ci[,2]<-pmin(ci[,2],1)
248  rownames(ci)<-parm
249  colnames(ci)<-format( c((1-level)/2, 1-(1-level)/2),3)
250  ci
251}
252
253
254predict.svycoxph<-function(object, newdata, se=FALSE,
255                           type=c("lp", "risk", "expected", "terms","curve"),
256                           ...){
257
258  type<-match.arg(type)
259  if(type!="curve") return(NextMethod())
260
261  design<-object$survey.design
262  response<-object$y
263
264  if (!is.null(attr(terms(object), "specials")$strata))
265    stop("Stratified models are not supported yet")
266
267  if (attr(response,"type")=="counting"){
268    time<-object$y[,2]
269    status<-object$y[,'status']
270    entry<-object$y[,1]
271  } else if (attr(response,'type')=="right"){
272    time<-object$y[,"time"]
273    status<-object$y[,"status"]
274    entry<-rep(-Inf,length(time))
275  } else stop("unsupported survival type")
276  if(is.null(object$na.action)){
277    design<-object$survey.design
278  } else {
279    design<-object$survey.design[-object$na.action,]
280  }
281
282  ff<-delete.response(terms(formula(object)))
283  zmf<-model.frame(ff, newdata)
284  z.pred<-model.matrix(ff, zmf)[,-1,drop=FALSE]
285
286##
287## The simple case first
288##
289  risk<-getS3method("predict","coxph")(object,type="risk",se.fit=FALSE)
290  if(se==FALSE){
291    tt<-c(time,entry)
292    ss<-c(status,rep(0,length(entry)))
293    ee<-c(rep(1,length(status)),rep(-1,length(entry)))
294    oo<-order(tt,-ee,-ss)
295    dN<-ss[oo]
296    w<-rep(weights(design),2)[oo]
297    risks<-rep(risk,2)
298    Y<-rev(cumsum(rev(risks[oo]*w*ee[oo])))
299    keep<-dN>0
300    s<-vector("list",nrow(z.pred))
301    beta<-coef(object)
302    h0<- cumsum( (w*dN/Y)[keep] )
303    for(i in 1:nrow(z.pred)){
304      zi<-z.pred[i,]-object$means
305      s[[i]]<-list(time=time[oo][keep],
306              surv=exp(-h0 * exp(sum(beta*zi))),
307              call=sys.call())
308    class(s[[i]])<-c("svykmcox","svykm")
309    }
310    names(s)<-rownames(newdata)
311    return(s)
312  }
313##
314## The hard case: curves with standard errors
315##
316  if(!inherits(design,"survey.design"))
317    stop("replicate-weight designs not supported yet")
318
319  keep<-which((status==1) & (weights(design)!=0))
320  y<-outer(time,time[keep],">=")*risk*outer(entry,time[keep],"<=")
321  dN<-diag(status)[,keep]
322  oo<-order(time[keep], -status[keep])
323  okeep<-keep[oo]
324  ntimes<-length(oo)
325  ttime<-time[okeep]
326  sstatus<-status[okeep]
327  totals<-svytotal(cbind(dN[,oo],y[,oo]), design)
328  rm(dN)
329
330  y<-coef(totals)[-(1:ntimes)]
331  dNbar<-coef(totals)[1:ntimes]
332  vtotals<-vcov(totals)
333  rm(totals)
334
335  h<-cumsum(dNbar/y)
336
337  dVn<- vtotals[(1:ntimes),(1:ntimes)]/outer(y,y)
338  dVy <- vtotals[-(1:ntimes),-(1:ntimes)]*outer(dNbar/y^2,dNbar/y^2)
339  dCVny<- -vtotals[(1:ntimes),-(1:ntimes)]*outer(1/y,dNbar/y^2)
340  dV<-dVn+dVy+dCVny+t(dCVny)
341
342  det<-suppressWarnings(coxph.detail(object))
343  ze<-sweep(as.matrix(det$means)[rep(1:length(det$time), det$nevent),,drop=FALSE],
344            2, object$means)
345  rm(det)
346
347  dH<-dNbar/y
348  h.ze<-dH*ze
349  varbeta<-vcov(object)
350
351  Vh<-numeric(ntimes)
352  Vh[1]<-dV[1,1]
353  for(i in 2:ntimes) Vh[i]<-Vh[i-1]+sum(dV[1:(i-1),i])+sum(dV[i,1:i])
354  dVb<-numeric(ntimes)
355  for(i in 1:ntimes) dVb[i]<-crossprod(h.ze[i,],varbeta%*%(h.ze[i,]))
356  Vb<-cumsum(dVb)
357  dCV<-matrix(nrow=ntimes,ncol=NCOL(ze))
358  for(i in 1:ntimes) dCV[i,] <- -varbeta%*%(h.ze[i,])
359  CV<-apply(dCV,2,cumsum)
360
361  V0<-Vh+Vb
362  s0<-exp(-h)
363
364  s<-vector("list",nrow(z.pred))
365  for(i in 1:nrow(z.pred)){
366    zi<-z.pred[i,]-object$means
367    riski<-exp(sum(zi*coef(object)))
368    Vz<-drop(crossprod(zi,varbeta%*%zi))*riski^2*h^2
369    CVz<-colSums(t(dCV)*zi)*riski^2*h
370
371    V<-V0*riski^2+Vz+CVz*2
372    s[[i]]<-list(time=ttime,surv=exp(-h*riski), varlog=V)
373    class(s[[i]])<-c("svykm.cox","svykm")
374  }
375  names(s)<-rownames(newdata)
376  scall<-sys.call()
377  scall[[1]]<-as.name(.Generic)
378  attr(s,"call")<-scall
379  class(s)<-c("svykmlist.cox","svykmlist")
380  return(s)
381}
382
383
384