benwhalley

clear all
set seed 123

* recruit subjects (if only this easy...)
set obs 2000

* assign to groups
gen mbct =  mod(_n, 2)

* depression baseline score

benwhalley

participant__username |         N
----------------------+----------
           Alex_Fowke |         5
        Carolyn_Hinds |         3
     Donna_Rutherford |         4
       Ellie_Campbell |         3
        Hannah_Wilson |         0
        Laura_Brummer |         5
          Sandy_Waite |         5
           Sara_Clark |         5

benwhalley

// simulated data
// y = depression, t=timepoint (1 to 3), therapist=therapistcluster, id=patient 

// 3 level model: observations nested in people in therapists
mixed y t || therapist: t, nocons || id:, stddev

// same model (???) using stata bayesian model fit
bayesmh y i.id i.therapist i.therapist#c.t, likelihood(normal({var_residual})) noconstant ///
		prior({y:_cons}, normal(0,100)) ///
		prior({var_therapist}, igamma(.001, .001)) ///

benwhalley

mysim <- function(){
  des <- expand.grid(
    congruence=c(1,0),
    cuedresp=c("L", "R", "None"),
    trial=1:5,
    block=1:10,
    participant=1:20
  )
  samp <- des %>% rowwise() %>% mutate(y=rbinom(1,1,.03*congruence))

benwhalley

# effect on total variance
library(dplyr)

N=10000
prob <-  .25

e_var = 1
u_var = .25*e_var  # see http://www.wolframalpha.com/input/?i=.1%3Db%2F%28a%2Bb%29
icc = u_var/(u_var+e_var)

benwhalley

# http://stackoverflow.com/questions/21477040/reshape2-multiple-results-of-aggregation-function

dcastMult <- function(data, formula, value.var = "value", 
                   funs = list("min" = min, "max" = max)) {
  require(reshape2)
  if (is.null(names(funs)) | any(names(funs) == "")) stop("funs must be named")
  Form <- formula(formula)
  LHS <- as.character(Form[[2]])
  if (length(LHS) > 1) LHS <- LHS[-1]
  temp <- lapply(seq_along(funs), function(Z) {

benwhalley

# http://stackoverflow.com/questions/21477040/reshape2-multiple-results-of-aggregation-function

dcastMult <- function(data, formula, value.var = "value", 
                   funs = list("min" = min, "max" = max)) {
  require(reshape2)
  if (is.null(names(funs)) | any(names(funs) == "")) stop("funs must be named")
  Form <- formula(formula)
  LHS <- as.character(Form[[2]])
  if (length(LHS) > 1) LHS <- LHS[-1]
  temp <- lapply(seq_along(funs), function(Z) {

benwhalley

# calculate ICC for lmer model object
icc <- function(m){
  vc <- as.data.frame((VarCorr(m)))
  l <- vc$vcov
  data_frame(grp=vc$grp, icc=sapply(l, function(x){x/sum(l)}))
}

benwhalley

# http://stackoverflow.com/questions/11693599/alternative-to-expand-grid-for-data-frames
expand.grid.df <-
  function(...)
    Reduce(function(...)
      merge(..., by = NULL), list(...))

benwhalley

# https://en.wikipedia.org/wiki/Akaike_information_criterion#How_to_apply_AIC_in_practice
compareAIC <- function(..., models) {
  mods = c(list(...), models)
  aics <- sapply(mods, FUN = AIC)
  round(exp((min(aics) - aics) / 2), 2)
}