Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(rethinking)
- library(brms)
- options(mc.cores = 4)
- # Import data
- data(Trolley)
- data <- Trolley
- # Subset some data to make process faster
- data <- subset(data, edu == "Bachelor's Degree" | edu == "Some College"| edu == "Master's Degree")
- data$edu <- factor(data$edu)
- levels(data$edu) <- c("bsc","msc","col")
- # Make order and response into fewer categories
- data$order <- round(data$order/6)
- data$response <- round(data$response/2.4)+1
- data$response <- ordered(data$response)
- # Make plot as wished for
- # Get usage for each dur per participant
- edus <- c("col","bsc","msc")
- cols <- brewer.pal(3, "Set2")
- par(mfrow=c(2,2))
- j <- 1
- for(j in 1:4){
- resp.trolley <- data %>%
- group_by(order, edu) %>%
- summarise(
- resp = mean(response==j))
- # plot raw
- test <- subset(resp.trolley, edu==edus[1])
- plot(jitter(test$order, 0.4),test$resp,lwd=3, pch=20, ylim=c(0,max(resp.trolley$resp)*1),
- col=cols[1],xlab="Duration", ylab = paste("Proportion Rating",j))
- lines(test$order,test$resp,
- col=cols[1],lwd=0.5, lty=2)
- if (j == 3) {
- # Add legend
- legend(x = "bottomright", # Position
- legend = c("Col", "Bsc","Msc"), # Legend texts
- lty = c(2), # Line types
- col = cols[1:3], # Line colors
- lwd = 1) # Line width
- }
- for(i in 2:3){
- test <- subset(resp.trolley, edu==edus[i])
- points(jitter(test$order, 0.4),test$resp,ylim=c(0,0.8),lwd=3, pch=20,
- col=cols[i])
- lines(jitter(test$order, 0.4),test$resp,ylim=c(0,0.8),
- col=cols[i],lwd=0.5, lty=2)}
- }
- # Make BRMS model
- data$response <- ordered(data$response)
- m1 <- brm(response ~ edu*mo(order) + (1|id), data= data,
- family = cumulative(threshold = "flexible"),
- chains = 2)
- # Plot
- conditional_effects(m1, categorical=T)
- # Trying to make for conditions
- conditions <- make_conditions(m1, "response")
- conditional_effects(m1,categorical=TRUE, conditions = conditions)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement