2

I'm modelling the salinity levels in four estuaries in relation to freshwater discharge from a large river (which impacts all 4 estuaries), where each salinity station is being included as a random effect (multiple sensors in each estuary).

My df looks like this: *Note = in reality, each station has dozens, if not hundreds, of observations each

estuary    station   salinity     discharge
A          s1        10.8         1185.18
A          s2        11.1         1551.91
A          s3        12.5         1749.98
B          s4        9.2          1100.65
B          s5        11.8         1692.67
C          s6        8.1          591.41
C          s7        8.7          742.88
D          s8        10.4         1098.1
D          s9        11.9         1174.6

I've fitted an exponential decay function to the data using the self-starting function SSasymp. THis model incorporates station as a random effect. Code below:

df$station <- as.factor(df$station)
df.g <- groupedData(salinity ~ discharge | station, data = df)
df.nlis <- nlsList(salinity ~ SSasymp(discharge, Asym, R0, lrc), data = df.g)

## Fit a nonlinear mixed model
df.nlme <- nlme(df.nlis, control=list(MaxIter = 200, msVerbose = T))

I would also like to run another model that incorporates estuary as a fixed effect, however I'm unsure how I could write the model. I've tried:

df.g <- groupedData(salinity ~ discharge + estuary | station, data = df)
df.nlis <- nlsList(salinity ~ SSasymp(discharge, Asym, R0, lrc) + estuary, data = df.g)
df.nlme <- nlme(df.nlis, control=list(MaxIter = 200, msVerbose = T))

But when it runs, I receive the message:

Error in names(start$fixed) <- names(cf) : attempt to set an attribute on 
NULL
Splash1199
  • 379
  • 3
  • 14
  • 1
    please share your data with `dput()` instead of a code block – Mike Sep 21 '21 at 13:19
  • 1
    IMO code blocks aren't bad -- it's easy enough to read them with `read.table(header=TRUE, text = "...")` but it would be useful to have an example with enough data points per group to allow the code to run – Ben Bolker Sep 22 '21 at 00:22
  • I didn't know that, thank you for the `read.table()` tip – Mike Sep 22 '21 at 02:18

1 Answers1

0

This doesn't work yet for my made-up example but it's where I would start:

df.nlis2 <- gnls(salinity ~ SSasymp(discharge, Asym, R0, lrc), data = df.g,
                 params = Asym + R0 + lrc ~ estuary,
                 start = rep(fixef(df.nlme), 4))
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453