4
$\begingroup$

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
$\endgroup$
2
  • 2
    $\begingroup$ Your responses look like they have approximately Normal errors. The sharp drop in mean response from $x=$ to $x\approx 10000$ is much steeper than a mere exponential. The very best approach to take would rely on an underlying scientific theory to suggest the proper shape of the curve. Barring that, you need to bring in your statistical objectives in the form of stipulating what kinds of estimation errors you can tolerate and how large they can be. Until you provide some of that information, our ability to help you will be seriously hampered. $\endgroup$
    – whuber
    Commented Nov 8, 2021 at 21:17
  • $\begingroup$ Thanks for your comments @whuber, to be honest, it is quite hard to have a scientific theory about this. Y is a made-up index that is supposed to measure for a city it's cultural closeness to the Sea (Where negative values actually mean negative closeness). and X is the distance to the sea in meters. The aim is to explore the residuals and check what happens culturally in those places that could explain the difference. I would say that we are more interested in having lower errors closer to the sea (low x values) than far away. $\endgroup$ Commented Nov 9, 2021 at 13:54

1 Answer 1

1
$\begingroup$

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?

example of two clusters

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.