Final project on regression class
Directed by Professor Kim, Sun Young
Thang, TB
19 December 2018
import data
dat=read.csv("C:/Users/BINH THANG/Dropbox/PhD/Lectures/regression/Final assigment/mydata.txt")
names(dat)## [1] "ID" "longitude" "latitude" "log.m_to_a1"
## [5] "log.m_to_a2" "log.m_to_a3" "log.m_to_coast" "log.m_to_comm"
## [9] "ll_a1_s00500" "ll_a3_s00500" "ll_a1_s01000" "ll_a3_s01000"
## [13] "pop_s01000" "pop_s01500" "pop_s02000" "age1c"
## [17] "gender1" "wtlb1" "htcm1" "waistcm1"
## [21] "hipcm1" "bmi1c" "cig1c" "sbp1c"
## [25] "pm25"Q1:
Method 1
t(summary(dat[,c("pm25" ,"longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000", "pop_s01000" , "pop_s01500" , "pop_s02000" )]))##
## pm25 Min. :16.01 1st Qu.:18.54 Median :20.98
## longitude Min. :-118.4 1st Qu.:-118.2 Median :-118.1
## latitude Min. :33.95 1st Qu.:34.03 Median :34.08
## log.m_to_a1 Min. :1.778 1st Qu.:2.821 Median :3.107
## log.m_to_a2 Min. :1.362 1st Qu.:3.065 Median :3.522
## log.m_to_a3 Min. :1.114 1st Qu.:2.124 Median :2.403
## log.m_to_coast Min. :3.801 1st Qu.:4.363 Median :4.398
## log.m_to_comm Min. :1.000 1st Qu.:2.107 Median :2.505
## ll_a1_s00500 Min. : 10.0 1st Qu.: 10.0 Median : 10.0
## ll_a3_s00500 Min. : 10.0 1st Qu.: 382.5 Median :1014.0
## ll_a1_s01000 Min. : 10 1st Qu.: 10 Median : 10
## ll_a3_s01000 Min. : 699 1st Qu.: 4006 Median : 4896
## pop_s01000 Min. : 3804 1st Qu.: 7763 Median :10146
## pop_s01500 Min. : 8808 1st Qu.:17403 Median :21625
## pop_s02000 Min. : 14499 1st Qu.: 31148 Median : 40059
##
## pm25 Mean :20.87 3rd Qu.:22.49 Max. :30.13
## longitude Mean :-118.1 3rd Qu.:-118.1 Max. :-118.0
## latitude Mean :34.07 3rd Qu.:34.12 Max. :34.19
## log.m_to_a1 Mean :3.043 3rd Qu.:3.291 Max. :3.611
## log.m_to_a2 Mean :3.364 3rd Qu.:3.699 Max. :4.033
## log.m_to_a3 Mean :2.340 3rd Qu.:2.639 Max. :2.972
## log.m_to_coast Mean :4.323 3rd Qu.:4.398 Max. :4.398
## log.m_to_comm Mean :2.288 3rd Qu.:2.664 Max. :2.995
## ll_a1_s00500 Mean : 277.7 3rd Qu.: 10.0 Max. :3310.0
## ll_a3_s00500 Mean :1342.8 3rd Qu.:2043.5 Max. :4095.0
## ll_a1_s01000 Mean :1709 3rd Qu.:2876 Max. :9709
## ll_a3_s01000 Mean : 5649 3rd Qu.: 6652 Max. :18287
## pop_s01000 Mean :11118 3rd Qu.:12294 Max. :34367
## pop_s01500 Mean :24739 3rd Qu.:27514 Max. :87118
## pop_s02000 Mean : 44081 3rd Qu.: 51618 Max. :155190method 2
# install packeage (stargazer)
library(stargazer)## Warning: package 'stargazer' was built under R version 3.4.4##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerdf <- data.frame(dat)
cols <- c("pm25" ,"longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000", "pop_s01000" , "pop_s01500" , "pop_s02000" )
stargazer(
df[, cols], type = "text",
summary.stat = c("min", "p25", "median", "p75", "max", "mean", "sd")
)##
## =================================================================================
## Statistic Min Pctl(25) Median Pctl(75) Max Mean St. Dev.
## ---------------------------------------------------------------------------------
## pm25 16.010 18.543 20.977 22.492 30.131 20.874 2.976
## longitude -118.433 -118.187 -118.116 -118.064 -117.956 -118.143 0.114
## latitude 33.950 34.034 34.077 34.122 34.188 34.073 0.062
## log.m_to_a1 1.778 2.821 3.107 3.291 3.611 3.043 0.387
## log.m_to_a2 1.362 3.065 3.522 3.699 4.033 3.364 0.531
## log.m_to_a3 1.114 2.124 2.403 2.639 2.972 2.340 0.414
## log.m_to_coast 3.801 4.363 4.398 4.398 4.398 4.323 0.154
## log.m_to_comm 1.000 2.107 2.505 2.664 2.995 2.288 0.569
## ll_a1_s00500 10 10 10 10 3,310 277.709 750.456
## ll_a3_s00500 10 382.5 1,014 2,043.5 4,095 1,342.764 1,168.999
## ll_a1_s01000 10 10 10 2,876 9,709 1,708.982 2,693.932
## ll_a3_s01000 699 4,006 4,896 6,651.5 18,287 5,649.473 2,976.603
## pop_s01000 3,804 7,763 10,146 12,294 34,367 11,117.730 5,301.950
## pop_s01500 8,808 17,403 21,625 27,513.5 87,118 24,738.730 12,384.190
## pop_s02000 14,499 31,147.5 40,059 51,617.5 155,190 44,081.440 22,521.920
## ---------------------------------------------------------------------------------Q2: correlation and VIF
Correlation
rsubset=subset(dat, select=c("pm25" ,"longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000", "pop_s01000" , "pop_s01500" , "pop_s02000"))
r <- cor(dat[,c("longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000", "pop_s01000" , "pop_s01500" , "pop_s02000")])
cor(dat[,c("pm25","longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000", "pop_s01000" , "pop_s01500" , "pop_s02000")])[1,] ## pm25 longitude latitude log.m_to_a1 log.m_to_a2
## 1.0000000 -0.1145770 -0.1632393 -0.6933055 -0.2819862
## log.m_to_a3 log.m_to_coast log.m_to_comm ll_a1_s00500 ll_a3_s00500
## -0.2574758 0.1433480 -0.3891957 0.3428126 0.6868936
## ll_a1_s01000 ll_a3_s01000 pop_s01000 pop_s01500 pop_s02000
## 0.3602977 0.6250312 0.2985561 0.4073528 0.4562132other method for r and p-value
# Hide upper triangle to get r (correlation)
library(xtable)## Warning: package 'xtable' was built under R version 3.4.4mcor<-round(cor(rsubset),2)
mcor## pm25 longitude latitude log.m_to_a1 log.m_to_a2
## pm25 1.00 -0.11 -0.16 -0.69 -0.28
## longitude -0.11 1.00 0.40 0.06 -0.10
## latitude -0.16 0.40 1.00 0.12 0.17
## log.m_to_a1 -0.69 0.06 0.12 1.00 -0.14
## log.m_to_a2 -0.28 -0.10 0.17 -0.14 1.00
## log.m_to_a3 -0.26 0.11 -0.04 -0.09 0.15
## log.m_to_coast 0.14 0.84 0.50 -0.11 -0.20
## log.m_to_comm -0.39 0.00 0.12 0.07 0.09
## ll_a1_s00500 0.34 0.07 -0.06 -0.70 0.09
## ll_a3_s00500 0.69 -0.25 -0.11 -0.21 0.04
## ll_a1_s01000 0.36 -0.08 -0.18 -0.72 0.15
## ll_a3_s01000 0.63 -0.30 -0.20 -0.28 0.07
## pop_s01000 0.30 -0.42 -0.42 -0.31 -0.02
## pop_s01500 0.41 -0.47 -0.38 -0.39 -0.05
## pop_s02000 0.46 -0.51 -0.38 -0.41 -0.03
## log.m_to_a3 log.m_to_coast log.m_to_comm ll_a1_s00500
## pm25 -0.26 0.14 -0.39 0.34
## longitude 0.11 0.84 0.00 0.07
## latitude -0.04 0.50 0.12 -0.06
## log.m_to_a1 -0.09 -0.11 0.07 -0.70
## log.m_to_a2 0.15 -0.20 0.09 0.09
## log.m_to_a3 1.00 0.01 0.25 0.13
## log.m_to_coast 0.01 1.00 0.00 0.12
## log.m_to_comm 0.25 0.00 1.00 0.09
## ll_a1_s00500 0.13 0.12 0.09 1.00
## ll_a3_s00500 -0.69 -0.12 -0.48 -0.06
## ll_a1_s01000 0.10 -0.01 0.01 0.64
## ll_a3_s01000 -0.45 -0.20 -0.49 0.01
## pop_s01000 -0.08 -0.33 -0.17 0.17
## pop_s01500 -0.05 -0.34 -0.21 0.21
## pop_s02000 -0.07 -0.36 -0.24 0.21
## ll_a3_s00500 ll_a1_s01000 ll_a3_s01000 pop_s01000
## pm25 0.69 0.36 0.63 0.30
## longitude -0.25 -0.08 -0.30 -0.42
## latitude -0.11 -0.18 -0.20 -0.42
## log.m_to_a1 -0.21 -0.72 -0.28 -0.31
## log.m_to_a2 0.04 0.15 0.07 -0.02
## log.m_to_a3 -0.69 0.10 -0.45 -0.08
## log.m_to_coast -0.12 -0.01 -0.20 -0.33
## log.m_to_comm -0.48 0.01 -0.49 -0.17
## ll_a1_s00500 -0.06 0.64 0.01 0.17
## ll_a3_s00500 1.00 0.01 0.80 0.23
## ll_a1_s01000 0.01 1.00 0.20 0.15
## ll_a3_s01000 0.80 0.20 1.00 0.23
## pop_s01000 0.23 0.15 0.23 1.00
## pop_s01500 0.26 0.21 0.31 0.95
## pop_s02000 0.32 0.25 0.43 0.88
## pop_s01500 pop_s02000
## pm25 0.41 0.46
## longitude -0.47 -0.51
## latitude -0.38 -0.38
## log.m_to_a1 -0.39 -0.41
## log.m_to_a2 -0.05 -0.03
## log.m_to_a3 -0.05 -0.07
## log.m_to_coast -0.34 -0.36
## log.m_to_comm -0.21 -0.24
## ll_a1_s00500 0.21 0.21
## ll_a3_s00500 0.26 0.32
## ll_a1_s01000 0.21 0.25
## ll_a3_s01000 0.31 0.43
## pop_s01000 0.95 0.88
## pop_s01500 1.00 0.98
## pop_s02000 0.98 1.00upper<-mcor
upper[upper.tri(mcor)]<-""
upper<-as.data.frame(upper)
upper## pm25 longitude latitude log.m_to_a1 log.m_to_a2
## pm25 1
## longitude -0.11 1
## latitude -0.16 0.4 1
## log.m_to_a1 -0.69 0.06 0.12 1
## log.m_to_a2 -0.28 -0.1 0.17 -0.14 1
## log.m_to_a3 -0.26 0.11 -0.04 -0.09 0.15
## log.m_to_coast 0.14 0.84 0.5 -0.11 -0.2
## log.m_to_comm -0.39 0 0.12 0.07 0.09
## ll_a1_s00500 0.34 0.07 -0.06 -0.7 0.09
## ll_a3_s00500 0.69 -0.25 -0.11 -0.21 0.04
## ll_a1_s01000 0.36 -0.08 -0.18 -0.72 0.15
## ll_a3_s01000 0.63 -0.3 -0.2 -0.28 0.07
## pop_s01000 0.3 -0.42 -0.42 -0.31 -0.02
## pop_s01500 0.41 -0.47 -0.38 -0.39 -0.05
## pop_s02000 0.46 -0.51 -0.38 -0.41 -0.03
## log.m_to_a3 log.m_to_coast log.m_to_comm ll_a1_s00500
## pm25
## longitude
## latitude
## log.m_to_a1
## log.m_to_a2
## log.m_to_a3 1
## log.m_to_coast 0.01 1
## log.m_to_comm 0.25 0 1
## ll_a1_s00500 0.13 0.12 0.09 1
## ll_a3_s00500 -0.69 -0.12 -0.48 -0.06
## ll_a1_s01000 0.1 -0.01 0.01 0.64
## ll_a3_s01000 -0.45 -0.2 -0.49 0.01
## pop_s01000 -0.08 -0.33 -0.17 0.17
## pop_s01500 -0.05 -0.34 -0.21 0.21
## pop_s02000 -0.07 -0.36 -0.24 0.21
## ll_a3_s00500 ll_a1_s01000 ll_a3_s01000 pop_s01000
## pm25
## longitude
## latitude
## log.m_to_a1
## log.m_to_a2
## log.m_to_a3
## log.m_to_coast
## log.m_to_comm
## ll_a1_s00500
## ll_a3_s00500 1
## ll_a1_s01000 0.01 1
## ll_a3_s01000 0.8 0.2 1
## pop_s01000 0.23 0.15 0.23 1
## pop_s01500 0.26 0.21 0.31 0.95
## pop_s02000 0.32 0.25 0.43 0.88
## pop_s01500 pop_s02000
## pm25
## longitude
## latitude
## log.m_to_a1
## log.m_to_a2
## log.m_to_a3
## log.m_to_coast
## log.m_to_comm
## ll_a1_s00500
## ll_a3_s00500
## ll_a1_s01000
## ll_a3_s01000
## pop_s01000
## pop_s01500 1
## pop_s02000 0.98 1### p value of correlation ( install.packages(Hmisc))
library("Hmisc", lib.loc="~/R/win-library/3.4")## Warning: package 'Hmisc' was built under R version 3.4.4## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Warning: package 'Formula' was built under R version 3.4.4## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:xtable':
##
## label, label<-## The following objects are masked from 'package:base':
##
## format.pval, unitsmycor <- rcorr(as.matrix(rsubset), type="pearson")
mycor$P## pm25 longitude latitude log.m_to_a1
## pm25 NA 4.048691e-01 0.2337257563 4.444893e-09
## longitude 4.048691e-01 NA 0.0025702353 6.423193e-01
## latitude 2.337258e-01 2.570235e-03 NA 3.672013e-01
## log.m_to_a1 4.444893e-09 6.423193e-01 0.3672013496 NA
## log.m_to_a2 3.699987e-02 4.457039e-01 0.2144774121 3.124008e-01
## log.m_to_a3 5.772698e-02 4.445118e-01 0.7678230071 5.131854e-01
## log.m_to_coast 2.964470e-01 6.661338e-16 0.0001135217 4.154258e-01
## log.m_to_comm 3.316106e-03 9.952432e-01 0.3717314298 6.335699e-01
## ll_a1_s00500 1.040339e-02 6.275054e-01 0.6434151003 3.558073e-09
## ll_a3_s00500 7.015851e-09 7.045431e-02 0.4275277935 1.215331e-01
## ll_a1_s01000 6.890866e-03 5.399584e-01 0.1969403129 6.733327e-10
## ll_a3_s01000 3.380541e-07 2.573284e-02 0.1369776351 3.508624e-02
## pop_s01000 2.682471e-02 1.402588e-03 0.0016138791 2.046790e-02
## pop_s01500 2.023946e-03 3.099413e-04 0.0043338831 3.066908e-03
## pop_s02000 4.642234e-04 5.834682e-05 0.0039525939 1.735563e-03
## log.m_to_a2 log.m_to_a3 log.m_to_coast log.m_to_comm
## pm25 0.03699987 5.772698e-02 2.964470e-01 0.0033161064
## longitude 0.44570390 4.445118e-01 6.661338e-16 0.9952431519
## latitude 0.21447741 7.678230e-01 1.135217e-04 0.3717314298
## log.m_to_a1 0.31240084 5.131854e-01 4.154258e-01 0.6335699061
## log.m_to_a2 NA 2.766970e-01 1.358197e-01 0.5023132389
## log.m_to_a3 0.27669698 NA 9.510113e-01 0.0685366196
## log.m_to_coast 0.13581973 9.510113e-01 NA 0.9924473367
## log.m_to_comm 0.50231324 6.853662e-02 9.924473e-01 NA
## ll_a1_s00500 0.50959961 3.603580e-01 3.925622e-01 0.5211026106
## ll_a3_s00500 0.79098884 5.831522e-09 3.961423e-01 0.0002475775
## ll_a1_s01000 0.26273816 4.879052e-01 9.496653e-01 0.9228166471
## ll_a3_s01000 0.63498939 5.206378e-04 1.429183e-01 0.0001421568
## pop_s01000 0.91080716 5.373354e-01 1.426281e-02 0.2175033268
## pop_s01500 0.69849811 7.402650e-01 1.160186e-02 0.1317993010
## pop_s02000 0.81498049 6.204343e-01 6.650400e-03 0.0775349853
## ll_a1_s00500 ll_a3_s00500 ll_a1_s01000 ll_a3_s01000
## pm25 1.040339e-02 7.015851e-09 6.890866e-03 3.380541e-07
## longitude 6.275054e-01 7.045431e-02 5.399584e-01 2.573284e-02
## latitude 6.434151e-01 4.275278e-01 1.969403e-01 1.369776e-01
## log.m_to_a1 3.558073e-09 1.215331e-01 6.733327e-10 3.508624e-02
## log.m_to_a2 5.095996e-01 7.909888e-01 2.627382e-01 6.349894e-01
## log.m_to_a3 3.603580e-01 5.831522e-09 4.879052e-01 5.206378e-04
## log.m_to_coast 3.925622e-01 3.961423e-01 9.496653e-01 1.429183e-01
## log.m_to_comm 5.211026e-01 2.475775e-04 9.228166e-01 1.421568e-04
## ll_a1_s00500 NA 6.849689e-01 1.889517e-07 9.698146e-01
## ll_a3_s00500 6.849689e-01 NA 9.654566e-01 1.605382e-13
## ll_a1_s01000 1.889517e-07 9.654566e-01 NA 1.516042e-01
## ll_a3_s01000 9.698146e-01 1.605382e-13 1.516042e-01 NA
## pop_s01000 2.202474e-01 9.292586e-02 2.729037e-01 9.741389e-02
## pop_s01500 1.162690e-01 5.709074e-02 1.231923e-01 2.108048e-02
## pop_s02000 1.303907e-01 1.746253e-02 6.224999e-02 1.023506e-03
## pop_s01000 pop_s01500 pop_s02000
## pm25 0.026824709 0.0020239459 4.642234e-04
## longitude 0.001402588 0.0003099413 5.834682e-05
## latitude 0.001613879 0.0043338831 3.952594e-03
## log.m_to_a1 0.020467900 0.0030669082 1.735563e-03
## log.m_to_a2 0.910807155 0.6984981082 8.149805e-01
## log.m_to_a3 0.537335355 0.7402649953 6.204343e-01
## log.m_to_coast 0.014262810 0.0116018612 6.650400e-03
## log.m_to_comm 0.217503327 0.1317993010 7.753499e-02
## ll_a1_s00500 0.220247417 0.1162690352 1.303907e-01
## ll_a3_s00500 0.092925860 0.0570907415 1.746253e-02
## ll_a1_s01000 0.272903715 0.1231922866 6.224999e-02
## ll_a3_s01000 0.097413886 0.0210804842 1.023506e-03
## pop_s01000 NA 0.0000000000 0.000000e+00
## pop_s01500 0.000000000 NA 0.000000e+00
## pop_s02000 0.000000000 0.0000000000 NAVIF
### VIF
m1 <- lm(pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 + log.m_to_coast + log.m_to_comm
+ ll_a1_s00500 + ll_a3_s00500 + ll_a1_s01000 + ll_a3_s01000 + pop_s01000 + pop_s01500 + pop_s02000, data=dat)
### install.packages(car)
library(car)## Warning: package 'car' was built under R version 3.4.4## Loading required package: carData## Warning: package 'carData' was built under R version 3.4.4vif(m1)## longitude latitude log.m_to_a1 log.m_to_a2 log.m_to_a3
## 5.441836 2.194325 4.490411 1.753038 2.881996
## log.m_to_coast log.m_to_comm ll_a1_s00500 ll_a3_s00500 ll_a1_s01000
## 6.067637 1.529404 2.412744 6.364412 2.826847
## ll_a3_s01000 pop_s01000 pop_s01500 pop_s02000
## 5.573764 25.665293 142.350519 74.295539Based on the result for correlation and VIF, we found that 3 variables pop_s02000, pop_s01500, and pop_s01000 have strong correlation with each others and VIF>10. Additionally, pop_s01000 variable are showing low correaltion with outcome (PM2.5). So we can consider to exclude this variable.
Q3. Comparision model
m2 <- lm(pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 + log.m_to_coast + log.m_to_comm
+ ll_a1_s00500 + ll_a3_s00500 + ll_a1_s01000 + ll_a3_s01000 + pop_s01500 + pop_s02000, data=dat)
a=step(m2, direction="backward", test="F") ## Start: AIC=-109.72
## pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + log.m_to_comm + ll_a1_s00500 + ll_a3_s00500 +
## ll_a1_s01000 + ll_a3_s01000 + pop_s01500 + pop_s02000
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - ll_a3_s01000 1 0.011 4.507 -111.590 0.0996 0.7538693
## - ll_a1_s01000 1 0.013 4.509 -111.565 0.1179 0.7330310
## - ll_a1_s00500 1 0.065 4.561 -110.935 0.5915 0.4462359
## - pop_s02000 1 0.071 4.568 -110.858 0.6504 0.4246327
## - pop_s01500 1 0.075 4.571 -110.816 0.6818 0.4137368
## - log.m_to_comm 1 0.166 4.663 -109.724 1.5175 0.2250204
## <none> 4.496 -109.723
## - latitude 1 0.376 4.872 -107.306 3.4288 0.0712781 .
## - longitude 1 1.482 5.978 -96.055 13.5140 0.0006804 ***
## - log.m_to_coast 1 2.745 7.241 -85.513 25.0310 1.11e-05 ***
## - log.m_to_a3 1 19.551 24.047 -19.502 178.2741 < 2.2e-16 ***
## - log.m_to_a1 1 27.865 32.361 -3.171 254.0821 < 2.2e-16 ***
## - ll_a3_s00500 1 50.336 54.832 25.832 458.9815 < 2.2e-16 ***
## - log.m_to_a2 1 57.133 61.629 32.259 520.9608 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-111.59
## pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + log.m_to_comm + ll_a1_s00500 + ll_a3_s00500 +
## ll_a1_s01000 + pop_s01500 + pop_s02000
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - ll_a1_s01000 1 0.009 4.517 -113.475 0.0881 0.7681326
## - pop_s02000 1 0.064 4.571 -112.818 0.5933 0.4454727
## - pop_s01500 1 0.067 4.574 -112.781 0.6223 0.4346140
## - ll_a1_s00500 1 0.070 4.578 -112.739 0.6544 0.4230901
## - log.m_to_comm 1 0.157 4.664 -111.710 1.4601 0.2336694
## <none> 4.507 -111.590
## - latitude 1 0.380 4.888 -109.134 3.5441 0.0666962 .
## - longitude 1 1.550 6.057 -97.334 14.4426 0.0004603 ***
## - log.m_to_coast 1 2.843 7.351 -86.690 26.4951 6.59e-06 ***
## - log.m_to_a3 1 19.610 24.117 -21.342 182.7291 < 2.2e-16 ***
## - log.m_to_a1 1 27.854 32.361 -5.170 259.5498 < 2.2e-16 ***
## - log.m_to_a2 1 57.125 61.632 30.262 532.3006 < 2.2e-16 ***
## - ll_a3_s00500 1 81.103 85.611 48.336 755.7368 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-113.47
## pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + log.m_to_comm + ll_a1_s00500 + ll_a3_s00500 +
## pop_s01500 + pop_s02000
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - pop_s02000 1 0.082 4.599 -114.481 0.7835 0.3810099
## - ll_a1_s00500 1 0.085 4.602 -114.444 0.8133 0.3721622
## - pop_s01500 1 0.090 4.607 -114.392 0.8549 0.3603226
## - log.m_to_comm 1 0.163 4.680 -113.523 1.5531 0.2194296
## <none> 4.517 -113.475
## - latitude 1 0.371 4.888 -111.134 3.5311 0.0670158 .
## - longitude 1 1.547 6.063 -99.278 14.7240 0.0004027 ***
## - log.m_to_coast 1 2.834 7.351 -88.689 26.9798 5.346e-06 ***
## - log.m_to_a3 1 20.518 25.034 -21.290 195.3284 < 2.2e-16 ***
## - log.m_to_a1 1 39.316 43.832 9.517 374.2871 < 2.2e-16 ***
## - log.m_to_a2 1 57.651 62.168 28.738 548.8422 < 2.2e-16 ***
## - ll_a3_s00500 1 88.978 93.495 51.181 847.0781 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-114.48
## pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + log.m_to_comm + ll_a1_s00500 + ll_a3_s00500 +
## pop_s01500
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - pop_s01500 1 0.008 4.607 -116.391 0.0721 0.7895570
## - ll_a1_s00500 1 0.084 4.683 -115.486 0.8039 0.3748235
## <none> 4.599 -114.481
## - log.m_to_comm 1 0.195 4.794 -114.198 1.8649 0.1790013
## - latitude 1 0.366 4.965 -112.272 3.4995 0.0680464 .
## - longitude 1 1.861 6.460 -97.795 17.8019 0.0001205 ***
## - log.m_to_coast 1 2.973 7.572 -89.055 28.4470 3.182e-06 ***
## - log.m_to_a3 1 21.098 25.697 -21.853 201.8452 < 2.2e-16 ***
## - log.m_to_a1 1 39.798 44.397 8.221 380.7510 < 2.2e-16 ***
## - log.m_to_a2 1 57.576 62.175 26.744 550.8348 < 2.2e-16 ***
## - ll_a3_s00500 1 93.269 97.868 51.696 892.3166 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-116.39
## pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + log.m_to_comm + ll_a1_s00500 + ll_a3_s00500
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - ll_a1_s00500 1 0.083 4.689 -117.412 0.8082 0.37344
## <none> 4.607 -116.391
## - log.m_to_comm 1 0.188 4.794 -116.197 1.8321 0.18264
## - latitude 1 0.389 4.996 -113.930 3.8025 0.05742 .
## - longitude 1 1.932 6.539 -99.129 18.8729 7.861e-05 ***
## - log.m_to_coast 1 2.967 7.574 -91.045 28.9862 2.538e-06 ***
## - log.m_to_a3 1 21.094 25.700 -23.846 206.0542 < 2.2e-16 ***
## - log.m_to_a1 1 41.768 46.375 8.619 408.0192 < 2.2e-16 ***
## - log.m_to_a2 1 58.174 62.781 25.278 568.2821 < 2.2e-16 ***
## - ll_a3_s00500 1 93.279 97.885 49.705 911.2010 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-117.41
## pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + log.m_to_comm + ll_a3_s00500
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 4.689 -117.412
## - log.m_to_comm 1 0.174 4.863 -117.409 1.7066 0.19792
## - latitude 1 0.391 5.080 -115.008 3.8352 0.05627 .
## - longitude 1 1.851 6.540 -101.115 18.1567 9.968e-05 ***
## - log.m_to_coast 1 2.886 7.576 -93.032 28.3126 2.972e-06 ***
## - log.m_to_a3 1 21.342 26.031 -25.142 209.3528 < 2.2e-16 ***
## - log.m_to_a2 1 58.223 62.912 23.393 571.1391 < 2.2e-16 ***
## - log.m_to_a1 1 92.210 96.899 47.149 904.5349 < 2.2e-16 ***
## - ll_a3_s00500 1 98.163 102.852 50.428 962.9286 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(a)##
## Call:
## lm(formula = pm25 ~ longitude + latitude + log.m_to_a1 + log.m_to_a2 +
## log.m_to_a3 + log.m_to_coast + log.m_to_comm + ll_a3_s00500,
## data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.85217 -0.17676 0.08793 0.20640 0.52089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.451e+02 9.567e+01 -4.652 2.79e-05 ***
## longitude -3.401e+00 7.981e-01 -4.261 9.97e-05 ***
## latitude 1.794e+00 9.159e-01 1.958 0.0563 .
## log.m_to_a1 -4.051e+00 1.347e-01 -30.075 < 2e-16 ***
## log.m_to_a2 -2.320e+00 9.710e-02 -23.899 < 2e-16 ***
## log.m_to_a3 2.348e+00 1.623e-01 14.469 < 2e-16 ***
## log.m_to_coast 3.529e+00 6.632e-01 5.321 2.97e-06 ***
## log.m_to_comm -1.182e-01 9.047e-02 -1.306 0.1979
## ll_a3_s00500 2.032e-03 6.550e-05 31.031 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3193 on 46 degrees of freedom
## Multiple R-squared: 0.9902, Adjusted R-squared: 0.9885
## F-statistic: 580.6 on 8 and 46 DF, p-value: < 2.2e-16rq2=subset(dat, select=c("pm25" ,"longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000" , "pop_s01500" , "pop_s02000"))
m0 <- lm(pm25 ~ 1, data=rq2)
b= step(m0, direction="forward", scope=list(lower = ~1, upper = m2), test="F")## Start: AIC=120.95
## pm25 ~ 1
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_a1 1 229.862 248.35 86.912 49.0551 4.445e-09 ***
## + ll_a3_s00500 1 225.630 252.58 87.841 47.3451 7.016e-09 ***
## + ll_a3_s01000 1 186.819 291.39 95.703 33.9799 3.381e-07 ***
## + pop_s02000 1 99.530 378.68 110.115 13.9302 0.0004642 ***
## + pop_s01500 1 79.352 398.86 112.970 10.5443 0.0020239 **
## + log.m_to_comm 1 72.436 405.77 113.915 9.4612 0.0033161 **
## + ll_a1_s01000 1 62.078 416.13 115.302 7.9065 0.0068909 **
## + ll_a1_s00500 1 56.199 422.01 116.073 7.0581 0.0104034 *
## + log.m_to_a2 1 38.025 440.18 118.392 4.5784 0.0369999 *
## + log.m_to_a3 1 31.702 446.51 119.177 3.7630 0.0577270 .
## <none> 478.21 120.949
## + latitude 1 12.743 465.47 121.464 1.4510 0.2337258
## + log.m_to_coast 1 9.827 468.38 121.807 1.1119 0.2964470
## + longitude 1 6.278 471.93 122.222 0.7050 0.4048691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=86.91
## pm25 ~ log.m_to_a1
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + ll_a3_s00500 1 146.187 102.16 40.056 74.4101 1.322e-11 ***
## + ll_a3_s01000 1 95.149 153.20 62.342 32.2962 6.069e-07 ***
## + log.m_to_a2 1 69.734 178.61 70.784 20.3020 3.777e-05 ***
## + log.m_to_comm 1 56.713 191.63 74.654 15.3892 0.0002578 ***
## + log.m_to_a3 1 49.342 199.00 76.730 12.8932 0.0007306 ***
## + ll_a1_s01000 1 18.702 229.64 84.606 4.2349 0.0446289 *
## + ll_a1_s00500 1 18.194 230.15 84.728 4.1107 0.0477519 *
## + pop_s02000 1 16.661 231.69 85.093 3.7395 0.0585905 .
## + pop_s01500 1 10.373 237.97 86.566 2.2666 0.1382427
## <none> 248.35 86.912
## + latitude 1 2.901 245.44 88.266 0.6147 0.4365777
## + longitude 1 2.365 245.98 88.386 0.5000 0.4826660
## + log.m_to_coast 1 2.089 246.26 88.448 0.4410 0.5095612
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=40.06
## pm25 ~ log.m_to_a1 + ll_a3_s00500
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_a2 1 71.260 30.900 -23.711 117.6120 7.553e-15 ***
## + log.m_to_coast 1 10.410 91.749 36.145 5.7867 0.01981 *
## + log.m_to_a3 1 6.707 95.452 38.321 3.5838 0.06403 .
## + log.m_to_comm 1 4.238 97.922 39.726 2.2071 0.14353
## <none> 102.160 40.056
## + ll_a1_s01000 1 3.150 99.009 40.334 1.6227 0.20848
## + longitude 1 1.907 100.252 41.020 0.9703 0.32925
## + pop_s02000 1 0.925 101.234 41.556 0.4662 0.49782
## + pop_s01500 1 0.776 101.384 41.637 0.3903 0.53495
## + ll_a1_s00500 1 0.639 101.521 41.711 0.3212 0.57339
## + latitude 1 0.451 101.709 41.813 0.2260 0.63653
## + ll_a3_s01000 1 0.074 102.085 42.016 0.0371 0.84793
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-23.71
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_a3 1 19.9277 10.973 -78.656 90.8068 7.938e-13 ***
## + log.m_to_coast 1 1.8770 29.023 -25.158 3.2337 0.07818 .
## + ll_a1_s01000 1 1.1825 29.718 -23.857 1.9896 0.16458
## <none> 30.900 -23.711
## + log.m_to_comm 1 1.0987 29.802 -23.703 1.8434 0.18065
## + latitude 1 0.9344 29.966 -23.400 1.5592 0.21760
## + ll_a1_s00500 1 0.7265 30.174 -23.020 1.2038 0.27782
## + ll_a3_s01000 1 0.3241 30.576 -22.291 0.5299 0.47004
## + longitude 1 0.3073 30.593 -22.261 0.5023 0.48180
## + pop_s01500 1 0.0230 30.877 -21.752 0.0372 0.84780
## + pop_s02000 1 0.0059 30.894 -21.722 0.0096 0.92224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-78.66
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_coast 1 3.6092 7.3634 -98.595 24.0174 1.086e-05 ***
## + latitude 1 3.2857 7.6869 -96.230 20.9445 3.249e-05 ***
## + longitude 1 0.7863 10.1862 -80.746 3.7826 0.05754 .
## + pop_s01500 1 0.4082 10.5644 -78.741 1.8933 0.17509
## <none> 10.9726 -78.656
## + pop_s02000 1 0.2910 10.6815 -78.135 1.3351 0.25351
## + log.m_to_comm 1 0.0659 10.9067 -76.988 0.2961 0.58882
## + ll_a1_s01000 1 0.0628 10.9098 -76.972 0.2819 0.59786
## + ll_a3_s01000 1 0.0549 10.9177 -76.932 0.2464 0.62181
## + ll_a1_s00500 1 0.0000 10.9725 -76.657 0.0002 0.98997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-98.59
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + longitude 1 2.13579 5.2276 -115.436 19.6109 5.47e-05 ***
## + latitude 1 0.77893 6.5845 -102.744 5.6783 0.02118 *
## <none> 7.3634 -98.595
## + pop_s02000 1 0.14546 7.2179 -97.692 0.9673 0.33028
## + pop_s01500 1 0.04460 7.3188 -96.929 0.2925 0.59111
## + ll_a3_s01000 1 0.01944 7.3440 -96.740 0.1271 0.72304
## + log.m_to_comm 1 0.01791 7.3455 -96.729 0.1171 0.73374
## + ll_a1_s01000 1 0.00781 7.3556 -96.653 0.0510 0.82230
## + ll_a1_s00500 1 0.00012 7.3633 -96.596 0.0008 0.97810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-115.44
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + longitude
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + latitude 1 0.36429 4.8633 -117.41 3.5206 0.06683 .
## <none> 5.2276 -115.44
## + log.m_to_comm 1 0.14729 5.0803 -115.01 1.3627 0.24896
## + ll_a1_s00500 1 0.07169 5.1559 -114.19 0.6535 0.42295
## + ll_a1_s01000 1 0.03235 5.1952 -113.78 0.2927 0.59106
## + ll_a3_s01000 1 0.02005 5.2076 -113.65 0.1809 0.67251
## + pop_s01500 1 0.00991 5.2177 -113.54 0.0893 0.76643
## + pop_s02000 1 0.00010 5.2275 -113.44 0.0009 0.97600
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-117.41
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + longitude + latitude
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_comm 1 0.173976 4.6893 -117.41 1.7066 0.1979
## <none> 4.8633 -117.41
## + ll_a1_s00500 1 0.069159 4.7942 -116.20 0.6636 0.4195
## + ll_a1_s01000 1 0.066168 4.7971 -116.16 0.6345 0.4298
## + ll_a3_s01000 1 0.032363 4.8309 -115.78 0.3082 0.5815
## + pop_s02000 1 0.006697 4.8566 -115.48 0.0634 0.8023
## + pop_s01500 1 0.000069 4.8632 -115.41 0.0006 0.9798
##
## Step: AIC=-117.41
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + longitude + latitude + log.m_to_comm
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 4.6893 -117.41
## + ll_a1_s00500 1 0.082735 4.6066 -116.39 0.8082 0.3734
## + ll_a1_s01000 1 0.054209 4.6351 -116.05 0.5263 0.4719
## + ll_a3_s01000 1 0.008237 4.6811 -115.51 0.0792 0.7797
## + pop_s01500 1 0.006249 4.6831 -115.49 0.0600 0.8075
## + pop_s02000 1 0.000002 4.6893 -115.41 0.0000 0.9965summary(b)##
## Call:
## lm(formula = pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 +
## log.m_to_a3 + log.m_to_coast + longitude + latitude + log.m_to_comm,
## data = rq2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.85217 -0.17676 0.08793 0.20640 0.52089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.451e+02 9.567e+01 -4.652 2.79e-05 ***
## log.m_to_a1 -4.051e+00 1.347e-01 -30.075 < 2e-16 ***
## ll_a3_s00500 2.032e-03 6.550e-05 31.031 < 2e-16 ***
## log.m_to_a2 -2.320e+00 9.710e-02 -23.899 < 2e-16 ***
## log.m_to_a3 2.348e+00 1.623e-01 14.469 < 2e-16 ***
## log.m_to_coast 3.529e+00 6.632e-01 5.321 2.97e-06 ***
## longitude -3.401e+00 7.981e-01 -4.261 9.97e-05 ***
## latitude 1.794e+00 9.159e-01 1.958 0.0563 .
## log.m_to_comm -1.182e-01 9.047e-02 -1.306 0.1979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3193 on 46 degrees of freedom
## Multiple R-squared: 0.9902, Adjusted R-squared: 0.9885
## F-statistic: 580.6 on 8 and 46 DF, p-value: < 2.2e-16c= step(m0, direction="both", scope=list(lower = ~1, upper = m2), test="F")## Start: AIC=120.95
## pm25 ~ 1
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_a1 1 229.862 248.35 86.912 49.0551 4.445e-09 ***
## + ll_a3_s00500 1 225.630 252.58 87.841 47.3451 7.016e-09 ***
## + ll_a3_s01000 1 186.819 291.39 95.703 33.9799 3.381e-07 ***
## + pop_s02000 1 99.530 378.68 110.115 13.9302 0.0004642 ***
## + pop_s01500 1 79.352 398.86 112.970 10.5443 0.0020239 **
## + log.m_to_comm 1 72.436 405.77 113.915 9.4612 0.0033161 **
## + ll_a1_s01000 1 62.078 416.13 115.302 7.9065 0.0068909 **
## + ll_a1_s00500 1 56.199 422.01 116.073 7.0581 0.0104034 *
## + log.m_to_a2 1 38.025 440.18 118.392 4.5784 0.0369999 *
## + log.m_to_a3 1 31.702 446.51 119.177 3.7630 0.0577270 .
## <none> 478.21 120.949
## + latitude 1 12.743 465.47 121.464 1.4510 0.2337258
## + log.m_to_coast 1 9.827 468.38 121.807 1.1119 0.2964470
## + longitude 1 6.278 471.93 122.222 0.7050 0.4048691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=86.91
## pm25 ~ log.m_to_a1
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + ll_a3_s00500 1 146.187 102.16 40.056 74.4101 1.322e-11 ***
## + ll_a3_s01000 1 95.149 153.20 62.342 32.2962 6.069e-07 ***
## + log.m_to_a2 1 69.734 178.61 70.784 20.3020 3.777e-05 ***
## + log.m_to_comm 1 56.713 191.63 74.654 15.3892 0.0002578 ***
## + log.m_to_a3 1 49.342 199.00 76.730 12.8932 0.0007306 ***
## + ll_a1_s01000 1 18.702 229.64 84.606 4.2349 0.0446289 *
## + ll_a1_s00500 1 18.194 230.15 84.728 4.1107 0.0477519 *
## + pop_s02000 1 16.661 231.69 85.093 3.7395 0.0585905 .
## + pop_s01500 1 10.373 237.97 86.566 2.2666 0.1382427
## <none> 248.35 86.912
## + latitude 1 2.901 245.45 88.266 0.6147 0.4365777
## + longitude 1 2.365 245.98 88.386 0.5000 0.4826660
## + log.m_to_coast 1 2.089 246.26 88.448 0.4410 0.5095612
## - log.m_to_a1 1 229.862 478.21 120.949 49.0551 4.445e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=40.06
## pm25 ~ log.m_to_a1 + ll_a3_s00500
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_a2 1 71.260 30.900 -23.711 117.6120 7.553e-15 ***
## + log.m_to_coast 1 10.410 91.749 36.145 5.7867 0.01981 *
## + log.m_to_a3 1 6.707 95.452 38.321 3.5838 0.06403 .
## + log.m_to_comm 1 4.238 97.922 39.726 2.2071 0.14353
## <none> 102.160 40.056
## + ll_a1_s01000 1 3.150 99.009 40.334 1.6227 0.20848
## + longitude 1 1.907 100.252 41.020 0.9703 0.32925
## + pop_s02000 1 0.925 101.234 41.556 0.4662 0.49782
## + pop_s01500 1 0.776 101.384 41.637 0.3903 0.53495
## + ll_a1_s00500 1 0.639 101.521 41.711 0.3212 0.57339
## + latitude 1 0.451 101.709 41.813 0.2260 0.63653
## + ll_a3_s01000 1 0.074 102.085 42.016 0.0371 0.84793
## - ll_a3_s00500 1 146.187 248.347 86.912 74.4101 1.322e-11 ***
## - log.m_to_a1 1 150.419 252.579 87.841 76.5642 8.473e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-23.71
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_a3 1 19.928 10.973 -78.656 90.8068 7.938e-13 ***
## + log.m_to_coast 1 1.877 29.023 -25.158 3.2337 0.07818 .
## + ll_a1_s01000 1 1.183 29.718 -23.857 1.9896 0.16458
## <none> 30.900 -23.711
## + log.m_to_comm 1 1.099 29.802 -23.703 1.8434 0.18065
## + latitude 1 0.934 29.966 -23.400 1.5592 0.21760
## + ll_a1_s00500 1 0.726 30.174 -23.020 1.2038 0.27782
## + ll_a3_s01000 1 0.324 30.576 -22.291 0.5299 0.47004
## + longitude 1 0.307 30.593 -22.261 0.5023 0.48180
## + pop_s01500 1 0.023 30.877 -21.752 0.0372 0.84780
## + pop_s02000 1 0.006 30.894 -21.722 0.0096 0.92224
## - log.m_to_a2 1 71.260 102.160 40.056 117.6120 7.553e-15 ***
## - ll_a3_s00500 1 147.712 178.612 70.784 243.7947 < 2.2e-16 ***
## - log.m_to_a1 1 176.518 207.418 79.007 291.3374 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-78.66
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_coast 1 3.609 7.363 -98.595 24.0174 1.086e-05 ***
## + latitude 1 3.286 7.687 -96.230 20.9445 3.249e-05 ***
## + longitude 1 0.786 10.186 -80.746 3.7826 0.05754 .
## + pop_s01500 1 0.408 10.564 -78.741 1.8933 0.17509
## <none> 10.973 -78.656
## + pop_s02000 1 0.291 10.682 -78.135 1.3351 0.25351
## + log.m_to_comm 1 0.066 10.907 -76.988 0.2961 0.58882
## + ll_a1_s01000 1 0.063 10.910 -76.972 0.2819 0.59786
## + ll_a3_s01000 1 0.055 10.918 -76.932 0.2464 0.62181
## + ll_a1_s00500 1 0.000 10.973 -76.657 0.0002 0.98997
## - log.m_to_a3 1 19.928 30.900 -23.711 90.8068 7.938e-13 ***
## - log.m_to_a2 1 84.480 95.452 38.321 384.9589 < 2.2e-16 ***
## - log.m_to_a1 1 125.991 136.964 58.181 574.1176 < 2.2e-16 ***
## - ll_a3_s00500 1 132.551 143.524 60.754 604.0120 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-98.59
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + longitude 1 2.136 5.228 -115.436 19.6109 5.470e-05 ***
## + latitude 1 0.779 6.584 -102.744 5.6783 0.02118 *
## <none> 7.363 -98.595
## + pop_s02000 1 0.145 7.218 -97.692 0.9673 0.33028
## + pop_s01500 1 0.045 7.319 -96.929 0.2925 0.59111
## + ll_a3_s01000 1 0.019 7.344 -96.740 0.1271 0.72304
## + log.m_to_comm 1 0.018 7.345 -96.729 0.1171 0.73374
## + ll_a1_s01000 1 0.008 7.356 -96.653 0.0510 0.82230
## + ll_a1_s00500 1 0.000 7.363 -96.596 0.0008 0.97810
## - log.m_to_coast 1 3.609 10.973 -78.656 24.0174 1.086e-05 ***
## - log.m_to_a3 1 21.660 29.023 -25.158 144.1361 3.323e-16 ***
## - log.m_to_a2 1 74.546 81.909 31.905 496.0662 < 2.2e-16 ***
## - log.m_to_a1 1 112.729 120.092 52.951 750.1570 < 2.2e-16 ***
## - ll_a3_s00500 1 136.100 143.463 62.731 905.6815 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-115.44
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + longitude
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + latitude 1 0.364 4.863 -117.409 3.5206 0.06683 .
## <none> 5.228 -115.436
## + log.m_to_comm 1 0.147 5.080 -115.008 1.3627 0.24896
## + ll_a1_s00500 1 0.072 5.156 -114.195 0.6535 0.42295
## + ll_a1_s01000 1 0.032 5.195 -113.777 0.2927 0.59106
## + ll_a3_s01000 1 0.020 5.208 -113.647 0.1809 0.67251
## + pop_s01500 1 0.010 5.218 -113.540 0.0893 0.76643
## + pop_s02000 1 0.000 5.227 -113.437 0.0009 0.97600
## - longitude 1 2.136 7.363 -98.595 19.6109 5.470e-05 ***
## - log.m_to_coast 1 4.959 10.186 -80.746 45.5304 1.793e-08 ***
## - log.m_to_a3 1 22.299 27.526 -26.071 204.7458 < 2.2e-16 ***
## - log.m_to_a2 1 68.013 73.241 27.753 624.4985 < 2.2e-16 ***
## - log.m_to_a1 1 95.738 100.965 45.409 879.0674 < 2.2e-16 ***
## - ll_a3_s00500 1 130.364 135.592 61.627 1197.0099 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-117.41
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + longitude + latitude
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + log.m_to_comm 1 0.174 4.689 -117.412 1.7066 0.1979198
## <none> 4.863 -117.409
## + ll_a1_s00500 1 0.069 4.794 -116.197 0.6636 0.4194967
## + ll_a1_s01000 1 0.066 4.797 -116.162 0.6345 0.4298068
## + ll_a3_s01000 1 0.032 4.831 -115.776 0.3082 0.5815020
## + pop_s02000 1 0.007 4.857 -115.485 0.0634 0.8022805
## - latitude 1 0.364 5.228 -115.436 3.5206 0.0668280 .
## + pop_s01500 1 0.000 4.863 -115.410 0.0006 0.9797855
## - longitude 1 1.721 6.584 -102.744 16.6336 0.0001739 ***
## - log.m_to_coast 1 2.799 7.662 -94.408 27.0469 4.255e-06 ***
## - log.m_to_a3 1 22.661 27.524 -24.074 219.0015 < 2.2e-16 ***
## - log.m_to_a2 1 60.149 65.012 23.198 581.2927 < 2.2e-16 ***
## - log.m_to_a1 1 92.047 96.910 45.155 889.5616 < 2.2e-16 ***
## - ll_a3_s00500 1 129.883 134.746 63.283 1255.2155 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=-117.41
## pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 + log.m_to_a3 +
## log.m_to_coast + longitude + latitude + log.m_to_comm
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 4.689 -117.412
## - log.m_to_comm 1 0.174 4.863 -117.409 1.7066 0.19792
## + ll_a1_s00500 1 0.083 4.607 -116.391 0.8082 0.37344
## + ll_a1_s01000 1 0.054 4.635 -116.052 0.5263 0.47193
## + ll_a3_s01000 1 0.008 4.681 -115.509 0.0792 0.77970
## + pop_s01500 1 0.006 4.683 -115.486 0.0600 0.80753
## + pop_s02000 1 0.000 4.689 -115.412 0.0000 0.99654
## - latitude 1 0.391 5.080 -115.008 3.8352 0.05627 .
## - longitude 1 1.851 6.540 -101.115 18.1567 9.968e-05 ***
## - log.m_to_coast 1 2.886 7.576 -93.032 28.3126 2.972e-06 ***
## - log.m_to_a3 1 21.342 26.031 -25.142 209.3528 < 2.2e-16 ***
## - log.m_to_a2 1 58.223 62.912 23.393 571.1391 < 2.2e-16 ***
## - log.m_to_a1 1 92.210 96.899 47.149 904.5349 < 2.2e-16 ***
## - ll_a3_s00500 1 98.163 102.852 50.428 962.9286 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(c)##
## Call:
## lm(formula = pm25 ~ log.m_to_a1 + ll_a3_s00500 + log.m_to_a2 +
## log.m_to_a3 + log.m_to_coast + longitude + latitude + log.m_to_comm,
## data = rq2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.85217 -0.17676 0.08793 0.20640 0.52089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.451e+02 9.567e+01 -4.652 2.79e-05 ***
## log.m_to_a1 -4.051e+00 1.347e-01 -30.075 < 2e-16 ***
## ll_a3_s00500 2.032e-03 6.550e-05 31.031 < 2e-16 ***
## log.m_to_a2 -2.320e+00 9.710e-02 -23.899 < 2e-16 ***
## log.m_to_a3 2.348e+00 1.623e-01 14.469 < 2e-16 ***
## log.m_to_coast 3.529e+00 6.632e-01 5.321 2.97e-06 ***
## longitude -3.401e+00 7.981e-01 -4.261 9.97e-05 ***
## latitude 1.794e+00 9.159e-01 1.958 0.0563 .
## log.m_to_comm -1.182e-01 9.047e-02 -1.306 0.1979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3193 on 46 degrees of freedom
## Multiple R-squared: 0.9902, Adjusted R-squared: 0.9885
## F-statistic: 580.6 on 8 and 46 DF, p-value: < 2.2e-16Interpretation: Using Stepwise approach in the multivariate regression model, we found that Sum of distances of interstate local roads within 500 meters (ll_a1_s00500), Log distance to the nearest local roads in meter (log.m_to_a3), Log distance to the nearest local roads in meter (log.m_to_a3), Log distance to the nearest coastline in meter (log.m_to_coast), latitude are positively significantly associated with increasing concentration of PM2.5 (p<0.01). While Log distance to the nearest interstate highway in meter (log.m_to_a1), Log distance to the nearest highway in meter (log.m_to_a2), longitude, and Log distance to the nearest commertial area in meter (log.m_to_comm) are negatively associated significance with concentration of PM2.5 (p<0.01).
Conclusion: Comparing with the stepwide model, we found that both backward and forward models are producing the similar results.
Q4
library(stargazer)
cols1 <- c("pm25" ,"longitude","latitude", "log.m_to_a1" , "log.m_to_a2" , "log.m_to_a3" , "log.m_to_coast", "log.m_to_comm", "ll_a1_s00500" , "ll_a3_s00500" ,
"ll_a1_s01000" , "ll_a3_s01000", "pop_s01000" , "pop_s01500" , "pop_s02000","age1c", "gender1", "wtlb1", "htcm1", "waistcm1", "hipcm1", "sbp1c" )
stargazer(
df[, cols1], type = "text",
summary.stat = c("min", "p25", "median", "p75", "max", "mean", "sd")
)##
## =================================================================================
## Statistic Min Pctl(25) Median Pctl(75) Max Mean St. Dev.
## ---------------------------------------------------------------------------------
## pm25 16.010 18.543 20.977 22.492 30.131 20.874 2.976
## longitude -118.433 -118.187 -118.116 -118.064 -117.956 -118.143 0.114
## latitude 33.950 34.034 34.077 34.122 34.188 34.073 0.062
## log.m_to_a1 1.778 2.821 3.107 3.291 3.611 3.043 0.387
## log.m_to_a2 1.362 3.065 3.522 3.699 4.033 3.364 0.531
## log.m_to_a3 1.114 2.124 2.403 2.639 2.972 2.340 0.414
## log.m_to_coast 3.801 4.363 4.398 4.398 4.398 4.323 0.154
## log.m_to_comm 1.000 2.107 2.505 2.664 2.995 2.288 0.569
## ll_a1_s00500 10 10 10 10 3,310 277.709 750.456
## ll_a3_s00500 10 382.5 1,014 2,043.5 4,095 1,342.764 1,168.999
## ll_a1_s01000 10 10 10 2,876 9,709 1,708.982 2,693.932
## ll_a3_s01000 699 4,006 4,896 6,651.5 18,287 5,649.473 2,976.603
## pop_s01000 3,804 7,763 10,146 12,294 34,367 11,117.730 5,301.950
## pop_s01500 8,808 17,403 21,625 27,513.5 87,118 24,738.730 12,384.190
## pop_s02000 14,499 31,147.5 40,059 51,617.5 155,190 44,081.440 22,521.920
## age1c 45 51.5 57 64.5 84 58.709 9.116
## gender1 0 0 0 1 1 0.436 0.501
## wtlb1 86.200 139.700 163.000 180.200 305.000 163.967 35.389
## htcm1 144.000 156.850 167.500 173.900 189.300 166.107 11.108
## waistcm1 62.000 88.150 95.200 103.900 140.000 95.296 13.858
## hipcm1 77.000 96.650 103.500 106.650 141.400 102.835 10.829
## sbp1c 92.500 106.000 115.000 129.000 202.500 119.527 20.556
## ---------------------------------------------------------------------------------Q5
hist(dat$sbp1c, breaks=20)
dat$sb1=dat$sbp1c-mean(dat$sbp1c)
qqnorm(dat$sb1)
qqline(dat$sb1)
dat$sb1=log10(dat$sbp1c+1)
qqnorm(dat$sb1)
qqline(dat$sb1)
However, we tried difference methods of transformation (log10(sbp1c) or sbp1c-mean(sbp1c), the distribution is remaining non-distribution. So in this case, we assume SBP are normal distribution.
Q6
library("olsrr", lib.loc="~/R/win-library/3.4")## Warning: package 'olsrr' was built under R version 3.4.4##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riverslibrary("kableExtra", lib.loc="~/R/win-library/3.4")## Warning: package 'kableExtra' was built under R version 3.4.4m_6 <- lm(sbp1c ~ pm25 + age1c + gender1 + bmi1c + cig1c, data=dat)
mod3_table <- ols_step_all_possible(m_6)
kable(mod3_table[c(7,20,28,31),1:8], digits=2) %>%
kable_styling(bootstrap_options = "striped", full_width = F)| mindex | n | predictors | rsquare | adjr | predrsq | cp | aic | |
|---|---|---|---|---|---|---|---|---|
| 6 | 7 | 2 | pm25 age1c | 0.27 | 0.25 | 0.20 | 9.45 | 477.95 |
| 18 | 20 | 3 | pm25 age1c cig1c | 0.27 | 0.23 | 0.16 | 11.45 | 479.95 |
| 30 | 28 | 4 | age1c gender1 bmi1c cig1c | 0.39 | 0.34 | 0.24 | 4.56 | 472.88 |
| 31 | 31 | 5 | pm25 age1c gender1 bmi1c cig1c | 0.39 | 0.33 | 0.23 | 6.00 | 474.25 |
Based on this results with r2, AIC, we can choose the model 3 is the best model.
Q7
m3 <- lm(sbp1c ~ pm25 + age1c + gender1, data=dat)
m4 <- lm(sbp1c ~ pm25 + age1c + gender1 + bmi1c, data=dat)
summary(m3)$coef ## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.3680828 22.8385120 2.9059723 5.404109e-03
## pm25 -0.7621631 0.8380563 -0.9094413 3.673952e-01
## age1c 1.1706728 0.2693639 4.3460640 6.621456e-05
## gender1 0.7784350 4.9760660 0.1564358 8.763071e-01summary(m4)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.4331863 24.9749708 1.0183470 3.134154e-01
## pm25 -0.5843847 0.7781024 -0.7510383 4.561500e-01
## age1c 1.1932540 0.2495127 4.7823382 1.569039e-05
## gender1 1.6198713 4.6154467 0.3509674 7.270873e-01
## bmi1c 1.3173080 0.4276323 3.0804688 3.355227e-03m5 <- lm(sbp1c ~ pm25 + age1c + gender1 + bmi1c + pm25:bmi1c, data=dat)
summary(m5)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -25.80777524 90.1967113 -0.2861277 7.759866e-01
## pm25 1.92916002 4.3207157 0.4464908 6.572098e-01
## age1c 1.22856862 0.2581485 4.7591540 1.757659e-05
## gender1 1.51193902 4.6493280 0.3251952 7.464179e-01
## bmi1c 3.16077939 3.1459668 1.0047084 3.199732e-01
## pm25:bmi1c -0.09432072 0.1594488 -0.5915424 5.568765e-01IF we look at the coefficient of PM2.5 in the model 1 and model 2, there are some changes of coefficient. We can consider BMI as a confounder. Std. Error of PM2.5 in model 1 (0.84) is larger than model 2 (0.78), so we also consider BMI a precision variable. In term of interaction (PM25:BMI), with t=-0.59 and p>0.05, we can conclude that there are no interaction PM25 and BMI in the model relationship between PM2.5 and SBP after adjustment for other variables such as Age, gender, and BMI.
Q8
m6 <- lm(sbp1c ~ pm25, data=dat)
m7 <- lm(sbp1c ~ pm25 + age1c + gender1 + wtlb1 + htcm1 + waistcm1 + hipcm1 + bmi1c + cig1c, data=dat)
step(m6, direction="both", scope=list(lower = m6, upper = m7), test="F")## Start: AIC=335.22
## sbp1c ~ pm25
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + age1c 1 6136.2 16549 319.87 19.2809 5.559e-05 ***
## + waistcm1 1 2981.1 19704 329.47 7.8673 0.007063 **
## + bmi1c 1 2372.1 20313 331.14 6.0724 0.017074 *
## + hipcm1 1 2011.0 20674 332.11 5.0582 0.028772 *
## + wtlb1 1 1473.9 21211 333.52 3.6132 0.062872 .
## <none> 22685 335.22
## + htcm1 1 199.3 22486 336.73 0.4608 0.500268
## + cig1c 1 79.2 22606 337.03 0.1821 0.671338
## + gender1 1 18.0 22667 337.17 0.0412 0.839911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=319.87
## sbp1c ~ pm25 + age1c
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + bmi1c 1 2612.2 13937 312.42 9.5590 0.003222 **
## + waistcm1 1 2240.1 14309 313.87 7.9841 0.006721 **
## + hipcm1 1 1964.5 14585 314.92 6.8696 0.011525 *
## + wtlb1 1 1805.4 14744 315.52 6.2451 0.015716 *
## <none> 16549 319.87
## + htcm1 1 123.6 16426 321.46 0.3838 0.538357
## + gender1 1 7.9 16541 321.84 0.0245 0.876307
## + cig1c 1 0.1 16549 321.87 0.0004 0.984641
## - age1c 1 6136.2 22685 335.22 19.2809 5.559e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=312.42
## sbp1c ~ pm25 + age1c + bmi1c
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 13937 312.42
## + gender1 1 34.2 13903 314.29 0.1232 0.727087
## + waistcm1 1 19.8 13917 314.34 0.0712 0.790686
## + hipcm1 1 16.9 13920 314.36 0.0608 0.806290
## + cig1c 1 13.4 13924 314.37 0.0482 0.827181
## + wtlb1 1 4.7 13932 314.40 0.0167 0.897636
## + htcm1 1 4.0 13933 314.41 0.0143 0.905333
## - bmi1c 1 2612.2 16549 319.87 9.5590 0.003222 **
## - age1c 1 6376.3 20313 331.14 23.3331 1.286e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## lm(formula = sbp1c ~ pm25 + age1c + bmi1c, data = dat)
##
## Coefficients:
## (Intercept) pm25 age1c bmi1c
## 25.3227 -0.5378 1.1947 1.3084Interpretation: 0.54 mmhg decrease of SBP corresponding to one micrometers increase of PM2.5 after adjusting for Age and BMI. But this decrease are not significant (p>0.05).
Q9
res <- m7$resid
lev <- hatvalues(m7)
r.jack <- rstudent(m7)
cook <- cooks.distance(m7)
par(mfrow=c(2,2))
plot(res, ylab="residual")
abline(h=0, lty=3)
plot(r.jack, ylab="jackknife residual")
abline(h=0, lty=3)
plot(cook, ylab="cook's distance")
abline(h=1, lty=3)
library(car)
(outs <- influencePlot(m7))
## StudRes Hat CookD
## 16 -1.938621 0.1433459 0.05925566
## 17 -1.094830 0.4772912 0.10896933
## 19 5.429844 0.1209000 0.24830597
## 31 1.034266 0.5508492 0.13098845# Q10
summary(m7)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 252.1691578 186.7627667 1.35021108 0.183700489
## pm25 -0.7763883 0.8288655 -0.93668786 0.353919189
## age1c 1.1532359 0.2805483 4.11065007 0.000165031
## gender1 5.4414702 9.5787000 0.56808024 0.572804354
## wtlb1 0.6129677 0.5587619 1.09701064 0.278472547
## htcm1 -1.4232702 1.1601966 -1.22674918 0.226296885
## waistcm1 0.1029533 0.3986637 0.25824592 0.797394605
## hipcm1 0.2409366 0.8110122 0.29708139 0.767771427
## bmi1c -3.1513040 3.8521861 -0.81805600 0.417635993
## cig1c -0.3843460 3.8830521 -0.09898039 0.921593454predict(m7, interval="confidence")## fit lwr upr
## 1 118.74753 108.66279 128.8323
## 2 135.38547 120.89167 149.8793
## 3 109.31460 88.93789 129.6913
## 4 122.29859 111.80981 132.7874
## 5 113.47082 100.98952 125.9521
## 6 119.26593 107.52468 131.0072
## 7 105.20456 89.98792 120.4212
## 8 106.48905 94.56992 118.4082
## 9 125.14448 109.97390 140.3150
## 10 121.15361 109.47750 132.8297
## 11 134.53472 120.14468 148.9248
## 12 119.63293 108.53250 130.7334
## 13 125.39301 114.33290 136.4531
## 14 126.45770 109.17702 143.7384
## 15 111.52562 100.14198 122.9093
## 16 140.55276 127.39503 153.7105
## 17 115.12784 91.11850 139.1372
## 18 127.49897 116.54123 138.4567
## 19 133.75750 121.67375 145.8412
## 20 100.45875 87.75757 113.1599
## 21 105.07081 87.13288 123.0087
## 22 143.39441 129.20865 157.5802
## 23 131.45029 119.21202 143.6886
## 24 120.24720 103.70459 136.7898
## 25 122.49233 107.20355 137.7811
## 26 104.96289 92.58936 117.3364
## 27 108.90148 91.84711 125.9559
## 28 149.36405 133.79226 164.9358
## 29 130.77317 113.66675 147.8796
## 30 111.80469 93.44199 130.1674
## 31 144.54915 118.75598 170.3423
## 32 126.80615 113.82299 139.7893
## 33 100.41992 84.69085 116.1490
## 34 93.12470 75.87788 110.3715
## 35 110.92397 98.91538 122.9326
## 36 103.77512 90.81454 116.7357
## 37 119.32726 105.88151 132.7730
## 38 101.70220 89.25893 114.1455
## 39 139.38335 123.98016 154.7865
## 40 117.91692 99.20381 136.6300
## 41 137.99631 123.06142 152.9312
## 42 114.11923 104.57152 123.6669
## 43 133.44397 117.29114 149.5968
## 44 114.43757 101.80307 127.0721
## 45 110.34207 99.69368 120.9905
## 46 104.34620 90.24895 118.4434
## 47 98.27743 87.39417 109.1607
## 48 119.38543 108.64595 130.1249
## 49 117.37691 106.44098 128.3128
## 50 136.71139 120.17671 153.2461
## 51 115.71987 98.86856 132.5712
## 52 115.55127 106.79754 124.3050
## 53 127.27858 104.99689 149.5603
## 54 104.19658 88.61646 119.7767
## 55 127.01272 110.73058 143.2948predict(m7, interval="prediction")## Warning in predict.lm(m7, interval = "prediction"): predictions on current data refer to _future_ responses## fit lwr upr
## 1 118.74753 82.56119 154.9339
## 2 135.38547 97.73152 173.0394
## 3 109.31460 69.02863 149.6006
## 4 122.29859 85.99758 158.5996
## 5 113.47082 76.54479 150.3968
## 6 119.26593 82.58343 155.9484
## 7 105.20456 67.26652 143.1426
## 8 106.48905 69.74922 143.2289
## 9 125.14448 87.22488 163.0641
## 10 121.15361 84.49191 157.8153
## 11 134.53472 96.92060 172.1488
## 12 119.63293 83.15049 156.1154
## 13 125.39301 88.92282 161.8632
## 14 126.45770 87.64570 165.2697
## 15 111.52562 74.95601 148.0952
## 16 140.55276 103.39263 177.7129
## 17 115.12784 72.88808 157.3676
## 18 127.49897 91.05969 163.9383
## 19 133.75750 96.96394 170.5511
## 20 100.45875 63.45782 137.4597
## 21 105.07081 65.96175 144.1799
## 22 143.39441 105.85796 180.9309
## 23 131.45029 94.60569 168.2949
## 24 120.24720 81.75815 158.7362
## 25 122.49233 84.52529 160.4594
## 26 104.96289 68.07314 141.8526
## 27 108.90148 70.18972 147.6132
## 28 149.36405 111.28216 187.4459
## 29 130.77317 92.03845 169.5079
## 30 111.80469 72.49899 151.1104
## 31 144.54915 101.27055 187.8277
## 32 126.80615 89.70748 163.9048
## 33 100.41992 62.27346 138.5664
## 34 93.12470 54.32777 131.9216
## 35 110.92397 74.15502 147.6929
## 36 103.77512 66.68434 140.8659
## 37 119.32726 82.06418 156.5903
## 38 101.70220 64.78900 138.6154
## 39 139.38335 101.37010 177.3966
## 40 117.91692 78.44630 157.3875
## 41 137.99631 100.17039 175.8222
## 42 114.11923 78.07886 150.1596
## 43 133.44397 95.12083 171.7671
## 44 114.43757 77.45947 151.4157
## 45 110.34207 73.99462 146.6895
## 46 104.34620 66.84311 141.8493
## 47 98.27743 61.86047 134.6944
## 48 119.38543 83.01119 155.7597
## 49 117.37691 80.94418 153.8096
## 50 136.71139 98.22575 175.1970
## 51 115.71987 77.09713 154.3426
## 52 115.55127 79.71307 151.3895
## 53 127.27858 85.99634 168.5608
## 54 104.19658 66.11129 142.2819
## 55 127.01272 88.63490 165.3905Interpretation If we conduct the experiment 100 times examning the level of PM2.5 exposure on the SBP, we are confident that there are 95 times which the SBP will fall into 95% CI. If we have a new level of exposure with pm2.5, the expected sbp will fall into 95% PI.
LS0tDQp0aXRsZTogIkZpbmFsIHByb2plY3Qgb24gcmVncmVzc2lvbiBjbGFzcyAiDQpzdWJ0aXRsZTogIkRpcmVjdGVkIGJ5IFByb2Zlc3NvciBLaW0sIFN1biBZb3VuZyINCmF1dGhvcjogIlRoYW5nLCBUQiINCmRhdGU6ICIxOSBEZWNlbWJlciAyMDE4Ig0Kb3V0cHV0OiANCiBodG1sX2RvY3VtZW50OiANCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlDQotLS0NCg0KDQojIyBpbXBvcnQgZGF0YQ0KDQpgYGB7cn0NCmRhdD1yZWFkLmNzdigiQzovVXNlcnMvQklOSCBUSEFORy9Ecm9wYm94L1BoRC9MZWN0dXJlcy9yZWdyZXNzaW9uL0ZpbmFsIGFzc2lnbWVudC9teWRhdGEudHh0IikNCg0KbmFtZXMoZGF0KQ0KYGBgDQoNCiMgUTE6IA0KIyMgTWV0aG9kIDENCg0KYGBge3J9DQoNCnQoc3VtbWFyeShkYXRbLGMoInBtMjUiICAsImxvbmdpdHVkZSIsImxhdGl0dWRlIiwgImxvZy5tX3RvX2ExIiAgLCAgImxvZy5tX3RvX2EyIiAgLCAgImxvZy5tX3RvX2EzIiAsICJsb2cubV90b19jb2FzdCIsICJsb2cubV90b19jb21tIiwgICJsbF9hMV9zMDA1MDAiICwgICJsbF9hM19zMDA1MDAiICwgDQogICAgICAgICAgICAgICAgICJsbF9hMV9zMDEwMDAiICwgICJsbF9hM19zMDEwMDAiLCAicG9wX3MwMTAwMCIgLCAicG9wX3MwMTUwMCIgLCAicG9wX3MwMjAwMCIgKV0pKQ0KDQoNCmBgYA0KIyMgbWV0aG9kIDINCg0KYGBge3J9DQojIGluc3RhbGwgcGFja2VhZ2UgKHN0YXJnYXplcikgDQoNCmxpYnJhcnkoc3RhcmdhemVyKQ0KDQpkZiA8LSBkYXRhLmZyYW1lKGRhdCkNCg0KY29scyA8LSBjKCJwbTI1IiAgLCJsb25naXR1ZGUiLCJsYXRpdHVkZSIsICJsb2cubV90b19hMSIgICwgICJsb2cubV90b19hMiIgICwgICJsb2cubV90b19hMyIgLCAibG9nLm1fdG9fY29hc3QiLCAibG9nLm1fdG9fY29tbSIsICAibGxfYTFfczAwNTAwIiAsICAibGxfYTNfczAwNTAwIiAsIA0KICAgICAgICAgICJsbF9hMV9zMDEwMDAiICwgICJsbF9hM19zMDEwMDAiLCAicG9wX3MwMTAwMCIgLCAicG9wX3MwMTUwMCIgLCAicG9wX3MwMjAwMCIgKQ0Kc3RhcmdhemVyKA0KICBkZlssIGNvbHNdLCB0eXBlID0gInRleHQiLCANCiAgc3VtbWFyeS5zdGF0ID0gYygibWluIiwgInAyNSIsICJtZWRpYW4iLCAicDc1IiwgIm1heCIsICJtZWFuIiwgInNkIikNCikNCmBgYA0KDQojIFEyOiBjb3JyZWxhdGlvbiBhbmQgVklGDQoNCiMjIENvcnJlbGF0aW9uDQoNCmBgYHtyfQ0KcnN1YnNldD1zdWJzZXQoZGF0LCBzZWxlY3Q9YygicG0yNSIgICwibG9uZ2l0dWRlIiwibGF0aXR1ZGUiLCAibG9nLm1fdG9fYTEiICAsICAibG9nLm1fdG9fYTIiICAsICAibG9nLm1fdG9fYTMiICwgImxvZy5tX3RvX2NvYXN0IiwgImxvZy5tX3RvX2NvbW0iLCAgImxsX2ExX3MwMDUwMCIgLCAgImxsX2EzX3MwMDUwMCIgLCANCiAgICAgICAgICAgICAgICAgICAgICAibGxfYTFfczAxMDAwIiAsICAibGxfYTNfczAxMDAwIiwgInBvcF9zMDEwMDAiICwgInBvcF9zMDE1MDAiICwgInBvcF9zMDIwMDAiKSkNCg0KDQpyIDwtIGNvcihkYXRbLGMoImxvbmdpdHVkZSIsImxhdGl0dWRlIiwgImxvZy5tX3RvX2ExIiAgLCAgImxvZy5tX3RvX2EyIiAgLCAgImxvZy5tX3RvX2EzIiAsICJsb2cubV90b19jb2FzdCIsICJsb2cubV90b19jb21tIiwgICJsbF9hMV9zMDA1MDAiICwgICJsbF9hM19zMDA1MDAiICwgDQogICAgICAgICAgICAgICAgImxsX2ExX3MwMTAwMCIgLCAgImxsX2EzX3MwMTAwMCIsICJwb3BfczAxMDAwIiAsICJwb3BfczAxNTAwIiAsICJwb3BfczAyMDAwIildKQ0KDQpjb3IoZGF0WyxjKCJwbTI1IiwibG9uZ2l0dWRlIiwibGF0aXR1ZGUiLCAibG9nLm1fdG9fYTEiICAsICAibG9nLm1fdG9fYTIiICAsICAibG9nLm1fdG9fYTMiICwgImxvZy5tX3RvX2NvYXN0IiwgImxvZy5tX3RvX2NvbW0iLCAgImxsX2ExX3MwMDUwMCIgLCAgImxsX2EzX3MwMDUwMCIgLCANCiAgICAgICAgICAgImxsX2ExX3MwMTAwMCIgLCAgImxsX2EzX3MwMTAwMCIsICJwb3BfczAxMDAwIiAsICJwb3BfczAxNTAwIiAsICJwb3BfczAyMDAwIildKVsxLF0gIA0KYGBgDQoNCg0KIyMjIG90aGVyIG1ldGhvZCBmb3IgciBhbmQgcC12YWx1ZQ0KDQpgYGB7cn0NCiMgSGlkZSB1cHBlciB0cmlhbmdsZSB0byBnZXQgciAoY29ycmVsYXRpb24pDQoNCmxpYnJhcnkoeHRhYmxlKQ0KbWNvcjwtcm91bmQoY29yKHJzdWJzZXQpLDIpDQoNCm1jb3INCg0KDQp1cHBlcjwtbWNvcg0KdXBwZXJbdXBwZXIudHJpKG1jb3IpXTwtIiINCnVwcGVyPC1hcy5kYXRhLmZyYW1lKHVwcGVyKQ0KdXBwZXINCg0KDQoNCg0KIyMjIHAgdmFsdWUgb2YgY29ycmVsYXRpb24gKCBpbnN0YWxsLnBhY2thZ2VzKEhtaXNjKSkgDQoNCmxpYnJhcnkoIkhtaXNjIiwgbGliLmxvYz0ifi9SL3dpbi1saWJyYXJ5LzMuNCIpDQoNCm15Y29yIDwtIHJjb3JyKGFzLm1hdHJpeChyc3Vic2V0KSwgdHlwZT0icGVhcnNvbiIpDQoNCm15Y29yJFANCg0KYGBgDQoNCg0KIyMgIFZJRg0KDQpgYGB7cn0NCiMjIyBWSUYNCg0KbTEgPC0gbG0ocG0yNSB+IGxvbmdpdHVkZSArIGxhdGl0dWRlICsgbG9nLm1fdG9fYTEgKyBsb2cubV90b19hMiArIGxvZy5tX3RvX2EzICsgbG9nLm1fdG9fY29hc3QgKyBsb2cubV90b19jb21tICANCiAgICAgICAgICsgbGxfYTFfczAwNTAwICsgbGxfYTNfczAwNTAwICsgbGxfYTFfczAxMDAwICsgbGxfYTNfczAxMDAwICsgcG9wX3MwMTAwMCArIHBvcF9zMDE1MDAgKyBwb3BfczAyMDAwLCBkYXRhPWRhdCkNCiMjIyBpbnN0YWxsLnBhY2thZ2VzKGNhcikNCmxpYnJhcnkoY2FyKQ0KDQp2aWYobTEpDQoNCmBgYA0KDQoNCkJhc2VkIG9uIHRoZSByZXN1bHQgZm9yIGNvcnJlbGF0aW9uIGFuZCBWSUYsIHdlIGZvdW5kIHRoYXQgMyB2YXJpYWJsZXMgcG9wX3MwMjAwMCwgcG9wX3MwMTUwMCwgYW5kIHBvcF9zMDEwMDAgaGF2ZSBzdHJvbmcgIGNvcnJlbGF0aW9uIHdpdGggZWFjaCBvdGhlcnMgYW5kIFZJRj4xMC4gIEFkZGl0aW9uYWxseSwgIHBvcF9zMDEwMDAgdmFyaWFibGUgYXJlIHNob3dpbmcgbG93IGNvcnJlYWx0aW9uIHdpdGggb3V0Y29tZSAoUE0yLjUpLiBTbyB3ZSBjYW4gY29uc2lkZXIgdG8gZXhjbHVkZSB0aGlzIHZhcmlhYmxlLiANCg0KDQojUTMuIENvbXBhcmlzaW9uIG1vZGVsDQoNCg0KDQpgYGB7cn0NCg0KbTIgPC0gbG0ocG0yNSB+IGxvbmdpdHVkZSArIGxhdGl0dWRlICsgbG9nLm1fdG9fYTEgKyBsb2cubV90b19hMiArIGxvZy5tX3RvX2EzICsgbG9nLm1fdG9fY29hc3QgKyBsb2cubV90b19jb21tICANCiAgICAgICAgICsgbGxfYTFfczAwNTAwICsgbGxfYTNfczAwNTAwICsgbGxfYTFfczAxMDAwICsgbGxfYTNfczAxMDAwICsgcG9wX3MwMTUwMCArIHBvcF9zMDIwMDAsIGRhdGE9ZGF0KQ0KDQoNCg0KYT1zdGVwKG0yLCBkaXJlY3Rpb249ImJhY2t3YXJkIiwgdGVzdD0iRiIpIA0KDQpzdW1tYXJ5KGEpDQoNCnJxMj1zdWJzZXQoZGF0LCBzZWxlY3Q9YygicG0yNSIgICwibG9uZ2l0dWRlIiwibGF0aXR1ZGUiLCAibG9nLm1fdG9fYTEiICAsICAibG9nLm1fdG9fYTIiICAsICAibG9nLm1fdG9fYTMiICwgImxvZy5tX3RvX2NvYXN0IiwgImxvZy5tX3RvX2NvbW0iLCAgImxsX2ExX3MwMDUwMCIgLCAgImxsX2EzX3MwMDUwMCIgLCANCiAgICAgICAgICAgICAgICAgICAgICAibGxfYTFfczAxMDAwIiAsICAibGxfYTNfczAxMDAwIiAsICJwb3BfczAxNTAwIiAsICJwb3BfczAyMDAwIikpDQoNCm0wIDwtIGxtKHBtMjUgfiAxLCBkYXRhPXJxMikNCg0KDQpiPSBzdGVwKG0wLCBkaXJlY3Rpb249ImZvcndhcmQiLCBzY29wZT1saXN0KGxvd2VyID0gfjEsIHVwcGVyID0gbTIpLCB0ZXN0PSJGIikNCnN1bW1hcnkoYikNCg0KYz0gc3RlcChtMCwgZGlyZWN0aW9uPSJib3RoIiwgc2NvcGU9bGlzdChsb3dlciA9IH4xLCB1cHBlciA9IG0yKSwgdGVzdD0iRiIpDQpzdW1tYXJ5KGMpDQpgYGANCg0KSW50ZXJwcmV0YXRpb246IFVzaW5nIFN0ZXB3aXNlIGFwcHJvYWNoIGluIHRoZSBtdWx0aXZhcmlhdGUgcmVncmVzc2lvbiBtb2RlbCwgd2UgZm91bmQgdGhhdCAgU3VtIG9mIGRpc3RhbmNlcyBvZiBpbnRlcnN0YXRlIGxvY2FsIHJvYWRzIHdpdGhpbiA1MDAgbWV0ZXJzIChsbF9hMV9zMDA1MDApLCBMb2cgZGlzdGFuY2UgdG8gdGhlIG5lYXJlc3QgbG9jYWwgcm9hZHMgaW4gbWV0ZXIgKGxvZy5tX3RvX2EzKSwgTG9nIGRpc3RhbmNlIHRvIHRoZSBuZWFyZXN0IGxvY2FsIHJvYWRzIGluIG1ldGVyIChsb2cubV90b19hMyksIExvZyBkaXN0YW5jZSB0byB0aGUgbmVhcmVzdCBjb2FzdGxpbmUgaW4gbWV0ZXIgKGxvZy5tX3RvX2NvYXN0KSwgbGF0aXR1ZGUgYXJlIHBvc2l0aXZlbHkgc2lnbmlmaWNhbnRseSBhc3NvY2lhdGVkIHdpdGggaW5jcmVhc2luZyBjb25jZW50cmF0aW9uIG9mIFBNMi41IChwPDAuMDEpLiBXaGlsZSBMb2cgZGlzdGFuY2UgdG8gdGhlIG5lYXJlc3QgaW50ZXJzdGF0ZSBoaWdod2F5IGluIG1ldGVyIChsb2cubV90b19hMSksIExvZyBkaXN0YW5jZSB0byB0aGUgbmVhcmVzdCBoaWdod2F5IGluIG1ldGVyIChsb2cubV90b19hMiksICBsb25naXR1ZGUsIGFuZCBMb2cgZGlzdGFuY2UgdG8gdGhlIG5lYXJlc3QgY29tbWVydGlhbCBhcmVhIGluIG1ldGVyIChsb2cubV90b19jb21tKSBhcmUgbmVnYXRpdmVseSBhc3NvY2lhdGVkIHNpZ25pZmljYW5jZSB3aXRoIGNvbmNlbnRyYXRpb24gb2YgUE0yLjUgKHA8MC4wMSkuDQoNCg0KQ29uY2x1c2lvbjogQ29tcGFyaW5nIHdpdGggdGhlIHN0ZXB3aWRlIG1vZGVsLCB3ZSBmb3VuZCB0aGF0IGJvdGggYmFja3dhcmQgYW5kIGZvcndhcmQgbW9kZWxzIGFyZSBwcm9kdWNpbmcgdGhlIHNpbWlsYXIgcmVzdWx0cy4gDQoNCg0KIyBRNA0KDQpgYGB7cn0NCmxpYnJhcnkoc3RhcmdhemVyKQ0KDQpjb2xzMSA8LSBjKCJwbTI1IiAgLCJsb25naXR1ZGUiLCJsYXRpdHVkZSIsICJsb2cubV90b19hMSIgICwgICJsb2cubV90b19hMiIgICwgICJsb2cubV90b19hMyIgLCAibG9nLm1fdG9fY29hc3QiLCAibG9nLm1fdG9fY29tbSIsICAibGxfYTFfczAwNTAwIiAsICAibGxfYTNfczAwNTAwIiAsIA0KICAgICAgICAgICJsbF9hMV9zMDEwMDAiICwgICJsbF9hM19zMDEwMDAiLCAicG9wX3MwMTAwMCIgLCAicG9wX3MwMTUwMCIgLCAicG9wX3MwMjAwMCIsImFnZTFjIiwgImdlbmRlcjEiLCAid3RsYjEiLCAiaHRjbTEiLCAid2Fpc3RjbTEiLCAiaGlwY20xIiwgInNicDFjIiApDQoNCnN0YXJnYXplcigNCiAgZGZbLCBjb2xzMV0sIHR5cGUgPSAidGV4dCIsIA0KICBzdW1tYXJ5LnN0YXQgPSBjKCJtaW4iLCAicDI1IiwgIm1lZGlhbiIsICJwNzUiLCAibWF4IiwgIm1lYW4iLCAic2QiKQ0KKQ0KDQoNCg0KYGBgDQoNCg0KIyBRNQ0KDQoNCmBgYHtyfQ0KaGlzdChkYXQkc2JwMWMsIGJyZWFrcz0yMCkNCg0KZGF0JHNiMT1kYXQkc2JwMWMtbWVhbihkYXQkc2JwMWMpDQoNCg0KcXFub3JtKGRhdCRzYjEpDQpxcWxpbmUoZGF0JHNiMSkNCg0KDQpkYXQkc2IxPWxvZzEwKGRhdCRzYnAxYysxKQ0KDQpxcW5vcm0oZGF0JHNiMSkNCnFxbGluZShkYXQkc2IxKQ0KDQpgYGANCg0KSG93ZXZlciwgd2UgdHJpZWQgZGlmZmVyZW5jZSBtZXRob2RzIG9mIHRyYW5zZm9ybWF0aW9uIChsb2cxMChzYnAxYykgb3Igc2JwMWMtbWVhbihzYnAxYyksIHRoZSBkaXN0cmlidXRpb24gaXMgcmVtYWluaW5nIG5vbi1kaXN0cmlidXRpb24uIFNvIGluIHRoaXMgY2FzZSwgd2UgYXNzdW1lIFNCUCBhcmUgbm9ybWFsIGRpc3RyaWJ1dGlvbi4gDQoNCg0KIyBRNg0KDQpgYGB7cn0NCmxpYnJhcnkoIm9sc3JyIiwgbGliLmxvYz0ifi9SL3dpbi1saWJyYXJ5LzMuNCIpDQoNCmxpYnJhcnkoImthYmxlRXh0cmEiLCBsaWIubG9jPSJ+L1Ivd2luLWxpYnJhcnkvMy40IikNCg0KbV82IDwtIGxtKHNicDFjIH4gcG0yNSArIGFnZTFjICsgZ2VuZGVyMSArIGJtaTFjICsgY2lnMWMsIGRhdGE9ZGF0KQ0KbW9kM190YWJsZSA8LSBvbHNfc3RlcF9hbGxfcG9zc2libGUobV82KQ0Ka2FibGUobW9kM190YWJsZVtjKDcsMjAsMjgsMzEpLDE6OF0sIGRpZ2l0cz0yKSAgJT4lDQogIGthYmxlX3N0eWxpbmcoYm9vdHN0cmFwX29wdGlvbnMgPSAic3RyaXBlZCIsIGZ1bGxfd2lkdGggPSBGKQ0KDQpgYGANCg0KQmFzZWQgb24gdGhpcyByZXN1bHRzIHdpdGggcjIsIEFJQywgd2UgY2FuIGNob29zZSB0aGUgbW9kZWwgMyBpcyB0aGUgYmVzdCBtb2RlbC4gDQoNCiMgUTcNCg0KDQpgYGB7cn0NCg0KbTMgPC0gbG0oc2JwMWMgfiBwbTI1ICsgYWdlMWMgKyBnZW5kZXIxLCBkYXRhPWRhdCkNCm00IDwtIGxtKHNicDFjIH4gcG0yNSArIGFnZTFjICsgZ2VuZGVyMSArICBibWkxYywgZGF0YT1kYXQpDQoNCnN1bW1hcnkobTMpJGNvZWYgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgDQpzdW1tYXJ5KG00KSRjb2VmDQoNCm01IDwtIGxtKHNicDFjIH4gcG0yNSArIGFnZTFjICsgZ2VuZGVyMSArIGJtaTFjICsgcG0yNTpibWkxYywgZGF0YT1kYXQpDQoNCnN1bW1hcnkobTUpJGNvZWYNCg0KYGBgDQoNCklGIHdlIGxvb2sgYXQgdGhlIGNvZWZmaWNpZW50IG9mIFBNMi41IGluIHRoZSBtb2RlbCAxIGFuZCBtb2RlbCAyLCB0aGVyZSBhcmUgc29tZSBjaGFuZ2VzIG9mIGNvZWZmaWNpZW50LiBXZSBjYW4gY29uc2lkZXIgQk1JIGFzIGEgY29uZm91bmRlci4gIFN0ZC4gRXJyb3Igb2YgUE0yLjUgaW4gbW9kZWwgMSAoMC44NCkgaXMgbGFyZ2VyIHRoYW4gbW9kZWwgMiAoMC43OCksIHNvIHdlIGFsc28gY29uc2lkZXIgQk1JIGEgcHJlY2lzaW9uIHZhcmlhYmxlLiBJbiB0ZXJtIG9mIGludGVyYWN0aW9uIChQTTI1OkJNSSksIHdpdGggdD0tMC41OSBhbmQgcD4wLjA1LCB3ZSBjYW4gY29uY2x1ZGUgdGhhdCB0aGVyZSBhcmUgbm8gaW50ZXJhY3Rpb24gUE0yNSBhbmQgQk1JIGluIHRoZSBtb2RlbCByZWxhdGlvbnNoaXAgYmV0d2VlbiBQTTIuNSBhbmQgU0JQIGFmdGVyIGFkanVzdG1lbnQgZm9yIG90aGVyIHZhcmlhYmxlcyBzdWNoIGFzIEFnZSwgZ2VuZGVyLCBhbmQgQk1JLiANCg0KDQojIFE4DQoNCmBgYHtyfQ0KDQptNiA8LSBsbShzYnAxYyB+IHBtMjUsIGRhdGE9ZGF0KQ0KbTcgPC0gbG0oc2JwMWMgfiBwbTI1ICsgYWdlMWMgKyBnZW5kZXIxICsgd3RsYjEgKyBodGNtMSArIHdhaXN0Y20xICsgaGlwY20xICsgYm1pMWMgKyBjaWcxYywgZGF0YT1kYXQpDQoNCnN0ZXAobTYsIGRpcmVjdGlvbj0iYm90aCIsIHNjb3BlPWxpc3QobG93ZXIgPSBtNiwgdXBwZXIgPSBtNyksIHRlc3Q9IkYiKQ0KYGBgDQoNCkludGVycHJldGF0aW9uOiAwLjU0IG1taGcgZGVjcmVhc2Ugb2YgU0JQIGNvcnJlc3BvbmRpbmcgdG8gb25lIG1pY3JvbWV0ZXJzIGluY3JlYXNlIG9mIFBNMi41IGFmdGVyIGFkanVzdGluZyBmb3IgQWdlIGFuZCBCTUkuIEJ1dCB0aGlzIGRlY3JlYXNlIGFyZSBub3Qgc2lnbmlmaWNhbnQgKHA+MC4wNSkuDQoNCiMgUTkNCg0KDQpgYGB7cn0NCg0KcmVzIDwtIG03JHJlc2lkDQpsZXYgPC0gaGF0dmFsdWVzKG03KQ0Kci5qYWNrIDwtIHJzdHVkZW50KG03KQ0KY29vayA8LSBjb29rcy5kaXN0YW5jZShtNykNCg0KcGFyKG1mcm93PWMoMiwyKSkNCnBsb3QocmVzLCB5bGFiPSJyZXNpZHVhbCIpDQphYmxpbmUoaD0wLCBsdHk9MykNCnBsb3Qoci5qYWNrLCB5bGFiPSJqYWNra25pZmUgcmVzaWR1YWwiKQ0KYWJsaW5lKGg9MCwgbHR5PTMpDQpwbG90KGNvb2ssIHlsYWI9ImNvb2sncyBkaXN0YW5jZSIpDQphYmxpbmUoaD0xLCBsdHk9MykNCg0KbGlicmFyeShjYXIpDQoNCihvdXRzIDwtIGluZmx1ZW5jZVBsb3QobTcpKQ0KDQoNCmBgYA0KDQoNCmBgYHtyfQ0KIyBRMTANCnN1bW1hcnkobTcpJGNvZWYNCg0KcHJlZGljdChtNywgaW50ZXJ2YWw9ImNvbmZpZGVuY2UiKQ0KcHJlZGljdChtNywgaW50ZXJ2YWw9InByZWRpY3Rpb24iKQ0KYGBgDQoNCkludGVycHJldGF0aW9uDQpJZiB3ZSBjb25kdWN0IHRoZSBleHBlcmltZW50IDEwMCB0aW1lcyBleGFtbmluZyB0aGUgbGV2ZWwgb2YgUE0yLjUgZXhwb3N1cmUgb24gdGhlIFNCUCwgd2UgYXJlIGNvbmZpZGVudCB0aGF0IHRoZXJlIGFyZSA5NSB0aW1lcyB3aGljaCB0aGUgU0JQIHdpbGwgZmFsbCBpbnRvIDk1JSBDSS4gSWYgd2UgaGF2ZSBhIG5ldyBsZXZlbCBvZiBleHBvc3VyZSB3aXRoIHBtMi41LCB0aGUgZXhwZWN0ZWQgc2JwIHdpbGwgZmFsbCBpbnRvIDk1JSBQSS4NCg0K