Estimating different asymptotes given unknown time

Modeling
1

Hello all,

I’m far afield here. We’ve got some old, old data as part of a class that we’re having students analyze. No one involved in the collection of the data are around, and the data has some… issues.

I have a “handle” on the data only in the sense that I understand the generating process well enough to simulate it. But I’m having a lot of trouble determining how to model this simulation!

The code to simulate a fair facsimile is:

##### Simulation
# number of trials
N <- 100

# Time spent is a random, unknown number that increases monotonically by day
byInd <- runif(N)
Day <- sample(1:4, 100, replace=T)
byDay <- c(0.01, 0.05, 0.5, 1)
tTime <- byInd * byDay[Day]

# Asymptote varies by treatment.
Start <- 10
Int <- -2
byx1 <- c(-4, -2)
byx2 <- c(-3, 0)

# will be the variables for the simulation & data
x1 <- sample(c(1,2), 100, replace=T)
x2 <- sample(c(1,2), 100, replace=T)

# input to simulation
simDat <- data.frame(time=tTime, K=Start+byx1[x1]+byx2[x2]+Int)

# logistic decline
consumed <- function(Time, K, Start=10, r=3)	{
	Num <- Start * K * exp(r * Time)
	Denom <- (K - Start) + ( Start * exp(r * Time))
	return(Num / Denom)
}

# generate the simulated output
out <- apply(simDat, 1, function(x) consumed(x[1], x[2]))

# fake dataset with output & treatment & habitat types + individual values
fakeData <- data.frame(y=out, x1=as.factor(x1), x2=as.factor(x2), Ind=as.factor(1:N), day=Day)

This is quite different from what I’m used to dealing with, and I’d welcome any suggestions on the best approach to model it. So far the best I’ve got is:

Model <- bf(y | resp_trunc(lb=0, ub=10) ~ x1 + x2 + (1 | Ind), sigma ~ day + (1 | Ind), family=skew_normal())
bpriors <- get_prior(Model, data=fakeData)
npriors <- c(prior(normal(0,1), class=b), prior(normal(7,3), class=Intercept))
fit <- brm(Model, data=fakeData, prior=npriors, iter=5e3)

Which gets the intercept, x1, and x2 values kind of close, and at least the right relative order, but produces lots of divergent transitions and other mayhem and is largely the result of guessing wildly at a model, and not a good approach.

2

Howdy!

If you can simulate it, then it seems like you could use that generating process as the model then? I tried to substitute in your definition of K from simDat into the consumed() function definition you created and make a non-linear model formula in brms. Unfortunately it doesn’t converge, but maybe it will give you some ideas. But other than the runif(N), I don’t see any stochastic part to this data generating process?

fit_nl <- brm(
  bf(y ~ (Start^2 + byx1 + byx2 + Int * exp(r * Time)) / (byx1 + byx2 + Int + Start * exp(r * Time)),
     Start ~ 1, byx1 ~ 0 + x1, byx2 ~ 0 + x2, Int ~ 1, r ~ 1, Time ~ 1 + (1|Ind) + (1|day),
     nl = TRUE),
  data = fakeData, family = gaussian(),
  prior = c(
    prior(normal(10, 0.25), nlpar = "Start", coef = "Intercept"),
    prior(normal(3, 0.25), nlpar = "r", coef = "Intercept"),
    prior(normal(-4, 0.5), nlpar = "byx1", coef = "x11"),
    prior(normal(-2, 0.5), nlpar = "byx1", coef = "x12"),
    prior(normal(-3, 0.5), nlpar = "byx2", coef = "x21"),
    prior(normal(0, 0.5), nlpar = "byx2", coef = "x22"),
    prior(normal(-2, 0.5), nlpar = "Int", coef = "Intercept"),
    prior(normal(0, 0.2), lb=0, nlpar = "Time"),
    prior(normal(0, 0.1), nlpar = "Time", coef = "Intercept", group = "day", class = sd),
    prior(normal(0, 0.1), nlpar = "Time", coef = "Intercept", group = "Ind", class = sd),
    prior(normal(0, 0.5), class = sigma)
  ),
  control = list(adapt_delta = 0.95), cores=4
)

fit_nl
1 Like