R/sprfunctions.R

In chaddlewick/spr: Simulation-based Parameter Recovery

#### SPR Functions ####

designdf = function(variablespecification){

if(is.null(observedvariables$participant) == TRUE){
  print("error: participant details missing")
} else {
  pptdesign = observedvariables$participant
  dataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
  dataset$participant = eval(parse(text = pptdesign))
}

observedvariablesnop = subset(observedvariables, names(observedvariables) != "participant")
variablenumber = 1


for(i in observedvariablesnop){

  tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
  variablename = names(observedvariablesnop)[variablenumber]
  tempnames = c("observation", variablename)

  # For repeat vs observed structures:
  if(unlist(strsplit(as.character(i), split = "\\("))[1] == "rep"){
    # For repeat-strucure variables:
    variabledesign = observedvariablesnop[variablenumber]
    tempdataset$tempname = eval(parse(text = variabledesign))
    names(tempdataset) = tempnames
    dataset = merge.data.frame(dataset, tempdataset, by ="observation")


  } else {
    # For observed-structure variables:
    variabledesign = observedvariablesnop[variablenumber]
    observedvariablenames = names(observedvariables)

    # Start matching process:
    matchresult = as.data.frame(matrix(ncol=0, nrow=0))
    for(i in 1:length(observedvariablenames)){
      tempmatchresult = grepl(observedvariablenames[i], variabledesign)
      names(tempmatchresult) = "matchtf"
      matchresult = rbind(matchresult, tempmatchresult)
    }
    names(matchresult) = "matchtf"
    matchresult$variableid = 1:lengths(matchresult)
    matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
    matchresult[,2] # the match number from observedvariablenames

    dfobservedvariables = as.data.frame(observedvariables)
    levelinformation = as.character(dfobservedvariables[1, matchresult[,2]])

    leveltotal = unlist(strsplit(levelinformation, split = "rep"))[2]
    leveltotal = unlist(strsplit(leveltotal, split = "[(]"))[2]
    leveltotal = unlist(strsplit(leveltotal, split = "1:"))[2]
    leveltotal = unlist(strsplit(leveltotal, split = ","))[1]
    leveltotal = as.numeric(leveltotal)

    levelname = colnames(dfobservedvariables[matchresult[,2]])

    tempobserved = as.data.frame(1:leveltotal)
    names(tempobserved) = levelname

    bits = unlist(strsplit(as.character(variabledesign), levelname))
    togetherbits = paste(bits[1], leveltotal, bits[2])
    tempobserved$tempname = eval(parse(text = togetherbits))

    names(tempobserved) = c(levelname, variablename)

    dataset = merge.data.frame(dataset ,tempobserved, by = levelname)

  }
  variablenumber = variablenumber +1
}

assign("specificationdf", dataset, envir=globalenv())
}




effectsdf = function(variablespecification, effectspecifications){

  if(is.null(observedvariables$participant) == TRUE){
    print("error: participant details missing")
  } else {
    pptdesign = observedvariables$participant
    dataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
    dataset$participant = eval(parse(text = pptdesign))
  }

  observedvariablesnop = subset(observedvariables, names(observedvariables) != "participant")
  variablenumber = 1


  for(i in observedvariablesnop){

    tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
    variablename = names(observedvariablesnop)[variablenumber]
    tempnames = c("observation", variablename)

    # For repeat vs observed structures:
    if(unlist(strsplit(as.character(i), split = "\\("))[1] == "rep"){
      # For repeat-strucure variables:
      variabledesign = observedvariablesnop[variablenumber]
      tempdataset$tempname = eval(parse(text = variabledesign))
      names(tempdataset) = tempnames
      dataset = merge.data.frame(dataset, tempdataset, by ="observation")

    } else {
      # For observed-structure variables:
      variabledesign = observedvariablesnop[variablenumber]
      observedvariablenames = names(observedvariables)

      # Start matching process:
      matchresult = as.data.frame(matrix(ncol=0, nrow=0))
      for(i in 1:length(observedvariablenames)){
        tempmatchresult = grepl(observedvariablenames[i], variabledesign)
        names(tempmatchresult) = "matchtf"
        matchresult = rbind(matchresult, tempmatchresult)
      }
      names(matchresult) = "matchtf"
      matchresult$variableid = 1:lengths(matchresult)
      matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
      matchresult[,2] # the match number from observedvariablenames

      dfobservedvariables = as.data.frame(observedvariables)
      levelinformation = as.character(dfobservedvariables[1, matchresult[,2]])

      leveltotal = unlist(strsplit(levelinformation, split = "rep"))[2]
      leveltotal = unlist(strsplit(leveltotal, split = "[(]"))[2]
      leveltotal = unlist(strsplit(leveltotal, split = "1:"))[2]
      leveltotal = unlist(strsplit(leveltotal, split = ","))[1]
      leveltotal = as.numeric(leveltotal)

      levelname = colnames(dfobservedvariables[matchresult[,2]])

      tempobserved = as.data.frame(1:leveltotal)
      names(tempobserved) = levelname

      bits = unlist(strsplit(as.character(variabledesign), levelname))
      togetherbits = paste(bits[1], leveltotal, bits[2])
      tempobserved$tempname = eval(parse(text = togetherbits))

      names(tempobserved) = c(levelname, variablename)

      dataset = merge.data.frame(dataset ,tempobserved, by = levelname)

    }
    variablenumber = variablenumber +1
  }

effectnumber = 1

for(i in effectvariables){

  tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
  effectname = names(effectvariables)[effectnumber]
  tempnames = c("observation", effectname)

  if(grepl("r", i) == "FALSE"){
    tempdataset$tempname = as.numeric(i)
    names(tempdataset) = tempnames
    dataset = cbind(dataset, tempdataset[2])
  } else {
    # Start matching process:
    matchresult = as.data.frame(matrix(ncol=0, nrow=0))
    for(j in 1:length(observedvariablenames)){
      tempmatchresult = grepl(observedvariablenames[j], i)
      names(tempmatchresult) = "matchtf"
      matchresult = rbind(matchresult, tempmatchresult)
    }
    names(matchresult) = "matchtf"
    matchresult$variableid = 1:lengths(matchresult)
    matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
    matchresult[,2] # the match number from observedvariablenames

    tempeffect = as.data.frame(unique(dataset[,observedvariablenames[matchresult[,2]]]))

    numberlevels = as.numeric(lengths(tempeffect))

    bits = unlist(strsplit(as.character(i), observedvariablenames[matchresult[,2]]))
    togetherbits = paste(bits[1], numberlevels, bits[2])
    tempeffect$tempname = eval(parse(text = togetherbits))

    names(tempeffect) = c(observedvariablenames[matchresult[,2]], effectname)

    dataset = merge.data.frame(dataset ,tempeffect, by = observedvariablenames[matchresult[,2]])

  }
  effectnumber = effectnumber +1
}
assign("specificationdf", dataset, envir=globalenv())
}



simdf = function(variablespecification, effectspecifications, outcomespecification){

    if(is.null(observedvariables$participant) == TRUE){
      print("error: participant details missing")
    } else {
      pptdesign = observedvariables$participant
      dataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
      dataset$participant = eval(parse(text = pptdesign))
    }

    observedvariablesnop = subset(observedvariables, names(observedvariables) != "participant")
    variablenumber = 1


    for(i in observedvariablesnop){

      tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
      variablename = names(observedvariablesnop)[variablenumber]
      tempnames = c("observation", variablename)

      # For repeat vs observed structures:
      if(unlist(strsplit(as.character(i), split = "\\("))[1] == "rep"){
        # For repeat-strucure variables:
        variabledesign = observedvariablesnop[variablenumber]
        tempdataset$tempname = eval(parse(text = variabledesign))
        names(tempdataset) = tempnames
        dataset = merge.data.frame(dataset, tempdataset, by ="observation")

      } else {
        # For observed-structure variables:
        variabledesign = observedvariablesnop[variablenumber]
        observedvariablenames = names(observedvariables)

        # Start matching process:
        matchresult = as.data.frame(matrix(ncol=0, nrow=0))
        for(i in 1:length(observedvariablenames)){
          tempmatchresult = grepl(observedvariablenames[i], variabledesign)
          names(tempmatchresult) = "matchtf"
          matchresult = rbind(matchresult, tempmatchresult)
        }
        names(matchresult) = "matchtf"
        matchresult$variableid = 1:lengths(matchresult)
        matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
        matchresult[,2] # the match number from observedvariablenames

        dfobservedvariables = as.data.frame(observedvariables)
        levelinformation = as.character(dfobservedvariables[1, matchresult[,2]])

        leveltotal = unlist(strsplit(levelinformation, split = "rep"))[2]
        leveltotal = unlist(strsplit(leveltotal, split = "[(]"))[2]
        leveltotal = unlist(strsplit(leveltotal, split = "1:"))[2]
        leveltotal = unlist(strsplit(leveltotal, split = ","))[1]
        leveltotal = as.numeric(leveltotal)

        levelname = colnames(dfobservedvariables[matchresult[,2]])

        tempobserved = as.data.frame(1:leveltotal)
        names(tempobserved) = levelname

        bits = unlist(strsplit(as.character(variabledesign), levelname))
        togetherbits = paste(bits[1], leveltotal, bits[2])
        tempobserved$tempname = eval(parse(text = togetherbits))

        names(tempobserved) = c(levelname, variablename)

        dataset = merge.data.frame(dataset ,tempobserved, by = levelname)

      }
      variablenumber = variablenumber +1
    }


    effectnumber = 1


    for(i in effectvariables){

      tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
      effectname = names(effectvariables)[effectnumber]
      tempnames = c("observation", effectname)

      if(grepl("r", i) == "FALSE"){
        tempdataset$tempname = as.numeric(i)
        names(tempdataset) = tempnames
        dataset = cbind(dataset, tempdataset[2])
      } else {
        # Start matching process:
        matchresult = as.data.frame(matrix(ncol=0, nrow=0))
        for(j in 1:length(observedvariablenames)){
          tempmatchresult = grepl(observedvariablenames[j], i)
          names(tempmatchresult) = "matchtf"
          matchresult = rbind(matchresult, tempmatchresult)
        }
        names(matchresult) = "matchtf"
        matchresult$variableid = 1:lengths(matchresult)
        matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
        matchresult[,2] # the match number from observedvariablenames

        tempeffect = as.data.frame(unique(dataset[,observedvariablenames[matchresult[,2]]]))

        numberlevels = as.numeric(lengths(tempeffect))

        bits = unlist(strsplit(as.character(i), observedvariablenames[matchresult[,2]]))
        togetherbits = paste(bits[1], numberlevels, bits[2])
        tempeffect$tempname = eval(parse(text = togetherbits))

        names(tempeffect) = c(observedvariablenames[matchresult[,2]], effectname)

        dataset = merge.data.frame(dataset ,tempeffect, by = observedvariablenames[matchresult[,2]])


      }
      effectnumber = effectnumber +1
    }

    noutcome = 1

  for(k in outcomegeneration){

    tempsplit.b = strsplit(as.character(k), split = "\\(")
    x = unlist(tempsplit.b)[1]
    found = find(x)


    if(length(found) == 0){
      # do nothing
    } else {

      if(grepl("package", found) == TRUE){

        tempsplit = strsplit(as.character(k), split = "\\(")
        tempsplit = strsplit(as.character(unlist(tempsplit)[2]), split = "\\)")
        tempsplit = strsplit(as.character(tempsplit), split = ",")



        for(l in unlist(tempsplit)){

            if(grepl("dataset", l) == TRUE){

            pname = unlist(strsplit(l, split = "\\$"))[2]
            matchresult = as.data.frame(grepl(pname, names(outcomegeneration)))

            names(matchresult) = "matchtf"

            matchresult$outcomeid = 1:lengths(matchresult)
            matchresult = subset(matchresult, matchresult$matchtf == "TRUE")

            if(is.null(outcomegeneration[matchresult[,2]]) == TRUE){
              print("error: missing outcome generation parameter")
            } else {
              tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))

              tempdataset$tempname = eval(parse(text = outcomegeneration[matchresult[,2]]))

              tempnames = c("observation", names(outcomegeneration)[matchresult[,2]])
              names(tempdataset) = tempnames

              dataset = cbind(dataset, tempdataset[2])

              }


            } else {
            # do nothing
          }
        }

        tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
        nobs = lengths(tempdataset)
        datagencommand = sub("observation", nobs, as.character(k))

        tempdataset$tempname = eval(parse(text = datagencommand))

        tempnames = c("observation", names(outcomegeneration[noutcome]))
        names(tempdataset) = tempnames

        dataset = cbind(dataset, tempdataset[2])

      } else {
        print("package for data generating function not found - check install")
      }
      #donothing
      #this outcomedistribution[n] will be used elsewhere
    }
    noutcome = noutcome +1
  }

    assign("simulateddf", dataset, envir=globalenv())
  }




simmodelfits = function(variablespecification, effectspecifications, outcomespecification, analyticmodel, ntimes){

   modelfits = as.data.frame(matrix(ncol=0, nrow=0))
   nmodel = 1

   for(o in 1:ntimes){

   if(is.null(observedvariables$participant) == TRUE){
     print("error: participant details missing")
   } else {
     pptdesign = observedvariables$participant
     dataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
     dataset$participant = eval(parse(text = pptdesign))
   }

   observedvariablesnop = subset(observedvariables, names(observedvariables) != "participant")
   variablenumber = 1


   for(i in observedvariablesnop){

     tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
     variablename = names(observedvariablesnop)[variablenumber]
     tempnames = c("observation", variablename)

     # For repeat vs observed structures:
     if(unlist(strsplit(as.character(i), split = "\\("))[1] == "rep"){
       # For repeat-strucure variables:
       variabledesign = observedvariablesnop[variablenumber]
       tempdataset$tempname = eval(parse(text = variabledesign))
       names(tempdataset) = tempnames
       dataset = merge.data.frame(dataset, tempdataset, by ="observation")

     } else {
       # For observed-structure variables:
       variabledesign = observedvariablesnop[variablenumber]
       observedvariablenames = names(observedvariables)

       # Start matching process:
       matchresult = as.data.frame(matrix(ncol=0, nrow=0))
       for(i in 1:length(observedvariablenames)){
         tempmatchresult = grepl(observedvariablenames[i], variabledesign)
         names(tempmatchresult) = "matchtf"
         matchresult = rbind(matchresult, tempmatchresult)
       }
       names(matchresult) = "matchtf"
       matchresult$variableid = 1:lengths(matchresult)
       matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
       matchresult[,2] # the match number from observedvariablenames


       dfobservedvariables = as.data.frame(observedvariables)
       levelinformation = as.character(dfobservedvariables[1, matchresult[,2]])

       leveltotal = unlist(strsplit(levelinformation, split = "rep"))[2]
       leveltotal = unlist(strsplit(leveltotal, split = "[(]"))[2]
       leveltotal = unlist(strsplit(leveltotal, split = "1:"))[2]
       leveltotal = unlist(strsplit(leveltotal, split = ","))[1]
       leveltotal = as.numeric(leveltotal)


       levelname = colnames(dfobservedvariables[matchresult[,2]])

       tempobserved = as.data.frame(1:leveltotal)
       names(tempobserved) = levelname

       bits = unlist(strsplit(as.character(variabledesign), levelname))
       togetherbits = paste(bits[1], leveltotal, bits[2])
       tempobserved$tempname = eval(parse(text = togetherbits))

       names(tempobserved) = c(levelname, variablename)


       dataset = merge.data.frame(dataset ,tempobserved, by = levelname)

     }
     variablenumber = variablenumber +1
   }


   effectnumber = 1

   for(i in effectvariables){

     tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
     effectname = names(effectvariables)[effectnumber]
     tempnames = c("observation", effectname)

     if(grepl("r", i) == "FALSE"){
       tempdataset$tempname = as.numeric(i)
       names(tempdataset) = tempnames
       dataset = cbind(dataset, tempdataset[2])
     } else {

       # Start matching process:
       matchresult = as.data.frame(matrix(ncol=0, nrow=0))
       for(j in 1:length(observedvariablenames)){
         tempmatchresult = grepl(observedvariablenames[j], i)
         names(tempmatchresult) = "matchtf"
         matchresult = rbind(matchresult, tempmatchresult)
       }
       names(matchresult) = "matchtf"
       matchresult$variableid = 1:lengths(matchresult)
       matchresult = subset(matchresult, matchresult$matchtf == "TRUE")
       matchresult[,2] # the match number from observedvariablenames


       tempeffect = as.data.frame(unique(dataset[,observedvariablenames[matchresult[,2]]]))

       numberlevels = as.numeric(lengths(tempeffect))

       bits = unlist(strsplit(as.character(i), observedvariablenames[matchresult[,2]]))
       togetherbits = paste(bits[1], numberlevels, bits[2])
       tempeffect$tempname = eval(parse(text = togetherbits))

       names(tempeffect) = c(observedvariablenames[matchresult[,2]], effectname)

       dataset = merge.data.frame(dataset ,tempeffect, by = observedvariablenames[matchresult[,2]])

     }
     effectnumber = effectnumber +1
   }

   noutcome = 1

   for(k in outcomegeneration){

     tempsplit.b = strsplit(as.character(k), split = "\\(")
     x = unlist(tempsplit.b)[1]
     found = find(x)


     if(length(found) == 0){
       # do nothing
     } else {

       if(grepl("package", found) == TRUE){

         tempsplit = strsplit(as.character(k), split = "\\(")
         tempsplit = strsplit(as.character(unlist(tempsplit)[2]), split = "\\)")
         tempsplit = strsplit(as.character(tempsplit), split = ",")


         for(l in unlist(tempsplit)){

             if(grepl("dataset", l) == TRUE){

             pname = unlist(strsplit(l, split = "\\$"))[2]
             matchresult = as.data.frame(grepl(pname, names(outcomegeneration)))

             names(matchresult) = "matchtf"

             matchresult$outcomeid = 1:lengths(matchresult)
             matchresult = subset(matchresult, matchresult$matchtf == "TRUE")

             if(is.null(outcomegeneration[matchresult[,2]]) == TRUE){
               print("error: missing outcome generation parameter")
             } else {
               tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))

               tempdataset$tempname = eval(parse(text = outcomegeneration[matchresult[,2]]))

               tempnames = c("observation", names(outcomegeneration)[matchresult[,2]])
               names(tempdataset) = tempnames

               dataset = cbind(dataset, tempdataset[2])


               }


              } else {
             # do nothing
           }
         }

         tempdataset = as.data.frame(cbind(observation = 1:length(eval(parse(text = pptdesign)))))
         nobs = lengths(tempdataset)
         datagencommand = sub("observation", nobs, as.character(k))

         tempdataset$tempname = eval(parse(text = datagencommand))

         tempnames = c("observation", names(outcomegeneration[noutcome]))
         names(tempdataset) = tempnames

         dataset = cbind(dataset, tempdataset[2])

       } else {
         print("package for data generating function not found - check install")
       }
       #donothing
       #this outcomedistribution[n] will be used elsewhere
     }
     noutcome = noutcome +1
   }

  tempmodel = eval(parse(text = analyticmodel))

  modelfunction = unlist(strsplit(analyticmodel, split = "\\("))[1]

  #### lm ####
  if(modelfunction == "lm"){

    # Estimates
    tempsummary = summary(tempmodel)

    tempcoefs = as.data.frame(tempsummary$coefficients)
    tempcoefs$parameter = c(row.names(tempcoefs))
    tempcoefs$nmodel = nmodel

    # Model fit (&sigma)
    tempfit = as.data.frame(matrix(0))

    tempfit$sigma = tempsummary$sigma
    tempfit$rsqr = tempsummary$r.squared
    tempfit$adjrsqr = tempsummary$adj.r.squared
    tempfit$Fstatistic = tempsummary$fstatistic[1]
    tempfit$pofFstatistic = pf(tempsummary$fstatistic[1],tempsummary$fstatistic[2],tempsummary$fstatistic[3],lower.tail=FALSE)

    tempsummary$convergence = "NA" #there is no index of this in lm?

    tempcoefs = cbind(tempcoefs, rep(tempfit[2:length(tempfit)]))



  } else {
    # do nothing
  }

  #### glm ####
  if(modelfunction == "glm"){

    # Estimates
    tempsummary = summary(tempmodel)

    tempcoefs = as.data.frame(tempsummary$coefficients)
    tempcoefs$parameter = c(row.names(tempcoefs))
    tempcoefs$nmodel = nmodel

    # Model fit
    tempfit = as.data.frame(matrix(0))

    tempfit$AIC = tempsummary$aic
    tempfit$logLik = logLik(tempmodel)
    tempfit$deviance = tempsummary$deviance
    tempfit$df.resid = tempsummary$df.residual
    tempfit$convergence = ifelse(tempmodel["converged"] == TRUE, 1, 0)

    tempcoefs = cbind(tempcoefs, rep(tempfit[2:length(tempfit)]))

  } else {
    # do nothing
  }

  #### lmer ####
  if(modelfunction == "lmer"){

    # Estimates
    tempsummary = summary(tempmodel)

    tempcoefs = as.data.frame(tempsummary$coefficients)
    tempcoefs$parameter = c(row.names(tempcoefs))

    # 'Random' effects
    varcorrraneffs = as.data.frame(VarCorr(tempmodel))
    ranefnrow = nrow(varcorrraneffs)
    ranefncol = ncol(varcorrraneffs)

    ranefdfpnames = as.data.frame(matrix(ncol=0, nrow=0))

    for(rownumber in 1:ranefnrow){
      templist = as.list(varcorrraneffs[rownumber, 1:(ranefncol-2)], sep = ".")
      tempranefdf = as.data.frame(do.call("paste", c(templist, sep = ".")))
      ranefdfpnames = rbind(ranefdfpnames, tempranefdf)
    }
    names(ranefdfpnames) = "parameter"

    ranefdf = as.data.frame(cbind("Estimate" = varcorrraneffs$vcov,
                                  "Std. Error" = varcorrraneffs$sdcor,
                                  "t value" = NA))
    ranefdf = cbind(ranefdf, ranefdfpnames)
    tempcoefs = rbind(tempcoefs, ranefdf)

    # Sigma
    sigmadf = as.data.frame(cbind("Estimate" = tempsummary$sigma,
                                  "Std. Error" = NA,
                                  "t value" = NA,
                                  "parameter" = "sigma"))

    tempcoefs = rbind(tempcoefs, sigmadf)

    tempcoefs$nmodel = nmodel

    # Model fit
    tempcoefs$AIC = AIC(tempmodel)
    tempcoefs$BIC = BIC(tempmodel)

    # Convergence
    tempcoefs$convergence = ifelse(length(tempsummary$fitMsgs) == 0, 1, 0)

    assign("tempcoefs", tempcoefs, envir=globalenv())
  } else {
    # do nothing
  }

  #### glmer ####
  if(modelfunction == "glmer"){

    # Estimates
    tempsummary = summary(tempmodel)
    
    tempcoefs = as.data.frame(tempsummary$coefficients)
    tempcoefs$parameter = c(row.names(tempcoefs))
    
    # 'Random' effects
    varcorrraneffs = as.data.frame(VarCorr(tempmodel))
    ranefnrow = nrow(varcorrraneffs)
    ranefncol = ncol(varcorrraneffs)
    
    ranefdfpnames = as.data.frame(matrix(ncol=0, nrow=0))
    
    for(rownumber in 1:ranefnrow){
      templist = as.list(varcorrraneffs[rownumber, 1:(ranefncol-2)], sep = ".")
      tempranefdf = as.data.frame(do.call("paste", c(templist, sep = ".")))
      ranefdfpnames = rbind(ranefdfpnames, tempranefdf)
    }
    names(ranefdfpnames) = "parameter"
    
    ranefdf = as.data.frame(cbind("Estimate" = varcorrraneffs$vcov,
                                  "Std. Error" = varcorrraneffs$sdcor,
                                  "z value" = NA,
                                  "Pr(>|z|)" = NA))
    ranefdf = cbind(ranefdf, ranefdfpnames)
    tempcoefs = rbind(tempcoefs, ranefdf)
    
    tempcoefs$nmodel = nmodel
    
    # Model fit
    tempcoefs$AIC = AIC(tempmodel)
    tempcoefs$BIC = BIC(tempmodel)
    
    # Convergence
    tempcoefs$convergence = ifelse(is.null(tempsummary$optinfo$conv$lme4$messages) == TRUE, 1, 0)
        
    assign("tempcoefs", tempcoefs, envir=globalenv())
    
  } else {
    # do nothing
  }

  #### clm ####
  if(modelfunction == "clm"){
    print("model function 'clm' not yet supported")
  } else {
    # do nothing
  }

  #### clmm ####
  if(modelfunction == "clmm"){
    print("model function 'clmm' not yet supported")
  } else {
    # do nothing
  }

  #### brm ####
  if(modelfunction == "brm"){
    tempsummary = summary(tempmodel)

    # find the family:- this may be generically formatted?
    familyname = as.character(tempmodel$family[1])

    # Estimates
    tempcoefs = as.data.frame(tempsummary$fixed)

    tempcoefs = rbind(tempcoefs, tempsummary$spec_pars)
    tempcoefs$parameter = c(row.names(tempcoefs))

    # 'Random' (group-level effects)
    if(is.null(tempsummary$group) == FALSE){
      for(m in 1:length(tempsummary$group)){
        tempraneff = as.data.frame(tempsummary$random[m])
        tempraneff$parameter = paste(tempsummary$group[m] ,c(row.names(tempraneff)), sep = ".")
        names(tempraneff) = colnames(tempcoefs)
        tempcoefs = rbind(tempcoefs, tempraneff)
      }
    } else {}

    tempcoefs$nmodel = nmodel

    ### Convergence ###
    # extract divergent transitions
    tempcoefs$divtrans = sum(subset(nuts_params(tempmodel), Parameter == "divergent__")$Value)


  } else {
    # do nothing
  }

  # Join model fits
  modelfits = rbind(modelfits, tempcoefs)
  nmodel = nmodel+1

  }
  assign("simulatedmodelfits", modelfits, envir=globalenv())
 }