box_1_ceiling_effect.Rmd

# Ceiling Effect Simulation

## Generate data

```{r}
set.seed(12345)
library(tidyverse)
n <- 500
ceiling <- tibble(
  Moderator = c(rep(0, n/5*4), rep(1, n/5)),
  pred = 1.4 * Moderator + 0.7 * rnorm(n),
  out = 1.1 * Moderator + 0.6 * pred + 0.5 * rnorm(n)
) %>% 
  mutate(
  Moderator = as.factor(Moderator)
  )
sd(ceiling$out)
```

## Plot latent variable
```{r}
library(ggplot2)
theme_set(ggthemes::theme_clean())
p1 <- ggplot(data = ceiling, aes(x = pred, y = out, color = Moderator)) +
  geom_point(alpha = .3) +
  geom_smooth(method = "lm", formula = y~x) +
  coord_cartesian(ylim = c(min(ceiling$out), max(ceiling$out))) +
  xlab("Predictor") +
  ylab("Outcome") +
  scale_color_viridis_d(guide = F, begin = 0.3, end = 0.8) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
p1

qplot(residuals(lm(out ~ pred*Moderator, ceiling)))
qplot(ceiling$out)
```

## Plot censored variable
```{r}
ceiling$cens <- ceiling$out
ceiling$cens[ceiling$cens > 2] <- 2
p2 <- ggplot(data = ceiling, aes(x = pred, y = cens, color = Moderator)) +
  geom_point(alpha = .3) +
  geom_smooth(method = "lm", formula = y~x) +
  coord_cartesian(ylim = c(min(ceiling$out), max(ceiling$out))) +
  xlab("Predictor") +
  ylab("Outcome with ceiling") +
  scale_color_viridis_d(begin = 0.3, end = 0.8) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
p2

reg <- lm(cens ~ pred*Moderator, ceiling)
summary(reg)
qplot(residuals(reg))
qplot(ceiling$cens)

```

## Plot ordinal variable
```{r}
ceiling$ordinal <- ceiling$out
ceiling$ordinal[ceiling$ordinal > 2] <- 2
ceiling$ordinal[ceiling$ordinal < -2] <- -2
ceiling$ordinal <- round(ceiling$ordinal)
p3 <- ggplot(data = ceiling, aes(x = pred, y = ordinal, color = Moderator)) +
  geom_point(alpha = .3) +
  geom_smooth(method = "lm", formula = y~x) +
  coord_cartesian(ylim = c(min(ceiling$out), max(ceiling$out))) +
  scale_color_viridis_d(guide = F, begin = 0.3, end = 0.8) +
  xlab("Predictor") +
  ylab("Outcome measured on Likert scale") +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
p3

reg <- lm(ordinal ~ pred*Moderator, ceiling)
summary(reg)
qplot(residuals(reg))
qplot(ceiling$ordinal)
```


## Multi-plot
```{r}
library(cowplot)
plot_grid(p1, p2 + theme(legend.title = element_text(size = 9), legend.background = element_blank(), legend.position = c(0.5,0.85), legend.justification = "right"), p3, labels = c("A", "B", "C"), nrow = 1, align = "hv", axis = "lb")
ggsave("plots/ceiling_multiplot.png", width = 7, height = 3)
ggsave("plots/ceiling_multiplot.pdf", width = 7, height = 3)
```

## Fit models with improper measurement models/distributional assumptions
```{r}
library(brms)

# original data: no interaction
summary(brm(out ~ pred*Moderator, data = ceiling,  cores = 4, backend = "cmdstanr", file = "models/ceiling/original_no_interaction"))

# with ceiling effect: interaction
summary(uncensored <-brm(cens ~ pred*Moderator, data = ceiling,  cores = 4, backend = "cmdstanr", file = "models/ceiling/ceiling_effect"))

# with ordinal variable: interaction
summary(brm(ordinal ~ pred*Moderator, data = ceiling, cores = 4, backend = "cmdstanr", file = "models/ceiling/ordinal_as_gaussian"))


qplot(ceiling$ordinal)
qplot(ceiling$cens)
```

## Fit proper models
Now let's set up a better model in brms.
First we set a variable that tells us which values are censored and in what way

```{r}
ceiling$censored <- "none"
ceiling$censored[ceiling$cens >= 2] <- "right"
table(ceiling$censored)
```

Run model
```{r}
censored <- brm(cens | cens(censored) ~ pred * Moderator, data = ceiling,  cores = 4, backend = "cmdstanr", file = "models/ceiling/censored")
summary(censored)
loo(uncensored, censored)
```


Now we treat the ordinal variable as ordinal instead of Gaussian. We're assuming a cumulative distribution with equidistant thresholds (because that's what we simulated).

```{r}
ceiling$ordinal <- factor(ceiling$ordinal, ordered = TRUE)
ordinal <- brm(ordinal ~ pred*Moderator, data = ceiling, family = cumulative(),  cores = 4, backend = "cmdstanr", file = "models/ceiling/cumulative")
summary(ordinal)
```