Hello,
I am looking at a nonlinear model (simple first order on “Treatment”,
top - (top-bottom)exp(-rateTreatment) where top is the asymptote and “bottom” is a constant lifting the curve away from 0). Think of this as similar to rate kinetics and the “Treatment” could just be time. In my case Treatment is a supplement given in an animal trial that results in this curve shape. Here is some synthetic data:
Synthetic Data
library(tidyverse)
library(brms)
library(ggmcmc)
set.seed(1)
# set up Treatment levels
Treatment <- c(0, 2.1, 4.5, 6.5)
#replicates at each Treatment level
rep = 1:50
# lets also have a grouping G, so 50 are by G1 and 50 by G2
Group <- c("A","B")
# set up first order curve
#curve parameters: Response ~ top - ((top-bottom)*exp(-rate * Treatment)
top = 14
bottom = 8
rate = 0.22
# this is set up as continuous Response in "Response"
SynthData <- expand_grid(Treatment, Group, rep) %>%
mutate(Response = top - (top-bottom)*exp(-rate*Treatment)) %>%
mutate(Response = ifelse(Group == "A",Response+0.5,Response )) %>% # +0.5 bias in A group
rowwise() %>%
mutate(Response = Response + rnorm(1,0,0.5)) %>%
ungroup() %>%
mutate(Control = "no")
So there is a grouping, perhaps one could think of the grouping factor as different animal houses resulting in a slight bias of around +0.5 for the “A” group. Secondly, I want to introduce a one level control. This could be an alternative formulation of the treatment added to the experiment at only one level to minimize costs of the trial:
# make a control group at a single Treatment level.
# This is just treatment 2.1 shifted up by Y +2
Control <- SynthData %>%
filter(Treatment == 2.1) %>%
filter(Group == "A") %>%
mutate(Control = "Yes") %>%
mutate(Response = Response+2)
SynthData <- rbind(SynthData, Control)
This lifts the response up by 2 units compared to the rest of the data.
Visualization:
ggplot(SynthData, aes(x=Treatment, y= Response))+
geom_point(aes(color = Group, shape = Control))+
theme_bw()
BRMS specification
# nonlinear, continuous scale for response
# Initial kicks
init_func <- function(chain_id=1) {
list ( b_top = array(14) ,
b_bottom = array(8.5) ,
b_rate = array(0.2),
b_Group = array(0.5),
b_Control = array(2))
}
# setting up chains
init_list <- list(
init_func(chain_id = 1),
init_func(chain_id = 2),
init_func(chain_id = 3),
init_func(chain_id = 4)
)
# NUTS configuration
control <- list(
adapt_engaged = TRUE,
adapt_delta = 0.9, #increased from default of 0.8
stepsize = 0.08, # 0.05 default
max_treedepth = 15 # was 15
)
prior1 <- prior(normal(14, 1), nlpar = "top") +
prior(normal(8.0, 1), nlpar = "bottom")+
prior(normal(0.2, 0.1), nlpar = "rate")+
prior(normal(0.5, 1), nlpar = "Group" )+
prior(normal(2, 1), nlpar = "Control" )
fit1 <- brm(bf(Response ~ top - ((top-bottom)*exp(-rate * Treatment)) + Group + Control, top + bottom + rate + Group + Control ~ 1, nl = TRUE),
data = SynthData,
prior = prior1,
#family = cumulative(link="logit", threshold="flexible"),
iter = 4000 , warmup = 2000, thin = 2,
chains = 4, cores=6,
init = init_list,
control = control,
#file = "./model/fit14"
)
# inspect results
summary(fit1)
Fitting results
Model seems to fit OK
Family: gaussian
Links: mu = identity; sigma = identity
Formula: Response ~ top - ((top - bottom) * exp(-rate * Treatment)) + Group + Control
top ~ 1
bottom ~ 1
rate ~ 1
Group ~ 1
Control ~ 1
Data: SynthData (Number of observations: 450)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 2;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
top_Intercept 12.81 0.63 11.60 14.03 1.00 3131 3358
bottom_Intercept 8.10 0.62 6.91 9.30 1.00 3160 3361
rate_Intercept 0.41 0.03 0.35 0.47 1.00 3734 3380
Group_Intercept -0.63 0.77 -2.18 0.84 1.00 3385 3329
Control_Intercept 0.89 0.77 -0.61 2.40 1.00 3427 3075
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.82 0.03 0.77 0.87 1.00 3921 3192
..but the control group is not at +2
Chains do all seem to converge and the model is taking my initial values:
# fit process visualization
ggs(fit1, burnin = TRUE) %>%
filter(Parameter == "b_Control_Intercept") %>%
mutate(chain = factor(Chain),
intercept = value) %>%
ggplot(aes(x = Iteration, y = intercept, color = chain)) +
geom_line()+
ggtitle("default settings")+
theme_bw()
Here is a visualization of the synthetic data overlayed with the predicted values from the model:
testPred <- predict(
fit1,
type = "response")
SynthData$Estimate <- testPred[,1]
#Visualize
ggplot(SynthData, aes(x=Treatment, y= Response))+
geom_point(aes(color = Group, shape = Control))+
geom_point(aes(x=Treatment, y=Estimate, color = Group, shape = Control), size = 6)+
theme_bw()
I think this must be a simple model specification problem when I moved over to running nonlinear models. My actual data is ordinal, but I wanted to make sure the non-linear and grouping is correct before I move onto a logit response dataset. Can anyone point me in the right direction?
Thanks!
