## a couple of modifications denoted by mwa-cf, made by mwa-cf (23 january 2007) ## Splus5 functions. Similar to Splus, with certain exceptions, ## such as integer vs dble format of numbers. Also, for linking ## Fortran 77 subroutines, Splus5 uses ``dyn.open", as below, ## whereas Splus uses ``dyn.load". ## The authors are not responsible for any consequences ## that may result in use of this file by any person. index.survival <- function(survival, censoring) { index.temp <- c(1:length(survival)) for(i in 1:length(survival)) { if(survival[i] < censoring[i]) { index.temp[i] <- 1 } else { index.temp[i] <- 0 } } index.temp } is.in <- function(x, r) { dummy <- 1 - 1 * ((x < r[1]) || (x > r[2])) dummy } index.time <- function(t0, jump) { i <- 1 l <- length(jump) if(l > 1) { if((jump[i]) > t0) { temp <- 0 } else { while(((jump[(i + 1)]) <= t0) && ((i + 2) <= l)) { i <- i + 1 } if((jump[(l - 1)]) < t0) { temp <- l } else { temp <- i } } } else { temp <- 1 - 1 * (t0 < jump[1]) } temp } cens.surv <- function(x, d) { n.total <- length(x) step <- cbind(sort(x), as.double(c(n.total:1))) fail.times <- sort(x[d == 1]) list(step = step, fail.times = fail.times) } dblsample.surv <- function(x.ct,Delta,Robs,S) { # # INPUTS. # # (all variables have equal length=the number of people # in the cohort.) # # 1. Robs is indicator: 1 if subject is not a study dropout. # 2. S is indicator: 1 if subject is a study dropout and is # double sampled. # 3,4 x.ct and Delta: (known only for Robs=1 or S=1) # x.ct=min(survival T, administrative censoring time C) # Delta= 1 if Tt) gets distinct values. # 4. probs: the estimator for pr(T>t) in Frangakis and Rubin (2001), # computed at t in ``failtimes". # 5. se : the pointwise standard error of the estimator, # computed, at t in failtimes, using the delta method on the standard # error on the corresponding cumulative hazard. # 6. lower: the lower point of the approximate 95% pointwise # confidence interval (CI) for Pr(T>t), evaluated at t=failtimes. # It is evaluated as the reciprocal of the exponentiated # upper limit for the CI of the hazard function. # 7. upper: the upper point of the CI in 6. #preliminary calcs. # dyn.open("routds.o") ##mwa-cf: commented data <- data.frame(cbind(x.ct, Delta, Robs, S)) ##mwa-cf: added dyn.load("routds.dll") ##mwa-cf: added n.total <- length(x.ct) n.nlfu <- as.integer(sum(Robs)) n.lfu <- as.integer(n.total - n.nlfu) p.lfu <- n.lfu/n.total n.sampled <- as.integer(sum(S)) data.0 <- data[data$Robs == 1, ] data.1 <- data[data$S == 1, ] x.ct.0<-x.ct[Robs==1] ##mwa-cf: changed from xx.ct to x.ct x.ct.1<-x.ct[S==1] ##mwa-cf: changed from xx.ct to x.ct Delta.0<-Delta[Robs==1] Delta.1<-Delta[S==1] sufstat.0 <- cens.surv(x.ct.0, Delta.0) sufstat.1 <- cens.surv(x.ct.1, Delta.1) n.failed.0 <- length(sufstat.0$fail.times) n.failed.1 <- length(sufstat.1$fail.times) n.failed <- as.integer(n.failed.0 + n.failed.1) #mat will have times,dsample indic. and survprobs if(n.failed > 0) { mat <- cbind(c(sufstat.0$fail.times, sufstat.1$fail.times), c( rep(0, n.failed.0), rep(1, n.failed.1))) mat <- mat[order(mat[, 1]), ] # Fortran function should start here estimator <- .Fortran("estimator", sufstat.0$step, n.nlfu, sufstat.1$step, n.sampled, p.lfu, n.total, mat[, 1:2], n.failed, numeric(n.nlfu), numeric(n.sampled), probs = numeric(n.failed), se = numeric(n.failed), upper = numeric(n.failed), lower = numeric(n.failed)) failtimes.d <- mat[, 1] lfustatus.d <- mat[, 2] probs.d <- estimator$probs se.d <- estimator$se lower.d <- estimator$lower upper.d <- estimator$upper } else { failtimes.d <- 1 lfustatus.d <- 1 probs.d <- 1 se.d <- 0 lower.d <- 1 upper.d <- 1 } list(lfustatus = lfustatus.d, size = n.failed, failtimes = failtimes.d, probs = probs.d, se = se.d, lower = lower.d, upper = upper.d) }