This the the R code to replicate the study on “The role of urban greening in limiting land take in case of population growth: evidence for Europe 2005-2018”, published in Land Use Policy (Vol.18 November 2025, 107726), available at: https://www.sciencedirect.com/science/article/pii/S0264837725002601.

This project has received funding from the European Union’s Horizon Europe Research and Innovation programme and the State Secretariat for Education, Research and Innovation from Switzerland under Grant Agreement No 101060423 as a part of ‘LAnd use and MAnagement modelling for SUStainable governance’ (LAMASUS) project.

Data

This paper uses merged dataset on grid level (2 km x 2 km) for almost whole Europe that includes land use (state and changes) and population (state and changes). Dataset can be downloaded from here: https://www.sugarsync.com/pf/D1836703_00572095_821512

Original data (before merging) are available here:

You can read the data using the following code. Please remember to set your working directory.

# read packages
library(sf)
library(ggplot2)

# set Working Directory 
setwd("C:/Users/Admin/Desktop/my_folder")

# read data into data.frame and sf class
dane<-read.csv("LUCAS_merged.csv", header=TRUE, sep=",", dec=".")
dane.sf<-st_as_sf(dane, coords=c("X_WGS84", "Y_WGS84"), crs=4326)

Overview of data

Codes for coloured matrix that illustrates the structure of land-use changes.

library(plot.matrix)        # for plot.matrix()

a1<-table(dane$change)
vec.dane<-paste0(dane$STR05, "to", dane$STR18a)
a2<-table(vec.dane)
a3<-matrix(a2, nrow=7, ncol=7, byrow=TRUE)      # matrix of nominal changes between clases
rownames(a3)<-paste0(rep("STR05_", times=7), 1:7)   
colnames(a3)<-paste0(rep("STR18_", times=7), 1:7)
a4<-a3/sum(a3)                      # matrix of relative changes between classes
str0<-a4

# Figure 1a - relative changes between classes 
par(mar=c(4,4,4,4))
plot(a4, breaks=c(0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.6), digits=3, col=brewer.pal(7, "BuPu"), cex=1.2, cex.axis=0.7, main="Structure of land use changes (2005-2018)", xlab="Land use structure in 2018", ylab="Land use structure in 2005") 

# Figure 1b - nominal changes between classes 
par(mar=c(4,4,4,4))
plot(a3, breaks=c(0, 1000, 5000, 10000, 50000, 100000, 250000, 500000), digits=0, col=brewer.pal(7, "BuPu"), cex=1.1, cex.axis=0.7, main="Structure of land use changes (2005-2018)", xlab="Land use structure in 2018", ylab="Land use structure in 2005")

Figure 1. Structure of land use changes in 2005-2018: a) in percent, b) in nominal values

Codes for circular flows of land-use classes

library(circlize)
library(reshape)
# we use a3 object - created for Fig.1 - full matrix of flows between land use classes
colnames(a3)<-c("18_1", "18_2", "18_3", "18_4", "18_5", "18_6", "18_7")
rownames(a3)<-c("05_1", "05_2", "05_3", "05_4", "05_5", "05_6", "05_7")
chordDiagram(a3)

Figure 2: Structure of flows of land use between years 2005-2018

Codes for point maps of changes from/to class 6 (artificial land cover). This code presents the visualisation of locations that undergo the changes from/to class 6. In grey are all available points in the dataset. In colours are the points of interest. First two figures are unconditional of population, another four figures refer to two factors: land use change and population change

# Changes from / to class 6 (indepedent of population changes)

# subsets – changes in land use - when class 6 disappears
sub.disappeared6<-dane[dane$change=="6to1" | dane$change=="6to2" | dane$change=="6to3" | dane$change=="6to4" | dane$change=="6to5" | dane$change=="6to7", ]

# subsets – changes in land use - when class 6 appears
sub.appeared6<-dane[dane$change=="1to6" | dane$change=="2to6" | dane$change=="3to6" | dane$change=="4to6" | dane$change=="5to6" | dane$change=="7to6", ]

# Plotting
# selecting a subset
selector<-sample(1:dim(dane)[1], 10000, replace=FALSE)

# Fig.3a
nn<-dim(sub.disappeared6)[1]
share<-round(nn/dim(dane)[1],4)*100
plot(dane[selector, 4:5], col="grey80", main="Class 6 (artificial land) disappeared")
points(sub.disappeared6[, 4:5], bg="red", pch=21, cex=0.01)
text(0, 69, paste0("location of ",nn, " points 
(",share,"% of data)"))

# Fig.3b
nn<-dim(sub.appeared6)[1]
share<-round(nn/dim(dane)[1],4)*100
plot(dane[selector, 4:5], col="grey80", main="Class 6 (artificial land) appeared")
points(sub.appeared6[, 4:5], bg="red", pch=21, cex=0.01)
text(0, 69, paste0("location of ",nn, " points 
(",share,"% of data)"))

Figure 3 (a & b): Location of land use changes of artificial land (class 6) between 2005 and 2018 

# conditional changes (on population)

# artificial cover disappears & population increases
cond1<-sub.disappeared6$popul.change.lab=="h.pop_growth_1.2-5" | sub.disappeared6$popul.change.lab=="g.pop_growth_from_1.05_to_1.2" | sub.disappeared6$popul.change.lab=="f.pop_growth_up_to_1.05" | sub.disappeared6$popul.change.lab=="appeared_new_inhabitants" | sub.disappeared6$popul.change.lab=="i.pop_growth_5_to_100_times" | sub.disappeared6$popul.change.lab=="j.pop_growth_100_times_or_more"

# artificial cover disappears & population decreases
cond2<-sub.disappeared6$popul.change.lab=="d.depop_between_5%_&_20%" | sub.disappeared6$popul.change.lab=="always_no_inhabitants" | sub.disappeared6$popul.change.lab=="c.depop_between_20%_&_50%" | sub.disappeared6$popul.change.lab=="e.depop_not_more_than_5%" | sub.disappeared6$popul.change.lab=="b.depop_more_than_50%" | sub.disappeared6$popul.change.lab=="all_inhabitants_disappeared"

# plotting of conditional changes
# Fig.3c
selector<-sample(1:dim(dane)[1], 10000, replace=FALSE)
nn<-dim(sub.disappeared6[cond1,])[1]
share<-round(nn/dim(dane)[1],4)*100
plot(dane[selector, 4:5], col="grey80", main="Class 6 (artificial land) disappeared / population inflow")
points(sub.disappeared6[cond1, 4:5], bg="red", pch=21, cex=0.01)
text(0, 69, paste0("location of ",nn, " points 
(",share,"% of data)"))

# Fig.3d
nn<-dim(sub.disappeared6[cond2,])[1]
share<-round(nn/dim(dane)[1],4)*100
plot(dane[selector, 4:5], col="grey80", main="Class 6 (artificial land) disappeared / population outflow")
points(sub.disappeared6[cond2, 4:5], bg="red", pch=21, cex=0.01)
text(0, 69, paste0("location of ",nn, " points 
(",share,"% of data)"))

Figure 3 (c & d): Location of land use changes of artificial land (class 6) between 2005 and 2018 

# conditional changes

# artificial cover appears & population increases
cond3<-sub.appeared6$popul.change.lab=="h.pop_growth_1.2-5" | sub.appeared6$popul.change.lab=="g.pop_growth_from_1.05_to_1.2" | sub.appeared6$popul.change.lab=="f.pop_growth_up_to_1.05" | sub.appeared6$popul.change.lab=="appeared_new_inhabitants" | sub.appeared6$popul.change.lab=="i.pop_growth_5_to_100_times" | sub.appeared6$popul.change.lab=="j.pop_growth_100_times_or_more"

# artificial cover appears & population decreases
cond4<-sub.appeared6$popul.change.lab=="d.depop_between_5%_&_20%" | sub.appeared6$popul.change.lab=="always_no_inhabitants" | sub.appeared6$popul.change.lab=="c.depop_between_20%_&_50%" | sub.appeared6$popul.change.lab=="e.depop_not_more_than_5%" | sub.appeared6$popul.change.lab=="b.depop_more_than_50%" | sub.appeared6$popul.change.lab=="all_inhabitants_disappeared"

# Fig.3e
nn<-dim(sub.appeared6[cond3,])[1]
share<-round(nn/dim(dane)[1],4)*100
plot(dane[selector, 4:5], col="grey80", main="Class 6 (artificial land) appeared / population inflow")
points(sub.appeared6[cond1, 4:5], bg="red", pch=21, cex=0.01)
text(0, 69, paste0("location of ",nn, " points 
(",share,"% of data)"))

# Fig.3f
nn<-dim(sub.appeared6[cond4,])[1]
share<-round(nn/dim(dane)[1],4)*100
plot(dane[selector, 4:5], col="grey80", main="Class 6 (artificial land) appeared / population outflow")
points(sub.appeared6[cond2, 4:5], bg="red", pch=21, cex=0.01)
text(0, 69, paste0("location of ",nn, " points 
(",share,"% of data)"))

Figure 3 (e & f): Location of land use changes of artificial land (class 6) between 2005 and 2018

Codes for interactive maps in mapview:: - it is to illustrate particular changes at grid level in the selected cities.

# necessary packages for interactive plotting and OpenStreetmap contours
library(sf)
library(mapview)
library(osmdata)
# Getting countour of Paris (France) from OpenStreetMap
paris_bbox <- getbb("Paris, France")  # Bounding box for Paris
paris_boundary <- opq(paris_bbox) %>%
  add_osm_feature(key = "boundary", value = "administrative") %>%
  osmdata_sf()
paris_sf <- paris_boundary$osm_multipolygons[1:4,]

# selection of data within the contour of Paris
data.paris.sf<-st_join(dane.sf, paris_sf, join=st_intersects)
uu<-which(is.na(data.paris.sf$osm_id)==FALSE)
uu   # IDs of observations within the contour of Paris
#  [1] 309914 314149 315700 320279 321102 323996 324926 325036 325924 326213
# [11] 326320 326550 326774 327162 327652 327975 328336 328486 328556 328696
# [21] 328751 329274 329678 330271 331638 332199 332661 333091 333158 333198
# [31] 333274 333340 333423 333535 333753 334289 334497 334796 334810 335357
# [41] 335733 336006 336092 336551 336856 337814 337864 337916 337922 338285
# [51] 338711 338935 338940 342160 343488 343726 343740 343805 343838 343905
# [61] 343936 344022 344045 344086 344089 344096 344102 344147 344187 344189
# [71] 344197 344207 344250 344295 344308 344313 344350 344371 344395 344403
# [81] 344509 345643 346470 346920 347332 347460 347974 348372 349169 350573
# [91] 351874 352416 352450 352673 353287 355822 355867 356686 356878 357751
#[101] 359741 361181 362300 363523 363902 365209 365701 366056 366070 366533
#[111] 366688 368551 370292 370600 371587 372253 372631 374371 375024 375105
#[121] 377237 377257 378344 378605 378606 378892 379474 380978 381872 382400
#[131] 385856 386948 387275 387718 388864 389257 390089 390997 391598 392151
#[141] 392808 392981 395224 395520 395706 396288 399346 399931 400305 401129
#[151] 401277 402011 402048 402786 402964 403278 403412 404781 405571 405801
#[161] 406306 406384 406446 406811 406964 407009 407454 407755 408712 408864
#[171] 408898 409268 410169 410567 411070 411715 412154 412348 412970 413319
#[181] 414314 414909 415858 416403 417616 417662 417815 418319 421875 422213
#[191] 868046 870153 878935 883047

paris.sf<-data.paris.sf[uu,] # subset of data by selecting proper rows

# creating 4 labels: greening, green, urbanisation, urban
paris.sf$greening<-ifelse(paris.sf$change=="6to1" |paris.sf$change=="6to2" |paris.sf$change=="6to3" | paris.sf$change=="6to4" | paris.sf$change=="6to5" | paris.sf$change=="6to7", "greening", ifelse(paris.sf$change=="1to6" |paris.sf$change=="2to6" |paris.sf$change=="3to6" | paris.sf$change=="4to6" | paris.sf$change=="5to6" | paris.sf$change=="7to6", "urbanisation", ifelse(paris.sf$change=="no_change_6", "urban", "green")))

# Plotting values of labels with mapview()
mapview(paris_sf, col.regions = "white", alpha.regions = 0.4, layer.name = "Paris") + mapview(paris.sf, zcol="greening", cex="greening")

Figure 4a: Land use changes in selected European cities - Paris 

# code for Warsaw, Poland
warsaw_bbox <- getbb("Warsaw, Poland")  # bounding box for Warsaw
warsaw_boundary <- opq(warsaw_bbox) %>%
  add_osm_feature(key = "boundary", value = "administrative") %>%
  osmdata_sf()
warsaw_sf <- warsaw_boundary$osm_multipolygons[3,]
mapview(warsaw_sf, col.regions = "white", alpha.regions = 0.4, layer.name = "Warsaw") 
# selection of observations within the contour
data.warsaw.sf<-st_join(dane.sf, warsaw_sf, join=st_intersects)
uu<-which(is.na(data.warsaw.sf$osm_id)==FALSE)
uu
warsaw.sf<-data.warsaw.sf[uu,]

# creating 4 labels: greening, green, urbanisation, urban
warsaw.sf$greening<-ifelse(warsaw.sf$change=="6to3" | warsaw.sf$change=="6to4" | warsaw.sf$change=="6to5", "greening", ifelse(warsaw.sf$change=="1to6" | warsaw.sf$change=="3to6" | warsaw.sf$change=="4to6", "urbanisation", ifelse(warsaw.sf$change=="no_change_6", "urban", "green")))

# Plotting values of labels with mapview()
mapview(warsaw_sf, col.regions = "white", alpha.regions = 0.4, layer.name = "Warsaw") + mapview(warsaw.sf, zcol="greening", cex="greening")
# code for Berlin, Germany
berlin_bbox <- getbb("Berlin, Germany")  # Bounding box for Berlin
berlin_boundary <- opq(berlin_bbox) %>%
  add_osm_feature(key = "boundary", value = "administrative") %>%
  osmdata_sf()
berlin_sf <- berlin_boundary$osm_multipolygons[c(57,60),]
mapview(berlin_sf, col.regions = "white", alpha.regions = 0.4, layer.name = "Berlin") 
# selection of observations within the contour
data.berlin.sf<-st_join(dane.sf, berlin_sf, join=st_intersects)
uu<-which(is.na(data.berlin.sf$osm_id)==FALSE)
uu
berlin.sf<-data.berlin.sf[uu,]

# creating 4 labels: greening, green, urbanisation, urban
berlin.sf$greening<-ifelse(berlin.sf$change=="6to1" |berlin.sf$change=="6to2" |berlin.sf$change=="6to3" | berlin.sf$change=="6to4" | berlin.sf$change=="6to5" | berlin.sf$change=="6to7", "greening", ifelse(berlin.sf$change=="1to6" |berlin.sf$change=="2to6" |berlin.sf$change=="3to6" | berlin.sf$change=="4to6" | berlin.sf$change=="5to6" | berlin.sf$change=="7to6", "urbanisation", ifelse(berlin.sf$change=="no_change_6", "urban", "green")))

# Rysowanie mapy z mapview
mapview(berlin_sf, col.regions = "white", alpha.regions = 0.4, layer.name = "Berlin") + mapview(berlin.sf, zcol="greening", cex="greening")

Figure 4 (b & c): Land use changes in selected European cities: b) Warsaw, c) Berlin

Association rules

The first step is preparing data for market basket analysis / association rules. Running market basket analysis requires selecting variables of interest in a separate object and reading this file in transaction format. Of course, all observations must be labels - with information on variable and its level.

# checking unique levels of variables that will be used in association rules

unique(dane$STR05_label)
[1] "STR05_1" "STR05_4" "STR05_3" "STR05_2" "STR05_6" "STR05_5" "STR05_7"

unique(dane$STR18a_label)
[1] "STR18a_2" "STR18a_1" "STR18a_3" "STR18a_5" "STR18a_4" "STR18a_6" "STR18a_7"

unique(dane$change)
 [1] "1to2"        "4to2"        "3to2"        "2to1"        "1to3"       
 [6] "no_change_1" "no_change_2" "1to5"        "1to4"        "3to4"       
[11] "no_change_6" "4to3"        "6to1"        "3to1"        "3to5"       
[16] "6to5"        "6to3"        "no_change_3" "3to6"        "6to4"       
[21] "no_change_5" "4to5"        "4to1"        "5to3"        "1to7"       
[26] "no_change_4" "5to4"        "1to6"        "7to5"        "5to1"       
[31] "7to3"        "4to6"        "6to2"        "5to6"        "no_change_7"
[36] "3to7"        "4to7"        "7to1"        "7to4"        "7to6"       
[41] "6to7"        "5to7"        "2to4"        "2to3"        "2to6"       
[46] "2to7"        "2to5"        "7to2"        "5to2"       

unique(dane$change_dummy)
[1] "change"    "no_change"

unique(dane$ELEV_class)
[1] "elev_0-250"     "elev_250-500"   "elev_<0"        "elev_500-1000" 
[5] "elev_1000-2000" "elev_>2000"    

unique(dane$popul.change.lab)
 [1] "g.pop_growth_from_1.05_to_1.2"  "h.pop_growth_1.2-5"            
 [3] "e.depop_not_more_than_5%"       "d.depop_between_5%_&_20%"      
 [5] "f.pop_growth_up_to_1.05"        "i.pop_growth_5_to_100_times"   
 [7] "c.depop_between_20%_&_50%"      "b.depop_more_than_50%"         
 [9] "j.pop_growth_100_times_or_more" "always_no_inhabitants"         
[11] "all_inhabitants_disappeared"    "appeared_new_inhabitants"      

unique(dane$dist_bord_label)
[1] "border_far_>25km"

unique(dane$dist_coast_label)
[1] "coast_far_>25km" 
# selected variables saved as a new file
write.csv(dane[, c("STR05_label", "change", "change_dummy", "ELEV_class", "popul.change.lab", "dist_bord_label", "dist_coast_label")], file="LUCAS_basket.csv")
# reading data as transactions - for association rules
library(arules)
library(arulesViz)

trans1<-read.transactions("LUCAS_basket.csv", format="basket", sep=",", skip=0) 
trans1
inspect(head(trans1))
sort(itemFrequency(trans1, type="relative"))

The second step is running the association rules - they are to be executed one after another. In the option of apriori() function you have to specify in each run the rhs option. In fact it should be each and every of 49 levels of ‘change’ variable.

# setting rules - change 'rhs' option
rules<-apriori(data=trans1, parameter=list(supp=0.001, conf=0.001), appearance=list(default="lhs", rhs="no_change_1"), control=list(verbose=F)) 

# cleaning the rules
rules.clean<-rules[!is.redundant(rules)]
rules.clean<-rules.clean[is.significant(rules.clean, trans1)]
rules.clean<-rules.clean[is.maximal(rules.clean)]

# displaying cleaned rules interactively
ruleExplorer(rules.clean)

Figure 1. Association rules in ruleExplorer in web browser

The next step is screening of rules related to class 6 - artificial land cover. The results of codes above can be summarised in as below - in Table 2.

Table 2: Association rules related to class 6 of land use

The last step is the visualisation of observations following given rules. This codes makes the subsets of data that make the rules from Tab.2. Those points are mapped in colours.

# rules reported in Table 2 were recoded into data subsets
# each condition 'cond' is one rule from Table 2 in the paper

cond1<-dane$popul.change.lab=="always_no_inhabitants" & dane$change=="6to4"
cond2<-dane$popul.change.lab=="c.depop_between_20%_&_50%" & dane$change=="6to4"
cond3<-dane$popul.change.lab=="h.pop_growth_1.2-5" & dane$change=="6to3"
cond4<-dane$popul.change.lab=="e.depop_not_more_than_5%" & dane$change=="no_change_6" & dane$ELEV_class=="elev_0-250"
cond5<-dane$popul.change.lab=="f.pop_growth_up_to_1.05" & dane$change=="no_change_6" & dane$ELEV_class=="elev_0-250"
cond6<-dane$popul.change.lab=="g.pop_growth_from_1.05_to_1.2" & dane$change=="no_change_6" & dane$ELEV_class=="elev_0-250"
cond7<-dane$popul.change.lab=="h.pop_growth_1.2-5" & dane$change=="no_change_6" & dane$ELEV_class=="elev_0-250"

# plots are having two major layers
# one is all points in grey as a background
# second are the points of interest
# for better plotting only a random sample of all points is visualised (grey one)
par(mar=c(2,2,2,2))                 # proper margins of figure
selector<-sample(1:dim(dane)[1], 20000, replace=FALSE)  # random subsample

# points fulfilling association rules were plotted below
# plots are in fact making locations of subsets

# Fig.5a
plot(dane[selector, 4:5], col="grey80", bg="grey80", pch=21, cex=0.9) # grey background
points(dane$X_WGS[cond7], dane$Y_WGS[cond7], bg="red", col="red", pch=21, cex=0.4)
points(dane$X_WGS[cond6], dane$Y_WGS[cond6], bg="orange", col="orange", pch=21, cex=0.3)
points(dane$X_WGS[cond5], dane$Y_WGS[cond5], bg="yellow", col="yellow", pch=21, cex=0.2)
points(dane$X_WGS[cond4], dane$Y_WGS[cond4], bg="blue", col="dodgerblue", pch=21, cex=0.15)

title(main="Strong association rules for population changes, 
no land use changes (class 6) of land elevated 0-250 m")
legend("topleft", c("depopulation not stronger than 0.95", "population growth up tp 1.05", "population growth 1.05-1.2", "population growth 1.2-5"), fill=c("dodgerblue", "yellow", "orange", "red"), bty="n")

# Fig.5b
plot(dane[selector, 4:5], col="grey80", bg="grey80", pch=21, cex=0.9)
points(dane$X_WGS[cond1], dane$Y_WGS[cond1], bg="green", col="green", pch=21, cex=0.4)
points(dane$X_WGS[cond2], dane$Y_WGS[cond2], bg="dodgerblue", col="dodgerblue", pch=21, cex=0.4)
points(dane$X_WGS[cond3], dane$Y_WGS[cond3], bg="red", col="red", pch=21, cex=0.4)

title(main="Strong association rules for population and land use changes")
legend("topleft", c("always no inhabitants, 6to4 change", "depopulation between 20% and 50%, 6to4 change", "population growth 1.2-5, 6to3 change"), fill=c("green", "dodgerblue", "red"), bty="n")

Figure 5: Geolocation of strong association rules involving built-up land (class 6) as RHS

Probit models

This is the code for three models: situation of no changes, appearing class 6 (artificial cover) = urbanisation, disapppearing class 6 (artificial cover) = greening.

# Read packages
library(DescTools)  # for PseudoR2() to get R2 McFadden in probits
library(texreg)     # for screenreg() to make synthetic table with results
####### preparing data for probit models ################

# dummy variable for 'change to 6'
dane$change_to_6<-0
dane$change_to_6[dane$change=="1to6" | dane$change=="2to6" | dane$change=="3to6" | dane$change=="4to6" | dane$change=="5to6" | dane$change=="7to6"]<-1

# dummy variable for 'change from 6'
dane$change_from_6<-0
dane$change_from_6[dane$change=="6to1" | dane$change=="6to2" | dane$change=="6to3" | dane$change=="6to4" | dane$change=="6to5" | dane$change=="6to7"]<-1

# dummy variable for 'no change in 6'
dane$always_6<-0
dane$always_6[dane$change=="no_change_6"]<-1

# rescaling variables (for nicer coefficients)
dane$TOT_P_2006_K<-dane$TOT_P_2006/1000
dane$DIST_COAST_K<-dane$DIST_COAST/1000
dane$popul.change.06.18_K<-dane$popul.change.06.18/1000

# factorising labels (technical issue)
dane$change_to_6.f<-as.factor(dane$change_to_6)
dane$change_dummy.f<-as.factor(dane$change_dummy)
dane$change_from_6.f<-as.factor(dane$change_from_6)
dane$always_6.f<-as.factor(dane$always_6)

######### equations to estimate and models ###################

# change to 6 - new grids covered with artificial  urban cover - urbanisation
eq1<-change_to_6.f ~ ELEV_class + NUTS0_13 + TOT_P_2006_K + DIST_COAST_K + popul.change.06.18_K + popul.change.lab
probit1<-glm(eq1, family=binomial(link="probit"), data=dane)
summary(probit1)
PseudoR2(probit1, which = NULL)

# change from 6 - disappearing grids covered with artificial / urban cover - de-urbanisation
eq2<-change_from_6.f ~ ELEV_class + NUTS0_13 + TOT_P_2006_K + DIST_COAST_K + popul.change.06.18_K + popul.change.lab
probit2<-glm(eq2, family=binomial(link="probit"), data=dane)
summary(probit2)
PseudoR2(probit2, which = NULL)

# no change in 6 - stable grids covered with artificial / urban cover
eq3<-always_6.f ~ ELEV_class + NUTS0_13 + TOT_P_2006_K + DIST_COAST_K + popul.change.06.18_K + popul.change.lab
probit3<-glm(eq3, family=binomial(link="probit"), data=dane)
summary(probit3)
PseudoR2(probit3, which = NULL)

####### synthetic summary of all models in a single table
screenreg(list(probit1, probit2, probit3), custom.header=list("Appearing 6"=1, "Disappearing 6"=2, "Always 6"=3), digits=3)