I have the following data set containing socioeconomic variables:

    > glimpse(df)
    Rows: 730,099
    Columns: 9
    $ id     <int> 25500, 25501, 25502, 25503, 25504, 25505, 25506, 25507, 25508, 25509, 25510, 25511, 25512, 25513, 25514, 25515, 25516, 255…
    $ Prov   <fct> Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al Hoceïma, Al…
    $ Age    <dbl> 45, 15, 65, 55, 35, 75, 45, 25, 55, 40, 35, 50, 70, 20, 30, 75, 50, 35, 50, 40, 35, 70, 70, 35, 35, 75, 55, 35, 25, 50, 30…
    $ Edu    <dbl> 5, 4, 0, 0, 0, 0, 0, 16, 0, 5, 4, 0, 0, 6, 0, 0, 0, 0, 0, 0, 9, 0, 0, 14, 4, 0, 0, 3, 5, 0, 9, 0, 0, 0, 15, 0, 0, 0, 14, 3…
    $ Kids   <int> 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 3, 4, 1, 0, 2, 0, 3, 2, 2, 1, 0, 5, 0, 2, 0, 0, 0, 1, 0, 2, 2, 1, 0, 5, 0, 0, 3, 0, 0, 1, 0,…
    $ Mil    <fct> Urban, Rural, Rural, Rural, Rural, Rural, Rural, Rural, Rural, Rural, Rural, Rural, Urban, Urban, Rural, Urban, Urban, Rur…
    $ DSize  <dbl+lbl> 2, 1, 4, 4, 1, 1, 3, 2, 3, 4, 1, 6, 3, 2, 5, 3, 4, 2, 4, 2, 3, 6, 3, 2, 4, 4, 3, 2, 2, 4, 1, 4, 6, 4, 5, 5, 4, 2, 5, 3…
    $ taille <dbl> 2, 5, 5, 7, 5, 1, 2, 1, 4, 1, 5, 9, 5, 1, 5, 3, 7, 4, 6, 7, 2, 14, 3, 4, 1, 2, 9, 3, 3, 12, 4, 7, 2, 9, 2, 6, 5, 2, 4, 3, …
    $ EP     <fct> 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,…

The variable `Prov` indicates the 72 provinces of the country, it looks like this:

    > levels(df$Prov)
     [1] "Agadir-Ida-Ou-Tanane"     "Al Haouz"                 "Al Hoceïma"               "Meknès"                   "Azilal"                  
     [6] "Béni Mellal"              "Benslimane"               "Berkane"                  "Berrechid"                "Boujdour"                
    [11] "Boulemane"                "Casablanca"               "Chefchaouen"              "Chichaoua"                "Chtouka-Ait Baha"        
    [16] "Driouch"                  "El Hajeb"                 "El Jadida"                "El Kelâa des Sraghna"     "Errachidia"              
    [21] "Essaouira"                "Fahs-Anjra"               "Fès"                      "Figuig"                   "Fquih Ben Salah"         
    [26] "Guelmim"                  "Guercif"                  "Ifrane"                   "Inezgane-Ait Melloul"     "Jerada"                  
    [31] "Kénitra"                  "Khémisset"                "Khénifra"                 "Khouribga"                "Laâyoune"                
    [36] "Larache"                  "Marrakech"                "Médiouna"                 "Midelt"                   "Mohammadia"              
    [41] "Nador"                    "Nouaceur"                 "Ouarzazate"               "Ouezzane"                 "Oujda-Angad"             
    [46] "Rabat"                    "Rehamna"                  "Safi"                     "Salé"                     "Sefrou"                  
    [51] "Settat"                   "Sidi Bennour"             "Sidi Ifni"                "Sidi Kacem"               "Sidi Slimane"            
    [56] "Skhirate-Témara"          "Tanger-Assilah"           "Taounate"                 "Taourirt"                 "Taroudannt"              
    [61] "Tata"                     "Taza"                     "Tétouan"                  "M'Diq-Fnideq"             "Tinghir"                 
    [66] "Tiznit"                   "Youssoufia"               "Zagora"                   "Moulay Yacoub"            "Tan-Tan / Assa-Zag"      
    [71] "Es-Semara / Tarfaya"      "Oued Ed Dahab / Aousserd"

and I have another shapefile containing the location polygons for each province, it looks like this:

    > map_3
    Simple feature collection with 72 features and 1 field
    Geometry type: POLYGON
    Dimension:     XY
    Bounding box:  xmin: -17.10496 ymin: 20.7715 xmax: -0.9987581 ymax: 35.92243
    Geodetic CRS:  WGS 84
    # A tibble: 72 × 2
       PROV                                                                                                     geometry
     * <chr>                                                                                               <POLYGON [°]>
     1 Agadir-Ida-Ou-Tanane ((-9.504189 30.3498, -9.503555 30.34989, -9.503493 30.34991, -9.502746 30.35016, -9.50251...
     2 Al Haouz             ((-7.344661 31.67321, -7.344881 31.67344, -7.345 31.6736, -7.346125 31.67503, -7.346162 3...
     3 Al Hoceïma           ((-3.926041 35.26096, -3.926289 35.26113, -3.926532 35.26142, -3.926776 35.26171, -3.9270...
     4 Azilal               ((-5.908287 32.30598, -5.908415 32.3063, -5.913961 32.32005, -5.915668 32.32325, -5.92306...
     5 Benslimane           ((-7.195987 33.18768, -7.191716 33.18865, -7.169873 33.19075, -7.164904 33.19066, -7.1557...
     6 Berkane              ((-2.313083 34.58991, -2.312776 34.58996, -2.312381 34.59001, -2.308861 34.59052, -2.3085...
     7 Berrechid            ((-7.250729 33.21856, -7.249378 33.22028, -7.246216 33.22429, -7.245882 33.22472, -7.2454...
     8 Boujdour             ((-12.54495 24.46466, -12.49897 24.47041, -12.4321 24.48361, -12.38604 24.49452, -12.3777...
     9 Boulemane            ((-4.202123 32.58081, -4.200447 32.58095, -4.183846 32.58234, -4.17291 32.58488, -4.15931...
    10 Casablanca           ((-7.669198 33.49706, -7.668895 33.49699, -7.66875 33.49705, -7.668693 33.49707, -7.66854...
    # ℹ 62 more rows
    # ℹ Use `print(n = ...)` to see more rows

When I want to run a simple logistic regression I simply run this command and it gives me the results I want:

    > model <- glm(EP~Age+Edu+Kids+Mil+DSize+taille,family = binomial(link = "logit"), data = df)
    > summary(model)
    
    Call:
    glm(formula = EP ~ Age + Edu + Kids + Mil + DSize + taille, family = binomial(link = "logit"), 
        data = df)
    
    Coefficients:
                  Estimate Std. Error  z value Pr(>|z|)    
    (Intercept) -0.0040682  0.0165088   -0.246    0.805    
    Age         -0.0060011  0.0002926  -20.511   <2e-16 ***
    Edu         -0.1165767  0.0010738 -108.562   <2e-16 ***
    Kids         0.0927567  0.0038183   24.293   <2e-16 ***
    MilRural     1.5594703  0.0078742  198.048   <2e-16 ***
    DSize       -0.4697581  0.0034187 -137.409   <2e-16 ***
    taille      -0.1368148  0.0026621  -51.393   <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    (Dispersion parameter for binomial family taken to be 1)
    
        Null deviance: 613335  on 719910  degrees of freedom
    Residual deviance: 497842  on 719904  degrees of freedom
      (10188 observations deleted due to missingness)
    AIC: 497856
    
    Number of Fisher Scoring iterations: 6

Now I want to run a Geographically Weighted Logistic Regression, and for that I checked the `GWModel` package manual, and found the function `ggwr.basic`, which has an option for logistic regression. 

    DM<-gw.dist(dp.locat=coordinates(londonhp))
    bw.f2 <- bw.ggwr(BATH2~FLOORSZ,data=londonhp, dMat=DM,family ="binomial")
    res.binomial<-ggwr.basic(BATH2~FLOORSZ, bw=bw.f2,data=londonhp, dMat=DM,
                  family ="binomial")

But in their example, the location variable is a point (X Y coordinates), whereas in my case the locations are polygons, also they have one observation per location (they use the `LondonHP` data set), whereas I have more than 730,000 observations scattered across 72 provinces. 

Therefore, how can I run a Geographically Weighted Logistic Regression on my data? I just gave the `GWModel` package as an example; do you know a different package or function that can get the job done?