R/mstFreq.R

Defines functions mst.perm g.total.mst g.within.mst rbd.rbd rbdP rbd mstFREQ.formula mstFREQ.default mstFREQ

Documented in mstFREQ

#############################################################################
############# MST main functions ################################################

#############################################################################
############# MST main functions ################################################
#' Analysis of Multisite Randomised Education Trials using Multilevel Model under a Frequentist Setting.
#'
#' \code{mstFREQ} performs analysis of multisite randomised education trials using a multilevel model under a frequentist setting.
#'
#' @export
#' @param formula the model to be analysed is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the independent variables.
#' @param random a string variable specifying the "clustering variable" as contained in the data. See example below.
#' @param intervention a string variable specifying the "intervention variable" as appearing in the formula and the data. See example below.
#' @param baseln A string variable allowing the user to specify the reference category for intervention variable. When not specified, the first level will be used as a reference.
#' @param nBoot number of bootstraps required to generate bootstrap confidence intervals.
#' @param nPerm number of permutations required to generate permutated p-value.
#' @param seed seed required for bootstrapping and permutation procedure, if not provided default seed will be used.
#' @param data data frame containing the data to be analysed.
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{Beta}: Estimates and confidence intervals for variables specified in the model.
#' \item \code{ES}: Conditional Hedge's g effect size (ES) and its 95% confidence intervals. If nBoot is not specified, 95% confidence intervals are based on standard errors. If nBoot is specified, they are non-parametric bootstrapped confidence intervals.
#' \item \code{covParm}: A list of variance decomposition into between cluster variance-covariance matrix (schools and school by intervention) and within cluster variance (Pupils). It also contains intra-cluster correlation (ICC).
#' \item \code{SchEffects}: A vector of the estimated deviation of each school from the intercept and intervention slope.
#' \item \code{Perm}: A "nPerm x 2w" matrix containing permutated effect sizes using residual variance and total variance. "w" denotes number of intervention. "w=1" for two arm trial and "w=2" for three arm trial excluding the control group. It is produced only when \code{nPerm} is specified.
#' \item \code{Bootstrap}: A "nBoot x 2w" matrix containing the bootstrapped effect sizes using residual variance (Within) and total variance (Total). "w" denotes number of intervention. "w=1" for two arm trial and "w=2" for three arm trial excluding the control group. It is only prduced when \code{nBoot} is specified.
#' \item \code{Unconditional}: A list of unconditional effect sizes, covParm, Perm and Bootstrap obtained based on variances from the unconditional model (model with only the intercept as a fixed effect).
#' }
#' @example inst/examples/mstExample.R
mstFREQ<- function(formula,random,intervention,baseln,nPerm,data,seed,nBoot)UseMethod("mstFREQ")

#'@export
mstFREQ.default<- function(formula,random,intervention,baseln,nPerm,data,seed,nBoot){stop("No correct formula input given.")}


#' @export
mstFREQ.formula <- function(formula,random,intervention,baseln,nPerm,data,seed,nBoot){

  data <- na.omit(data[ ,unique(c(all.vars(formula),random, intervention))])
  data <- data[order(data.frame(data)[,which(colnames(data)==random)],data[,which(colnames(data)==intervention)]),]
  trt <- data[,which(colnames(data)==intervention)]
  #trt <- as.factor(trt)
  if(missing(baseln)){trt <- as.factor(trt)}
  if(!missing(baseln)){trt <- relevel(as.factor(trt),baseln)}

  tmp2 <- which(colnames(data)==random)
  cluster2 <- data[,tmp2]

  chk <- sum(rowSums(table(cluster2,trt)!=0)>1)
  if(chk ==0){stop("This is not a MST design")}

  if(!missing(nPerm) & !missing(nBoot)){stop("Either nPerm or nBoot must be specified")}
  if(missing(nPerm)){nPerm <-0}
  if(missing(nBoot)){nBoot <-0}


  tmp3 <- which(colnames(data)==intervention)
  if( missing(baseln)){data[,tmp3] <- as.factor(data[,tmp3])}
  if(!missing(baseln)){data[,tmp3] <- relevel(as.factor(data[,tmp3]),baseln)}
  mf <- model.frame(formula=formula, data=data)
  mf <- mf[order(cluster2),]
  cluster <- cluster2[order(cluster2)]
  trt <- trt[order(cluster2)]
  fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
  tmp <- colnames(fixedDesignMatrix )
  tmp[1]  <- "Intercept"
  colnames(fixedDesignMatrix)<- tmp
  posttest <- model.response(mf)
  intervention <- intervention
  btp  <- which(tmp %in% paste0(intervention, unique(data[, intervention]) ))


  if(length(tmp2)!= 1){stop("Cluster variable misspecified")}
  if(length(tmp3)!= 1){stop("Intervention variable misspecified")}


  output <- rbd(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt,cluster=cluster,btp=btp)

  if(nPerm > 0){

    if(nPerm<999){stop("nPerm must be greater than 1000")}
    output$Perm <- mst.perm(formula,data,trt,intervention,nPerm,random,cluster,btp,seed)
    output$Conditional$Perm <- round(data.frame(output$Perm$Conditional),2)
    output$Unconditional$Perm <- round(data.frame(output$Perm$Unconditional),2)
  }

  if(nBoot>0){

    if(nBoot<999){stop("nBoot must be greater than 1000")}

    tid <- c(1:nrow(fixedDesignMatrix))

    #set.seed(1020252)
    if(!missing(seed)){set.seed(seed)}
    bootSamples <- NULL

    for(ii in 1:length(unique(cluster))){
      selID <- tid[cluster==unique(cluster)[ii]]
      if(length(selID)>0){
        selID2<- sapply(c(1:nBoot),function(x)selID [sample(1:length(selID), length(selID),replace=TRUE)])

        bootSamples <- rbind(bootSamples ,selID2)
      }

    }

    bootResults <- apply(bootSamples ,2,function(bt)rbd.rbd(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt,cluster=cluster,bt=bt, btp=btp))

    bootES <- bootCompile(output=output,trt=trt,bootResults=bootResults,intervention=intervention)
    output$ES <- bootES
    output$Bootstrap <- bootResults
    output$Conditional$Bootstrap <- data.frame(t(sapply(1:nBoot, function(i)unlist(bootResults[[i]]$Conditional))))
    output$Unconditional$Bootstrap <- data.frame(t(sapply(1:nBoot, function(i)unlist(bootResults[[i]]$Unconditional))))
    Bnames <- gsub("Estimate1","Within", gsub("Estimate2", "Total",names(output$Conditional$Bootstrap)))
    names(output$Conditional$Bootstrap) <- Bnames
    names(output$Unconditional$Bootstrap) <- Bnames
  }
  output1 <- list()
  output1$Beta   <- output$Beta
  output1$covParm <- output$covParm$Conditional
  output1$ES     <- output$ES$Conditional
  output1$SchEffects     <- output$SchEffects
  if(nPerm > 0){output1$permES    <- output$Conditional$Perm}
  if(nBoot > 0){output1$Bootstrap <- output$Conditional$Bootstrap}

  output1$Unconditional$ES     <- output$ES$Unconditional
  output1$Unconditional$covParm <- output$covParm$Unconditional
  if(nPerm > 0){output1$Unconditional$permES <- output$Unconditional$Perm}
  if(nBoot > 0){output1$Unconditional$Bootstrap <-  output$Unconditional$Bootstrap}


  output1$Method <- "MLM"
  output1$Function <- "mstFREQ"
  class(output1) <- "eefAnalytics"
  return(output1)
}



###########################################################################################

## - internal

rbd <- function(posttest,fixedDesignMatrix,intervention,trt,cluster,btp){

  freqFit <- lmer(posttest~ fixedDesignMatrix-1+(1+trt|cluster))
  np<- row.names(summary(freqFit)$coef)
  cit <- confint(freqFit,np)
  betaB <- data.frame(cbind(summary(freqFit)$coefficients[,1],cit))
  row.names(betaB)<- colnames(fixedDesignMatrix)
  colnames(betaB) <- c("Estimate","95% LB ","95% UB")
  betaB <- betaB
  var.B2<- as.matrix(VarCorr(freqFit)$cluster)
  var.B2_1<-as.data.frame(VarCorr(freqFit)$cluster)
  var.B2_2<-as.data.frame(VarCorr(freqFit))
  var.B3 <- diag(var.B2)[-1]
  var.sch <- var.B2[1,1]
  vcov.schTrt <- t(as.matrix(na.omit(var.B2_2$vcov[var.B2_2$var1=="(Intercept)"][-1])))
  var.E <- var.B3
  var.W<- summary(freqFit)$sigma^2
  N <- as.numeric(length(summary(freqFit)$res))
  N.t <- as.matrix(summary(trt))[-1];
  var.tt <- var.sch+var.W+sum(N.t/N*(var.B3+2*vcov.schTrt))
  ICC1 <- sum(var.B2)/var.tt
  sigmaBE1 <- list(round(var.B2_1,2),round(var.W,2),round(var.tt,2),round(ICC1,2))
  names(sigmaBE1)<-c('School:trt', 'Pupils', 'Total', 'ICC')
  freqFit1<- lmer(posttest~1+(1|cluster))
  var.B21<- as.matrix(VarCorr(freqFit1)$cluster)
  var.B21_1<-as.data.frame(VarCorr(freqFit1)$cluster)
  var.sch1 <- var.B21[1,1]
  var.W1<- summary(freqFit1)$sigma^2
  var.tt1 <- var.sch1+var.W1
  ICC2 <- sum(var.B21)/var.tt1
  sigmaBE2 <- list(round(var.B21_1,2),round(var.W1,2),round(var.tt1,2),round(ICC2,2))
  names(sigmaBE2)<-c('School', 'Pupils', 'Total', 'ICC')
  ICC <-(c(conditional=ICC1,unconditional=ICC2))
  sigmaBE <- list(sigmaBE1,sigmaBE2)
  names(sigmaBE)<-c('Conditional','Unconditional')
  sigmaBE <- sigmaBE
  schRand <- data.frame(unique(cluster),ranef(freqFit)$cluster)
  names(schRand)[1]<- "Schools"
  names(schRand)[2]<- "Intercept"
  #names(schRand)<- c("Schools","Intercept", paste0('trt', 1:(ncol(ranef(freqFit)$cluster)-1)))
  #btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)

  output2 <- list()
  for( i in 1:length(btp)){

    beta <- betaB[btp[i],1]
    group <- fixedDesignMatrix[,btp[i]]
    var.e<- var.E[i]
    esWithin <- g.within.mst(var.w=var.W, var.e=var.e, beta=beta, group=group, schoolID=cluster)
    esTotal <- g.total.mst(var.w=var.W, var.e=var.e, var.tt=var.tt, beta=beta, group=group, schoolID=cluster)
    outputc <- data.frame(rbind(esWithin,esTotal))
    colnames(outputc) <- c("Estimate","95% LB","95% UB")
    rownames(outputc) <- c("Within","Total")
    output2[[i]] <- round(outputc,2)
  }
  names(output2) <- row.names(betaB)[btp]
  output3<-list()
  for( i in 1:length(btp)){

    beta <- betaB[btp[i],1]
    group <- fixedDesignMatrix[,btp[i]]
    esWithin1 <- g.within(var.w=var.W1, beta=beta, icc=ICC2, group=group, schoolID=cluster)
    esTotal1 <- g.total(var.tt=var.tt1, beta=beta, icc=ICC2, group=group, schoolID=cluster)
    outputu <- data.frame(rbind(esWithin1,esTotal1))
    colnames(outputu) <- c("Estimate","95% LB","95% UB")
    rownames(outputu) <- c("Within","Total")
    output3[[i]] <- round(outputu,2)
  }
  names(output3) <- row.names(betaB)[btp]
  output2.3<-list(Conditional=output2,Unconditional=output3)
  output <- list(Beta=round(betaB,2),covParm=sigmaBE,ES=output2.3,SchEffects=round(schRand,2))

  return(output)
}


## - internal

rbdP <- function(posttest,fixedDesignMatrix,intervention,trt,cluster,btp){

  freqFit <- lmer(posttest~ fixedDesignMatrix-1+(1+trt|cluster))
  np<- row.names(summary(freqFit)$coef)
  cit <- confint(freqFit,np)
  betaB <- data.frame(cbind(summary(freqFit)$coefficients[,1],cit))
  row.names(betaB)<- colnames(fixedDesignMatrix)
  #colnames(betaB) <- c("Estimate","95% LB ","95% UB")
  betaB <- betaB
  var.B2<- as.matrix(VarCorr(freqFit)$cluster)
  var.B2_1<-as.data.frame(VarCorr(freqFit)$cluster)
  var.B2_2<-as.data.frame(VarCorr(freqFit))
  var.B3 <- diag(var.B2)[-1]
  var.sch <- var.B2[1,1]
  vcov.schTrt <- t(as.matrix(na.omit(var.B2_2$vcov[var.B2_2$var1=="(Intercept)"][-1])))
  var.E <- var.B3
  var.W<- summary(freqFit)$sigma^2
  N <- as.numeric(length(summary(freqFit)$res))
  N.t <- as.matrix(summary(trt))[-1];
  var.tt <- var.sch+var.W+sum(N.t/N*(var.B3+2*vcov.schTrt))
  ICC1 <- sum(var.B2)/var.tt
  sigmaBE1 <- list(round(var.B2_1,2),round(var.W,2),round(var.tt,2),round(ICC1,2))
  names(sigmaBE1)<-c('School:trt', 'Pupils', 'Total', 'ICC')
  freqFit1<- lmer(posttest~1+(1|cluster))
  var.B21<- as.matrix(VarCorr(freqFit1)$cluster)
  var.B21_1<-as.data.frame(VarCorr(freqFit1)$cluster)
  var.sch1 <- var.B21[1,1]
  var.W1<- summary(freqFit1)$sigma^2
  var.tt1 <- var.sch1+var.W1
  ICC2 <- sum(var.B21)/var.tt1
  sigmaBE2 <- list(round(var.B21_1,2),round(var.W1,2),round(var.tt1,2),round(ICC2,2))
  names(sigmaBE2)<-c('School', 'Pupils', 'Total', 'ICC')
  ICC <-(c(conditional=ICC1,unconditional=ICC2))
  sigmaBE <- list(sigmaBE1,sigmaBE2)
  names(sigmaBE)<-c('Conditional','Unconditional')
  sigmaBE <- sigmaBE
  schRand <- data.frame(unique(cluster),ranef(freqFit)$cluster)
  names(schRand)<- c("Schools","Intercept", paste0('trt', 1:(ncol(ranef(freqFit)$cluster)-1)))
  #btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)

  output2 <- list()
  for( i in 1:length(btp)){

    beta <- betaB[btp[i],1]
    group <- fixedDesignMatrix[,btp[i]]
    var.e<- var.E[i]
    esWithin <- g.within.mst(var.w=var.W, var.e=var.e, beta=beta, group=group, schoolID=cluster)
    esTotal <- g.total.mst(var.w=var.W, var.e=var.e, var.tt=var.tt, beta=beta, group=group, schoolID=cluster)
    outputc <- round(data.frame(rbind(esWithin,esTotal)),2)
    colnames(outputc) <- c("Estimate","95% LB","95% UB")
    rownames(outputc) <- c("Within","Total")
    output2[[i]] <- round(outputc,2)
  }
  names(output2) <- row.names(betaB)[btp]
  output3<-list()
  for( i in 1:length(btp)){

    beta <- betaB[btp[i],1]
    group <- fixedDesignMatrix[,btp[i]]
    esWithin1 <- g.within(var.w=var.W1, beta=beta, icc=ICC2, group=group, schoolID=cluster)
    esTotal1 <- g.total(var.tt=var.tt1, beta=beta, icc=ICC2, group=group, schoolID=cluster)
    outputu <- data.frame(rbind(esWithin1,esTotal1))
    colnames(outputu) <- c("Estimate","95% LB","95% UB")
    rownames(outputu) <- c("Within","Total")
    output3[[i]] <- round(outputu,2)
  }
  names(output3) <- row.names(betaB)[btp]
  output2.3<-list(Conditional=output2,Unconditional=output3)
  output <- list(ES=output2.3)

  return(output)
}

## - internal

rbd.rbd <- function(posttest,fixedDesignMatrix,intervention,trt,cluster,bt,btp){

  posttest2 <- posttest[bt]
  fixedDesignMatrix2 <- fixedDesignMatrix[bt,]
  trt2 <- trt[bt]
  cluster2 <- cluster[bt]


  freqFit <- try(lmer(posttest2~ fixedDesignMatrix2-1+(1+trt2|cluster2)),silent=TRUE)
  output2 <- NULL
  if(!is(freqFit, "try-error")){

    betaB <- data.frame(summary(freqFit)$coefficients[,1])
    row.names(betaB)<- colnames(fixedDesignMatrix)
    betaB <- betaB
    var.B2<- as.matrix(VarCorr(freqFit)$cluster)
    var.B2_1<-as.data.frame(VarCorr(freqFit)$cluster)
    var.B2_2<-as.data.frame(VarCorr(freqFit))
    var.B3 <- diag(var.B2)[-1]
    var.sch <- var.B2[1,1]
    vcov.schTrt <- t(as.matrix(na.omit(var.B2_2$vcov[var.B2_2$var1=="(Intercept)"][-1])))
    var.E <- var.B3
    var.W<- summary(freqFit)$sigma^2
    N <- as.numeric(length(summary(freqFit)$res))
    N.t <- as.matrix(summary(trt))[-1];
    var.tt <- var.sch+var.W+sum(N.t/N*(var.B3+2*vcov.schTrt))
    ICC1 <- sum(var.B2)/var.tt
    sigmaBE1 <- list(round(var.B2_1,2),round(var.W,2),round(var.tt,2),round(ICC1,2))
    names(sigmaBE1)<-c('School:trt', 'Pupils', 'Total', 'ICC')
    freqFit1<- lmer(posttest~1+(1|cluster2))
    var.B21<- as.matrix(VarCorr(freqFit1)$cluster)
    var.B21_1<-as.data.frame(VarCorr(freqFit1)$cluster)
    var.sch1 <- var.B21[1,1]
    var.W1<- summary(freqFit1)$sigma^2
    var.tt1 <- var.sch1+var.W1
    ICC2 <- sum(var.B21)/var.tt1
    sigmaBE2 <- list(round(var.B21_1,2),round(var.W1,2),round(var.tt1,2),round(ICC2,2))
    names(sigmaBE2)<-c('School', 'Pupils', 'Total', 'ICC')
    ICC <-(c(conditional=ICC1,unconditional=ICC2))
    sigmaBE <- list(sigmaBE1,sigmaBE2)
    names(sigmaBE)<-c('Conditional','Unconditional')
    sigmaBE <- sigmaBE
    schRand <- data.frame(unique(cluster),ranef(freqFit)$cluster)
    names(schRand)<- c("Schools","Intercept", paste0('trt', 1:(ncol(ranef(freqFit)$cluster)-1)))
    #btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)


    output3 <- list()
    for( i in 1:length(btp)){

      beta <- betaB[btp[i],1]
      group <- fixedDesignMatrix[,btp[i]]
      var.e<- var.E[i]
      esWithin <- beta/sqrt(var.W)
      esTotal <- beta/sqrt(var.tt)

      outputc <- data.frame(rbind(esWithin,esTotal))
      names(outputc) <- c("Estimate")
      rownames(outputc) <- c("Within","Total")
      output3[[i]] <- round(outputc,2)
    }
    names(output3) <- row.names(betaB)[btp]
    output4<-list()
    for( i in 1:length(btp)){

      beta <- betaB[btp[i],1]
      group <- fixedDesignMatrix[,btp[i]]
      esWithin1 <- beta/sqrt(var.W1)
      esTotal1 <- beta/sqrt(var.tt1)

      outputu <- data.frame(rbind(esWithin1,esTotal1))
      names(outputu) <- c("Estimate")
      rownames(outputu) <- c("Within","Total")
      output4[[i]] <- round(outputu,2)
    }
    names(output4) <- row.names(betaB)[btp]
    output2<-list(Conditional=output3,Unconditional=output4)

  }
  return(output2)
}



## - internal

g.within.mst <- function(var.w, var.e, beta, group, schoolID){
  t <- group; id <- schoolID
  d.w <- (beta/sqrt(var.w))
  n.itc <- as.data.frame.matrix(table(id,t))
  n.it <- n.itc[n.itc$`0`!=0 & n.itc$`1`!=0, "1"]
  n.ic <- n.itc[n.itc$`0`!=0 & n.itc$`1`!=0, "0"]
  M <- length(unique(id))
  N.t <- sum(table(id[t==1])); N.c <- sum(table(id[t==0]))
  N <- (N.t + N.c)
  nin<-(n.it*n.ic)/(n.it+n.ic)
  ni <- n.it+n.ic
  vterm1 <- 1/(sum(var.w/(var.e+var.w/nin)))
  vterm2 <- ((d.w^2)/((2*N-4*M)))
  se <- sqrt(vterm1+vterm2)
  LB <- (d.w-1.96*se); UB <- (d.w+1.96*se)
  output <- data.frame(d.w, LB, UB)
  names(output) <- c("g", "LB", "UB")
  return(output)
}

## - internal

g.total.mst <- function(var.w, var.e, var.tt, beta, group, schoolID){
  t <- group; id <- schoolID
  n.itc <- as.data.frame.matrix(table(id,t))
  n.it <- n.itc[n.itc$`0`!=0 & n.itc$`1`!=0, "1"]
  n.ic <- n.itc[n.itc$`0`!=0 & n.itc$`1`!=0, "0"]
  M <- length(unique(id))
  N.t <- sum(table(id[t==1])); N.c <- sum(table(id[t==0]))
  N <- (N.t + N.c)
  nin<-(n.it*n.ic)/(n.it+n.ic)
  ni <- n.it+n.ic
  d.t <- (beta/sqrt(var.tt))
  vterm1 <- 1/(sum(var.tt/(var.e+var.w/nin)))
  vterm2 <- ((d.t^2)/((2*N-4*M)))
  se <- sqrt(vterm1+vterm2)
  LB <- (d.t-1.96*se); UB <- (d.t+1.96*se)
  output <- data.frame(d.t, LB, UB)
  names(output)<- c("g", "LB", "UB")
  return(output)
}


## - internal

mst.perm <- function(formula,data,trt,intervention,nPerm,random,cluster,btp,seed){

  data2 <- data
  g <- matrix(NA,nPerm,2*(length(unique(trt))-1))
  g.unc <- matrix(NA,nPerm,2*(length(unique(trt))-1))

  for(i in 1:nPerm){
    #set.seed(12890*i+1)
    if(!missing(seed)){set.seed(seed*i+1)}
    tryCatch({data2[,which(colnames(data)==intervention)]<-unlist(tapply(trt,cluster,function(x)sample(x)))
    data3 <- data2[order(data2[,which(colnames(data2)==random)],data2[,which(colnames(data2)==intervention)]),]
    cluster = data3[,which(colnames(data3)==random)]
    mf <- model.frame(formula=formula, data=data3)
    fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data3)))
    tmp <- colnames(fixedDesignMatrix )
    tmp[1]  <- "Intercept"
    colnames(fixedDesignMatrix)<- tmp
    posttest <- model.response(mf)
    intervention <- intervention
    trt2 <- data3[,which(colnames(data3)==intervention)]
    p2CRTFREQ <-rbdP(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt2,cluster=cluster, btp=btp)

    chkppp <- data.frame(cond=unlist(p2CRTFREQ$ES$Conditional),uncond=unlist(p2CRTFREQ$ES$Unconditional))
    chkppp2 <- c(seq(1,6*(length(unique(trt))-1),6),seq(2,6*(length(unique(trt))-1),6))
    chkppp3 <- chkppp2[order(chkppp2)]
    g[i,]  <-  chkppp[chkppp3,"cond"]
    g.unc[i,]  <-  chkppp[chkppp3,"uncond"]}, error=function(e){})
  }
  ntpp <- rep(names(p2CRTFREQ$ES$Conditional),2)
  ntpp <- ntpp[order(ntpp)]
  wt <- rep(c("Within","Total"),length(names(p2CRTFREQ$ES$Conditional)))
  colnames(g) <- paste(ntpp ,wt,sep="")
  colnames(g.unc) <- paste(ntpp ,wt,sep="")
  g1 <- list(Conditional=g,Unconditional=g.unc)
  return(g1)
}

Try the eefAnalytics package in your browser

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

eefAnalytics documentation built on March 16, 2021, 5:08 p.m.