draft si_gravity.R

[R-icntay]

Hello @Robinlovelace , @adamdennett, @Nowosad,

I've had a look at #26 and here's a draft version of how I approached it. I was hoping to get your thoughts on the following:

• The function's dependency on dplyr
• Whether a poisson model should be included to calibrate $\mu$, $\alpha$, $\beta$. Probably do this when a user doesn't provide values for these arguments?

Thank you!

Here is the function in action:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
source("https://raw.githubusercontent.com/R-icntay/simodels/main/R/si_gravity.R")

# Import data
mdatasub <- read.csv("https://raw.githubusercontent.com/R-icntay/sim_australia/main/mdatasub.csv")

# Run a classic gravity model
mdatasub = si_gravity(
  O_i = vi1_origpop,
  D_j = wj3_destmedinc, 
  C_ij = dist,
  flows = Flow,
  data = mdatasub,
  beta = -2
) %>% 
  mutate(T_ij = round(T_ij))

# View estimated flows T_ij
mdatasub %>% 
  glimpse()
#> Rows: 210
#> Columns: 18
#> $ Origin            <chr> "Greater Sydney", "Greater Sydney", "Greater Sydney"~
#> $ Orig_code         <chr> "1GSYD", "1GSYD", "1GSYD", "1GSYD", "1GSYD", "1GSYD"~
#> $ Destination       <chr> "Rest of NSW", "Greater Melbourne", "Rest of Vic", "~
#> $ Dest_code         <chr> "1RNSW", "2GMEL", "2RVIC", "3GBRI", "3RQLD", "4GADE"~
#> $ Flow              <int> 91031, 22601, 4416, 22888, 27445, 5817, 795, 10574, ~
#> $ vi1_origpop       <int> 4391673, 4391673, 4391673, 4391673, 4391673, 4391673~
#> $ wj1_destpop       <int> 2512952, 3999981, 1345717, 2065998, 2253723, 1225235~
#> $ vi2_origunemp     <dbl> 5.74, 5.74, 5.74, 5.74, 5.74, 5.74, 5.74, 5.74, 5.74~
#> $ wj2_destunemp     <dbl> 6.12, 5.47, 5.17, 5.86, 6.22, 5.78, 5.45, 4.76, 4.42~
#> $ vi3_origmedinc    <dbl> 780.64, 780.64, 780.64, 780.64, 780.64, 780.64, 780.~
#> $ wj3_destmedinc    <dbl> 509.97, 407.95, 506.58, 767.08, 446.48, 445.53, 522.~
#> $ vi4_origpctrent   <dbl> 31.77, 31.77, 31.77, 31.77, 31.77, 31.77, 31.77, 31.~
#> $ wj4_destpctrent   <dbl> 27.20, 27.34, 24.08, 33.19, 32.57, 28.27, 26.17, 27.~
#> $ flow_no_intra     <int> 91031, 22601, 4416, 22888, 27445, 5817, 795, 10574, ~
#> $ offset            <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
#> $ dist              <dbl> 391.4379, 682.7450, 685.8484, 707.9081, 1386.4854, 1~
#> $ unconstrainedEst1 <int> 45500, 11964, 14723, 20926, 3175, 4923, 3342, 960, 1~
#> $ T_ij              <dbl> 45500, 11964, 14723, 20926, 3175, 4923, 3342, 960, 1~

^{Created on 2022-08-13 by the reprex package (v2.0.1)}

[Robinlovelace]

Hi @R-icntay Many thanks for this, I'm on holiday and flying to Italy tomorrow, have added to the TODO list (on awesome app LogSeq ; ) and will get back with feedback asap.

[Robinlovelace]

Looking like a solid start, many thanks! I have made a few comments and suggested adding some Roxygen tags in there. Keen to see examples of this before merging.

[Robinlovelace]

Two minor issues with this:

1. Use of both equals and arrow assignment in same script file, please use = throughout as that's the package style
2. Use of pipe operator in function calls, not sure about this but if it works I'm happy to give it a spin

[Robinlovelace]

Can you try adding some examples? You can use the "Insert Roxygen Skeleton" button, in the Code tab if I remember, to do this I think.

[R-icntay]

Looking like a solid start, many thanks! I have made a few comments and suggested adding some Roxygen tags in there. Keen to see examples of this before merging.

Thank you Robin! So my main thoughts on this is that the function produces results from uncalibrated values of $\alpha \ , \beta$ etc. How useful would this be in a real world scenario? Tagging in @adamdennett. Should we also include a calibration model inside the function?

Thank you.