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