Model_risk_covid19_individual_github.Rmd

---
title: "Ecological Risk Model"
author: "Nicolas Ayala, Toni Monleon. BIOST3 U. of Barcelona ; Jaume Canela. Hospital Clinic of Barcelona"
date: "04/04/2022"
output: html_document
---
```{r setup, include=FALSE}
#knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(echo = F, message=FALSE, warning=FALSE)
```

## Representation of data in Spain and Catalonia for the article
Data until 22-7-2020
```{r}
library(BDSbiost3)
```


```{r}
 library(readr)
datos_ccaas <- read_csv("datos_ccaas.csv")
#datos de catalunya
cat<- datos_ccaas[ which(datos_ccaas$ccaa_iso=='CT'), ]
plot(cat$num_casos~cat$fecha,type="l",col="blue",xlab="date", ylab="daily cases (COVID19)")
 
library(BDSbiost3)
 #Predict evolution of cases in Spain with an "ad-hoc" model ()
pred.Casos.diarios<-covid.cases.cumcases3(X=seq(1:168), Y=cat$num_casos, max.cases.free=T, max.cases = max(cat$num_casos), lim.days=365,iterations = 20, label.ind="CATALONIA: PREDICTION COVID-19 (DAY1: 31/01/2020)", duration_pandemic=200)

#Number of total cases
sum((cat$num_casos))

sum((cat$num_casos))/7609499 #per inhabitants

100000*sum((cat$num_casos))/7609499 #per 100.000 inhabitants


```

##Times Series Plot
```{r}
library(dplyr)
library(ggplot2)
library(RColorBrewer)

Sys.setlocale("LC_ALL", "English")
theme_update(plot.title = element_text(hjust = 0.5))
ggplot(data=cat, aes(y=cat$num_casos, x=cat$fecha )) +
  geom_line(color="steelblue")+
   theme_bw()+
   ylim(0, 2000)+
   ggtitle (" COVID-19 daily cases from February to July in 2020")+
   xlab("2020")+
   ylab("COVID-19 Cases")+
   scale_x_date(date_labels = "%B")+
   theme(axis.text.x = element_text(angle=45, hjust = 1, ))+
   theme_minimal()

max(cat$num_casos, na.rm = TRUE)
```

## Analysis of a COVID19 risk model
Data management based on incidence data, IDESCAT data and air pollution data:
```{r }

#Data on incidence rates and density
library(readxl)                         
model_risk<- read_excel("dadesCanela_modelrisk (1).xlsx", sheet = "JUNTO_TODO", skip = 1)

colnames(model_risk)= c("health_area","positive_cases","suspected_cases","raw_rate_10k","standard_rate_10k","insured_people ","municipality","population","surface_km2", "density", "rank_number_inhabitant","municipality_1","NOX_reduction_january2020","NOX_reduction_february2020","NOX_reduction_march2020","NOX_reduction_april2020","NOX_reduction_may2020","health_area_1","RFDB (milers d'euros)","RFDB per habitant (milers d'euros)","RFDB per habitant (índex Catalunya=100)")

head(model_risk)
```

## Exploratory analysis of population density, wealth and pollution reduction (indirect mobility)
Below is an exploratory analysis of the regressors and the response variable and different predictive models of the risk of coronavirus.
Some descriptive analyzes of the rate of COVID-19 in municipalities of Catalonia and density and population, population and reduction of pollution by NOx during 2017 - 2020 are presented.

# Mobility using indirect variables (densidad y poblacion por municipio)
Data obtained from IDESCAT
```{r}
library(Rmisc)
library(gplots)
library(ggplot2)
library(DescTools)
#DESCRIPTIVE OF NO2 BY YEAR
res.RMISC.Q <- summarySE(model_risk, measurevar="density", groupvars=c("municipality"),na.rm = T)
res.RMISC.Q 

```

# Mobility using indirect variables
Data obtained from IDESCAT.
https://www.idescat.cat/pub/?id=aec&tema=macro&n=82&lang=es 

income_euros :  RFDB (milers d'euros)
income_euros_inhab: RFDB per habitant (milers d'euros)
income_index_100: RFDB per habitant (índex Catalunya=100)

```{r}
library(Rmisc)
library(gplots)
library(ggplot2)
library(DescTools)

model_risk$income_euros<- as.numeric(model_risk$`RFDB (milers d'euros)`)
model_risk$income_euros_inhab<- as.numeric(model_risk$`RFDB per habitant (milers d'euros)`)
model_risk$income_index_100<- as.numeric(model_risk$`RFDB per habitant (índex Catalunya=100)`)


#Descriptive analysis for family incomes in thousands of euros (income_euros).
res.RMISC.Q <- summarySE(model_risk, measurevar="income_euros", groupvars=c("municipality"),na.rm = T)
res.RMISC.Q #descriptive


```

#Descriptive analysis for family incomes in thousands of euros per inhabitants (income_euros_inhab).
```{r}
##Descriptive analysis for family incomes in thousands of euros per inhabitants (income_euros_inhab). 
res.RMISC.Q <- summarySE(model_risk, measurevar="income_euros_inhab", groupvars=c("municipality"),na.rm = T)
res.RMISC.Q #descriptive analysis
```

#Descriptive analysis for family incomes as index (index Catalonia = 100)
```{r}
##Descriptive analysis for family incomes as index (index Catalonia = 100)
res.RMISC.Q <- summarySE(model_risk, measurevar="income_index_100", groupvars=c("municipality"),na.rm = T)
res.RMISC.Q #descriptiva
```

# Mobility using indirect variables (contamination in the control points)
Nox will be studied in aggregates from different control points in Catalonia by year (2007-2020), month and municipality.
See source: 
https://analisi.transparenciacatalunya.cat/Medi-Ambient/Dades-d-immissi-dels-punts-de-mesurament-de-la-Xar/uy6k-2s8r/data

```{r}
## open data frame "DATOS_AGREGADOS_2017_2020.Rda" whit variables as NOx per month, year and Catalonia municipality
load("DATOS_AGREGADOS_2017_2020.Rda")
DATOS_AGREGADOS_2017_2020.1<-DATOS_AGREGADOS_2017_2020

#Delete infinite data
is.na(DATOS_AGREGADOS_2017_2020.1) <- sapply(DATOS_AGREGADOS_2017_2020.1, is.infinite)
DATOS_AGREGADOS_2017_2020.1$MUNICIPI<-DATOS_AGREGADOS_2017_2020.1$Group.3

#Cols variables
colnames(DATOS_AGREGADOS_2017_2020.1) <- c("Group.1", "Group.2", "Group.3", "year", "month", "municipality", "NO.MAX")
head(DATOS_AGREGADOS_2017_2020.1)
```

Summary of NOx by year
```{r}
library(Rmisc)
library(gplots)
library(ggplot2)
library(DescTools)

#Descriptive analyses of  NOx by year.
res.RMISC.Q <- summarySE(DATOS_AGREGADOS_2017_2020.1, measurevar="NO.MAX", groupvars=c("year"),na.rm = T)
res.RMISC.Q 
plotmeans(NO.MAX~ year, data=DATOS_AGREGADOS_2017_2020.1, frame = T,xlab = "YEAR OF ANALYSIS", ylab="maximum NOx") 
```

Nox data per month and year.
```{r}

#NOx per month and year.

res.RMISC.Q1 <- summarySE(DATOS_AGREGADOS_2017_2020.1, measurevar="NO.MAX", groupvars=c("year","month"),na.rm = T)
res.RMISC.Q1 #descriptiva
plotmeans(NO.MAX~ month, data=DATOS_AGREGADOS_2017_2020.1, frame = T,xlab = "MONTH OF ANALYSIS", ylab="maximum NOx") 
```


Nox data per municipality, month and year:
```{r}
#Nox data per municipality, month and year
res.RMISC.Q <- summarySE(DATOS_AGREGADOS_2017_2020.1, measurevar="NO.MAX", groupvars=c("municipality", "year","month"),na.rm = T)
AGGREGATENOX.MUN.YEAR.MONTH<-res.RMISC.Q #descriptive
plotmeans(NO.MAX~ municipality , data=DATOS_AGREGADOS_2017_2020.1, frame = T,xlab = "CITY", ylab="maximum NOx") #PLOT DE LAB

```

Calculation of aggregated data
```{r}
#Calculation of aggregate data of NOX per municipality, year and month.
#AGGREGATENOX.MUN.YEAR.MONTH

MAM_2020 <- AGGREGATENOX.MUN.YEAR.MONTH[ which(AGGREGATENOX.MUN.YEAR.MONTH$year=="2020"), c(1,2,3,5)]
MAM_2019 <- AGGREGATENOX.MUN.YEAR.MONTH[ which(AGGREGATENOX.MUN.YEAR.MONTH$year=="2019"), c(1,2,3,5) ]
MAM_2018 <- AGGREGATENOX.MUN.YEAR.MONTH[ which(AGGREGATENOX.MUN.YEAR.MONTH$year=="2018"), c(1,2,3,5) ]
MAM_2017 <- AGGREGATENOX.MUN.YEAR.MONTH[ which(AGGREGATENOX.MUN.YEAR.MONTH$year=="2017"), c(1,2,3,5) ]

#Aggregate value 2020-2017
total.agregat1<- merge(MAM_2020,MAM_2019, by=c("municipality","month"))
total.agregat2<- merge(total.agregat1,MAM_2018, by=c("municipality","month"))
total.agregat3<- merge(total.agregat2,MAM_2017, by=c("municipality","month"))


#names(total.agregat3 )
colnames( total.agregat3)<- c("municipality", "month",      "year.x",    "NO.MAX.2020", "year.y", "NO.MAX.2019", "year.x1",   "NO.MAX.2018", "year.y1",    "NO.MAX.2017")
#percentage difference between the average per month of 2017-2019 compared to 2020.

for(i in 1:dim(total.agregat3)[1]){
   #i<-9
   total.agregat3$media2017_2019[i] <- 
      mean(c(total.agregat3$NO.MAX.2019[i] ,
                total.agregat3$NO.MAX.2018[i] , 
                total.agregat3$NO.MAX.2017[i]),na.rm = T) 
   
   total.agregat3$per.var2020_otros[i]<-100*(total.agregat3$media2017_2019[i]
                                      -total.agregat3$NO.MAX.2020[i])/total.agregat3$media2017_2019[i]
}
```

Percentage reduction (%) of Nox in the months of the year 2017-2020:
```{r}

#percent nox reduction
#Descriptive analysis of NOx per month, year and municipality
res.1 <- summarySE(total.agregat3, measurevar="per.var2020_otros", groupvars=c("municipality","month"),na.rm = T)
res.1

plotmeans(per.var2020_otros~month, data=total.agregat3, frame = T,xlab = "% REDUCCION NOX 2017-19 VS 2020 CATALUNYA") #PLOT DE LAB
```
Nox percentage reduction graph by municipality
```{r}
library(gplots)
plotmeans(per.var2020_otros~ municipality, data=total.agregat3, frame = T, main = "% NOX REDUCTION 2017-19 VS 2020 IN CATALONIA", xlab="CITY") 
```

## Plots of the relationship between population density, pollution and COVID-19
Raw rate vs density inhabitants/km2
```{r pressure, echo=FALSE}
plot(model_risk$raw_rate_10k ~ model_risk$density, data=model_risk, xlab="Density inhabitants/km2", ylab="raw rate ")


```
Standard rate vs Density.
```{r}
plot (model_risk$standard_rate_10k ~ model_risk$density, data=model_risk, xlab="Density inhab. /km2", ylab="Standard rate ")
```

Positive Cases of COVID-19 and Density
```{r}
plot(model_risk$positive_cases ~ model_risk$density, data=model_risk, xlab="Density inhab./km2", ylab="Positive Cases COVID-19")
```


Standard rate and NOx reduction (%) in March 2020 compared to 2017-2019
```{r}
plot(model_risk$standard_rate_10k ~ model_risk$NOX_reduction_march2020, data=model_risk, col="red", xlab="NOx reduction (%)", ylab="Standard rate COVID-19", xlim=c(0, 60))
```


Standard rate and NOx reduction (%) in April 2020 compared to 2017-2019
```{r}
plot(model_risk$standard_rate_10k ~ model_risk$NOX_reduction_april2020, data=model_risk, col="red", xlab="Nox reduction", ylab="standard_rate_10k COVID19", xlim=c(0, 80))
```

Standard rate and NOx reduction (%) in May 2020 compared to 2017-2019
```{r}
plot(model_risk$standard_rate_10k ~ model_risk$NOX_reduction_may2020, data=model_risk, col="red", xlab="Nox reduction", ylab="standard_rate_10k COVID19", xlim=c(0, 110))
```

Standard rate of COVID19 and  NOx reduction (%) in May 2020 compared to May 2017-2019
```{r}
plot(model_risk$standard_rate_10k ~ model_risk$NOX_reduction_may2020, data=model_risk, col="red", xlab="NOx reduction in may 2020", ylab="standard_rate_10k COVID19", xlim=c(0, 120))
```


## Modeling of the risk of COVID-19 in Catalonia in aggregate data.
For this modeling we will follow what is indicated by Morgensternn (1982). H Morgenstern. Uses of Ecologic Analysis in Epidemiologic Research. Am J Public Health, 1982 Dec;72(12):1336-44.

Abstract:Despite the widespread use of ecologic analysis in epidemiologic research and health planning, little attention has been given by health scientists and practitioners to the methodological aspects of this approach. This paper reviews the major types of ecologic study designs, the analytic methods appropriate for each, the limitations of ecologic data for making causal inferences and what can be done to minimize these problems, and the relative advantages of ecologic analysis. Numerous examples are provided to illustrate the important principles and methods. A careful distinction is made between ecologic studies that generate or test etiologic hypotheses and those that evaluate the impact of intervention programs or policies (given adequate knowledge of disease etiology). Failure to recognize this difference in the conduct of ecologic studies can lead to results that are not very informative or that are misinterpreted by others.

Where it is indicated how to calculate the relative ecological risk of an event (in this case COVID) with respect to a continuous risk factor (eg mobility, population density) continuous using a regression function. It is a very interesting approach for aggregated data as in this case.

## Linear regression and Relative Risk.
The incidence-type model ~f(mobilitat, densitat) seeks to determine the relative ecological risk (Morgensen, 1982) that a person has the disease and the population risk.
The level of aggregation in this case is: municipality/portion present in the system of the "Generalitat de Catalunya".
The dependent variables (Y) to predict are:
-standard_rate_10k: Standard rate (10.000 inhab) of COVID-19
-raw_rate_10_k: Raw rate of COVID-19
-positive_cases: Positive cases of COVID-19

Indepent variables(Xi) are:
-density: Population density in inhabitants per km2
-population: Population of each area of Catalonia
-surface_km2: Surface in km2
-income_euros: Family income in thousands of euros (Renta Familiar Disponible RFDB)
-income_euros_inhab: RFDB per inhabitant
-income_index_100: RFDB index (Catalonia index = 100)
-nox_reduction_march_2020: Nox reduction in March 2020 compared to March 2017-2019
-nox_reduction_april_2020: Nox reduction in April 2020 compared to April 2017-2019
-nox_reduction_may_2020: Nox reduction in May 2020 compared to May 2017-2019

Multiple regression models with regressors (density, population):
```{r}
# raw rate of COVID19
res1<-lm(raw_rate_10k ~ density + surface_km2 + population , data=model_risk)
summary(res1)
```

Standard rate and population + density + surface
```{r}
# taxa_ESTANDARITZADA
res2<-lm(standard_rate_10k ~ surface_km2 + population , data=model_risk)
summary(res2)

```

Positive cases of COVID-19 and population
```{r}
#casos_positius
res3<-lm(positive_cases~ density + surface_km2 + population, data=model_risk)
summary(res3)
```

Standard rate of COVID-19 and density
```{r}
# taxa_ESTANDARITZADA
res2<-lm(standard_rate_10k ~ density , data=model_risk)
summary(res2)

```

**Relative Ecologic Risk Morguenson 1982 ** 
Standard rate of COVID-19 compared to Density
```{r}
# Standard rate of COVID-19
res2<-lm(standard_rate_10k~ density , data=model_risk)
#summary(res2)
print("Relative Ecologic Risk ()")
(res2$coefficients[1] + res2$coefficients[2])/res2$coefficients[1] 


################################################################################
#compute RRE and its CI95 using resampling
#Function RRE
RRE<-function(X,Y,num.itera=100,nsize=10){
   DM<-data.frame(X=X, Y=Y)
   res2<-lm(Y ~ X , data=DM)

   #Relative ecologic risk (mean and CI95)
   RRE_mean<-(res2$coefficients[1] + res2$coefficients[2])/res2$coefficients[1] 
   print("========================================================")
   cat("Relative ecologic risk (RRE) [estimation]= ",RRE_mean,sep="")
   print("========================================================")
   
   #bootstrap Ci95 percent
   #length.Y<-num.itera
   my.array.RRE <- array(NA, dim=c(num.itera,1))
   for (i in 1:num.itera){
      DM.sample<-DM[sample(nrow(DM), nsize,replace = T), ]
      
      res3<-lm(Y ~ X , data=DM.sample)
      RRE_mean<-(res3$coefficients[1] + res3$coefficients[2])/res3$coefficients[1] 
      my.array.RRE[i,1] <- RRE_mean #valor real
      #print(RRE_mean)
  
   }

   #ci95
  library(Rmisc)
  res.REE.CI<-CI(my.array.RRE[,1], ci=0.95)   
  print("========================================================")
  #Relative ecologic risk (mean and CI95)
  cat("Relative ecologic risk (RRE) CI95 interval [UL, mean, LL]= ",res.REE.CI,sep=",")
  print("========================================================")
  print("")
   
  return(res.REE.CI)
}

final.RRE<-RRE(X=model_risk$density ,Y=model_risk$standard_rate_10k,num.itera=1000,nsize=200)

final.RRE

```

Pollution as an indirect measure of mobility
Multiple regression models with regressors (pollution reduction):

Standard rate of COVID19 and pollution reduction
```{r}
#change data format
model_risk$NOX_reduction_january2020<- as.numeric(model_risk$NOX_reduction_january2020)
model_risk$NOX_reduction_february2020<- as.numeric(model_risk$NOX_reduction_february2020)
model_risk$NOX_reduction_march2020<- as.numeric(model_risk$NOX_reduction_march2020)
model_risk$NOX_reduction_april2020<- as.numeric(model_risk$NOX_reduction_april2020)
model_risk$NOX_reduction_may2020<- as.numeric(model_risk$NOX_reduction_may2020)
```

Regression models for pollution reduction
standard_rate_10k vs pullution reduction (March 2020) 
**Relative Ecologic Risk Morguenson 1982 ** 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_march2020 , data=model_risk)
summary(res1)

print("**Relative Ecologic Risk (RRE)")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 


#using resamplen to compute CI95
final.RRE<-RRE(X=model_risk$NOX_reduction_march2020,Y=model_risk$standard_rate_10k,
               num.itera=1000,nsize=200)

final.RRE
```

Standard rate of COVID19 and pollution reduction (April 2020) and **Ecological Relative Risk Morguenson 1982 **:
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_april2020 , data=model_risk)
summary(res1)

print("Relative Ecologic Risk")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 

#using resamplen to compute CI95
final.RRE<-RRE(X=model_risk$NOX_reduction_april2020,Y=model_risk$standard_rate_10k,
               num.itera=1000,nsize=200)

final.RRE
```
Standard rate of COVID19 and pollution reduction (May 2020) and **Ecological Relative Risk Morguenson 1982 **:
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_may2020 , data=model_risk)
summary(res1)

print("Relative Ecologic Risk")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 

#using resamplen to compute CI95
final.RRE<-RRE(X=model_risk$NOX_reduction_may2020,Y=model_risk$standard_rate_10k,
               num.itera=100,nsize=200)
final.RRE
```



Standard rate of COVID19 vs Pollution reduction (%) using March, April and May:
```{r}
res1<-lm(positive_cases ~ NOX_reduction_march2020 + NOX_reduction_april2020 + NOX_reduction_may2020 , data=model_risk)
summary(res1)

```

## ###################################################################################
## UNIVARIATE MODELS OF COVID-19 RISK ---Y=standard_rate_10k: --------------
## ###################################################################################
The level of aggregation in this case is: municipality/portion present in the system of the "Generalitat de Catalunya"
Below are multiple regression models with various regressors (density, population):
# COVID-19 risk modeling (Y= standard_rate_10k, X=  density)

Standard rate of COVID-19 (standard_rate_10k):
```{r}
res1<-lm(standard_rate_10k ~ density , data=model_risk)
summary(res1)
print("Relative Ecologic Risk---------------------------------:")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 
```

# COVID-19 risk modeling  (Y= standard_rate_10k, X=  income_euros)
Familiar Income (RFDB (milers d'euros)) 	
```{r}
res1<-lm(standard_rate_10k ~ income_euros , data=model_risk)
summary(res1)
print("Relative Ecologic Risk---------------------------------:")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 
```

# COVID-19 risk modeling (Y= standard_rate_10k, X=  income_euros_inhab)
Familiar Income per inhabitant (RFDB per habitant (milers d'euros))	
```{r}
res1<-lm(standard_rate_10k ~ income_euros_inhab , data=model_risk)
print(res1)
summary(res1)

print("Relative Ecologic Risk---------------------------------:")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 

final.RRE<-RRE(X=model_risk$income_euros_inhab,Y=model_risk$standard_rate_10k,
               num.itera=1000,nsize=200)
final.RRE
```

# COVID-19 risk modeling (Y= standard_rate_10k, X=  income_euros_inhab)
Familiar Income per Inhabitante as Index (RFDB per habitant (índex Catalunya=100))
```{r}
res1<-lm(standard_rate_10k ~ income_index_100 , data=model_risk)
summary(res1)
print("Relative Ecologic Risk---------------------------------:")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 
```

# COVID-19 risk modeling  (Y= standard_rate_10k, X=  pollution reduction in April 2020)
X= Pollution reduction in march (AVERAGE PERCENTAGE OF REDUCTION IN AUTOMATIC SEASONS IN APRIL 2020 COMPARED TO APRIL 2017-2019)
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_april2020 , data=model_risk)
summary(res1)
print("Relative Ecologic Risk---------------------------------:")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 
```

# COVID-19 risk modeling  (Y= standard_rate_10k, X=  pollution reduction in May 2020)
X= Pollution reduction in March (AVERAGE PERCENTAGE OF REDUCTION IN AUTOMATIC SEASONS  IN APRIL 2020 COMPARED TO MAY 2017-2019)
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_may2020 , data=model_risk)
summary(res1)
print("Relative Ecologic Risk---------------------------------:")
(res1$coefficients[1] + res1$coefficients[2])/res1$coefficients[1] 
```

## ###################################################################################
## MULTIVARIATE MODELS OF COVID-19 RISK ---Y=standard_rate_10k: -----------------------
## ###################################################################################

# COVID-19 risk modeling with interaction (Y= standard_rate_10k, X= pollution reduction(%) in March, April, May, income_euros_inhab, density)
Linear regression with no interaction for March 2020
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_march2020 + density  + income_euros_inhab, data=model_risk)
summary(res1)
```

Linear regression with  interaction for March 2020 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_march2020 + density + NOX_reduction_march2020* density , data=model_risk)
summary(res1)
```

Linear regression with no interaction for April 2020 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_april2020 + density  , data=model_risk)
summary(res1)
```

Linear regression with interaction for April 2020 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_april2020 + density + NOX_reduction_april2020* density , data=model_risk)
summary(res1)
```

Linear regression with no interaction for May 2020 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_may2020 + density  , data=model_risk)
summary(res1)
```

Linear regression with interaction for May 2020 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_may2020 + density + NOX_reduction_may2020* density , data=model_risk)
summary(res1)
```

Linear regression with interaction for May 2020 and income 
```{r}
res1<-lm(standard_rate_10k ~ NOX_reduction_may2020 + density + income_euros_inhab +
            NOX_reduction_may2020*income_euros_inhab
            , data=model_risk)
summary(res1)
```