Correct GLM or NLS to model exponential model with response variable with positive and negative values

Question

I have been struggling to find the right way to model this dataset, this is a Data Frame with the dataset:

Example_Data <- structure(list(Y = c(
  -0.646528881822999, -0.5172231054584, -0.510079139913394,
  -0.402518490005774, -0.385693032482023, -0.372908622908623, -0.328852788959172,
  -0.323765626091208, -0.323360507610136, -0.322561692126909, -0.319093286835222,
  -0.31441821757193, -0.312787095395791, -0.304006133793368, -0.289417944590358,
  -0.283638477186864, -0.273148888533504, -0.269428962035966, -0.265198174045359,
  -0.262437948152234, -0.258498917035502, -0.254570747217806, -0.246911498081545,
  -0.234405021134879, -0.229087757591321, -0.226531875386837, -0.203956533853441,
  -0.20265079420009, -0.200511322132944, -0.195576484733111, -0.189092318124576,
  -0.186287982898152, -0.183940138485593, -0.183864161559329, -0.173351884195258,
  -0.168372627947096, -0.167440878378378, -0.163389817897733, -0.15827027027027,
  -0.157639711424572, -0.157263782065054, -0.156517276770441, -0.156517276770441,
  -0.156449763926399, -0.155584392339359, -0.151310009818614, -0.151252169600793,
  -0.149030386318522, -0.149004876153113, -0.148193136957182, -0.148082213950477,
  -0.147996452050506, -0.146546546546547, -0.143490725539683, -0.141565536914374,
  -0.13960331299041, -0.137387387387387, -0.137047487863915, -0.133134480375396,
  -0.131891891891892, -0.128133314661812, -0.122547634742757, -0.116427881027247,
  -0.115021998742929, -0.113050193050193, -0.110937105329629, -0.110678510678511,
  -0.109706749263126, -0.105046639559914, -0.100467797463142, -0.100425298394842,
  -0.0995410504844467, -0.0974570629743044, -0.0940293906073374,
  -0.0891161431701972, -0.084186313973548, -0.0791351351351351,
  -0.0785070785070785, -0.0781967738489478, -0.075945427193028,
  -0.0742123805581136, -0.0742123805581136, -0.0704049244974512,
  -0.052992277992278, -0.052476880063087, -0.0440618792511888,
  -0.0432905553255663, -0.0416500711237553, -0.039449219907844,
  -0.0393315780989737, -0.0381435419995831, -0.036036036036036,
  -0.0309121621621622, -0.0291903855164645, -0.0275924460024878,
  -0.0154681655854447, 0, 0, 0, 0.0367046081331796, 0.0666120666120666,
  0.0746910194762956, 0.0873456237694648, 0.0958516656191075, 0.0960377853581737,
  0.101125288857631, 0.10250665172945, 0.103399566814201, 0.106136179097552,
  0.110318496935597, 0.116603833696172, 0.125531622993552, 0.147769914026505,
  0.158059524237953, 0.169736943807487, 0.173461879045635, 0.207160039620773,
  0.213713713713714, 0.222289705435773, 0.234474474474474, 0.250163757013072,
  0.285171658144631, 0.287848997247346, 0.324577702702703, 0.3377032793818,
  0.367338160441609, 0.378598468175015, 0.399096502842901, 0.415241281836799,
  0.437396580253723, 0.525422983959569, 0.554944846669952, 0.567365370357889,
  0.590560709963695, 0.594765651727677, 0.608324902442549, 0.616999293055631,
  0.689630807277866, 0.706113379153532, 0.719933907706833, 0.727344992050874,
  0.80750137892995, 0.820986978881716, 0.827458637344569, 0.961454912856782,
  0.980114059013142, 0.985539044579635, 1
), X = c(
  73376.3934991615,
  108989.76203822, 80956.8990471898, 76099.2985812133, 76900.623531451,
  98260.7092009951, 5157.57790114048, 30820.6795387321, 90177.3769949584,
  75777.9251254972, 54144.4989988793, 88978.3064788961, 27888.23288993,
  56461.8353689865, 110448.861738318, 112173.820609583, 113866.481090548,
  88193.5496923359, 86624.8100834148, 104085.434778549, 85741.3583741874,
  41186.6833252131, 65284.5646726949, 85149.9525214302, 85559.9150389397,
  61195.0493263041, 73264.2724621648, 32068.2542638019, 50869.4024020612,
  107091.769124836, 82257.9544395188, 9849.60340487854, 99159.7652222794,
  32303.6208808293, 69944.3311966181, 90770.2172968661, 12471.3741701687,
  95130.3127734323, 90038.1538169062, 90718.7440528133, 71250.124670451,
  102463.634147458, 89278.7280055428, 92136.0067933825, 73703.8477614749,
  63523.8780758179, 74230.2280241016, 101628.622525462, 61167.5732580155,
  60638.6442549899, 90108.1047746354, 73212.7788821735, 9096.17727006117,
  102435.432907824, 106029.219757365, 85754.9862715359, 6477.95629541751,
  91342.5519621444, 98035.8912685832, 104090.254520417, 66597.7900871458,
  22715.0411128533, 52856.0764092185, 68351.0264782874, 94902.9575756876,
  37534.2825114403, 140021.073330153, 92888.6648850076, 73356.3839464727,
  17218.2818121902, 90718.7440528133, 66146.0292524013, 98169.2858363378,
  101260.526266798, 103114.827080279, 107334.571092804, 50497.6004070263,
  26118.1914947915, 75551.3027817467, 112555.168771616, 95222.2180390191,
  95222.2180390191, 99970.5161385444, 60124.3745720035, 84533.5577261424,
  88797.6323826909, 85979.6084009005, 40395.2474581588, 29622.485242982,
  34021.650581241, 104188.944439581, 97363.0590907177, 41991.8378331259,
  90983.4443890642, 37402.2319634528, 68992.367022227, 20453.7899292494,
  96289.402156546, 62203.7060150694, 25262.275303754, 43726.5886744504,
  17130.9939045931, 22884.9895481025, 7145.63899258642, 72058.1892781364,
  3227.33385134628, 51966.2649320912, 5360.16800260167, 6214.23467758991,
  54768.4704168219, 8099.89207587154, 17175.2636997988, 6198.3000754135,
  28297.3261202508, 45073.8418558188, 6124.85470448711, 4732.46793860604,
  85195.2162860544, 10431.1362666056, 44405.1338499369, 29838.2693636617,
  63630.8310422971, 5499.65480547125, 4992.30879089013, 8817.58452161529,
  5801.14057475208, 5921.92115270202, 8528.42240081873, 7063.86021329663,
  10278.4337738267, 6029.81388864381, 7040.69161676212, 57556.0389612851,
  6753.0933741021, 5665.54685867964, 5317.56598804313, 6897.71522436381,
  7805.38044965523, 6174.40416180664, 97386.2972828469, 6154.75085293128,
  5413.95743794229, 8126.04433315438, 5774.28385069807, 6322.22812627604,
  3479.24351961031, 4991.36126383291, 5385.39494240332
)), row.names = c(
  NA,
  -148L
), class = c("tbl_df", "tbl", "data.frame"))

For this problem I will use the following packages

# For data wrangling and graphs
library(tidyverse)
# For makin NLS models easier
library(minpack.lm)
# For doing some test with negative binomial distributions
library(MASS)
#> 
#> Attaching package: 'MASS'
#> The following object is masked from 'package:dplyr':
#> 
#>     select

I will show a lot of what I have done, but I am looking for the best way on how to model this data

First lets see a the structure of the dataset:

As seen in the following graph

the relationship is clearly exponential (at least to my eyes), however there are some complications. Namely the response variable has positive and negative values (which complicates log transformations and exponentials), and also the values are not intergers, which complicates poisson and similar models.

If we look at the distribution of Y, we can see that is more or less normally distributed

This is what I have tried

So I have tried 3 approximations, a NLS using the minpack.lm, a GLM using quasipoisson family and a negative binomial, here are my issues and ideas

NLS

First thing I just tried to use an NLS, which seems reasonable to me, for that I used the minpack.lm package, here is the fitted model:

FitNLS <- nlsLM(Y ~ a * X^b, data = Example_Data, start = list(a = 0.3, b = -2.32999e-05))

Which gives me the following model:

summary(FitNLS)
#> 
#> Formula: Y ~ a * X^b
#> 
#> Parameters:
#>     Estimate Std. Error t value Pr(>|t|)    
#> a  7.292e+04  2.027e+05   0.360     0.72    
#> b -1.402e+00  3.242e-01  -4.324 2.82e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.2654 on 146 degrees of freedom
#> 
#> Number of iterations to convergence: 42 
#> Achieved convergence tolerance: 1.49e-08

And the following graph of data vs prediction

This approach has some clear issues, as for example that this model will never make a prediction bellow zero, and in this case this might be quite important

GLM quasipoisson

In order to be able to make this model, I need to transform the Y variable so that is only has positive values, so in order to do that I tried this:

MinY <- abs(min(Example_Data$Y)) + (abs(min(Example_Data$Y)) * 0.01)

Example_Data <- Example_Data %>%
  mutate(Positive_Y = Y + MinY)

This will set my values to go all positive but at the same time some of the values will exceed 1, with this I can try the following:

FitGLM <- glm(Positive_Y ~ I(log(X)), data = Example_Data, family = quasipoisson)

This gives me the following model

summary(FitGLM)
#> 
#> Call:
#> glm(formula = Positive_Y ~ I(log(X)), family = quasipoisson, 
#>     data = Example_Data)
#> 
#> Deviance Residuals: 
#>      Min        1Q    Median        3Q       Max  
#> -0.99025  -0.15909  -0.00601   0.13789   1.03932  
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  2.48135    0.23382   10.61   <2e-16 ***
#> I(log(X))   -0.27896    0.02293  -12.16   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for quasipoisson family taken to be 0.07508724)
#> 
#>     Null deviance: 21.923  on 147  degrees of freedom
#> Residual deviance: 11.213  on 146  degrees of freedom
#> AIC: NA
#> 
#> Number of Fisher Scoring iterations: 4

And backtransforming everything I can get the following reponse:

Although this looks quite good and I can get negative values, I am not quite convinced, I choose quasipoisson family because of the fractional data, but this might not be the best option, I guess I could transform this values from 0 to 1 and go quasibinomial with a log link, but that does not sound right either. I have tried also Negative Binomial GLM, but I think I am making my point clear. Any help on the best way of modelling this dataset would be greatly apprecieated

Session info

Standard output and standard error

-- nothing to show --

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.1.2 (2021-11-01)
#>  os       Ubuntu 18.04.6 LTS          
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       Europe/Copenhagen           
#>  date     2021-11-08                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version date       lib source        
#>  assertthat    0.2.1   2019-03-21 [1] CRAN (R 4.1.0)
#>  backports     1.2.1   2020-12-09 [1] CRAN (R 4.1.0)
#>  broom         0.7.9   2021-07-27 [1] CRAN (R 4.1.1)
#>  cellranger    1.1.0   2016-07-27 [1] CRAN (R 4.1.0)
#>  cli           3.0.1   2021-07-17 [1] CRAN (R 4.1.0)
#>  codetools     0.2-18  2020-11-04 [1] CRAN (R 4.1.0)
#>  colorspace    2.0-2   2021-06-24 [1] CRAN (R 4.1.0)
#>  crayon        1.4.1   2021-02-08 [1] CRAN (R 4.1.0)
#>  curl          4.3.2   2021-06-23 [1] CRAN (R 4.1.0)
#>  DBI           1.1.1   2021-01-15 [1] CRAN (R 4.1.0)
#>  dbplyr        2.1.1   2021-04-06 [1] CRAN (R 4.1.0)
#>  digest        0.6.28  2021-09-23 [1] CRAN (R 4.1.1)
#>  dplyr       * 1.0.7   2021-06-18 [1] CRAN (R 4.1.0)
#>  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.1.0)
#>  evaluate      0.14    2019-05-28 [1] CRAN (R 4.1.0)
#>  fansi         0.5.0   2021-05-25 [1] CRAN (R 4.1.0)
#>  farver        2.1.0   2021-02-28 [1] CRAN (R 4.1.0)
#>  fastmap       1.1.0   2021-01-25 [1] CRAN (R 4.1.0)
#>  forcats     * 0.5.1   2021-01-27 [1] CRAN (R 4.1.0)
#>  fs            1.5.0   2020-07-31 [1] CRAN (R 4.1.0)
#>  generics      0.1.0   2020-10-31 [1] CRAN (R 4.1.0)
#>  ggplot2     * 3.3.5   2021-06-25 [1] CRAN (R 4.1.0)
#>  glue          1.4.2   2020-08-27 [1] CRAN (R 4.1.0)
#>  gtable        0.3.0   2019-03-25 [1] CRAN (R 4.1.0)
#>  haven         2.4.3   2021-08-04 [1] CRAN (R 4.1.1)
#>  highr         0.9     2021-04-16 [1] CRAN (R 4.1.0)
#>  hms           1.1.1   2021-09-26 [1] CRAN (R 4.1.1)
#>  htmltools     0.5.2   2021-08-25 [1] CRAN (R 4.1.1)
#>  httr          1.4.2   2020-07-20 [1] CRAN (R 4.1.0)
#>  jsonlite      1.7.2   2020-12-09 [1] CRAN (R 4.1.0)
#>  knitr         1.36    2021-09-29 [1] CRAN (R 4.1.1)
#>  labeling      0.4.2   2020-10-20 [1] CRAN (R 4.1.0)
#>  lifecycle     1.0.1   2021-09-24 [1] CRAN (R 4.1.1)
#>  lubridate     1.7.10  2021-02-26 [1] CRAN (R 4.1.0)
#>  magrittr      2.0.1   2020-11-17 [1] CRAN (R 4.1.0)
#>  MASS        * 7.3-54  2021-05-03 [1] CRAN (R 4.1.0)
#>  mime          0.12    2021-09-28 [1] CRAN (R 4.1.1)
#>  minpack.lm  * 1.2-1   2016-11-20 [1] CRAN (R 4.1.0)
#>  modelr        0.1.8   2020-05-19 [1] CRAN (R 4.1.0)
#>  munsell       0.5.0   2018-06-12 [1] CRAN (R 4.1.0)
#>  pillar        1.6.3   2021-09-26 [1] CRAN (R 4.1.1)
#>  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.1.0)
#>  purrr       * 0.3.4   2020-04-17 [1] CRAN (R 4.1.0)
#>  R.cache       0.15.0  2021-04-30 [1] CRAN (R 4.1.1)
#>  R.methodsS3   1.8.1   2020-08-26 [1] CRAN (R 4.1.1)
#>  R.oo          1.24.0  2020-08-26 [1] CRAN (R 4.1.1)
#>  R.utils       2.11.0  2021-09-26 [1] CRAN (R 4.1.1)
#>  R6            2.5.1   2021-08-19 [1] CRAN (R 4.1.1)
#>  Rcpp          1.0.7   2021-07-07 [1] CRAN (R 4.1.0)
#>  readr       * 2.0.1   2021-08-10 [1] CRAN (R 4.1.1)
#>  readxl        1.3.1   2019-03-13 [1] CRAN (R 4.1.0)
#>  rematch2      2.1.2   2020-05-01 [1] CRAN (R 4.1.0)
#>  reprex        2.0.0   2021-04-02 [1] CRAN (R 4.1.0)
#>  rlang         0.4.11  2021-04-30 [1] CRAN (R 4.1.0)
#>  rmarkdown     2.11    2021-09-14 [1] CRAN (R 4.1.1)
#>  rstudioapi    0.13    2020-11-12 [1] CRAN (R 4.1.0)
#>  rvest         1.0.0   2021-03-09 [1] CRAN (R 4.1.0)
#>  scales        1.1.1   2020-05-11 [1] CRAN (R 4.1.0)
#>  sessioninfo   1.1.1   2018-11-05 [1] CRAN (R 4.1.0)
#>  stringi       1.7.5   2021-10-04 [1] CRAN (R 4.1.1)
#>  stringr     * 1.4.0   2019-02-10 [1] CRAN (R 4.1.0)
#>  styler        1.6.2   2021-09-23 [1] CRAN (R 4.1.2)
#>  tibble      * 3.1.5   2021-09-30 [1] CRAN (R 4.1.1)
#>  tidyr       * 1.1.4   2021-09-27 [1] CRAN (R 4.1.1)
#>  tidyselect    1.1.1   2021-04-30 [1] CRAN (R 4.1.0)
#>  tidyverse   * 1.3.1   2021-04-15 [1] CRAN (R 4.1.0)
#>  tzdb          0.1.2   2021-07-20 [1] CRAN (R 4.1.1)
#>  utf8          1.2.2   2021-07-24 [1] CRAN (R 4.1.1)
#>  vctrs         0.3.8   2021-04-29 [1] CRAN (R 4.1.0)
#>  withr         2.4.2   2021-04-18 [1] CRAN (R 4.1.0)
#>  xfun          0.26    2021-09-14 [1] CRAN (R 4.1.1)
#>  xml2          1.3.2   2020-04-23 [1] CRAN (R 4.1.0)
#>  yaml          2.2.1   2020-02-01 [1] CRAN (R 4.1.0)
#> 
#> [1] /home/derek/R/x86_64-pc-linux-gnu-library/4.1
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library

Answer 1

I would not try and fit a line through the data points.

The image below, with the x-axis on a log scale, shows that the variation in the y-axis does not strongly relate to the position on the x-axis.

Instead of some smooth curve with 'y' some function of 'x' the distribution seems to be more like a sort of step function. You have very different behavior of the points below and above 15 km.

You might plot some sort of sinusoid curve through the points. But it won't be well determined with a regression algorithm.

• In the image below I have manually plotted a curve. In the region below 15 km (the blue points) you can see that the points mostly differentiate in the horizontal direction. In the vertical direction, there are large variations.
• If you would fit a curve through these points, then the optimum, that minimizes the vertical distance, would be a curve through the middle, around the value of 0.5. Probably not what you had expected (you might be thinking of a curve that goes to 1).

In my opinion, the best would be to just make the plot and potentially describe some clusters. But fitting a curve may have little value. Especially since, as you mentioned in the comments, you have no theory behind the curve. So what would be the point of putting a curve into the diagram?