R/within_between_func.R

In duncan-clark/Blolog: Bayesian Inference for LOLOG Network Models

#This function calculate the within between perm variation for Bayesian Lolog fit:

#' 
#' @export


within_between_func <- function(sample#samples from each chain - list same length as number of perms.
){
  
  #get number of perms and parameters:
  perms <- length(sample)
  lengths <- sapply(sample,length)
  params <- length(sample[[1]][[1]])
  
  means_within <- lapply(sample,function(x){
    tmp <- do.call(rbind,x)
    tmp <- apply(tmp,2,mean)
    return(tmp)
  })
  
  sd_within <- lapply(sample,function(x){
    tmp <- do.call(rbind,x)
    tmp <- apply(tmp,2,sd)
    return(tmp)
  })
  
  ss_between <- (perms-1)*apply(do.call(rbind,means_within),2,function(x){sd(sqrt(lengths)*x)**2})
  ss_within <- Reduce("+",mapply(sd_within,lengths,FUN = function(x,length){(x**2)*(length-1)},SIMPLIFY  = F))
  ss_tot <- (sum(lengths)-1)*apply(do.call(rbind,do.call(c,sample)),2,sd)**2
  
  
  if(sum(round(ss_between + ss_within,4)!=round(ss_tot,4))!=0){stop("The between and within SS do not sum to the total SS")}
  
  return(list(ss_tot = ss_tot,
              ss_between = ss_between,
              ss_within = ss_within))
  
}