0

I am not familiar with nonlinear regression and would appreciate some help with running an exponential decay model in R. Please see the graph for how the data looks like. My hunch is that an exponential model might be a good choice. I have one fixed effect and one random effect. y ~ x + (1|random factor). How to get the starting values for the exponential model (please assume that I know nothing about nonlinear regression) in R? How do I subsequently run a nonlinear model with these starting values? Could anyone please help me with the logic as well as the R code?

As I am not familiar with nonlinear regression, I haven't been able to attempt it in R.

raw plot

Revathe Thillaikumar
  • 1
  • 1

1 Answers1

0

The correct syntax will depend on your experimental design and model but I hope to give you a general idea on how to get started.

We begin by generating some data that should match the type of data you are working with. You had mentioned a fixed factor and a random one. Here, the fixed factor is represented by the variable treatment and the random factor is represented by the variable grouping_factor.

library(nlraa)
library(nlme)
library(ggplot2)

## Setting this seed should allow you to reach the same result as me
set.seed(3232333)

example_data <- expand.grid(treatment = c("A", "B"),
                            grouping_factor = c('1', '2', '3'),
                            replication = c(1, 2, 3),
                            xvar = 1:15)

The next step is to create some "observations". Here, we use an exponential function y=a∗exp(c∗x) and some random noise to create some data. Also, we add a constant to treatment A just to create some treatment differences.

example_data$y <- ave(example_data$xvar, example_data[, c('treatment', 'replication', 'grouping_factor')],
                        FUN = function(x) {expf(x = x, 
                                                a = 10, 
                                                c = -0.3) + rnorm(1, 0, 0.6)})

example_data$y[example_data$treatment == 'A'] <-  example_data$y[example_data$treatment == 'A'] + 0.8

All right, now we start fitting the model.

## Create a grouped data frame
exampleG <- groupedData(y ~ xvar|grouping_factor, data = example_data)


## Fit a separate model to each groupped level
fitL <- nlsList(y ~ SSexpf(xvar, a, c), data = exampleG)

## Grab the coefficients of the general model
fxf <- fixed.effects(fit1)

## Add treatment as a fixed effect. Also, use the coeffients from the previous
## regression model as starting values.
fit2 <- update(fit1, fixed = a + c ~ treatment,
               start = c(fxf[1], 0,
                         fxf[2], 0))

Looking at the model output, it will give you information like the following:

Nonlinear mixed-effects model fit by maximum likelihood
  Model: y ~ SSexpf(xvar, a, c) 
  Data: exampleG 
       AIC      BIC    logLik
  475.8632 504.6506 -229.9316

Random effects:
 Formula: list(a ~ 1, c ~ 1)
 Level: grouping_factor
 Structure: General positive-definite, Log-Cholesky parametrization
              StdDev       Corr  
a.(Intercept) 3.254827e-04 a.(In)
c.(Intercept) 1.248580e-06 0     
Residual      5.670317e-01       

Fixed effects:  a + c ~ treatment 
                  Value Std.Error  DF   t-value p-value
a.(Intercept)  9.634383 0.2189967 264  43.99329  0.0000
a.treatmentB   0.353342 0.3621573 264   0.97566  0.3301
c.(Intercept) -0.204848 0.0060642 264 -33.77976  0.0000
c.treatmentB  -0.092138 0.0120463 264  -7.64867  0.0000
 Correlation: 
              a.(In) a.trtB c.(In)
a.treatmentB  -0.605              
c.(Intercept) -0.785  0.475       
c.treatmentB   0.395 -0.792 -0.503

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.93208903 -0.34340037  0.04767133  0.78924247  1.95516431 

Number of Observations: 270
Number of Groups: 3

Then, if you wanted to visualize the model fit, you could do the following.

## Here we store the model predictions for visualization purposes
predictionsDf <- cbind(example_data,
                       predict_nlme(fit2, interval = 'conf'))

## Here we make a graph to check it out
ggplot()+
  geom_ribbon(data = predictionsDf, 
              aes( x = xvar , ymin = Q2.5, ymax = Q97.5, fill = treatment),
              color = NA, alpha = 0.3)+
  geom_point(data = example_data, aes( x = xvar, y = y, col = treatment))+
geom_line(data = predictionsDf, aes(x = xvar, y = Estimate, col = treatment), size = 1.1)

This shows the model fit.

clsantos
  • 26
  • 2