RAISE.Rmd

---
title: "Exploring Air Quality"
author: "Daniela"
date: "5/26/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Import Data

The data came from the Air Quality Historical Data Platform, an open data framework which gathers data from worldwide EPA monitering programs.

```{r cars}
#airqualitybeijing.csv <- #"/home/CAMPUS/mwl04747/github/China_air_pollution_COVID/beijing-air-quality.csv"
airqualitybeijing.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/beijing-air-quality.csv"

AQbeijing <- read.csv(airqualitybeijing.csv)

str(AQbeijing)

```

## Scatter Plot

You can also embed plots, for example:

```{r pressure, echo=FALSE}
plot(AQbeijing$pm10)
```

Dates aren't correct -- Index...

## Next Steps

- Fix Dates

 - as.Date()
 
```{r fixdates}
# Fix the dates!
AQbeijing$NewDate = as.Date(AQbeijing$date)
str(AQbeijing)
head(AQbeijing)
```


```{r}
plot(pm25 ~ NewDate, data=AQbeijing)

```


 
- Select Data Before and After Lockdown. file is city_size_lockdowndate.csv
- Different for each city? Yes 
- 80 days

#merge files, or extract by city, date. before/after by date. subset(). how can I extract a window of dates by factor. create different data frames by each window. new column where: if date before lockdown, -1, if date after lockdown +1, get rid of zeros. create binary columns.

  - subset()
  
  
```{r subset}
# Dates for Beijing: "2020-02-10" Start

AQbeijingLockDown = as.Date("2020-02-10"); str(AQbeijingLockDown)

AQbeijingBefore = subset(AQbeijing, 
        subset=(NewDate > AQbeijingLockDown-80 &
                NewDate < AQbeijingLockDown))
AQbeijingAfter = subset(AQbeijing, 
        subset=(NewDate >= AQbeijingLockDown &
                NewDate < AQbeijingLockDown + 80))

plot(pm25~NewDate, AQbeijingBefore)
plot(pm25~NewDate, AQbeijingAfter)

```


- t-test Before and After
  - t.test()
  
  
```{r ttest}
#run ttest
t.test(AQbeijingBefore$pm25, AQbeijingAfter$pm25)


```

```{r ttest}
#extract mean_before, mean_after, p-value
ttest = t.test(AQbeijingBefore$pm25, AQbeijingAfter$pm25)
ttest$estimate
ttest$p.value


#export whole ttest
capture.output(ttest, file = "Welch_output.csv", append = TRUE)

#export just city, estimate, pvalue
estimate = ttest$estimate
pvalue = ttest$p.value
city = "Beijing"
capture.output(city, estimate, pvalue, file = "output.csv", append = TRUE)


```

#repeating import data, scatterplot, fix dates, welch 2-sample t test, output for the other cities.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Import Data

The data came from the Air Quality Historical Data Platform, an open data framework which gathers data from worldwide EPA monitering programs.

```{r cars}
#airqualitybeijing.csv <- #"/home/CAMPUS/mwl04747/github/China_air_pollution_COVID/beijing-air-quality.csv"


xianning.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/e/xianning-air-quality.csv"

AQ <- read.csv(xianning.csv)

str(AQ)

```

## Scatter Plot

You can also embed plots, for example:

```{r pressure, echo=FALSE}
plot(AQ$pm25)
```

Dates aren't correct -- Index...

## Next Steps

- Fix Dates

 - as.Date()
 
```{r fixdates}
# Fix the dates
AQ$NewDate = as.Date(AQ$date)
str(AQ)
head(AQ)
```


```{r}
plot(pm25 ~ NewDate, data=AQ)

```


 
- Select Data Before and After Lockdown
- Different for each city? Yes 
- 80 days

  - subset()
  
  
```{r subset}
# Dates for city: "yyyy-mm-dd" Start Lockdown

AQLockDown = as.Date("2020-01-24"); str(AQLockDown)

AQxianningBefore = subset(AQ, 
        subset=(NewDate > AQLockDown-80 &
                NewDate < AQLockDown))
AQxianningAfter = subset(AQ, 
        subset=(NewDate >= AQLockDown &
                NewDate < AQLockDown + 80))

plot(pm25~NewDate, AQxianningBefore) 
plot(pm25~NewDate, AQxianningAfter)

```

- t-test Before and After
  - t.test()
  
  
```{r ttest}
#run ttest
t.test(AQxianningBefore$pm25, AQxianningAfter$pm25)


```

```{r ttest}
#extract mean_before, mean_after, p-value
ttest = t.test(AQxianningBefore$pm25, AQxianningAfter$pm25)
ttest$estimate
ttest$p.value


#export whole ttest
capture.output(ttest, file = "Welch_output.csv", append = TRUE)

#export just city, estimate, pvalue
estimate = ttest$estimate
pvalue = ttest$p.value
city = "Xianning"
capture.output(city, estimate, pvalue, file = "output.csv", append = TRUE)

```


- BACI "Before and Ater Control-Impact"
  - aov() 
  
  
visualization:

- time series -- with abreak B/A
- box plot B/A



```{r}
# cityname, date, pm25: two files: one before, one after lockdown starts

#renamed cities in the file tbl from id, using sed at the command line 
#sed 's/"\/home\/CAMPUS\/dapa2019\/R\/China_air_pollution_COVID\/data\/baotou-air-quality\.csv"/baotou/g' tbl.csv>tbl2.csv
#renamed cities are in tbl_28cities.csv and tbl_2cities.csv. 
#tbl_28cities.csv includes 28 cities, excludes two problematic files (tianjin and xianning) which are in tbl_2cities.csv

library(dplyr)
library(readr)
files <- list.files(path = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data", pattern = "*.csv", full.names = T)

tbl.csv <- sapply(files, read_csv, simplify=FALSE) %>% 
bind_rows(.id = "id")

#write tbl as csv
write.csv(tbl, "tbl.csv")
head(tbl.csv)

#trying to include xianning and tianjin, which have no O3 values.
#opened tbl_2cities in BBEdit, found that o3 included TRUE (likely an unspecified trace amount) and NA, changed TRUE to NA and renamed to tbl_2cities_editedbetter.csv
#merge 28 and 2:
cities2 <- read.csv("/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/tbl/tbl_2cities_editedbetter.csv")
cities28  <- read.csv("/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/tbl/tbl_28cities.csv")

str(cities2)

cities30 <- rbind(cities2, cities28)
write.csv(cities30, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/tbl/cities30.csv")
head(cities30)

str(cities30) #date column is formatted as yyyy/mm/dd

```

## Next Steps

 
```{r fixdates}
# Fix the dates! Change from <factor> (date column) to <date> (newdate column).
cities30$newdate = as.Date(cities30$date) #newdate format is yyyy-mm-dd
str(cities30)
head(cities30)


```{r}
plot(pm25 ~ newdate, data=cities30, pch=20, main="Seasonal PM 2.5 Patterns",
  xlab="Time", ylab="Daily PM 2.5 AQI", col="darkblue") #points instead of circles, rotate y-axis, change labels from newdate etc. to be more readable, points--> grey. instert a moving average regression line. 
#look into AIC
help(blues9)
str(cities30)
```




The goal of the next section is to split the cities30.csv into four new csv files:
the 80 days leading up to the date of lockdown in 2020
the 80 days after and including the date of lockdown in 2020
the same two ranges of dates but for 2019



```{r}
#this file includes city name and the date of lockdown formatted yyyy-mm-dd
city_size_lockdowndate.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/city_size_lockdowndate.csv"
lockdowndate <- read.csv(city_size_lockdowndate.csv)

str(lockdowndate)


#before newdate > "2020-02-10"-80 & newdate < "2020-02-10")}

#after newdate >= "2020-02-10" & newdate < "2020-02-10"+80)}

#worked with a duplicate file (same as cities30) for experimenting
cities30_duplicate = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/tbl/cities30_duplicate.csv"
duplicate <- read.csv(cities30_duplicate)
newdate <- cities30$newdate

#CITY2020BEFORE GUIDELINE

beijing_beforelockdown2020 <- subset(cities30, subset=c(id=="beijing" & newdate > as.Date("2020-02-10")-80 & newdate < "2020-02-10"))
str(beijing_beforelockdown2020)
#write.csv(beijing_beforelockdown2020, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/beforelockdown2020/beijing_beforelockdown2020.csv")

#results can be accessed as .csv files or as data frames.

all_beforelockdown2020 <- do.call("rbind", list(beijing_beforelockdown2020, wuhan_beforelockdown2020, xianning_beforelockdown2020, tangshan_beforelockdown2020, wuxi_beforelockdown2020, nanjing_beforelockdown2020, hangzhou_beforelockdown2020, xuzhou_beforelockdown2020, changzhou_beforelockdown2020, dalian_beforelockdown2020, nanchang_beforelockdown2020, benxi_beforelockdown2020, dandong_beforelockdown2020, fuxin_beforelockdown2020, hefei_beforelockdown2020, huludao_beforelockdown2020, panjin_beforelockdown2020, yangzhou_beforelockdown2020, jinzhou_beforelockdown2020, tianjin_beforelockdown2020, yibin_beforelockdown2020, ganzhou_beforelockdown2020, guangzhou_beforelockdown2020, chuzhou_beforelockdown2020, suzhou_beforelockdown2020, fangchenggang_beforelockdown2020, cangzhou_beforelockdown2020, shanghai_beforelockdown2020, baotou_beforelockdown2020, chifeng_beforelockdown2020))

write.csv(all_beforelockdown2020, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_beforelockdown2020.csv")

#CITY2020AFTER GUIDELINE
beijing_afterlockdown2020 <- subset(cities30, subset=c(id=="beijing" & newdate >= "2020-02-10" & newdate < as.Date("2020-02-10")+80))
head(beijing_afterlockdown2020)
write.csv(beijing_afterlockdown2020, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/afterlockdown2020/beijing_afterlockdown2020.csv")

all_afterlockdown2020 <- do.call("rbind", list(beijing_afterlockdown2020, wuhan_afterlockdown2020, xianning_afterlockdown2020, tangshan_afterlockdown2020, wuxi_afterlockdown2020, nanjing_afterlockdown2020, hangzhou_afterlockdown2020, xuzhou_afterlockdown2020, changzhou_afterlockdown2020, dalian_afterlockdown2020, nanchang_afterlockdown2020, benxi_afterlockdown2020, dandong_afterlockdown2020, fuxin_afterlockdown2020, hefei_afterlockdown2020, huludao_afterlockdown2020, panjin_afterlockdown2020, yangzhou_afterlockdown2020, jinzhou_afterlockdown2020, tianjin_afterlockdown2020, yibin_afterlockdown2020, ganzhou_afterlockdown2020, guangzhou_afterlockdown2020, chuzhou_afterlockdown2020, suzhou_afterlockdown2020, fangchenggang_afterlockdown2020, cangzhou_afterlockdown2020, shanghai_afterlockdown2020, baotou_afterlockdown2020, chifeng_afterlockdown2020))

write.csv(all_afterlockdown2020, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_afterlockdown2020.csv")

head(all_afterlockdown2020)

```
repeating the above for 2019
```{r}

#CITY2019BEFORE GUIDELINE

chifeng_beforelockdown2019 <- subset(cities30, subset=c(id=="chifeng" & newdate > as.Date("2019-02-13")-80 & newdate < "2019-02-13"))
head(chifeng_beforelockdown2019)
write.csv(chifeng_beforelockdown2019, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/beforelockdown2019/chifeng_beforelockdown2019.csv")

all_beforelockdown2019 <- do.call("rbind", list(beijing_beforelockdown2019, wuhan_beforelockdown2019, xianning_beforelockdown2019, tangshan_beforelockdown2019, wuxi_beforelockdown2019, nanjing_beforelockdown2019, hangzhou_beforelockdown2019, xuzhou_beforelockdown2019, changzhou_beforelockdown2019, dalian_beforelockdown2019, nanchang_beforelockdown2019, benxi_beforelockdown2019, dandong_beforelockdown2019, fuxin_beforelockdown2019, hefei_beforelockdown2019, huludao_beforelockdown2019, panjin_beforelockdown2019, yangzhou_beforelockdown2019, jinzhou_beforelockdown2019, tianjin_beforelockdown2019, yibin_beforelockdown2019, ganzhou_beforelockdown2019, guangzhou_beforelockdown2019, chuzhou_beforelockdown2019, suzhou_beforelockdown2019, fangchenggang_beforelockdown2019, cangzhou_beforelockdown2019, shanghai_beforelockdown2019, baotou_beforelockdown2019, chifeng_beforelockdown2019))

write.csv(all_beforelockdown2019, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_beforelockdown2019.csv")


#CITY2019AFTER GUIDELINE
chifeng_afterlockdown2019 <- subset(cities30, subset=c(id=="chifeng" & newdate >= "2019-02-13" & newdate < as.Date("2019-02-13")+80))
head(chifeng_afterlockdown2019)
write.csv(chifeng_afterlockdown2019, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/afterlockdown2019/chifeng_afterlockdown2019.csv")

all_afterlockdown2019 <- do.call("rbind", list(beijing_afterlockdown2019, wuhan_afterlockdown2019, xianning_afterlockdown2019, tangshan_afterlockdown2019, wuxi_afterlockdown2019, nanjing_afterlockdown2019, hangzhou_afterlockdown2019, xuzhou_afterlockdown2019, changzhou_afterlockdown2019, dalian_afterlockdown2019, nanchang_afterlockdown2019, benxi_afterlockdown2019, dandong_afterlockdown2019, fuxin_afterlockdown2019, hefei_afterlockdown2019, huludao_afterlockdown2019, panjin_afterlockdown2019, yangzhou_afterlockdown2019, jinzhou_afterlockdown2019, tianjin_afterlockdown2019, yibin_afterlockdown2019, ganzhou_afterlockdown2019, guangzhou_afterlockdown2019, chuzhou_afterlockdown2019, suzhou_afterlockdown2019, fangchenggang_afterlockdown2019, cangzhou_afterlockdown2019, shanghai_afterlockdown2019, baotou_afterlockdown2019, chifeng_afterlockdown2019))

write.csv(all_afterlockdown2019, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_afterlockdown2019.csv")


all_2019 <- rbind(all_beforelockdown2019, all_afterlockdown2019)
all_2020 <- rbind(all_beforelockdown2020, all_afterlockdown2020)

```
 
#plotting all cities data
#scatterplots
#boxplots

```{r}


library(RColorBrewer)
display.brewer.all()

#points instead of circles, rotate y-axis, change labels from newdate etc. to be more readable, points--> grey. instert a moving average regression line. 

ggplot(all_afterlockdown2020, aes(newdate, pm25), color="firebrick")

#type-c
library(lattice)
xyplot(pm25 ~ newdate, data=all_afterlockdown2020, 
       auto.key=list(space="right"), main="Pollution After Lockdown", xlab="Date", ylab="PM 2.5", 
       jitter.x=TRUE, jitter.y=TRUE, type=c("p","smooth"), col=brewer.pal(3, "YlOrRd"), pch = 16, cex=.4)  
#cannot add groups=pm25. how to impose gradient?

install.packages("tidyverse")

g <- ggplot(all_afterlockdown2020, aes(newdate, pm25)) 
g + scale_colour_gradient(low = "yellow", high = "red", space = "Lab", na.value = "grey50", guide = "colourbar", aesthetics = "colour") #no result
g + scale_colour_brewer(type = "seq", palette = "YlOrRd", direction = 1, aesthetics = "colour")  #no result

#for LOESS
library(dplyr)
  all_afterlockdown2020 %>%
  mutate(loess = predict(loess(pm25 ~ newdate, data = all_afterlockdown2020))) %>%
  ggplot(aes(newdate, pm25)) +
  geom_point(color = "grey50") +
  geom_smooth(aes(pm25 = loess)) #not working
  


#ggplot(all_afterlockdown2020, aes(newdate, pm25)) + 
#geom_point(alpha=.4, size=.8, aes(colour=factor(pm25))) +
#geom_smooth(method = "loess", se = FALSE, colour="black") + theme(legend.position='none', panel.background = element_blank()) + scale_colour_brewer(type = "seq", palette = "YlOrRd", direction = 1, aesthetics = "colour") #n too large, allowed maximum for palette YlOrRd is 9 Returning the palette you asked for with that many colors

#ggplot(all_afterlockdown2020, aes(newdate, pm25)) + 
#geom_point(alpha=.4, size=.8, aes(colour=factor(pm25))) +
#geom_smooth(method = "loess", se = FALSE, colour="black") + theme(legend.position='none', panel.background = element_blank()) + scale_color_gradient2(low = "yellow", mid = "orange", high = "red", midpoint = mean(pm25)) #error discrete values applied to a continuous scale


ggplot(all_afterlockdown2020, aes(newdate, pm25), jitter.x=TRUE, jitter.y=TRUE, auto.key=list(space="right")) + labs(x="Date", y="PM 2.5", title="PM 2.5 Pollution Levels After 2020 Lockdowns") +
geom_point(alpha=.4, size=.8, aes(colour=factor(pm25))) +
geom_smooth(method = "loess", se = FALSE, colour="black") + theme(legend.position='none', panel.background = element_blank()) #WORKS but how to change the gradient color? I'm hoping to use the brewer color palatte YlOrRd


#overall patterns are more obvious here:
ggplot(all_2020, aes(newdate, pm25), jitter.x=TRUE, jitter.y=TRUE, auto.key=list(space="right")) + labs(x="Date", y="PM 2.5 AQI", title="PM 2.5 during 80 days before and after 2020 lockdown") +
geom_point(alpha=.4, size=.8, aes(colour=factor(pm25))) + 
geom_smooth(method = "loess", se = FALSE, colour="black") + theme(legend.position='none', panel.background = element_blank()) #hope to graph all2019 and all2020 with a "divider" marking lockdowndate.

ggplot(all_2019, aes(newdate, pm25), jitter.x=TRUE, jitter.y=TRUE, auto.key=list(space="right")) + labs(x="Date", y="PM 2.5 AQI", title="PM 2.5 during 80 days before and after 2019 lockdown") +
geom_point(alpha=.4, size=.8, aes(colour=factor(pm25))) + 
geom_smooth(method = "loess", se = FALSE, colour="black") + theme(legend.position='none', panel.background = element_blank())

plot(y=all_2020$pm25, x=all_2020$newdate)

plot(all_beforelockdown2020$newdate, all_beforelockdown2020$pm25, type="l",col="red", xlim=as.Date(c("2019-10-05","2020-06-01")))
lines(all_afterlockdown2020$newdate, all_afterlockdown2020$pm25,col="green") #fix date range.



#make grpah for each pollutatnt type
#fix boxplots to work with means
#find real population sizes (even if old) preferably from same source.
#ask email questions!


#learn overall format for a mult. regression
#proximity to factories
#coastal vs inland
#general geography (eg. elevation)
#duplicate_sorting = subset(duplicate, subset=c(id == "beijing"))
#duplicate_sorting2 = subset(duplicate_sorting, subset=c(newdate < "2020-02-10"))
#duplicate_sorting3 = subset(duplicate_sorting2, subset=c(newdate > as.Date("2020-02-10")-80))

```
```{r boxplots}

#neither work right now.

geom_boxplot(outlier.colour="black", outlier.shape=16,
             outlier.size=2, notch=FALSE)

boxplot(pm25~newdate,data=all_afterlockdown2020, main="Pollution After 2020 Lockdown",
   xlab="PM 2.5", ylab="Date")

#next: boxplot of the city means before, boxplot of the city means after.

#NEW idea: given file with id, newdate, before vs after, pollutant levels

all_beforelockdown2020$timegroup = "before_2020"
all_afterlockdown2020$timegroup = "after_2020"

head(all_afterlockdown2020)

all_2020_timegroup <- rbind(all_beforelockdown2020, all_afterlockdown2020)
write.csv(all_2020_timegroup, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_2020_timegroup.csv")
str(all_2020_timegroup) #created timegroup as chr

boxplot(pm25~timegroup,data=all_2020_timegroup, main="PM 2.5 Before and After 2020 Lockdown",
   xlab="time groupings ranging 80 days each", ylab="PM 2.5") #works, but need to reoirent x axis

#all four timegroups: 
all_beforelockdown2019$timegroup = "before_2019"
all_afterlockdown2019$timegroup = "after_2019"
all_2019_timegroup <- rbind(all_beforelockdown2019, all_afterlockdown2019)
write.csv(all_2019_timegroup, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_2019_timegroup.csv")
str(all_2019_timegroup)

all_four_timegroup <- rbind(all_2019_timegroup, all_2020_timegroup)
write.csv(all_four_timegroup, "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/data/all_four_timegroup.csv")
str(all_four_timegroup)

boxplot(pm25~timegroup,data=all_four_timegroup, main="PM 2.5 Before and After Lockdown in 2020 and 2019", xlab="time groupings ranging 80 days each", ylab="PM 2.5") #works, and can use this format for all pollutants, but first need to reoirent x axis. also.. why am I getting so many outliers???
#boxplot for all data, not just the means.


#boxplot for just the means #works but only for pm25 in 2020.
cities_means_pm25.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/cities_means_pm25_revised.csv"
cities_means_pm25 <- read.csv(cities_means_pm25.csv)
boxplot(pm25~before_or_after,data=cities_means_pm25, main="Pollution Means Before and After 2020 Lockdowns", xlab="city", ylab="mean pm2.5")
head(cities_means_pm25)


```
creating boxplots of means for each pollutant.
first: create one csv with mean data from all cities.

```{r}

#find the mean for each pollutant within each time frame for each city.
#created: 
#means_allcities_allpollutants_beforelockdown2019.csv 
#means_allcities_allpollutants_afterlockdown2019.csv 
#means_allcities_allpollutants_beforelockdown2020.csv 
#means_allcities_allpollutants_afterlockdown2020.csv 

description = "yibin_afterlockdown2020"

capture.output(description, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)

citymean <- mean(yibin_afterlockdown2020$pm25, trim=0, na.rm=TRUE)
pollutiontype = "pm25"
capture.output(pollutiontype, citymean, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)

citymean <- mean(yibin_afterlockdown2020$pm10, trim=0, na.rm=TRUE)
pollutiontype = "pm10"
capture.output(pollutiontype, citymean, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)

citymean <- mean(yibin_afterlockdown2020$o3, trim=0, na.rm=TRUE)
pollutiontype = "o3"
capture.output(pollutiontype, citymean, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)

citymean <- mean(yibin_afterlockdown2020$no2, trim=0, na.rm=TRUE)
pollutiontype = "no2"
capture.output(pollutiontype, citymean, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)

citymean <- mean(yibin_afterlockdown2020$so2, trim=0, na.rm=TRUE)
pollutiontype = "so2"
capture.output(pollutiontype, citymean, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)

citymean <- mean(yibin_afterlockdown2020$co, trim=0, na.rm=TRUE)
pollutiontype = "co"
capture.output(pollutiontype, citymean, file= "means_allcities_allpollutants_afterlockdown2020.csv", append = TRUE)


#edited together manually in Excel

allmeans.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/means_allcities_allpollutants_alltimeframes.csv"

allmeans <- read.csv(allmeans.csv)

str(allmeans)

#boxplot showing all cities, able to specify individual pollutants, using the mean of the pollutant for each city, grouped by all four time frames.

boxplot(so2~timeframe,data=allmeans, main="Mean Sulfur Dioxide Before and After Lockdown", xlab="80 Day Time Frames", ylab="Mean SO2 AQI", names=c("2019 before","2019 after","2020 before","2020 after"), col=c('powderblue', 'powderblue','mistyrose','mistyrose'))
allmeans$timeframe <- factor(allmeans$timeframe , levels=c("2019_before", "2019_after", "2020_before", "2020_after"))
 
#legend("topright", legend = c("80 days before the 2019 equivalent of the date lockdowns started", "80 days after the 2019 equivalent of the date lockdowns started", "80 days before the date lockdowns started in 2020", "80 days after the date lockdowns started in 2020"))
#trouble including legend; it blocks view of graph.

#creates tables summarizing the data shown in the boxplots
summary()

```


Running t-tests

```{r regression}

#how much of the variation in each pollutant is accounted for by time frame (by year? by before/after?), by population?

#hypothesis: pollutant levels decreased significantly after lockdown in 2020, as compared with before lockdown, controlling for seasonal effects.

#to test what's going wrong, I created regression2020only.csv file.
#also created a file taht includes comparing beforevsafter in 2020 and 2020vs2019 for "after" timeframe
regression.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/regression.csv"
regression <- read.csv(regression.csv)

plot <- lm(pm25 ~ timeframe+population, regression) #only three time frames showingup, and the tiny values make no sense
summary(plot)
anova(plot) #why are DF = 1 for both? 

#some of the variation in the pollutants can be explained by city population size
#pollution decreases significantly during the 2020 lockdowns


# x = time, x' = population, x" = year, y = pm25
# data frame is: city, [year, before/after], population, proximity to nearest factory, each pollutant.
#doing * instead of + can show interaction.
#look at covariance, anova, assumptions


```

```{r ttest}
 
#extract mean_before, mean_after, p-value
ttest = t.test(all_beforelockdown2020$pm25, all_afterlockdown2020$pm25)
ttest$estimate
ttest$p.value

#export whole ttest
capture.output(ttest, file = "Welch_output_allcities.csv", append = TRUE)

#export just description, estimate, pvalue
estimate = ttest$estimate
pvalue = ttest$p.value
description = "all_beforelockdown2020 vs all_afterlockdown2020"
capture.output(description, estimate, pvalue, file = "output_allcities.csv", append = TRUE)

```

run BACI analysis to see how % decreased

```{r BACI}
#group into 3 city size categories --> run an ANOVA
#pollutant as a function of after - before and city size. BACI.
#aov delta citysize
#before = WINTER 
#after = SPRING


#RESTRUCTURE FILE eg. pm25_before_2020, _after_2020, ditto for 2019.
# delta = 2019(winter-spring) - 2020(winter-spring)

regression_withdelta.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/regression_withdelta.csv"
regression_withdelta <- read.csv(regression_withdelta.csv)

str(regression_withdelta)

pm25.aov <- aov(pm25_2020after ~ pm25_2020minus2019 + population_discrete, data=regression_withdelta) #broke population into three categories
summary(pm25.aov) #run for every pollutant
#plot(pm25_2020after ~ pm25_2020minus2019, data=regression_withdelta)
#hist(regression_withdelta$pm25_2020afterminus2020before)
#hist(regression_withdelta$pm25_2019afterminus2019before)
#hist(regression_withdelta$pm25_2020minus2019)

#the anova is ancova? population is categorical, and have about the same number of cities in each category. analysis of covariance. see if slope and p-value are significant based on *
#presentation: tell a narrative, showing results thorughout.



#alternative approach: run two regressions per pollutant, one for each yearly data set. tehn compare experimentla and control results.
#expect population to explain more of the 2020 data set.

regression2020only.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/regression2020only.csv"
regression2020only <- read.csv(regression2020only.csv)
str(regression_withdelta)

pm25.aov <- aov(pm25 ~ timeframe + population, data=regression2020only) #timeframe is significant
summary(pm25.aov)

#alt approach edited for population to be discrete


regression2020only.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/regression2020only.csv"
regression2020only <- read.csv(regression2020only.csv)
str(regression2020only)

pm25.aov <- aov(pm25 ~ timeframe + population_discrete, data=regression2020only) #timeframe is significant
summary(pm25.aov)

#boxplots for 3 city sizes
#scatterplot color coded for covid
```

```{r}

#now trying a hierarchical multiple regression
#(pm25_2020after ~ population + pm25before + pm25_2019after, dataset)

hierarchical_reg.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/hierarchical_reg.csv"
hierarchical_reg <- read.csv(hierarchical_reg.csv)
str(hierarchical_reg)

var1 <- lm(formula = pm25_after2020 ~ population_30cities, data = hierarchical_reg) #plot of residuals shows this might not be a linear model
var2 <- lm(formula = pm25_after2020 ~ population_30cities + pm25_after2019, data = hierarchical_reg)
var3 <- lm(formula = pm25_after2020 ~ population_30cities + pm25_after2019 + pm25_before2020, data = hierarchical_reg)
anova(var1,var2,var3)

#reformatted
hierarchical_reg_reformatted.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/hierarchical_reg_reformatted.csv"
hierarchical_reg_reformatted <- read.csv(hierarchical_reg_reformatted.csv)
str(hierarchical_reg_reformatted)

var1 <- lm(formula = pm25_all2020 ~ population_discrete, data = hierarchical_reg_reformatted) #plot of residuals shows this might not be a linear model
plot(var1) #see constant leverage
var2 <- lm(formula = pm25_all2020 ~ population_discrete + pm25_all2019, data = hierarchical_reg_reformatted)
var3 <- lm(formula = pm25_all2020 ~ population_discrete + pm25_all2019 + beforevsafter, data = hierarchical_reg_reformatted)
anova(var1,var2,var3)

#revised to include industrial GRP
hierarchical_reg_grp.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/hierarchical_reg_discrete_withfactories.csv"
hierarchical_reg_grp <- read.csv(hierarchical_reg_grp.csv)
str(hierarchical_reg_grp)

var_initial <- lm(formula = pm25_all2020 ~ 1, data = hierarchical_reg_grp)
var1 <- lm(formula = pm25_all2020 ~ population_30cities, data = hierarchical_reg_grp) #plot of residuals shows this might not be a linear model
plot(var1) #see constant leverage
var2 <- lm(formula = pm25_all2020 ~ population_30cities + industry, data = hierarchical_reg_grp)
var3 <- lm(formula = pm25_all2020 ~ population_30cities + industry + pm25_all2019, data = hierarchical_reg_grp)
var4 <- lm(formula = pm25_all2020 ~ population_30cities + industry + pm25_all2019 + beforevsafter, data = hierarchical_reg_grp)
anova(var_initial, var1,var2,var3, var4)

summary(var2)

#changed the dependent variable to be pm25_after2020 and reformatted the csv
hierarchical_reg_grp_2.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/hierarchical_reg_discrete_withfactories_reformatteddependentvar.csv"
hierarchical_reg_grp_2 <- read.csv(hierarchical_reg_grp_2.csv)
str(hierarchical_reg_grp_2)

var_initial <- lm(formula = pm25_2020after ~ 1, data = hierarchical_reg_grp_2)
var1 <- lm(formula = pm25_2020after ~ population_30cities, data = hierarchical_reg_grp_2) #plot of residuals shows this might not be a linear model
var2 <- lm(formula = pm25_2020after ~ population_30cities + industry, data = hierarchical_reg_grp_2)
var3 <- lm(formula = pm25_2020after ~ population_30cities + industry + pm25_2019after, data = hierarchical_reg_grp_2)
var4 <- lm(formula = pm25_2020after ~ population_30cities + industry + pm25_2019after + pm25_2020before, data = hierarchical_reg_grp_2)
anova(var_initial, var1,var2,var3, var4)

summary(var4)

#summary(lm(formula = pm25_2020after ~ pm25_2020before, data = hierarchical_reg_grp_2))


#new method: AIC

AIC(var_initial,var1,var2,var3,var4)
#high evidence for model 3 and 4. amount of evidence supporting 4 is slightly greater than 3.
#penalizes for number of parameters

hierarchical_reg_grp.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/hierarchical_reg_discrete_withfactories.csv"
hierarchical_reg_grp <- read.csv(hierarchical_reg_grp.csv)
str(hierarchical_reg_grp)

var_initial <- lm(formula = pm25_all2020 ~ 1, data = hierarchical_reg_grp)
var1 <- lm(formula = pm25_all2020 ~ population_30cities, data = hierarchical_reg_grp) #plot of residuals shows this might not be a linear model
#plot(var1) #see constant leverage
var2 <- lm(formula = pm25_all2020 ~ population_30cities + industry, data = hierarchical_reg_grp)
AIC(var_initial, var1,var2)

summary(var2)

#or...
hierarchical_reg_grp_2.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/hierarchical_reg_discrete_withfactories_reformatteddependentvar.csv"
hierarchical_reg_grp_2 <- read.csv(hierarchical_reg_grp_2.csv)
str(hierarchical_reg_grp_2)

var_initial <- lm(formula = pm25_2020after ~ 1, data = hierarchical_reg_grp_2)
var1 <- lm(formula = pm25_2020after ~ population_30cities, data = hierarchical_reg_grp_2) #plot of residuals shows this might not be a linear model
var2 <- lm(formula = pm25_2020after ~ population_30cities + industry, data = hierarchical_reg_grp_2)
AIC(var_initial, var1,var2)

summary(var2)

help(AIC)



```

#delta approach results
#only reports change not decrease!
                   Df Sum Sq Mean Sq F value Pr(>F)
pm25_2020minus2019  1     66   66.17   0.242  0.627
population          1    251  251.36   0.918  0.346
Residuals          27   7390  273.70               
90 observations deleted due to missingness

#delta approach results with 3 discrete population groupings
                    Df Sum Sq Mean Sq F value Pr(>F)
pm25_2020minus2019   1    388   388.0   1.436  0.242
population_discrete  2    295   147.6   0.546  0.586
Residuals           26   7024   270.2               
90 observations deleted due to missingness


#alternative approach results
#2019only
            Df Sum Sq Mean Sq F value Pr(>F)  
timeframe    1   2900  2899.5   6.652 0.0125 *
population   1    875   875.5   2.008 0.1619  
Residuals   57  24847   435.9  

#2020only
            Df Sum Sq Mean Sq F value   Pr(>F)    
timeframe    1  13941   13941  28.597 1.64e-06 ***
population   1   1719    1719   3.527   0.0655 .  
Residuals   57  27788     488                     


#hierarchical multiple regression. 
#(degrees of freedom here are confusing)
Analysis of Variance Table 

Model 1: pm25_after2020 ~ population_30cities
Model 2: pm25_after2020 ~ population_30cities + pm25_after2019
Model 3: pm25_after2020 ~ population_30cities + pm25_after2019 + pm25_before2020
  Res.Df    RSS Df Sum of Sq       F    Pr(>F)    
1     28 7445.3                                   
2     27 2121.1  1    5324.2 88.9942 7.037e-10 ***
3     26 1555.5  1     565.6  9.4537  0.004906 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

#reformatted. seems more effective!
Analysis of Variance Table

Model 1: pm25_all2020 ~ population_30cities
Model 2: pm25_all2020 ~ population_30cities + pm25_all2019
Model 3: pm25_all2020 ~ population_30cities + pm25_all2019 + beforevsafter
  Res.Df   RSS Df Sum of Sq       F    Pr(>F)    
1     58 41729                                   
2     57 15830  1   25899.3 132.352 2.236e-16 ***
3     56 10958  1    4871.7  24.896 6.218e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#reformatted, and making population discrete makes little/no difference.
Analysis of Variance Table

Model 1: pm25_all2020 ~ population_discrete
Model 2: pm25_all2020 ~ population_discrete + pm25_all2019
Model 3: pm25_all2020 ~ population_discrete + pm25_all2019 + beforevsafter
  Res.Df   RSS Df Sum of Sq       F    Pr(>F)    
1     57 41131                                   
2     56 15026  1     26105 141.863 < 2.2e-16 ***
3     55 10121  1      4905  26.655 3.457e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#CHANGE to include GRP

#including industrial gross regional product where the region is each city's province or autonomous region
Analysis of Variance Table

Model 1: pm25_all2020 ~ 1
Model 2: pm25_all2020 ~ population_discrete
Model 3: pm25_all2020 ~ population_discrete + industry
Model 4: pm25_all2020 ~ population_discrete + industry + pm25_all2019
Model 5: pm25_all2020 ~ population_discrete + industry + pm25_all2019 + 
    beforevsafter
  Res.Df   RSS Df Sum of Sq        F    Pr(>F)    
1     59 43449                                    
2     57 41131  2    2317.7   7.5306  0.001304 ** 
3     56 41065  1      66.3   0.4307  0.514426    
4     55 12698  1   28366.8 184.3378 < 2.2e-16 ***
5     54  8310  1    4388.3  28.5167   1.9e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#summary statistics: questions:
#why are there 54 DF???
#why is medium population significant, and big nonexistent in the data?
Call:
lm(formula = pm25_all2020 ~ population_discrete + industry + 
    pm25_all2019 + beforevsafter, data = hierarchical_reg_grp)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.142  -7.517  -0.605   5.147  66.137 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -0.1059525  9.9437266  -0.011  0.99154    
population_discretemedium  9.6812602  3.9366970   2.459  0.01716 *  
population_discretesmall   4.0572823  4.1974786   0.967  0.33805    
industry                  -0.0005519  0.0001609  -3.431  0.00116 ** 
pm25_all2019               0.8896606  0.0804611  11.057 1.73e-15 ***
beforevsafterbefore       18.1173713  3.3927008   5.340 1.90e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.41 on 54 degrees of freedom
Multiple R-squared:  0.8087,	Adjusted R-squared:  0.791 
F-statistic: 45.67 on 5 and 54 DF,  p-value: < 2.2e-16

#same as above with continuous not discrete population
Analysis of Variance Table

Model 1: pm25_all2020 ~ 1
Model 2: pm25_all2020 ~ population_30cities
Model 3: pm25_all2020 ~ population_30cities + industry
Model 4: pm25_all2020 ~ population_30cities + industry + pm25_all2019
Model 5: pm25_all2020 ~ population_30cities + industry + pm25_all2019 + 
    beforevsafter
  Res.Df   RSS Df Sum of Sq        F    Pr(>F)    
1     59 43449                                    
2     58 41729  1    1719.5  10.4843  0.002043 ** 
3     57 41484  1     245.1   1.4946  0.226721    
4     56 13451  1   28033.2 170.9278 < 2.2e-16 ***
5     55  9020  1    4430.7  27.0156 3.051e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#summary statistics
Call:
lm(formula = pm25_all2020 ~ population_30cities + industry + 
    pm25_all2019 + beforevsafter, data = hierarchical_reg_grp)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.091  -6.211  -1.306   5.058  70.779 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          6.354e+00  1.065e+01   0.597  0.55320    
population_30cities -3.157e-07  2.684e-07  -1.176  0.24458    
industry            -5.250e-04  1.527e-04  -3.438  0.00113 ** 
pm25_all2019         8.833e-01  8.312e-02  10.627 6.04e-15 ***
beforevsafterbefore  1.821e+01  3.503e+00   5.198 3.05e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.81 on 55 degrees of freedom
Multiple R-squared:  0.7924,	Adjusted R-squared:  0.7773 
F-statistic: 52.48 on 4 and 55 DF,  p-value: < 2.2e-16

#multiple r-squared values from the summary() for each individual var#. 
#Compared r squared values to check for multicollinearity and none found.
var1 Multiple R-squared:  0.03958 #population
var2 Multiple R-squared:  0.04522 #industrial GRP
var3 Multiple R-squared:  0.6904 #expected seasonal changes
var4 Multiple R-squared:  0.7924 #lockdown



#reformatted to change the dependent variable (unsure if I should even have done this)
Analysis of Variance Table

Model 1: pm25_2020after ~ 1
Model 2: pm25_2020after ~ population_30cities
Model 3: pm25_2020after ~ population_30cities + industry
Model 4: pm25_2020after ~ population_30cities + industry + pm25_2019after
Model 5: pm25_2020after ~ population_30cities + industry + pm25_2019after + pm25_2020before
  Res.Df    RSS Df Sum of Sq        F    Pr(>F)    
1     29 7707.5                                    
2     28 7445.3  1     262.2   5.4270   0.02820 *  
3     27 7391.8  1      53.5   1.1080   0.30260    
4     26 1461.4  1    5930.4 122.7348 3.903e-11 ***
5     25 1208.0  1     253.4   5.2438   0.03074 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#with summary(var#) hierarchical method
#more meaningful statistic
var1 Multiple R-squared:  0.03402 #population
var2 Multiple R-squared:  0.04097 #industrial GRP
var3 Multiple R-squared:  0.8104 #compared to the same months in 2019
var4 Multiple R-squared:  0.8433 #compared to before lockdown in 2020

#with summary(lm(formula = pm25_2020after ~ #one_variable, data = hierarchical_reg_grp_2)) 
#less meaningful statistic
population_30cities Multiple R-squared:  0.03402
industry            Multiple R-squared:  0.006941
pm25_2019after      Multiple R-squared:  0.7133
pm25_2020before     Multiple R-squared:  0.5757


#using AIC
var         df AIC
var_initial	2	255.5991		
var1	3	256.5606		
var2	4	258.3441		
var3	5	211.7139		
var4	6	208.0015	
#the smallest wins!
#high evidence for model 3 and 4. Amount of evidence supporting 4 is slightly greater than 3.
#penalizes for number of parameters

when redone, with dependent var = all 2020
var1	3	568.9495	#population	
var2	4	570.5961	#industrial grp

when redone, with dependent var = 2020after
var1	3	256.5606	#population
var2	4	258.3441	#industrial grp
#this just means we can trust the baci data more.
```{r}

baci.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/baci.csv"
baci <- read.csv(baci.csv)
str(baci)

attach(baci)

baci[Event=="before_lockdown",2]
mean(baci[Event=="before_lockdown",2]) #for C_2019 137.9755

baci[Event=="after_lockdown",2]
mean(baci[Event=="after_lockdown",2]) #for C_2019 124.0723

baci[Event=="before_lockdown",3]
mean(baci[Event=="before_lockdown",3]) #for I_2020 137.097

baci[Event=="after_lockdown",3]
mean(baci[Event=="after_lockdown",3]) #for I_2020 106.6104


baci2.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/baci_updated.csv"
baci2 <- read.csv(baci2.csv)
str(baci2)
baci2$controlvsimpact = as.factor(baci2$X2019vs2020)

baci.aov <- aov(pm2.5 ~ controlvsimpact + beforevsafter + controlvsimpact:beforevsafter, data=baci2); summary(baci.aov) 

#boxplots
boxplot(pm2.5 ~ controlvsimpact + beforevsafter, data=baci2)
```

#redoing BACI with a new data set including all the pollutants' means.

```{r}

baci_allmeans.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/baci_allmeans.csv"
baci_allmeans <- read.csv(baci_allmeans.csv)
str(baci_allmeans)

attach(baci_allmeans)

baci_allmeans[Event=="before_lockdown",7]
mean(baci_allmeans[Event=="before_lockdown",7], na.rm=TRUE) #for C_2019 137.9755

baci_allmeans[Event=="after_lockdown",7]
mean(baci_allmeans[Event=="after_lockdown",7], na.rm=TRUE) #for C_2019 124.0723

baci_allmeans[Event=="before_lockdown",13]
mean(baci_allmeans[Event=="before_lockdown",13], na.rm=TRUE) #for I_2020 137.097

baci_allmeans[Event=="after_lockdown",13]
mean(baci_allmeans[Event=="after_lockdown",13], na.rm=TRUE) #for I_2020 106.6104


baci2_allmeans.csv = "/home/CAMPUS/dapa2019/R/China_air_pollution_COVID/means/baci2_allmeans.csv"
baci2_allmeans <- read.csv(baci2_allmeans.csv)
str(baci2_allmeans)
baci2_allmeans$controlvsimpact = as.factor(baci2_allmeans$X2019vs2020)

baci.aov <- aov(co ~ beforevsafter + controlvsimpact + beforevsafter:controlvsimpact, data=baci2_allmeans); summary(baci.aov) 


```

#pm2.5

baci[Event=="before_lockdown",2]
mean(baci[Event=="before_lockdown",2]) #for C_2019 137.9755

baci[Event=="after_lockdown",2]
mean(baci[Event=="after_lockdown",2]) #for C_2019 124.0723

baci[Event=="before_lockdown",3]
mean(baci[Event=="before_lockdown",3]) #for I_2020 137.097

baci[Event=="after_lockdown",3]
mean(baci[Event=="after_lockdown",3]) #for I_2020 106.6104

means:
beforelockdown 2019 137.9755
afterlockdown 2019 124.0723
beforelockdown 2020 137.097
afterlockdown 2020 106.6104

                               Df Sum Sq Mean Sq F value   Pr(>F)    
beforevsafter                   1  14778   14778  31.039 1.66e-07 *** #seasonal effect
controlvsimpact                 1   2523    2523   5.299   0.0231 *  #covid year
beforevsafter:controlvsimpact   1   2063    2063   4.332   0.0396 *  #there is an interaction between the years and the before and after. both contribute to the diff in the means. the relationship is nonlinear. when you add control vs impact there is a non additive effect.
Residuals                     116  55230     476                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


#moving on to other pollutants

#pm10
beforelockdown 2019 [1] 77.62988
afterlockdown 2019 [1] 67.50539
beforelockdown 2020 [1] 70.61246
afterlockdown 2020 [1] 55.9869

                               Df Sum Sq Mean Sq F value   Pr(>F)    
beforevsafter                   1   4594    4594  14.125 0.000269 ***
controlvsimpact                 1   2577    2577   7.922 0.005738 ** 
beforevsafter:controlvsimpact   1    152     152   0.467 0.495664    
Residuals                     116  37731     325                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#o3
beforelockdown 2019 [1] 23.58492
afterlockdown 2019 [1] 39.49967
beforelockdown 2020 [1] 25.28063
afterlockdown 2020 [1] 39.29786

                              Df Sum Sq Mean Sq F value Pr(>F)    
beforevsafter                  1   5685    5685 216.937 <2e-16 ***
controlvsimpact                1     14      14   0.525  0.471    
beforevsafter:controlvsimpact  1     22      22   0.846  0.360    
Residuals                     96   2516      26                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
20 observations deleted due to missingness

#no2
beforelockdown 2019 [1] 23.76107
afterlockdown 2019 [1] 19.95249
beforelockdown 2020 [1] 22.61039
afterlockdown 2020 [1] 14.92706

                               Df Sum Sq Mean Sq F value   Pr(>F)    
beforevsafter                   1    990   990.5  36.047 2.26e-08 ***
controlvsimpact                 1    286   286.1  10.412  0.00163 ** 
beforevsafter:controlvsimpact   1    113   112.6   4.098  0.04523 *  
Residuals                     116   3187    27.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#so2
beforelockdown 2019 [1] 13.00346
afterlockdown 2019 [1] 9.460548
beforelockdown 2020 [1] 11.07439
afterlockdown 2020 [1] 7.28235

                               Df Sum Sq Mean Sq F value Pr(>F)  
beforevsafter                   1    404   403.5   5.124 0.0255 *
controlvsimpact                 1    127   126.5   1.607 0.2075  
beforevsafter:controlvsimpact   1      0     0.5   0.006 0.9388  
Residuals                     116   9134    78.7                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#co
beforelockdown 2019 [1] 12.11218
afterlockdown 2019 [1] 9.254269
beforelockdown 2020 [1] 11.72867
afterlockdown 2020 [1] 7.674415

Show in New WindowClear OutputExpand/Collapse Output
'data.frame':	2328 obs. of  7 variables:
 $ date: Factor w/ 2328 levels "2014/1/1","2014/1/10",..: 2317 2322 2323 2324 2325 2326 2327 2328 2307 2308 ...
 $ pm25: int  189 149 143 109 150 160 160 107 100 72 ...
 $ pm10: int  65 58 49 48 53 61 26 42 52 101 ...
 $ o3  : int  101 38 21 33 92 36 40 26 35 42 ...
 $ no2 : int  8 12 11 12 15 16 8 11 9 4 ...
 $ so2 : int  1 2 2 NA 1 4 NA NA NA NA ...
 $ co  : int  5 10 6 5 6 12 8 5 4 1 ...
Show in New WindowClear OutputExpand/Collapse Output
 
 
X
<int>
id
<fctr>
date
<fctr>
pm25
<int>
pm10
<int>
o3
<int>
no2
<int>
so2
<int>
co
<int>
1	1	tianjin	2020/7/2	94	53	NA	8	3	9	
2	2	tianjin	2020/7/3	141	NA	NA	NA	NA	NA	
3	3	tianjin	2020/6/2	84	78	NA	13	1	6	
4	4	tianjin	2020/6/3	102	180	NA	16	2	7	
5	5	tianjin	2020/6/4	145	101	NA	18	5	13	
6	6	tianjin	2020/6/5	97	38	NA	14	6	7	
6 rows | 1-10 of 10 columns
Show in New WindowClear OutputExpand/Collapse Output

Show in New WindowClear OutputExpand/Collapse Output
 
 
X
<int>
id
<fctr>
date
<fctr>
pm25
<int>
pm10
<int>
o3
<int>
no2
<int>
so2
<int>
co
<int>
17008	12301	chifeng	2019/1/2	95	64	18	25	29	13	
17009	12302	chifeng	2019/1/3	121	84	21	21	26	12	
17010	12303	chifeng	2019/1/4	149	29	27	11	8	5	
17011	12304	chifeng	2019/1/5	52	40	22	18	20	8	
17012	12305	chifeng	2019/1/6	77	87	6	26	45	18	
17013	12306	chifeng	2019/1/7	151	58	23	17	16	11	
6 rows | 1-10 of 10 columns
Show in New WindowClear OutputExpand/Collapse Output
Error: unexpected symbol in:
"
all_2019"
Show in New WindowClear OutputExpand/Collapse Output
	
`geom_smooth()` using formula 'y ~ x'
	
Removed 30 rows containing non-finite values (stat_smooth).
Removed 30 rows containing missing values (geom_point).

Show in New WindowClear OutputExpand/Collapse Output

Show in New WindowClear OutputExpand/Collapse Output
'data.frame':	120 obs. of  8 variables:
 $ city     : Factor w/ 30 levels "baotou","beijing",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ pm25     : num  150 113 148 161 137 ...
 $ pm10     : num  107.6 72.4 92.8 101.5 69.4 ...
 $ o3       : num  24.6 17.5 24.1 22.1 19.1 ...
 $ no2      : num  25.8 23 22.9 26.1 27.3 ...
 $ so2      : num  27.48 4.69 23.32 16.31 5.27 ...
 $ co       : num  21.63 9.18 21.04 12.46 9 ...
 $ timeframe: Factor w/ 4 levels "2019_after","2019_before",..: 2 2 2 2 2 2 2 2 2 2 ...
Show in New WindowClear OutputExpand/Collapse Output

Call:
lm(formula = pm25 ~ timeframe + population, data = regression)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.440 -10.142   1.848  10.683  85.181 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.275e+02  4.177e+00  30.532  < 2e-16 ***
timeframe2019_before  1.390e+01  5.528e+00   2.515  0.01328 *  
timeframe2020_after  -1.746e+01  5.528e+00  -3.159  0.00202 ** 
timeframe2020_before  1.302e+01  5.528e+00   2.356  0.02015 *  
population           -7.313e-07  3.116e-07  -2.347  0.02064 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.41 on 115 degrees of freedom
Multiple R-squared:  0.2934,	Adjusted R-squared:  0.2689 
F-statistic: 11.94 on 4 and 115 DF,  p-value: 3.789e-08

Show in New WindowClear OutputExpand/Collapse Output

Call:
lm(formula = pm25_all2020 ~ population_30cities + industry, data = hierarchical_reg_grp)

Residuals:
    Min      1Q  Median      3Q     Max 
-54.611 -15.740  -2.356  13.122  98.320 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          1.286e+02  6.410e+00  20.064   <2e-16 ***
population_30cities -8.536e-07  5.553e-07  -1.537    0.130    
industry            -1.825e-04  3.145e-04  -0.580    0.564    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.98 on 57 degrees of freedom
Multiple R-squared:  0.04522,	Adjusted R-squared:  0.01172 
F-statistic:  1.35 on 2 and 57 DF,  p-value: 0.2675

                               Df Sum Sq Mean Sq F value   Pr(>F)    
beforevsafter                   1  358.3   358.3  23.329 4.21e-06 ***
controlvsimpact                 1   28.9    28.9   1.882    0.173    
beforevsafter:controlvsimpact   1   10.7    10.7   0.699    0.405    
Residuals                     116 1781.8    15.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1