1hzdrate.sshzd <- ## Evaluate hazard estimate 2function (object,x,se=FALSE,include=c(object$terms$labels,object$lab.p)) { 3 if (!any(class(object)=="sshzd")) 4 stop("gss error in hzdrate.sshzd: not a sshzd object") 5 if (dim(object$mf)[2]==1&is.vector(x)) { 6 x <- data.frame(x,stringsAsFactors=TRUE) 7 colnames(x) <- colnames(object$mf) 8 } 9 if (!is.null(object$d)) s <- matrix(0,dim(x)[1],length(object$d)) 10 r <- matrix(0,dim(x)[1],length(object$id.basis)) 11 for (label in include) { 12 if (label=="1") { 13 iphi <- object$terms[[label]]$iphi 14 s[,iphi] <- rep(1,dim(x)[1]) 15 next 16 } 17 if (label%in%object$lab.p) next 18 xx <- object$mf[object$id.basis,object$terms[[label]]$vlist] 19 x.new <- x[,object$terms[[label]]$vlist] 20 nphi <- object$terms[[label]]$nphi 21 nrk <- object$terms[[label]]$nrk 22 if (nphi) { 23 iphi <- object$terms[[label]]$iphi 24 phi <- object$terms[[label]]$phi 25 for (i in 1:nphi) { 26 s[,iphi+(i-1)] <- phi$fun(x.new,nu=i,env=phi$env) 27 } 28 } 29 if (nrk) { 30 irk <- object$terms[[label]]$irk 31 rk <- object$terms[[label]]$rk 32 for (i in 1:nrk) { 33 r <- r + 10^object$theta[irk+(i-1)]* 34 rk$fun(x.new,xx,nu=i,env=rk$env,out=TRUE) 35 } 36 } 37 } 38 if (!is.null(object$partial)) { 39 vars.p <- as.character(attr(object$partial$mt,"variables"))[-1] 40 facs.p <- attr(object$partial$mt,"factors") 41 vlist <- vars.p[as.logical(apply(facs.p,1,sum))] 42 for (lab in object$lab.p) { 43 if (lab%in%include) { 44 vlist.wk <- vars.p[as.logical(facs.p[,lab])] 45 vlist <- vlist[!(vlist%in%vlist.wk)] 46 } 47 } 48 if (length(vlist)) { 49 for (lab in vlist) x[[lab]] <- 0 50 } 51 matx.p <- model.matrix(object$partial$mt,x)[,-1,drop=FALSE] 52 matx.p <- sweep(matx.p,2,object$partial$center) 53 matx.p <- sweep(matx.p,2,object$partial$scale,"/") 54 nu <- length(object$d)-dim(matx.p)[2] 55 for (label in object$lab.p) { 56 nu <- nu+1 57 if (label%in%include) s[,nu] <- matx.p[,label] 58 } 59 } 60 if (is.null(object$random)) rs <- cbind(r,s) 61 else { 62 nz <- length(object$b) 63 rs <- cbind(r,matrix(0,dim(x)[1],nz),s) 64 } 65 if (!se) as.vector(exp(rs%*%c(object$c,object$b,object$d))) 66 else { 67 fit <- as.vector(exp(rs%*%c(object$c,object$b,object$d))) 68 se.fit <- .Fortran("hzdaux2", 69 as.double(object$se.aux$v), as.integer(dim(rs)[2]), 70 as.integer(object$se.aux$jpvt), 71 as.double(t(rs)), as.integer(dim(rs)[1]), 72 se=double(dim(rs)[1]), PACKAGE="gss")[["se"]] 73 list(fit=fit,se.fit=se.fit) 74 } 75} 76 77hzdcurve.sshzd <- ## Evaluate hazard curve for plotting 78function (object,time,covariates=NULL,se=FALSE) { 79 tname <- object$tname 80 xnames <- object$xnames 81 if (!any(class(object)=="sshzd")) 82 stop("gss error in hzdcurve.sshzd: not a sshzd object") 83 if (length(xnames)&&(!all(xnames%in%names(covariates)))) 84 stop("gss error in hzdcurve.sshzd: missing covariates") 85 mn <- min(object$tdomain) 86 mx <- max(object$tdomain) 87 if ((min(time)<mn)|(max(time)>mx)) 88 stop("gss error in hzdcurve.sshzd: time range beyond the domain") 89 if (length(xnames)) { 90 xx <- covariates[,xnames,drop=FALSE] 91 xy <- data.frame(matrix(0,length(time),length(xnames)+1)) 92 names(xy) <- c(tname,xnames) 93 xy[,tname] <- time 94 } 95 else xx <- NULL 96 if (!se) { 97 if (is.null(xx)) 98 zz <- hzdrate.sshzd(object,time) 99 else { 100 zz <- NULL 101 for (i in 1:dim(xx)[1]) { 102 xy[,xnames] <- xx[rep(i,length(time)),] 103 zz <- cbind(zz,hzdrate.sshzd(object,xy)) 104 } 105 zz <- zz[,,drop=TRUE] 106 } 107 zz 108 } 109 else { 110 if (is.null(xx)) 111 zz <- hzdrate.sshzd(object,time,TRUE) 112 else { 113 fit <- se.fit <- NULL 114 for (i in 1:dim(xx)[1]) { 115 xy[,xnames] <- xx[rep(i,length(time)),] 116 wk <- hzdrate.sshzd(object,xy,TRUE) 117 fit <- cbind(fit,wk$fit) 118 se.fit <- cbind(se.fit,wk$se.fit) 119 } 120 zz <- list(fit=fit[,,drop=TRUE],se.fit=se.fit[,,drop=TRUE]) 121 } 122 zz 123 } 124} 125 126survexp.sshzd <- ## Compute expected survival 127function(object,time,covariates=NULL,start=0) { 128 tname <- object$tname 129 xnames <- object$xnames 130 ## Check inputs 131 if (!any(class(object)=="sshzd")) 132 stop("gss error in survexp.sshzd: not a sshzd object") 133 if (length(xnames)&&(!all(xnames%in%names(covariates)))) 134 stop("gss error in survexp.sshzd: missing covariates") 135 lmt <- cbind(start,time) 136 if (any(lmt[,1]>lmt[,2])) 137 stop("gss error in survexp.sshzd: start after follow-up time") 138 nt <- dim(lmt)[1] 139 if (is.null(covariates)) ncov <- 1 140 else ncov <- dim(covariates)[1] 141 mn <- min(object$tdomain) 142 mx <- max(object$tdomain) 143 if ((min(start)<mn)|(max(time)>mx)) 144 stop("gss error in survexp.sshzd: time range beyond the domain") 145 ## Calculate 146 qd.hize <- 200 147 qd <- gauss.quad(2*qd.hize,c(mn,mx)) 148 gap <- diff(qd$pt) 149 g.wk <- gap[qd.hize]/2 150 for (i in 1:(qd.hize-2)) g.wk <- c(g.wk,gap[qd.hize+i]-g.wk[i]) 151 g.wk <- 2*g.wk 152 g.wk <- c(g.wk,(mx-mn)/2-sum(g.wk)) 153 gap[qd.hize:1] <- gap[qd.hize+(1:qd.hize)] <- g.wk 154 brk <- cumsum(c(mn,gap)) 155 if (is.null(covariates)) { 156 zz <- NULL 157 d.qd <- hzdcurve.sshzd(object,qd$pt) 158 for (i in 1:nt) { 159 ind <- (1:(2*qd.hize))[(qd$pt<lmt[i,2])&(qd$pt>lmt[i,1])] 160 if (length(ind)) { 161 wk <- sum(d.qd[ind]*qd$wt[ind]) 162 id.mx <- max(ind) 163 if (lmt[i,2]<brk[id.mx+1]) 164 wk <- wk-d.qd[id.mx]*qd$wt[id.mx]*(brk[id.mx+1]-lmt[i,2])/gap[id.mx] 165 else wk <- wk+d.qd[id.mx+1]*qd$wt[id.mx+1]*(lmt[i,2]-brk[id.mx+1])/gap[id.mx+1] 166 id.mn <- min(ind) 167 if (lmt[i,1]<brk[id.mn]) 168 wk <- wk+d.qd[id.mn-1]*qd$wt[id.mn-1]*(brk[id.mn]-lmt[i,1])/gap[id.mn-1] 169 else wk <- wk-d.qd[id.mn]*qd$wt[id.mn]*(lmt[i,1]-brk[id.mn])/gap[id.mn] 170 } 171 else { 172 if (lmt[i,1]<=qd$pt[1]) 173 wk <- d.qd[1]*qd$wt[1]*(lmt[i,2]-lmt[i,1])/gap[1] 174 if (lmt[i,1]>=qd$pt[2*qd.hize]) 175 wk <- d.qd[2*qd.hize]*qd$wt[2*qd.hize]*(lmt[i,2]-lmt[i,1])/gap[2*qd.hize] 176 if ((lmt[i,1]>qd$pt[1])&(lmt[i,1]<qd$pt[2*qd.hize])) { 177 i.wk <- min((1:(2*qd.hize))[qd$pt>lmt[i,1]]) 178 if (brk[i.wk]<=lmt[i,1]) 179 wk <- d.qd[i.wk]*qd$wt[i.wk]*(lmt[i,2]-lmt[i,1])/gap[i.wk] 180 if (brk[i.wk]>=lmt[i,2]) 181 wk <- d.qd[i.wk-1]*qd$wt[i.wk-1]*(lmt[i,2]-lmt[i,1])/gap[i.wk-1] 182 if ((brk[i.wk]<lmt[i,2])&(brk[i.wk]>lmt[i,1])) 183 wk <- d.qd[i.wk]*qd$wt[i.wk]*(lmt[i,2]-brk[i.wk])/gap[i.wk]+ 184 d.qd[i.wk-1]*qd$wt[i.wk-1]*(brk[i.wk]-lmt[i,1])/gap[i.wk-1] 185 } 186 } 187 zz <- c(zz,wk) 188 } 189 } 190 else { 191 zz <- NULL 192 for (j in 1:ncov) { 193 zz.wk <- NULL 194 d.qd <- hzdcurve.sshzd(object,qd$pt,covariates[j,,drop=FALSE]) 195 for (i in 1:nt) { 196 ind <- (1:(2*qd.hize))[(qd$pt<lmt[i,2])&(qd$pt>lmt[i,1])] 197 if (length(ind)) { 198 wk <- sum(d.qd[ind]*qd$wt[ind]) 199 id.mx <- max(ind) 200 if (lmt[i,2]<=brk[id.mx+1]) 201 wk <- wk-d.qd[id.mx]*qd$wt[id.mx]*(brk[id.mx+1]-lmt[i,2])/gap[id.mx] 202 else wk <- wk+d.qd[id.mx+1]*qd$wt[id.mx+1]*(lmt[i,2]-brk[id.mx+1])/gap[id.mx+1] 203 id.mn <- min(ind) 204 if (lmt[i,1]<brk[id.mn]) 205 wk <- wk+d.qd[id.mn-1]*qd$wt[id.mn-1]*(brk[id.mn]-lmt[i,1])/gap[id.mn-1] 206 else wk <- wk-d.qd[id.mn]*qd$wt[id.mn]*(lmt[i,1]-brk[id.mn])/gap[id.mn] 207 } 208 else { 209 if (lmt[i,1]<=qd$pt[1]) 210 wk <- d.qd[1]*qd$wt[1]*(lmt[i,2]-lmt[i,1])/gap[1] 211 if (lmt[i,1]>=qd$pt[2*qd.hize]) 212 wk <- d.qd[2*qd.hize]*qd$wt[2*qd.hize]*(lmt[i,2]-lmt[i,1])/gap[2*qd.hize] 213 if ((lmt[i,1]>qd$pt[1])&(lmt[i,1]<qd$pt[2*qd.hize])) { 214 i.wk <- min((1:(2*qd.hize))[qd$pt>lmt[i,1]]) 215 if (brk[i.wk]<=lmt[i,1]) 216 wk <- d.qd[i.wk]*qd$wt[i.wk]*(lmt[i,2]-lmt[i,1])/gap[i.wk] 217 if (brk[i.wk]>=lmt[i,2]) 218 wk <- d.qd[i.wk-1]*qd$wt[i.wk-1]*(lmt[i,2]-lmt[i,1])/gap[i.wk-1] 219 if ((brk[i.wk]<lmt[i,2])&(brk[i.wk]>lmt[i,1])) 220 wk <- d.qd[i.wk]*qd$wt[i.wk]*(lmt[i,2]-brk[i.wk])/gap[i.wk]+ 221 d.qd[i.wk-1]*qd$wt[i.wk-1]*(brk[i.wk]-lmt[i,1])/gap[i.wk-1] 222 } 223 } 224 zz.wk <- c(zz.wk,wk) 225 } 226 zz <- cbind(zz,as.vector(zz.wk)) 227 } 228 if (ncov==1) zz <- as.vector(zz) 229 } 230 exp(-zz) 231} 232