Piethon
Piethon

Reputation: 247

Can't place constraints on nlme model coefficients

I'm running the code below to fit a decaying exponential function to some data, using the nlme package in R:

    reprex_df <- structure(list(PARTICIPANT_ID = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 10L, 10L, 10L, 10L), levels = c("s068133728", "s110540230", 
"s232476227", "s284908738", "s425227476", "s441171402", "s569306179", 
"s584875578", "s670756648", "s763441377"), class = "factor"), 
    GUILT_COMB = c(7, 1, 2, 1, 2, 8, 8, 2, 3, 2, 2, 2, 1, 2, 
    8, 6, 7, 8, 6, 2, 3, 2, 3, 5, 3, 3, 1, 10, 9, 7, 6, 8, 7, 
    7, 3, 0, 0, 0, 0, 0, 3, 1, 4, 10, 3, 0, 2, 1, 0, 0, 0, 4, 
    1, 4, 3, 2, 1, 1, 4, 0, 0, 0, 0, 0, 1, 2, 1, 0, 3, 0, 0, 
    1, 0, 0, 2, 2, 2, 0, 0, 3, 3, 0, 0, 0, 0, 3, 1, 0, 0, 0, 
    0, 0, 2, 1, 5, 6, 3, 6, 6, 5, 6, 4, 3, 4, 6, 6, 10, 5, 1, 
    3, 2, 2, 1, 4, 4, 6, 3, 4, 6, 4, 6, 6, 5, 3, 2, 2, 6, 4, 
    3, 5, 5, 5, 6, 4, 3, 8, 3, 2, 3, 2, 2, 2, 2, 2, 3, 3, 3, 
    2, 6, 2, 2, 2, 2, 2, 3, 2, 2, 1, 2, 2, 2, 2, 2, 3, 2, 2, 
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2, 2, 2, 2, 3, 
    2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 2, 3, 2, 1, 2, 2, 2, 2, 2, 
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 4, 2, 3, 2, NA, 2, 3, 
    3, 2, 2, 5, 7, 2, 6, 5, 8, 7, 10, 4, 5, 3, 5, 5, 8, 5, 5, 
    3, 3, 3, 3, 7, 6, 5, 10, 7, 7, 7, 7, 7, 7, 7, 7, 8, 10, 7, 
    10, 8, 7, 7, 6, 8, 8, 7, 8, 7, 8, 7, 7, 8, 7, 7, 7, 7, 8, 
    7, 8, 9, 8, 8, 8, 8, 5, 10, 8, 8, 8, 8, 7, 4, 2, 2, 2, 0, 
    0, 1, 0, 0, 2, 2, 6, 3, 2, 2, 0, 2), streaklength = c(0, 
    5.19972222222222, 6.78305555555556, 10.0163888888889, 14.6497222222222, 
    15.8830555555556, 0, 1.24972222222222, 3.88305555555556, 
    7.73305555555556, 19.0497222222222, 22.4830555555556, 25.1497222222222, 
    28.6997222222222, 0, 5.91638888888889, 8.63305555555556, 
    0, 8.34972222222222, 10.1997222222222, 14.0163888888889, 
    16.6997222222222, 19.9997222222222, 32.0330555555556, 37.2997222222222, 
    40.5830555555556, 43.4330555555556, 0, 0.05, 0.116666666666667, 
    0.0833333333333333, 0.433333333333333, 8.53333333333333, 
    10.4, 14.1166666666667, 34.2166666666667, 44.5, 56.2, 58.5666666666667, 
    61.3833333333333, 0.316666666666667, 0.15, 0, 0, 0, 4.68277777777778, 
    0, 3.68277777777778, 6.58277777777778, 18.1494444444444, 
    21.3661111111111, 0, 0.399444444444444, 3.54944444444444, 
    6.61611111111111, 9.53277777777778, 22.7994444444444, 25.1327777777778, 
    0, 1.98277777777778, 14.8994444444444, 17.9827777777778, 
    21.6827777777778, 24.6827777777778, 27.2994444444444, 39.6494444444444, 
    42.9327777777778, 47.3327777777778, 0, 0.532777777777778, 
    11.7494444444444, 15.7494444444444, 19.3161111111111, 21.8661111111111, 
    25.9494444444444, 0, 0, 27.8661111111111, 30.1494444444444, 
    0, 0, 13.1494444444444, 15.8661111111111, 22.3994444444444, 
    72.5327777777778, 0, 1.66611111111111, 19.7661111111111, 
    23.5827777777778, 25.9327777777778, 28.5161111111111, 31.0494444444444, 
    0, 0, 0, 0, 0.299444444444444, 3.78277777777778, 6.58277777777778, 
    8.46611111111111, 0, 7.78277777777778, 11.6494444444444, 
    13.9994444444444, 20.0661111111111, 0, 9.51555555555556, 
    0, 7.26638888888889, 10.8163888888889, 13.9497222222222, 
    29.0497222222222, 31.7830555555556, 1.13305555555556, 9.91638888888889, 
    0, 2.11638888888889, 6.66638888888889, 0, 14.0830555555556, 
    14.0663888888889, 39.3997222222222, 0, 8.21638888888889, 
    29.8830555555556, 32.9497222222222, 0, 2.61638888888889, 
    6.99972222222222, 9.94972222222222, 33.8830555555556, 0, 
    0, 20.4830555555556, 0, 0, 2.24944444444444, 4.84944444444444, 
    0, 1.18277777777778, 4.16611111111111, 7.03277777777778, 
    0, 1.24944444444444, 3.33277777777778, 6.49944444444444, 
    0, 6.99944444444444, 6.99944444444444, 2.18277777777778, 
    14.2161111111111, 17.7661111111111, 20.6994444444444, 0, 
    0.682777777777778, 3.74944444444444, 6.99944444444444, 0, 
    0.649444444444444, 3.66611111111111, 7.53277777777778, 0, 
    2.26611111111111, 4.73277777777778, 9.01611111111111, 19.5161111111111, 
    0, 0.966111111111111, 3.69944444444444, 7.78277777777778, 
    0, 2.91611111111111, 5.49944444444444, 8.69944444444444, 
    0, 1.61611111111111, 5.31611111111111, 7.81611111111111, 
    19.2827777777778, 22.9327777777778, 0, 1.76611111111111, 
    4.88277777777778, 7.53277777777778, 0, 3.66611111111111, 
    8.01611111111111, 18.5827777777778, 22.1827777777778, 0, 
    0.599444444444444, 2.99944444444444, 5.74944444444444, 0, 
    1.06611111111111, 6.89944444444444, 0, 1.93277777777778, 
    3.81611111111111, 0, 2.29944444444444, 5.41611111111111, 
    8.34944444444444, 0, 1.74944444444444, 4.31611111111111, 
    7.98277777777778, 0, 0, 2.08277777777778, 17.4661111111111, 
    20.3661111111111, 23.0327777777778, 25.5994444444444, 29.2827777777778, 
    0, 7.01611111111111, 0.382777777777778, 3.49944444444444, 
    0, 2.39944444444444, 7.59944444444444, 0, 1.58277777777778, 
    7.34944444444444, 0, 7.06611111111111, 7.03277777777778, 
    16.2827777777778, 0, 0.749444444444444, 7.03277777777778, 
    2.23277777777778, 0, 1.71611111111111, 7.04944444444444, 
    23.4161111111111, 26.1161111111111, 0, 7.01611111111111, 
    0.532777777777778, 1.43277777777778, 18.0661111111111, 21.3494444444444, 
    23.4661111111111, 25.9327777777778, 0, 7.06611111111111, 
    0.166111111111111, 7.03277777777778, 16.2494444444444, 19.2661111111111, 
    20.0661111111111, 7.04944444444444, 1.99944444444444, 0, 
    4.03277777777778, 6.11611111111111, 7.04944444444444, 0.116111111111111, 
    2.54944444444444, 0, 7.06611111111111, 10.7661111111111, 
    13.5161111111111, 35.2327777777778, 7.06611111111111, 0, 
    0.849444444444444, 6.99944444444444, 16.7161111111111, 0, 
    0.966111111111111, 7.24944444444444, 24.3161111111111, 0, 
    2.66611111111111, 4.24944444444444, 8.51611111111111, 0, 
    1.03277777777778, 3.11611111111111, 4.74944444444444, 0, 
    3.06611111111111, 5.66611111111111, 8.64944444444444, 20.2494444444444, 
    23.1161111111111, 25.5327777777778, 29.5494444444444, 31.6827777777778, 
    0, 0, 5.71472222222222, 8.24805555555556, 10.2980555555556, 
    25.3313888888889, 28.6147222222222, 30.8480555555556, 34.3313888888889, 
    37.9480555555556, 49.7313888888889, 52.4980555555556, 54.8313888888889, 
    0, 8.71472222222222, 13.7313888888889, 36.7480555555556, 
    60.4813888888889, 0)), row.names = c(NA, -311L), class = c("tbl_df", 
"tbl", "data.frame"))

model <- nlme(
    GUILT_COMB ~ a*(2-exp(b*streaklength))+c,
    data=reprex_df %>% filter(!is.na(streaklength) & !is.na(GUILT_COMB)),
    fixed = a + b + c ~ 1,
    random = c ~ 1,
    groups = ~ PARTICIPANT_ID,
    start = c(a = mean(reprex_df$GUILT_COMB, na.rm=TRUE), b = -0.01, c = mean(reprex_df$GUILT_COMB, na.rm=TRUE)),
    control = nlmeControl(opt='nlminb', msMaxIter=200, maxIter=500, tolerance=1, lower=c(0, -Inf, 0), upper=c(10, 0, 10))
)

But this returns model coefficients as follows:

> model$coefficients
$fixed
          a           b           c 
2.878675927 0.008859636 2.470470658 

Yet I've specified in the model call that the upper limit for b should be 0, using upper=c(10, 0, 10). Why is it exceeding this limit? In similar models, a and c also exceed their limits. So it appears that nlme is not recognizing the upper and lower arguments. But I believe I've followed the documentation correctly (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/nlminb.html)

Any ideas?

Many thanks

Upvotes: 0

Views: 47

Answers (0)

Related Questions