R/runModel.R

Defines functions runModel

Documented in runModel

runModel <- function(outdata, outpriors, link="logit", quantiles = c(0.025, 0.5, 0.975), verbose=FALSE, num.threads = 1){
  if(requireNamespace("INLA", quietly = TRUE)){
    if (!(sum(search()=="package:INLA"))==1){
      stop("INLA need to be loaded! \n
Please use the following command to load INLA,\n
library(INLA) \n")
    }
    
    model.type = outdata$model.type
    N = dim(outdata$internaldata)[1]
    varnames = names(outdata$internaldata)
    nnv = c("studynames", "Y", "id", "Ntrials") 
    nnvn = unlist(lapply(nnv, function(x) which(varnames==x)))
    varnames = varnames[-nnvn]
    
    
    if(outpriors$wishart.flag){
      data_temp = outdata$internaldata
      data <- rbind(data_temp[seq(1,N,by=2),], data_temp[seq(2,N,by=2),])
      data$id <- 1:N
    }else{
      data = outdata$internaldata
    }
    #data = outdata$internaldata
    
    if(model.type==1){names.fitted = c("Se","Sp")}
    if(model.type==2){names.fitted = c("Se","1-Sp")}
    if(model.type==3){names.fitted = c("1-Se","Sp")}
    if(model.type==4){names.fitted = c("1-Se","1-Sp")}
    
    
    if(!outdata$modality.flag){ # no modality
      if(!outdata$covariates.flag){ # no covariates
        
        names.summarized.fitted = paste("mean(",link,".",names.fitted,")",sep="")
        
        lc1.ind  = agrep("mu", varnames, max.distance=0)
        lc2.ind  = agrep("nu", varnames, max.distance=0)
        
        lc1text = paste("INLA::inla.make.lincomb(",paste(varnames[lc1.ind], "=1",sep="", collapse = ", "),")",sep="")
        lc2text = paste("INLA::inla.make.lincomb(",paste(varnames[lc2.ind], "=1",sep="", collapse = ", "),")",sep="")
        
        lc1 = eval(parse(text=lc1text))
        names(lc1) = names.summarized.fitted[1]
        lc2 = eval(parse(text=lc2text))
        names(lc2) = names.summarized.fitted[2]
        lc = c(lc1, lc2) 
      }else{ # has covariates
        studynames = outdata$originaldata$studynames
        
        mu.ind = which(data[,"mu"]==1)
        nu.ind = which(data[,"nu"]==1)
        
        lc1.ind  = c(agrep("mu", varnames, max.distance=0), agrep("alpha", varnames, max.distance=0))
        lc2.ind  = c(agrep("nu", varnames, max.distance=0), agrep("beta", varnames, max.distance=0))
        lc1text = paste("INLA::inla.make.lincombs(",paste(varnames[lc1.ind], "=", data[mu.ind,varnames[lc1.ind]],sep="", collapse = ", "),")",sep="")
        lc2text = paste("INLA::inla.make.lincombs(",paste(varnames[lc2.ind], "=", data[nu.ind,varnames[lc2.ind]],sep="", collapse = ", "),")",sep="")
        lc1 = eval(parse(text=lc1text))
        names(lc1) = paste("mean(",link,".",names.fitted[1],".",1:(0.5*N),")",sep="")
        lc2 = eval(parse(text=lc2text))
        names(lc2) = paste("mean(",link,".",names.fitted[2],".",1:(0.5*N),")",sep="")
        lc = c(lc1, lc2)   
      }
    }else{ # has modality
      um = as.character(unique(outdata$originaldata[,outdata$modality.setting]))
      if(!outdata$covariates.flag){ # no covariates
        lc1 = c()
        lc2 = c()
        for(i in 1:length(um)){
          umname = um[i]
          umname = paste(unlist(strsplit(basename(umname), "[-]")), collapse=".")
          names.summarized.fitted = paste("mean(",link,".",names.fitted,".",umname,")",sep="")
          lc1.ind = agrep(paste("mu.",umname,sep=""), varnames, max.distance=0)
          lc2.ind = agrep(paste("nu.",umname,sep=""), varnames, max.distance=0)
          lc1text_um = paste("INLA::inla.make.lincomb(",paste(varnames[lc1.ind], "=1", sep="", collapse = ", "),")",sep="")
          lc2text_um = paste("INLA::inla.make.lincomb(",paste(varnames[lc2.ind], "=1", sep="", collapse = ", "),")",sep="")
          lc1_um = eval(parse(text=lc1text_um))
          names(lc1_um) = names.summarized.fitted[1]
          lc1 = append(lc1, lc1_um)
          lc2_um = eval(parse(text=lc2text_um))
          names(lc2_um) = names.summarized.fitted[2]
          lc2 = append(lc2, lc2_um)
        }
        lc = c(lc1, lc2) 
      }else{ # has covariates
        lc1 = c()
        lc2 = c()
        for(i in 1:length(um)){
          umname = um[i]
          umname = paste(unlist(strsplit(basename(umname), "[-]")), collapse=".")
          mu.ind = which(data[,paste("mu.",umname,sep="")]==1)
          nu.ind = which(data[,paste("nu.",umname,sep="")]==1)
          
          lc1.ind  = c(agrep(paste("mu.",umname,sep=""), varnames, max.distance=0), agrep("alpha", varnames, max.distance=0))
          lc2.ind  = c(agrep(paste("nu.",umname,sep=""), varnames, max.distance=0), agrep("beta", varnames, max.distance=0))
          lc1text = paste("INLA::inla.make.lincombs(",paste(varnames[lc1.ind], "=", data[mu.ind,varnames[lc1.ind]],sep="", collapse = ", "),")",sep="")
          lc2text = paste("INLA::inla.make.lincombs(",paste(varnames[lc2.ind], "=", data[nu.ind,varnames[lc2.ind]],sep="", collapse = ", "),")",sep="")
          lc1_um = eval(parse(text=lc1text))
          names(lc1_um) = paste("mean(",link,".",names.fitted[1],".",umname,".",1:length(mu.ind),")",sep="")
          lc1 = append(lc1, lc1_um)
          lc2_um = eval(parse(text=lc2text))
          names(lc2_um) = paste("mean(",link,".",names.fitted[2],".",umname,".",1:length(nu.ind),")",sep="")
          lc2 = append(lc2, lc2_um)
        }
        lc = c(lc1, lc2) 
      }
    }
    
    quantiles = c(quantiles, 0.025, 0.5, 0.975)
    quantiles = sort(unique(quantiles))
    
    Ntrials = data$Ntrials
    
    if(!outpriors$wishart.flag){
      prec1 = outpriors$prec1
      prec2 = outpriors$prec2
      cor = outpriors$cor
      
      fm <- as.formula(paste("Y ~ f(id, model=\"2diid\", hyper=list(prec1 = prec1,prec2 = prec2,cor = cor), n=N) + ", 
                             paste(varnames, collapse= " + "), "-1"))
      
      model = INLA::inla(fm, family="binomial", data=data, Ntrials=Ntrials, quantiles=quantiles,
                         verbose=verbose, lincomb = lc, num.threads = num.threads,
                         control.inla=list(lincomb.derived.correlation.matrix = TRUE),
                         control.predictor=list(compute=TRUE),
                         control.family = list(link = link),
                         control.compute = list(dic = TRUE, waic = TRUE, cpo=TRUE, mlik = TRUE,config=TRUE, return.marginals.predictor=TRUE))
      
      if(!model$ok){
        stop("Something wrong while running model with data! Please set verbose=TRUE to check!!!!")
      }
    }else{
      prec1 = outpriors$prec1
      
      fm <- as.formula(paste("Y ~ f(id, model=\"iid2d\", hyper=list(prec1 = prec1), n=N) + ", 
                             paste(varnames, collapse= " + "), "-1"))
      
      
      
      model = INLA::inla(fm, family="binomial", data=data, Ntrials=Ntrials, quantiles=quantiles,
                         verbose=verbose, lincomb = lc,  num.threads = num.threads,
                         control.inla=list(lincomb.derived.correlation.matrix = TRUE),
                         control.predictor=list(compute=TRUE),
                         control.family = list(link = link),
                         control.compute = list(dic = TRUE, waic = TRUE, cpo=TRUE, mlik = TRUE, config=TRUE, return.marginals.predictor=TRUE))
      
      if(!model$ok){
        stop("Something wrong while running model with data! Please set verbose=TRUE to check!!!!")
      } 
    }
    
    model$link = link
    model$quantiles = quantiles
    model$verbose = verbose
    model$inherent.outdata = outdata
    model$inherent.priors = outpriors
    
    if(link=="logit"){
      model$linkfunc = function(x){return(log(x/(1-x)))}
      model$inv.linkfunc = function(x){return(exp(x)/(1+exp(x)))}
      model$c.inv.linkfunc = function(x){return(1-exp(x)/(1+exp(x)))}
    }
    if(link=="probit"){
      model$linkfunc = function(x){return(qnorm(x))}
      model$inv.linkfunc = function(x){return(pnorm(x))}
      model$c.inv.linkfunc = function(x){return(1-pnorm(x))}
    }
    if(link=="cloglog"){
      model$linkfunc = function(x){return(log(-log(1-x)))}
      model$inv.linkfunc = function(x){return(1-exp(-exp(x)))}
      model$c.inv.linkfunc = function(x){return(exp(-exp(x)))}
    }
    
    return(model)
  }else{
    stop("INLA need to be installed and loaded!\n
Please use the following command to install and load INLA,\n
install.packages(\"INLA\", repos=c(getOption(\"repos\"), INLA=\"https://inla.r-inla-download.org/R/testing\"), dep=TRUE) \n
library(INLA) \n")
  }
}

Try the meta4diag package in your browser

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

meta4diag documentation built on Nov. 29, 2023, 3:01 a.m.