This is an example of how to visualize effect modification using R package {effects}.
People living in a greener neighborhood, measured by NDVI, reported higher life satisfaction if they had little financial worries. For people with high financial worries, greener neighborhood was associated with lower life satisfaction. See the original publication for details.
### Not to run ###
library(betareg)
library(effects)
covco <- readRDS('data/covco_greenspace.rds')
# transform life satisfaction score into a proportion
covco$satisf_now_prop <- (covco$satisf_now+0.01)/10.1
covco$satisf_past_prop <- (covco$satisf_past+0.01)/10.1
# interaction with financial worries
model5 <- satisf_now_prop ~ ndvi_buff300 * finwor_3cat + satisf_past_prop +
agecat + sex + hh_income + kanton + urbanity + mon_entry_2cat
fit5 <- betareg(model5, data = covco)
png('results/ndvi_finwor.png', res = 300, height = 2000, width = 2000)
plot(predictorEffect('ndvi_buff300', fit5),
lines=list(multiline=TRUE),
confint=list(style='auto', alpha=0.1),
axes=list(x=list(ndvi_buff300=list(lab='NDVI 300m buffer')),
y=list(transform=function(x) (exp(x)/(1+exp(x)))*10.1-0.01,
lab='Life satisfaction')),
lattice=list(key.args=list(title='Financial worries')),
main = '', rug=FALSE, ylim=c(6.7,8.2))
dev.off()