Selecting segments with trees

# This script simulates data and estimates how many segments it should be cut into
# using the tree package

# Load library
library(brms)
require(tree)

# Clear work space
rm(list=ls())

# True rates (Unkown in reality)
trueRates <- c(0.76,0.65,0.58,0.63,0.79,0.67,0.65,0.70)
category <- as.factor(1:8)

### Simulate data
# Generate data for 1 subject
subject <- 1
indvar <- rnorm(1,0,0.02) # Variance how good or bad is this individual?
results <- cbind(subject,(rbinom(25, 1, trueRates[1]+indvar)),category[1])
for(i in 2:8){results <- rbind(results,cbind(subject,(rbinom(25, 1, trueRates[i]+indvar)),category[i]))}
# Loop for rest
for(s in 2:20){
  indvar <- rnorm(1,0,0.02)
  for(i in 1:8){results <- rbind(results,cbind(s,(rbinom(25, 1, trueRates[i]+indvar)),category[i]))}}

# Tweak and name Data frame
results <- data.frame(results); names(results) <- c("Subject","Cor","Category")
results$Subject <- as.factor(results$Subject)
results$Category <- as.factor(results$Category)

# # Make model that takes Subject as random intercept (It will take some processing time)
# Model8Categories <- brm(Cor ~ Category + (1|Subject), family=bernoulli(link = "logit"), data = results,
#                  chains = 2,iter = 2000, warmup = 500)
# saveplot <- marginal_effects(Model8Categories)
#
# # Retrieve estimates
# round(saveplot[["Category"]][["estimate__"]],2)-trueRates
# # Fairly close estimates

# Use Tree to segment driven by data (But not subject as random effect)
meanpercat <- aggregate(Cor ~ Subject + Category, results, mean) # Find mean for each participant, each condition
tree.groups <-tree(Cor~ Category, data=meanpercat)
cv.groups <- cv.tree(tree.groups,FUN=prune.tree)

# How many segments are best according to CV?
round(sum(cv.groups$size[which(cv.groups$dev == min(cv.groups$dev))]/(length(which(cv.groups$dev == min(cv.groups$dev))))))
# Running this a few times gives variance from 2-5 segments being best