R/custom.measure.r

Defines functions custom.measure

Documented in custom.measure

custom.measure <- function(X,fun=NULL,cond=NULL,verbose=TRUE,with.se=FALSE,...){
if(is.null(fun)){stop("'fun' need to be specified")}
if(is.function(fun)==FALSE){stop("'fun' needs to be a function")}
if(is.null(X$simulated.data)){stop("Set 'ret=TRUE' in gformula() to work with custom measures.")}
if(verbose==TRUE){if(X$setup$i.type=="custom"){cat("Note: you chose custom interventions. \n 'custom.measure' asssumes that the intervention values at the first time point (time=1) are not identical.\n\n")}}
if(X$setup$measure=="custom"){stop("You can't apply 'custom measure' on an X onto which you already applied 'custom measure'.")}
if(X$setup$B==0 & with.se==T){stop("Standard error can only be calculated if B>0")}
if(X$setup$B>0){if(X$setup$B!=(length(X$simulated.data)-1)){X$setup$B<- (length(X$simulated.data)-1) }}
sdata <- X$simulated.data
odata <- X$observed.data
results.table <- X$results
if(length(grep("mult",X$setup$fams))>0 & X$setup$i.type=="natural"){
if(is.null(sum(sapply(odata[[1]],is.factor))>0)==FALSE){if(sum(sapply(odata[[1]],is.factor))>0){
for(b in 1:(X$setup$B+1)){odata[[b]][,sapply(odata[[b]],is.factor)] <- apply(sapply(( subset(odata[[b]],select=sapply(odata[[b]],is.factor))),as.character),2,as.numeric)}
}}
}
if(X$setup$i.type=="individual"){
  custom <- any(unlist(lapply(lapply(apply(X$results[, grep("^a", names(X$results)), drop=F],1,unique),na.omit),length))>1 )
  if(custom==TRUE){X$setup$i.type<-"custom"}else{X$setup$i.type<-"standard"}
  if(custom==TRUE){X$setup$abar <- do.call(rbind,lapply(lapply(X$setup$abar,colnames),as.numeric))}else{
                   X$setup$abar <- do.call(rbind,lapply(lapply(X$setup$abar,colnames),as.numeric))[,1]
  }
}

if(X$setup$i.type=="standard"){
for(i in 1:length(sdata[[1]])){
if(is.null(cond)==FALSE){sel <- with(as.data.frame(sdata[[1]][[i]]), eval(parse(text=cond),envir=environment()))}else{sel<-rep(TRUE,nrow(sdata[[1]][[i]]))}
results.table[results.table$a1==X$setup$abar[i],"psi"] <- apply(as.matrix(sdata[[1]][[i]][sel,X$setup$Ynodes],ncol=length(X$setup$Ynodes)),2,fun,...)
if(X$setup$B>0){
sbdata <- lapply(sdata, '[[', i); bfs <- matrix(NA,ncol=length(X$setup$Ynodes),nrow=X$setup$B)
for(b in 2:(X$setup$B+1)){bfs[b-1,] <- apply(data.matrix(sbdata[[b]][sel,X$setup$Ynodes]),2,fun,...) } 
results.table[results.table$a1==X$setup$abar[i],c("l95","u95")] <- t(apply(bfs,2,quantile,probs=c(0.025,0.975)))
if(with.se==TRUE){
  if(i==1){results.table$se<-NA}
  results.table[results.table$a1==X$setup$abar[i],c("se")] <- apply(bfs,2,sd)} 
}
}}

if(X$setup$i.type=="custom"){
for(i in 1:length(sdata[[1]])){
if(is.null(cond)==FALSE){sel <- with(as.data.frame(sdata[[1]][[i]]), eval(parse(text=cond),envir=environment()))}else{sel<-rep(TRUE,nrow(sdata[[1]][[i]]))}
results.table[results.table$a1==X$setup$abar[i,1],"psi"] <- apply(as.matrix(sdata[[1]][[i]][sel,X$setup$Ynodes],ncol=length(X$setup$Ynodes)),2,fun,...)
if(X$setup$B>0){
sbdata <- lapply(sdata, '[[', i); bfs <- matrix(NA,ncol=length(X$setup$Ynodes),nrow=X$setup$B)
for(b in 2:(X$setup$B+1)){bfs[b-1,] <- apply(data.matrix(sbdata[[b]][sel,X$setup$Ynodes]),2,fun,...) } 
results.table[results.table$a1==X$setup$abar[i],c("l95","u95")] <- t(apply(bfs,2,quantile,probs=c(0.025,0.975)))
if(with.se==TRUE){
  if(i==1){results.table$se<-NA}
  results.table[results.table$a1==X$setup$abar[i],c("se")] <- apply(bfs,2,sd)}  
}
}}


if(X$setup$i.type=="natural"){
if(is.null(cond)==FALSE){sel_nat <- with(as.data.frame(sdata[[1]][[1]]), eval(parse(text=cond),envir=environment()))}else{sel_nat<-rep(TRUE,nrow(sdata[[1]][[1]]))}
if(is.null(cond)==FALSE){sel_obs <- with(as.data.frame(odata[[1]]), eval(parse(text=cond),envir=environment()))}else{sel_obs<-rep(TRUE,nrow(odata[[1]]))}
results.table[results.table$a1=="natural","psi"] <- apply(as.matrix(sdata[[1]][[1]][sel_nat,X$setup$Ynodes],ncol=length(X$setup$Ynodes)),2,fun,...)
results.table[results.table$a1=="observed","psi"] <- try(apply(data.frame(odata[[1]][sel_obs,X$setup$Ynodes]),2,fun,...))
  if(is.null(X$setup$Lnodes)==FALSE){
  if(verbose==TRUE){cat("Note: specified function and condition will also be applied to counterfactual Lnodes. \n\n")}
      lblocks <- rep(NA,length(X$setup$Ynodes))
    for(i in 1:(length(X$setup$Ynodes))){
    blocks <- make.interval(which(colnames(odata[[1]])%in%X$setup$Ynodes))
    lblocks[i] <- sum(which(colnames(odata[[1]])%in%X$setup$Lnodes)%in%blocks[[i]])
    }
   if(is.wholenumber(length(X$setup$Lnodes)/length(X$setup$Ynodes)) | X$setup$n.t==1){
   results.table[results.table$a1=="natural",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t))]  <- matrix(apply(matrix(sdata[[1]][[1]][sel_nat,X$setup$Lnodes],ncol=length(X$setup$Lnodes)),2,fun,...),ncol=(length(X$setup$Lnodes)/X$setup$n.t),byrow=T)
   results.table[results.table$a1=="observed",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t))] <- matrix(try(apply(data.frame(odata[[1]][sel_obs,X$setup$Lnodes]),2,fun,...)),ncol=(length(X$setup$Lnodes)/X$setup$n.t),byrow=T)
   results.table[results.table$a1=="difference",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t))] <- results.table[results.table$a1=="natural",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t))] - 
                                                                                                            results.table[results.table$a1=="observed",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t))]                                                                                                       
   }else{
    if(length(lblocks)>1 & length(unique(lblocks[-1]))==1){
        results.table[results.table$a1=="natural" & results.table$time!=1,paste0("L_",1:max(lblocks))] <- matrix(apply(matrix(sdata[[1]][[1]][sel_nat,X$setup$Lnodes],ncol=length(X$setup$Lnodes)),2,fun,...),ncol=max(lblocks),byrow=T)
        results.table[results.table$a1=="observed"& results.table$time!=1,paste0("L_",1:max(lblocks))] <- matrix(try(apply(data.frame(odata[[1]][sel_obs,X$setup$Lnodes]),2,fun,...)),ncol=max(lblocks),byrow=T)
        if(verbose==TRUE){cat("Note: it is not entirely clear to which time points your Lnodes belong to. I guessed what could make sense, but please check.\n\n")}                                                                                          
   }else{
      if(max(lblocks,na.rm=T)==1){
        updat.index <- unlist(lapply(apply(t(outer(which(colnames(odata[[1]])%in%X$setup$Lnodes), which(colnames(odata[[1]])%in%X$setup$Ynodes), "<")),1,which),max))
        updt.Lnodes <- apply(as.matrix(sdata[[1]][[1]][sel_nat,X$setup$Lnodes],ncol=length(X$setup$Lnodes)),2,fun,...)[updat.index]
        updt.Lnodes2 <-  try(apply(data.frame(odata[[1]][sel_obs,X$setup$Lnodes]),2,fun,...))[updat.index]
        results.table[results.table$a1=="natural",paste0("L_",1)] <- updt.Lnodes
        results.table[results.table$a1=="observed",paste0("L_",1)] <- updt.Lnodes2 
   }else{if(verbose==TRUE){cat("Note: your number of Lnodes are not a multiple of your number of Ynodes.\n This is fine, but you seem to have more than one Lnode per time point, and maybe unequally distributed. \n It is difficult for me to guess which Lnodes belong 'together', and at which time point. \n Thus, no natural course values for your Lnodes are calculated. \n Use 'ret=TRUE' and calculate manually.\n")}}
   }}
   }
results.table[results.table$a1=="difference","psi"] <- results.table[results.table$a1=="natural","psi"] - results.table[results.table$a1=="observed","psi"]
  
    
if(X$setup$B>0){
sbdata <- lapply(sdata, '[[', 1); bfs <- ofs <- matrix(NA,ncol=length(c(X$setup$Ynodes,X$setup$Lnodes)),nrow=X$setup$B,dimnames=list(NULL,c(X$setup$Ynodes,X$setup$Lnodes)))
for(b in 2:(X$setup$B+1)){
    if(is.null(cond)==FALSE){sel <- with(as.data.frame(sbdata[[b]]), eval(parse(text=cond)))}else{sel<-rep(TRUE,nrow(sbdata[[b]]))}
    bfs[b-1,] <- apply(data.matrix(sbdata[[b]][sel,c(X$setup$Ynodes,X$setup$Lnodes)]),2,fun,...)  }
for(b in 2:(X$setup$B+1)){
    if(is.null(cond)==FALSE){sel <- with(as.data.frame(odata[[b]]), eval(parse(text=cond)))}else{sel<-rep(TRUE,nrow(odata[[b]]))}
    ofs[b-1,] <- apply(data.matrix( odata[[b]][sel,c(X$setup$Ynodes,X$setup$Lnodes)]),2,fun,...) 
}
dfs <- bfs-ofs
results.table[results.table$a1=="natural",c("l95","u95")] <- t(apply(subset(bfs,select=X$setup$Ynodes),2,quantile,probs=c(0.025,0.975)))
results.table[results.table$a1=="observed",c("l95","u95")] <- t(apply(subset(ofs,select=X$setup$Ynodes),2,quantile,probs=c(0.025,0.975)))
results.table[results.table$a1=="difference",c("l95","u95")] <-    t(apply(subset(dfs,select=X$setup$Ynodes),2,quantile,probs=c(0.025,0.975)))
  if(with.se==TRUE){
  results.table$se<-NA
  results.table[results.table$a1=="natural",c("se")] <- apply(subset(bfs,select=X$setup$Ynodes),2,sd)
  results.table[results.table$a1=="observed",c("se")] <- apply(subset(ofs,select=X$setup$Ynodes),2,sd)
  results.table[results.table$a1=="difference",c("se")] <-    apply(subset(dfs,select=X$setup$Ynodes),2,sd)
  }
  if(is.null(X$setup$Lnodes)==FALSE){
    if(is.wholenumber(length(X$setup$Lnodes)/length(X$setup$Ynodes)) | X$setup$n.t==1){
    results.table[results.table$a1=="natural",grep(":l95",colnames(results.table))] <-  matrix(apply(bfs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=length(X$setup$Lnodes)/X$setup$n.t,byrow=T)
    results.table[results.table$a1=="natural",grep(":u95",colnames(results.table))] <-  matrix(apply(bfs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=length(X$setup$Lnodes)/X$setup$n.t,byrow=T)
    results.table[results.table$a1=="observed",grep(":l95",colnames(results.table))] <-  matrix(apply(ofs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=length(X$setup$Lnodes)/X$setup$n.t,byrow=T)
    results.table[results.table$a1=="observed",grep(":u95",colnames(results.table))] <-  matrix(apply(ofs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=length(X$setup$Lnodes)/X$setup$n.t,byrow=T)
    results.table[results.table$a1=="difference",grep(":l95",colnames(results.table))] <-  matrix(apply(dfs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=length(X$setup$Lnodes)/X$setup$n.t,byrow=T)
    results.table[results.table$a1=="difference",grep(":u95",colnames(results.table))] <-  matrix(apply(dfs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=length(X$setup$Lnodes)/X$setup$n.t,byrow=T)
        if(with.se==TRUE){
            results.table[results.table$a1=="natural",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t),":se")]    <-  matrix(apply(bfs[,X$setup$Lnodes],2,sd),ncol=(length(X$setup$Lnodes)/X$setup$n.t),byrow=T)
            results.table[results.table$a1=="observed",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t),":se")]   <-  matrix(apply(ofs[,X$setup$Lnodes],2,sd),ncol=(length(X$setup$Lnodes)/X$setup$n.t),byrow=T)
            results.table[results.table$a1=="difference",paste0("L_",1:(length(X$setup$Lnodes)/X$setup$n.t),":se")] <-  matrix(apply(dfs[,X$setup$Lnodes],2,sd),ncol=(length(X$setup$Lnodes)/X$setup$n.t),byrow=T)}
    }else{
      if(length(lblocks)>1 & length(unique(lblocks[-1]))==1){
        results.table[results.table$a1=="natural" & results.table$time!=1,grep(":l95",colnames(results.table))] <-  matrix(apply(bfs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=max(lblocks),byrow=T)
        results.table[results.table$a1=="natural" & results.table$time!=1,grep(":u95",colnames(results.table))] <-  matrix(apply(bfs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=max(lblocks),byrow=T)
        results.table[results.table$a1=="observed" & results.table$time!=1,grep(":l95",colnames(results.table))] <-  matrix(apply(ofs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=max(lblocks),byrow=T)
        results.table[results.table$a1=="observed" & results.table$time!=1,grep(":u95",colnames(results.table))] <-  matrix(apply(ofs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=max(lblocks),byrow=T)
        results.table[results.table$a1=="difference" & results.table$time!=1,grep(":l95",colnames(results.table))] <-  matrix(apply(dfs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=max(lblocks),byrow=T)
        results.table[results.table$a1=="difference" & results.table$time!=1,grep(":u95",colnames(results.table))] <-  matrix(apply(dfs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=max(lblocks),byrow=T)
          if(with.se==TRUE){
            results.table[results.table$a1=="natural" & results.table$time!=1,paste0("L_",1:max(lblocks),":se")]      <-  matrix(apply(bfs[,X$setup$Lnodes],2,sd),ncol=max(lblocks),byrow=T)
            results.table[results.table$a1=="observed" & results.table$time!=1,paste0("L_",1:max(lblocks),":se")]     <-  matrix(apply(ofs[,X$setup$Lnodes],2,sd),ncol=max(lblocks),byrow=T)
            results.table[results.table$a1=="difference" & results.table$time!=1,paste0("L_",1:max(lblocks),":se")]   <-  matrix(apply(dfs[,X$setup$Lnodes],2,sd),ncol=max(lblocks),byrow=T)}
      }else{
        if(max(lblocks,na.rm=T)==1){
          results.table[results.table$a1=="natural", paste(paste0("L_",1),"l95",sep=":")] <-
          matrix(apply(bfs[,X$setup$Lnodes],2,quantile,probs=c(0.025)),ncol=1,byrow=T)[updat.index]
          results.table[results.table$a1=="natural", paste(paste0("L_",1),"u95",sep=":")] <-
          matrix(apply(bfs[,X$setup$Lnodes],2,quantile,probs=c(0.975)),ncol=1,byrow=T)[updat.index]
          results.table[results.table$a1=="observed", c(paste(paste0("L_",1),"l95",sep=":"))] <-
          apply(ofs[,X$setup$Lnodes],2,quantile,probs=c(0.025))[updat.index]
          results.table[results.table$a1=="observed", c(paste(paste0("L_",1),"u95",sep=":"))] <-
          apply(ofs[,X$setup$Lnodes],2,quantile,probs=c(0.975))[updat.index]
          col.order <- c(1:grep("psi",colnames(results.table)),which(colnames(results.table)%in%c("l95","u95")))
          col.order <- c(col.order,sort(grep(paste0("L_",1),colnames(results.table))))
          if(with.se==TRUE){col.order<-c(col.order, which(colnames(results.table)%in%c("se")))}
          results.table <- results.table[,col.order]
          results.table[results.table$a1=="difference", paste(paste0("L_",1),"l95",sep=":")] <-
          apply(bfs[,X$setup$Lnodes]-ofs[,X$setup$Lnodes],2,quantile,probs=c(0.025))[updat.index]
          results.table[results.table$a1=="difference", paste(paste0("L_",1),"u95",sep=":")] <-
          apply(bfs[,X$setup$Lnodes]-ofs[,X$setup$Lnodes],2,quantile,probs=c(0.975))[updat.index]
          if(with.se==TRUE){
            results.table[results.table$a1=="natural",paste0("L_",1,":se")]      <-  matrix(apply(bfs[,X$setup$Lnodes],2,sd),ncol=1,byrow=T)[updat.index]
            results.table[results.table$a1=="observed",paste0("L_",1,":se")]     <-  matrix(apply(ofs[,X$setup$Lnodes],2,sd),ncol=1,byrow=T)[updat.index]
            results.table[results.table$a1=="difference",paste0("L_",1,":se")]   <-  matrix(apply(dfs[,X$setup$Lnodes],2,sd),ncol=1,byrow=T)[updat.index]}
         }else{
            if(verbose==TRUE){cat("Note: your number of Lnodes are not a multiple of your number of Ynodes.\n This is fine, but you seem to have more than one Lnode per time point, and maybe unequally distributed. \n It is difficult for me to guess which Lnodes belong 'together', and at which time point. \n Thus, no natural course values for your Lnodes are calculated. \n Use 'ret=TRUE' and calculate manually.\n")}
              }}
    }
    }
}                                                                                                            


}

custom.setup <- X$setup
custom.setup$measure<-"custom"
bsa <- ofa <- dfa <- NULL; if(with.se==TRUE){bsa<-bfs
                if(X$setup$i.type=="natural"){ofa<-ofs; dfa<-dfs}}

results <- list(results=results.table,diagnostics=X$diagnostics,setup=custom.setup,b.results = bsa, b.results2 = list(ofa,dfa))
class(results) <- "gformula"
results
}

Try the CICI package in your browser

Any scripts or data that you put into this service are public.

CICI documentation built on June 8, 2025, 11:01 a.m.