Goal:

Make a table for mvmeta output for the meta-regression section in the manuscript. And figures if applicable.

Load in R packages

Load in R packages

# lpath <- "/home/ad.abt.local/layc/R/library" # assign library for user w/o admin priveleges on ACE3
# .libPaths(lpath)
library(here) # convenience for filepaths dynamic
## here() starts at C:/Users/LayC/Projects/cira3-mortality
here <- here::here # to avoid overwrite

if(grepl("C://Users/", here())){
  lpath =  paste("C://Users/", Sys.info(  )["login"], "/R.lib", sep = "")
.libPaths(lpath)}

#library(chron)
library(tictoc) # for system time
library(readxl) # for reading in excel
library(readr)  # for csv
library(survival) # from original vodonos code. Loaded for compatibility
library(splines) # for splines
library(MASS)  # for mvmeta, likely
#library(nlme) 
library(mgcv) # from original vodonos code. Loaded for compatibility
## Loading required package: nlme
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
library(dlnm) # from original vodonos code. Loaded for compatibility
## This is dlnm 2.3.9. For details: help(dlnm) and vignette('dlnmOverview').
library(tsModel) # from original vodonos code. Loaded for compatibility
## Time Series Modeling for Air Pollution and Health (0.6)
library(mvmeta) # from original vodonos code. Loaded for compatibility
## This is mvmeta 1.0.3. For an overview type: help('mvmeta-package').
library(mixmeta) # UPDATED MVMETA package with better random specification syntax
## This is mixmeta 1.0.7. For an overview type: help('mixmeta-package').
library(data.table) # fread for faster read (read.table takes >15 min)
library(tidyverse)
## -- Attaching packages ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1     v dplyr   0.8.4
## v tibble  2.1.3     v stringr 1.4.0
## v tidyr   1.0.2     v forcats 0.4.0
## v purrr   0.3.3
## -- Conflicts -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::between()   masks data.table::between()
## x dplyr::collapse()  masks nlme::collapse()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks data.table::first()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks data.table::last()
## x dplyr::select()    masks MASS::select()
## x purrr::transpose() masks data.table::transpose()
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
select <- dplyr::select # to re-assign select to specify dplyr as is more commonly used than raster and MASS functions in script below.
# Temporary fix due to changes in nest/unnest syntax (deprecated arguments)
nest <- nest_legacy
unnest <- unnest_legacy
#library(RcppRoll) # rolling averages from vodonos? Likey not needed
library(purrr) # for iterative processing.
library(gridExtra) # for grid.arrange
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
source(here::here("/code/functions/01_cluster_load_cp.r"))# load cluster functions
rm("cluster_load_cp"); rm("cluster_load_mv") # remove extra functions used for meta analysis 

source(here::here("code/functions/03_lag_calc.R"))

source(here::here("code/functions/06_rr_calc.R"))

source(here::here("code/functions/metareg/000_ftab_mvmeta.R"))

if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/"), recursive = TRUE)}

ymat_get <- function(tp = 2, dframe = cluster_2_parms, ymat_name = "ymat_1"){
  # function to extract ymat from dframe
  cs = dframe %>% filter(time_period == tp ) %>% 
    pluck("parms" ,1,  ymat_name )    #cs[!rowSums( is.na(cs)),]
  
}

Slist_get <- function(tp =2, dframe = cluster_2_parms, Slist_name = "Slist_1"){
  # functio to get slist from dframe
  cs = dframe %>% filter(time_period == tp ) %>% 
    pluck("parms" , 1,  Slist_name ) 
  
}
cluster_fits <-  seq(1,9)  %>%  
  map_dfr(cluster_load, 
          nsim = 1000, 
          out_folder = "output/vodonos_out/2019-06-13_mvcl_4_time/", 
          time_periods = c(1,2,3,4)
          )  
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9

Examine

Look at some of the mvmeta output from the cluster specific meta-smoothing models

cluster_fits  %>% pluck("mv_out", 1, "ymat_1" ) # ymat of all cities
##                       spl1         spl2         spl3        spl4       spl5
## allentown      0.223362177  0.198929909  0.302113223  0.37267640 0.64065875
## annandale      1.592067947  1.102159027  1.398238321  1.41230645 1.82157873
## atlantic city  0.077264947 -0.016694798 -0.002552391  0.07002456 0.15926094
## baltimore     -0.066703992 -0.039379716 -0.007150994  0.07048845 0.21609647
## barnstable     0.199860389  0.089040430  0.117558485  0.29077362 1.46001543
## bergen        -0.026331468  0.048608556  0.113890584  0.19651228 0.33974409
## boston         0.029881558  0.056607030  0.048018196  0.22324269 0.64192642
## carlisle      -0.200703738 -0.053618663  0.004072488 -0.05877374 0.10789410
## dover          0.794429794  0.797074626  0.841675695  0.67448954 0.85258381
## elizabeth      0.191301233  0.210001140  0.279130764  0.34673129 0.46687332
## gettysburg    -0.068329262  0.027458093 -0.051612631  0.16795089 0.50103005
## harrisburg     0.365601390  0.292608332  0.379424365  0.28845810 0.40461245
## hartford       0.088259436  0.075639814  0.125614587  0.08145279 0.38426543
## jersy city     0.005533399  0.033438390  0.056586700  0.15170943 0.52639732
## lancast       -0.058229108  0.005920053  0.037262813  0.13721602 0.24517864
## marlboro      -0.197370303 -0.061205581 -0.104535652  0.13385565 0.48303328
## middles        0.381306724  0.306319053  0.345796044  0.33511468 0.55882694
## monmouth      -0.076003285  0.005408390 -0.041782831  0.07157760 0.10381079
## nassau        -0.165873679 -0.122767045 -0.122366542 -0.01200808 0.01801486
## newark         0.031663413  0.050541883  0.076475016  0.16638079 0.46967362
## newburgh       0.023559619  0.026587694  0.056810913  0.15977627 0.79748784
## newhaven      -0.015451022 -0.022268768 -0.034749404  0.11310331 0.48635866
## newlond       -0.010379816  0.001305222 -0.016449195  0.20137046 0.30852881
## newyork        0.112892307  0.141771793  0.002081847  0.24007191 0.84848270
## philly         0.080621509  0.083594127  0.095244735  0.20285817 0.36605050
## plymouth      -0.027413608  0.012622576  0.054512148  0.13879260 0.39423790
## reading        0.044844118  0.103148246  0.119939733  0.35204257 0.94036644
## richmond      -0.056445197 -0.037093835 -0.008144682  0.05923435 0.10311170
## rockville      0.295507678  0.333940778  0.376909419  0.47534103 0.27696733
## essex          0.055368475  0.034925302  0.061694096  0.14939194 0.78228920
## springfieldma -0.076787507 -0.021978272  0.022010201  0.13568199 0.58115590
## stamford      -0.040104706  0.040993705  0.018454321  0.04744804 0.42693811
## trenton       -0.081119457 -0.037277705 -0.022051622  0.06720402 0.23277161
## wdc            0.111375764  0.157426972  0.188700594  0.29083444 0.31139213
## wilmington    -0.080215580  0.006804586  0.085853112  0.08159035 0.26840127
## york                    NA           NA           NA          NA         NA
cluster_fits  %>% pluck("mv_out", 1, "Slist_1", "allentown") # specific vcov
##            btmean0b1  btmean0b2  btmean0b3  btmean0b4  btmean0b5
## btmean0b1 0.01627959 0.01156877 0.01318139 0.01247519 0.01307927
## btmean0b2 0.01156877 0.01132552 0.01125925 0.01130633 0.01120745
## btmean0b3 0.01318139 0.01125925 0.01468001 0.01303859 0.01457783
## btmean0b4 0.01247519 0.01130633 0.01303859 0.01588292 0.01249334
## btmean0b5 0.01307927 0.01120745 0.01457783 0.01249334 0.04486401

Cluster 1

Pulls in data used for mvmeta fit. Refits models fitted by Vodonos for QA and to identify missing city-specific data. Visual check against vodonos output in meta-regression_Cluster1.md. In cluster 1, winter and summer mean temp have r of 0.96. Avoided including both in models. Population variables are also tightly linked (0.91), but were left in for consistency with Vodonos results.

# Year 1: had data corrections, may change if data obtained from Harvard.
if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster1/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster1/"), recursive = T)}
source(here("code/scripts/cluster_1_meta_reg_data_correct.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std     totalpop    popover65     air_cond     time_con      totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  -----------  -----------  -----------  -----------  -----------  ------------
## summer_mean      1.0000000   -0.1604216     0.9557826    0.0873320   -0.0200521   -0.0839655    0.5021941    0.1770252   -0.0200521    -0.0839655
## summer_std      -0.1604216    1.0000000    -0.3379536    0.4091625   -0.0633439   -0.0544270   -0.0625313   -0.0778238   -0.0633439    -0.0544270
## winter_mean      0.9557826   -0.3379536     1.0000000   -0.0875072   -0.0284958   -0.0888064    0.4374341    0.1207286   -0.0284958    -0.0888064
## winter_std       0.0873320    0.4091625    -0.0875072    1.0000000   -0.0184793   -0.1527144   -0.1937093   -0.3628200   -0.0184793    -0.1527144
## totalpop        -0.0200521   -0.0633439    -0.0284958   -0.0184793    1.0000000    0.9114431    0.1262878    0.0712323    1.0000000     0.9114431
## popover65       -0.0839655   -0.0544270    -0.0888064   -0.1527144    0.9114431    1.0000000    0.1490393    0.2495949    0.9114431     1.0000000
## air_cond         0.5021941   -0.0625313     0.4374341   -0.1937093    0.1262878    0.1490393    1.0000000    0.7371860    0.1262878     0.1490393
## time_con         0.1770252   -0.0778238     0.1207286   -0.3628200    0.0712323    0.2495949    0.7371860    1.0000000    0.0712323     0.2495949
## totpop1         -0.0200521   -0.0633439    -0.0284958   -0.0184793    1.0000000    0.9114431    0.1262878    0.0712323    1.0000000     0.9114431
## popover65_1     -0.0839655   -0.0544270    -0.0888064   -0.1527144    0.9114431    1.0000000    0.1490393    0.2495949    0.9114431     1.0000000
meta_reg_tab_cluster1  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster1/cluster_1_meta_reg.csv"))
meta_reg_tab_cluster1
##         model_name cluster   lag
## 1         int_lag0       1 lag 0
## 2        time_lag0       1 lag 0
## 3        mod1_lag0       1 lag 0
## 4        mod2_lag0       1 lag 0
## 5  mod_ac_int_lag0       1 lag 0
## 6       modac_lag0       1 lag 0
## 7         int_ma15       1  MA15
## 8        time_ma15       1  MA15
## 9        mod1_ma15       1  MA15
## 10       mod2_ma15       1  MA15
## 11       mod3_ma15       1  MA15
##                                                                                          mod
## 1                                                                                (Intercept)
## 2                                                        (Intercept) + time1 + time2 + time4
## 3                  (Intercept) + summer_mean + time1 + time2 + time4 + totpop1 + popover65_1
## 4                                (Intercept) + time1 + time2 + time4 + totpop1 + popover65_1
## 5        (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 6                               (Intercept) + time1 + time2 + time4 + air_cond + summer_mean
## 7                                                                                (Intercept)
## 8                                                        (Intercept) + time1 + time2 + time4
## 9                  (Intercept) + time1 + time2 + time4 + winter_mean + totpop1 + popover65_1
## 10                               (Intercept) + time1 + time2 + time4 + totpop1 + popover65_1
## 11 (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##            ref         Q df_q          Q_p       I2       AIC       BIC
## 1              1196.9373  705 0.000000e+00 41.09967 -1629.921 -1538.616
## 2  (Intercept)  963.2552  690 2.370126e-11 28.36789 -1682.015 -1522.231
## 3  (Intercept)  851.1906  675 4.298623e-06 20.69931 -1728.296 -1500.032
## 4  (Intercept)  949.4542  680 3.185829e-11 28.37990 -1670.272 -1464.835
## 5  (Intercept)  794.8092  675 9.560312e-04 15.07396 -1770.727 -1542.464
## 6  (Intercept)  849.7423  680 9.130792e-06 19.97574 -1738.166 -1532.729
## 7              1108.1146  846 2.802937e-09 23.65411 -2096.715 -1968.530
## 8  (Intercept) 1008.5797  828 1.545127e-05 17.90436 -2114.704 -1901.063
## 9  (Intercept)  903.9894  810 1.169714e-02 10.39718 -2138.207 -1839.109
## 10 (Intercept)  949.8648  816 7.772441e-04 14.09303 -2113.013 -1842.401
## 11 (Intercept)  912.3799  810 6.944024e-03 11.22120 -2129.805 -1830.707

Cluster 2 metareg

Note that in cluster 2, total and 65+ populations had pearson r of 0.98; removed the total population from models.

if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster2/"))){dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster2/"))}
source(here("code/scripts/cluster_2_meta_reg.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std    totalpop    popover65    air_cond     time_con     totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  ----------  -----------  ----------  -----------  ----------  ------------
## summer_mean      1.0000000    0.3114383     0.4975557    0.3773406   0.1164384    0.1242576   0.5140942    0.1031048   0.1164384     0.1242576
## summer_std       0.3114383    1.0000000    -0.4136432    0.6228724   0.0271125    0.0026857   0.1209993   -0.0993570   0.0271125     0.0026857
## winter_mean      0.4975557   -0.4136432     1.0000000   -0.5157627   0.0482538    0.0659956   0.2700248    0.2087518   0.0482538     0.0659956
## winter_std       0.3773406    0.6228724    -0.5157627    1.0000000   0.0136148   -0.0097478   0.0714751   -0.2538976   0.0136148    -0.0097478
## totalpop         0.1164384    0.0271125     0.0482538    0.0136148   1.0000000    0.9882329   0.0893624    0.0197545   1.0000000     0.9882329
## popover65        0.1242576    0.0026857     0.0659956   -0.0097478   0.9882329    1.0000000   0.1263842    0.0793820   0.9882329     1.0000000
## air_cond         0.5140942    0.1209993     0.2700248    0.0714751   0.0893624    0.1263842   1.0000000    0.7729171   0.0893624     0.1263842
## time_con         0.1031048   -0.0993570     0.2087518   -0.2538976   0.0197545    0.0793820   0.7729171    1.0000000   0.0197545     0.0793820
## totpop1          0.1164384    0.0271125     0.0482538    0.0136148   1.0000000    0.9882329   0.0893624    0.0197545   1.0000000     0.9882329
## popover65_1      0.1242576    0.0026857     0.0659956   -0.0097478   0.9882329    1.0000000   0.1263842    0.0793820   0.9882329     1.0000000
##       model_name cluster   lag
## 1       int_lag0       1 lag 0
## 2      time_lag0       1 lag 0
## 4 modac_int_lag0       1 lag 0
##                                                                                   mod
## 1                                                                         (Intercept)
## 2                                                 (Intercept) + time1 + time2 + time4
## 4 (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
##           ref         Q df_q          Q_p       I2       AIC       BIC
## 1             1222.8547  828 0.000000e+00 32.28959 -1417.068 -1350.934
## 2 (Intercept) 1054.2999  816 2.942610e-08 22.60267 -1474.335 -1351.515
## 4 (Intercept)  909.3318  804 5.598028e-03 11.58343 -1550.106 -1370.600
##   model_name cluster  lag
## 1   int_ma15       1 MA15
## 2  time_ma15       1 MA15
## 3  mod1_ma15       1 MA15
## 5  mod3_ma15       1 MA15
##                                                                                         mod
## 1                                                                               (Intercept)
## 2                                                       (Intercept) + time1 + time2 + time4
## 3             (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + popover65_1
## 5 (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##           ref         Q df_q          Q_p       I2       AIC       BIC
## 1             1108.1146  846 2.802937e-09 23.65411 -2096.715 -1968.530
## 2 (Intercept) 1008.5797  828 1.545127e-05 17.90436 -2114.704 -1901.063
## 3 (Intercept)  913.3065  810 6.543176e-03 11.31126 -2130.044 -1830.946
## 5 (Intercept)  912.3799  810 6.944024e-03 11.22120 -2129.805 -1830.707
meta_reg_tab_cluster2  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster2/cluster_2_metareg.csv"))
meta_reg_tab_cluster2  
##        model_name cluster   lag
## 1        int_lag0       1 lag 0
## 2       time_lag0       1 lag 0
## 3      modac_lag0       1 lag 0
## 4  modac_int_lag0       1 lag 0
## 5       mod1_lag0       1 lag 0
## 6       mod2_lag0       1 lag 0
## 7        int_ma15       1  MA15
## 8       time_ma15       1  MA15
## 9       mod1_ma15       1  MA15
## 10      mod2_ma15       1  MA15
## 11      mod3_ma15       1  MA15
##                                                                                          mod
## 1                                                                                (Intercept)
## 2                                                        (Intercept) + time1 + time2 + time4
## 3                               (Intercept) + time1 + time2 + time4 + air_cond + summer_mean
## 4        (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 5                            (Intercept) + summer_mean + time1 + time2 + time4 + popover65_1
## 6                                          (Intercept) + time1 + time2 + time4 + popover65_1
## 7                                                                                (Intercept)
## 8                                                        (Intercept) + time1 + time2 + time4
## 9              (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + popover65_1
## 10                                         (Intercept) + time1 + time2 + time4 + popover65_1
## 11 (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##            ref         Q df_q          Q_p       I2       AIC       BIC
## 1              1222.8547  828 0.000000e+00 32.28959 -1417.068 -1350.934
## 2  (Intercept) 1054.2999  816 2.942610e-08 22.60267 -1474.335 -1351.515
## 3  (Intercept)  946.5510  808 5.125053e-04 14.63746 -1526.392 -1365.782
## 4  (Intercept)  909.3318  804 5.598028e-03 11.58343 -1550.106 -1370.600
## 5  (Intercept)  952.6739  808 3.126156e-04 15.18609 -1526.245 -1365.635
## 6  (Intercept) 1040.2301  812 8.851209e-08 21.94035 -1473.302 -1331.587
## 7              1108.1146  846 2.802937e-09 23.65411 -2096.715 -1968.530
## 8  (Intercept) 1008.5797  828 1.545127e-05 17.90436 -2114.704 -1901.063
## 9  (Intercept)  913.3065  810 6.543176e-03 11.31126 -2130.044 -1830.946
## 10 (Intercept)  968.5495  822 2.955231e-04 15.13082 -2110.686 -1868.559
## 11 (Intercept)  912.3799  810 6.944024e-03 11.22120 -2129.805 -1830.707

Cluster 3

This cluster had a frame-shift issue similar to Cluster 1. St. Louis was removed from the meta-regression, but the data had Wichita removed instead. Removed both cities.

# Year 1 missing the last city in the time period [1:26,] 
# Vodonos lag0 - mod 1 summer_mean+as.factor(time)+totpop1+popover65_1
# Vodonos lag04 mod1 summer_mean+winter_mean+as.factor(time)+totpop1+popover65_1
if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster3/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster3/"), recursive = T)}
source(here("code/scripts/cluster_3_meta_reg.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std     totalpop    popover65     air_cond     time_rev      totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  -----------  -----------  -----------  -----------  -----------  ------------
## summer_mean      1.0000000    0.3332959     0.8261337    0.1633614   -0.3777460   -0.3941669    0.6116407    0.1508590   -0.3777460    -0.3941669
## summer_std       0.3332959    1.0000000    -0.1411438    0.4852300   -0.2931388   -0.3464019    0.1852764   -0.1525042   -0.2931388    -0.3464019
## winter_mean      0.8261337   -0.1411438     1.0000000   -0.2544158   -0.2619416   -0.2377128    0.5307216    0.2852607   -0.2619416    -0.2377128
## winter_std       0.1633614    0.4852300    -0.2544158    1.0000000   -0.2754959   -0.3960677   -0.1888755   -0.4961840   -0.2754959    -0.3960677
## totalpop        -0.3777460   -0.2931388    -0.2619416   -0.2754959    1.0000000    0.9504208   -0.2083284   -0.0036072    1.0000000     0.9504208
## popover65       -0.3941669   -0.3464019    -0.2377128   -0.3960677    0.9504208    1.0000000   -0.1380769    0.1325274    0.9504208     1.0000000
## air_cond         0.6116407    0.1852764     0.5307216   -0.1888755   -0.2083284   -0.1380769    1.0000000    0.7708572   -0.2083284    -0.1380769
## time_rev         0.1508590   -0.1525042     0.2852607   -0.4961840   -0.0036072    0.1325274    0.7708572    1.0000000   -0.0036072     0.1325274
## totpop1         -0.3777460   -0.2931388    -0.2619416   -0.2754959    1.0000000    0.9504208   -0.2083284   -0.0036072    1.0000000     0.9504208
## popover65_1     -0.3941669   -0.3464019    -0.2377128   -0.3960677    0.9504208    1.0000000   -0.1380769    0.1325274    0.9504208     1.0000000
meta_reg_tab_cluster_3  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster3/cluster_3_metareg.csv"))
meta_reg_tab_cluster_3 
##         model_name cluster   lag
## 1         int_lag0       1 lag 0
## 2        time_lag0       1 lag 0
## 3        mod1_lag0       1 lag 0
## 4        mod2_lag0       1 lag 0
## 5  mod_ac_int_lag0       1 lag 0
## 6       modac_lag0       1 lag 0
## 7         int_ma15       1  MA15
## 8        time_ma15       1  MA15
## 9        mod1_ma15       1  MA15
## 10       mod2_ma15       1  MA15
## 11       mod3_ma15       1  MA15
##                                                                                          mod
## 1                                                                                (Intercept)
## 2                                                        (Intercept) + time1 + time2 + time4
## 3                            (Intercept) + summer_mean + time1 + time2 + time4 + popover65_1
## 4                                          (Intercept) + time1 + time2 + time4 + popover65_1
## 5        (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 6                               (Intercept) + time1 + time2 + time4 + air_cond + summer_mean
## 7                                                                                (Intercept)
## 8                                                        (Intercept) + time1 + time2 + time4
## 9              (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + popover65_1
## 10                                         (Intercept) + time1 + time2 + time4 + popover65_1
## 11 (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##            ref        Q df_q          Q_p       I2       AIC       BIC
## 1              565.4460  372 3.589863e-10 34.21122 -799.3973 -744.3831
## 2  (Intercept) 474.6819  360 4.546195e-05 24.15975 -815.6369 -713.4676
## 3  (Intercept) 410.6064  352 1.693118e-02 14.27314 -823.3892 -689.7832
## 4  (Intercept) 417.4161  356 1.369693e-02 14.71340 -826.4175 -708.5298
## 5  (Intercept) 414.9380  348 7.854455e-03 16.13205 -815.0832 -665.7588
## 6  (Intercept) 422.5535  352 5.794518e-03 16.69694 -818.0682 -684.4622
## 7              391.7059  372 2.313483e-01  5.03078 -846.9324 -791.9182
## 8  (Intercept) 357.5769  360 5.261694e-01  0.00000 -850.9162 -748.7469
## 9  (Intercept) 344.8606  348 5.374712e-01  0.00000 -838.3713 -689.0469
## 10 (Intercept) 351.7207  356 5.540928e-01  0.00000 -848.2381 -730.3504
## 11 (Intercept) 342.1529  348 5.783125e-01  0.00000 -840.9580 -691.6336

Cluster 4

#vodonos lag0 mod1 summer_mean + winter_mean + as.factor(time) 
#vodonos lag05 mod1 summer_mean+winter_mean+ as.factor(time)+totpop1+popover65_1+totpop1+air_cond
if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster4/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster4/"), recursive = T)}
source(here("code/scripts/cluster_4_meta_reg.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std    totalpop   popover65     air_cond     time_con     totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  ----------  ----------  -----------  -----------  ----------  ------------
## summer_mean      1.0000000   -0.3035588     0.8367802   -0.0342598   0.5439749   0.4941585    0.2955859    0.1241268   0.5439749     0.4941585
## summer_std      -0.3035588    1.0000000    -0.7404784    0.7548768   0.1813262   0.2362585    0.0326036   -0.1698272   0.1813262     0.2362585
## winter_mean      0.8367802   -0.7404784     1.0000000   -0.4691200   0.2335434   0.1589356    0.1140392    0.1209720   0.2335434     0.1589356
## winter_std      -0.0342598    0.7548768    -0.4691200    1.0000000   0.2438651   0.2501911   -0.1521836   -0.4196718   0.2438651     0.2501911
## totalpop         0.5439749    0.1813262     0.2335434    0.2438651   1.0000000   0.9696580    0.3350578    0.1500953   1.0000000     0.9696580
## popover65        0.4941585    0.2362585     0.1589356    0.2501911   0.9696580   1.0000000    0.4330929    0.2764142   0.9696580     1.0000000
## air_cond         0.2955859    0.0326036     0.1140392   -0.1521836   0.3350578   0.4330929    1.0000000    0.8187204   0.3350578     0.4330929
## time_con         0.1241268   -0.1698272     0.1209720   -0.4196718   0.1500953   0.2764142    0.8187204    1.0000000   0.1500953     0.2764142
## totpop1          0.5439749    0.1813262     0.2335434    0.2438651   1.0000000   0.9696580    0.3350578    0.1500953   1.0000000     0.9696580
## popover65_1      0.4941585    0.2362585     0.1589356    0.2501911   0.9696580   1.0000000    0.4330929    0.2764142   0.9696580     1.0000000
##   model_name cluster   lag
## 1   int_lag0       1 lag 0
## 2  time_lag0       1 lag 0
## 3 modac_lag0       1 lag 0
##                                                            mod         ref
## 1                                                  (Intercept)            
## 2                          (Intercept) + time1 + time2 + time4 (Intercept)
## 3 (Intercept) + time1 + time2 + time4 + air_cond + summer_mean (Intercept)
##          Q df_q          Q_p       I2       AIC       BIC
## 1 577.1849  420 5.122327e-07 27.23302 -836.0124 -779.3161
## 2 500.7304  408 1.138837e-03 18.51902 -864.5751 -759.2820
## 3 477.3856  400 4.656503e-03 16.21030 -862.4814 -724.7905
##   model_name cluster  lag
## 1   int_ma15       1 MA15
## 2  time_ma15       1 MA15
## 3  mod1_ma15       1 MA15
##                                                                                        mod
## 1                                                                              (Intercept)
## 2                                                      (Intercept) + time1 + time2 + time4
## 3 (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond + popover65_1
##           ref        Q df_q         Q_p       I2       AIC       BIC
## 1             512.5305  420 0.001319493 18.05366 -843.3958 -786.6995
## 2 (Intercept) 473.6110  408 0.013651237 13.85335 -855.0696 -749.7765
## 3 (Intercept) 445.5598  392 0.031753332 12.02078 -848.1147 -678.0259
meta_reg_tab_cluster_4  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster4/cluster_4_metareg.csv"))
meta_reg_tab_cluster_4
##        model_name cluster   lag
## 1        int_lag0       1 lag 0
## 2       time_lag0       1 lag 0
## 3      modac_lag0       1 lag 0
## 4  modac_int_lag0       1 lag 0
## 5       mod1_lag0       1 lag 0
## 6       mod2_lag0       1 lag 0
## 7       mod3_lag0       1 lag 0
## 8        int_ma15       1  MA15
## 9       time_ma15       1  MA15
## 10      mod1_ma15       1  MA15
## 11      mod2_ma15       1  MA15
## 12      mod3_ma15       1  MA15
##                                                                                          mod
## 1                                                                                (Intercept)
## 2                                                        (Intercept) + time1 + time2 + time4
## 3                               (Intercept) + time1 + time2 + time4 + air_cond + summer_mean
## 4        (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 5       (Intercept) + air_cond + summer_mean + time1 + time2 + time4 + totpop1 + popover65_1
## 6                     (Intercept) + time1 + time2 + time4 + air_cond + totpop1 + popover65_1
## 7                                (Intercept) + time1 + time2 + time4 + totpop1 + popover65_1
## 8                                                                                (Intercept)
## 9                                                        (Intercept) + time1 + time2 + time4
## 10  (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond + popover65_1
## 11                              (Intercept) + time1 + time2 + time4 + air_cond + popover65_1
## 12 (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##            ref        Q df_q          Q_p       I2       AIC       BIC
## 1              577.1849  420 5.122327e-07 27.23302 -836.0124 -779.3161
## 2  (Intercept) 500.7304  408 1.138837e-03 18.51902 -864.5751 -759.2820
## 3  (Intercept) 477.3856  400 4.656503e-03 16.21030 -862.4814 -724.7905
## 4  (Intercept) 473.9355  396 4.277516e-03 16.44433 -856.9782 -703.0883
## 5  (Intercept) 469.2884  392 4.387395e-03 16.46927 -852.3012 -682.2124
## 6  (Intercept) 477.1715  396 3.155248e-03 17.01096 -855.9226 -702.0327
## 7  (Intercept) 488.2034  400 1.651568e-03 18.06694 -858.8667 -721.1758
## 8              512.5305  420 1.319493e-03 18.05366 -843.3958 -786.6995
## 9  (Intercept) 473.6110  408 1.365124e-02 13.85335 -855.0696 -749.7765
## 10 (Intercept) 445.5598  392 3.175333e-02 12.02078 -848.1147 -678.0259
## 11 (Intercept) 460.2860  400 1.988300e-02 13.09751 -850.9217 -713.2308
## 12 (Intercept) 453.9418  396 2.331430e-02 12.76415 -848.2626 -694.3728

Cluster 5 was fitted with different numbers of knots. It appears that for time-period 2 in lag 0 and time period 4 in MA15, the meta-analysis results cut the length of the spline matrix to 4 knots rather than 5. These results are thus unreliable. Leaving out. The crosspred results from the meta-smoothing stage are still acceptable, but this cluster cannot be included in the meta-analysis.

Cluster 6

if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster6/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster6/"), recursive = T)}
source(here("code/scripts/cluster_6_meta_reg.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std    totalpop    popover65     air_cond     time_con     totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  ----------  -----------  -----------  -----------  ----------  ------------
## summer_mean      1.0000000   -0.2441263     0.6311149   -0.1060338   0.4714084    0.5206866    0.3145885    0.3328939   0.4714084     0.5206866
## summer_std      -0.2441263    1.0000000    -0.8641236    0.8133605   0.0763300    0.0021514    0.0128317   -0.1562215   0.0763300     0.0021514
## winter_mean      0.6311149   -0.8641236     1.0000000   -0.6526725   0.1462948    0.2140265    0.0307875    0.1549982   0.1462948     0.2140265
## winter_std      -0.1060338    0.8133605    -0.6526725    1.0000000   0.0434398   -0.0950771   -0.2803974   -0.4932011   0.0434398    -0.0950771
## totalpop         0.4714084    0.0763300     0.1462948    0.0434398   1.0000000    0.9623808    0.2219028    0.1542825   1.0000000     0.9623808
## popover65        0.5206866    0.0021514     0.2140265   -0.0950771   0.9623808    1.0000000    0.3116030    0.3067505   0.9623808     1.0000000
## air_cond         0.3145885    0.0128317     0.0307875   -0.2803974   0.2219028    0.3116030    1.0000000    0.8423861   0.2219028     0.3116030
## time_con         0.3328939   -0.1562215     0.1549982   -0.4932011   0.1542825    0.3067505    0.8423861    1.0000000   0.1542825     0.3067505
## totpop1          0.4714084    0.0763300     0.1462948    0.0434398   1.0000000    0.9623808    0.2219028    0.1542825   1.0000000     0.9623808
## popover65_1      0.5206866    0.0021514     0.2140265   -0.0950771   0.9623808    1.0000000    0.3116030    0.3067505   0.9623808     1.0000000
##   model_name cluster   lag
## 1   int_lag0       1 lag 0
## 2  time_lag0       1 lag 0
## 5  mod1_lag0       1 lag 0
##                                                                          mod
## 1                                                                (Intercept)
## 2                                        (Intercept) + time1 + time2 + time4
## 5 (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond
##           ref        Q df_q        Q_p        I2       AIC       BIC
## 1             332.4330  295 0.06583941 11.260301 -553.8002 -479.7246
## 2 (Intercept) 300.9885  280 0.18574372  6.973187 -547.7063 -418.0739
## 5 (Intercept) 278.3906  265 0.27392769  4.810012 -539.5465 -354.3574
##   model_name cluster  lag
## 1   int_ma15       1 MA15
## 2  time_ma15       1 MA15
## 3  moda_ma15       1 MA15
##                                                                                        mod
## 1                                                                              (Intercept)
## 2                                                      (Intercept) + time1 + time2 + time4
## 3 (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond + popover65_1
##           ref        Q df_q        Q_p        I2       AIC       BIC
## 1             283.7322  236 0.01810864 16.822977 -432.8219 -384.0930
## 2 (Intercept) 256.5406  224 0.06682955 12.684397 -424.5546 -334.0580
## 3 (Intercept) 227.7070  208 0.16609123  8.654551 -417.4109 -271.2241
meta_reg_tab_cluster_6  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster6/cluster_6_metareg.csv"))
meta_reg_tab_cluster_6
##        model_name cluster   lag
## 1        int_lag0       1 lag 0
## 2       time_lag0       1 lag 0
## 3      modac_lag0       1 lag 0
## 4  modac_int_lag0       1 lag 0
## 5       mod1_lag0       1 lag 0
## 6       mod2_lag0       1 lag 0
## 7       mod3_lag0       1 lag 0
## 8        int_ma15       1  MA15
## 9       time_ma15       1  MA15
## 10      moda_ma15       1  MA15
## 11      mod1_ma15       1  MA15
## 12      mod2_ma15       1  MA15
## 13      mod3_ma15       1  MA15
##                                                                                          mod
## 1                                                                                (Intercept)
## 2                                                        (Intercept) + time1 + time2 + time4
## 3                               (Intercept) + time1 + time2 + time4 + air_cond + summer_mean
## 4        (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 5                 (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond
## 6                     (Intercept) + time1 + time2 + time4 + air_cond + totpop1 + popover65_1
## 7                                (Intercept) + time1 + time2 + time4 + totpop1 + popover65_1
## 8                                                                                (Intercept)
## 9                                                        (Intercept) + time1 + time2 + time4
## 10  (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond + popover65_1
## 11                (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond
## 12                                            (Intercept) + time1 + time2 + time4 + air_cond
## 13 (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##            ref        Q df_q        Q_p        I2       AIC       BIC
## 1              332.4330  295 0.06583941 11.260301 -553.8002 -479.7246
## 2  (Intercept) 300.9885  280 0.18574372  6.973187 -547.7063 -418.0739
## 3  (Intercept) 291.5469  270 0.17562688  7.390533 -536.3533 -369.6831
## 4  (Intercept) 285.3658  265 0.18627170  7.136738 -532.7325 -347.5434
## 5  (Intercept) 278.3906  265 0.27392769  4.810012 -539.5465 -354.3574
## 6  (Intercept) 288.0457  265 0.15813082  8.000703 -528.0336 -342.8445
## 7  (Intercept) 295.3076  270 0.13890279  8.569909 -532.5744 -365.9041
## 8              283.7322  236 0.01810864 16.822977 -432.8219 -384.0930
## 9  (Intercept) 256.5406  224 0.06682955 12.684397 -424.5546 -334.0580
## 10 (Intercept) 227.7070  208 0.16609123  8.654551 -417.4109 -271.2241
## 11 (Intercept) 234.9262  212 0.13399160  9.758912 -419.9512 -287.6869
## 12 (Intercept) 244.6650  220 0.12179011 10.081136 -428.3016 -323.8824
## 13 (Intercept) 234.6059  212 0.13715607  9.635708 -418.2535 -285.9892
# vodonos mod1 lag0 = summer_mean+winter_mean+ as.factor(time)+air_cond
# vodonos mod2 lag0 = as.factor(time)+totpop1+popover65_1+air_cond
# vodonos mod1 lag05 = summer_mean+winter_mean+ as.factor(time)+air_cond

Cluster 7

# vodonos lag 0  mod1 summer_mean+ as.factor(time)+air_cond
# vodonos lag05 mod1 summer_mean+winter_mean+ as.factor(time)+air_cond
if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster7/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster7/"), recursive = T)}
source(here("code/scripts/cluster_7_meta_reg.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std     totalpop    popover65     air_cond     time_con      totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  -----------  -----------  -----------  -----------  -----------  ------------
## summer_mean      1.0000000    0.2281936     0.3987825    0.1759789    0.1774208    0.1109428    0.1683549    0.3447052    0.1774208     0.1109428
## summer_std       0.2281936    1.0000000    -0.7658609    0.9370874   -0.3084551   -0.4413713   -0.2600212   -0.0402018   -0.3084551    -0.4413713
## winter_mean      0.3987825   -0.7658609     1.0000000   -0.7493665    0.3371723    0.4173951    0.2532085    0.1197161    0.3371723     0.4173951
## winter_std       0.1759789    0.9370874    -0.7493665    1.0000000   -0.3801667   -0.5208813   -0.4884183   -0.2407328   -0.3801667    -0.5208813
## totalpop         0.1774208   -0.3084551     0.3371723   -0.3801667    1.0000000    0.8906829    0.4449854    0.4609526    1.0000000     0.8906829
## popover65        0.1109428   -0.4413713     0.4173951   -0.5208813    0.8906829    1.0000000    0.4998701    0.4344477    0.8906829     1.0000000
## air_cond         0.1683549   -0.2600212     0.2532085   -0.4884183    0.4449854    0.4998701    1.0000000    0.8427727    0.4449854     0.4998701
## time_con         0.3447052   -0.0402018     0.1197161   -0.2407328    0.4609526    0.4344477    0.8427727    1.0000000    0.4609526     0.4344477
## totpop1          0.1774208   -0.3084551     0.3371723   -0.3801667    1.0000000    0.8906829    0.4449854    0.4609526    1.0000000     0.8906829
## popover65_1      0.1109428   -0.4413713     0.4173951   -0.5208813    0.8906829    1.0000000    0.4998701    0.4344477    0.8906829     1.0000000
##       model_name cluster   lag
## 1       int_lag0       1 lag 0
## 2      time_lag0       1 lag 0
## 5      mod1_lag0       1 lag 0
## 4 modac_int_lag0       1 lag 0
##                                                                                   mod
## 1                                                                         (Intercept)
## 2                                                 (Intercept) + time1 + time2 + time4
## 5                        (Intercept) + summer_mean + time1 + time2 + time4 + air_cond
## 4 (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
##           ref        Q df_q        Q_p        I2       AIC       BIC
## 1             275.2446  236 0.04047038 14.258079 -565.9020 -517.1731
## 2 (Intercept) 242.8344  224 0.18476715  7.756062 -568.6033 -478.1067
## 5 (Intercept) 223.9754  216 0.34050737  3.560847 -570.4603 -452.1186
## 4 (Intercept) 219.7155  212 0.34364581  3.511590 -566.4482 -434.1839
##   model_name cluster  lag
## 1   int_ma15       1 MA15
## 2  time_ma15       1 MA15
## 4  mod1_ma15       1 MA15
##                                                               mod         ref
## 1                                                     (Intercept)            
## 2                             (Intercept) + time1 + time2 + time4 (Intercept)
## 4 (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean (Intercept)
##          Q df_q         Q_p       I2       AIC       BIC
## 1 298.8246  236 0.003478418 21.02391 -538.2779 -489.5489
## 2 266.4864  224 0.027187196 15.94319 -535.6081 -445.1115
## 4 242.0625  216 0.107756183 10.76686 -539.9236 -421.5819
meta_reg_tab_cluster_7  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster7/cluster_7_metareg.csv"))
meta_reg_tab_cluster_7
##          model_name cluster   lag
## 1          int_lag0       1 lag 0
## 2         time_lag0       1 lag 0
## 3  mod_add_all_lag0       1 lag 0
## 4    modac_int_lag0       1 lag 0
## 5         mod1_lag0       1 lag 0
## 6         mod2_lag0       1 lag 0
## 7         mod3_lag0       1 lag 0
## 8          int_ma15       1  MA15
## 9         time_ma15       1  MA15
## 10        moda_ma15       1  MA15
## 11        mod1_ma15       1  MA15
## 12        mod2_ma15       1  MA15
## 13        mod3_ma15       1  MA15
##                                                                                                   mod
## 1                                                                                         (Intercept)
## 2                                                                 (Intercept) + time1 + time2 + time4
## 3                (Intercept) + air_cond + summer_mean + time1 + time2 + time4 + totpop1 + popover65_1
## 4                 (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 5                                        (Intercept) + summer_mean + time1 + time2 + time4 + air_cond
## 6                                                      (Intercept) + time1 + time2 + time4 + air_cond
## 7                                         (Intercept) + time1 + time2 + time4 + totpop1 + popover65_1
## 8                                                                                         (Intercept)
## 9                                                                 (Intercept) + time1 + time2 + time4
## 10 (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + air_cond + totpop1 + popover65_1
## 11                                    (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean
## 12                                                     (Intercept) + time1 + time2 + time4 + air_cond
## 13                (Intercept) + time1 + time2 + time4 + air_cond + popover65_1 + air_cond:popover65_1
##            ref        Q df_q         Q_p        I2       AIC       BIC
## 1              275.2446  236 0.040470384 14.258079 -565.9020 -517.1731
## 2  (Intercept) 242.8344  224 0.184767148  7.756062 -568.6033 -478.1067
## 3  (Intercept) 217.9350  208 0.304282858  4.558691 -560.4659 -414.2791
## 4  (Intercept) 219.7155  212 0.343645813  3.511590 -566.4482 -434.1839
## 5  (Intercept) 223.9754  216 0.340507369  3.560847 -570.4603 -452.1186
## 6  (Intercept) 237.3567  220 0.200952685  7.312500 -565.6956 -461.2764
## 7  (Intercept) 235.4792  216 0.173029211  8.272139 -559.9172 -441.5755
## 8              298.8246  236 0.003478418 21.023914 -538.2779 -489.5489
## 9  (Intercept) 266.4864  224 0.027187196 15.943193 -535.6081 -445.1115
## 10 (Intercept) 236.2389  204 0.060382169 13.646728 -520.8527 -360.7433
## 11 (Intercept) 242.0625  216 0.107756183 10.766862 -539.9236 -421.5819
## 12 (Intercept) 265.3413  220 0.019700280 17.087915 -528.4953 -424.0761
## 13 (Intercept) 250.4450  212 0.036189344 15.350672 -524.8139 -392.5496

Cluster 8 # Has the same problem as cluster 5, and cannot be included in the meta-regression

Cluster 9

# vodonos lag0 mod1 summer_mean+winter_mean+ as.factor(time)
# vodonos lag5 mod1summer_mean+winter_mean+as.factor(time)
if(!dir.exists(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster9/"))){ dir.create(here("output/y4-final-manuscript-output/2020-03-03-Manu-Meta-reg/Cluster9/"), recursive = T)}
source(here("code/scripts/cluster_9_meta_reg.r"))
## Parsed with column specification:
## cols(
##   city = col_character(),
##   summer_mean = col_double(),
##   summer_std = col_double(),
##   winter_mean = col_double(),
##   winter_std = col_double(),
##   time = col_double(),
##   totalpop = col_double(),
##   popover65 = col_double(),
##   fips_sm = col_double(),
##   air_cond = col_double()
## )
## 
## 
##                summer_mean   summer_std   winter_mean   winter_std     totalpop    popover65     air_cond     time_con      totpop1   popover65_1
## ------------  ------------  -----------  ------------  -----------  -----------  -----------  -----------  -----------  -----------  ------------
## summer_mean      1.0000000    0.1627083     0.5568355    0.3320857    0.2044237    0.0733911    0.2656498    0.1961129    0.2044237     0.0733911
## summer_std       0.1627083    1.0000000    -0.4086842    0.5305974    0.0577095    0.0030144   -0.0858308    0.1328253    0.0577095     0.0030144
## winter_mean      0.5568355   -0.4086842     1.0000000   -0.4888799    0.2722466    0.2014908    0.3210037    0.1699737    0.2722466     0.2014908
## winter_std       0.3320857    0.5305974    -0.4888799    1.0000000   -0.1127992   -0.1947994   -0.1107561   -0.1613095   -0.1127992    -0.1947994
## totalpop         0.2044237    0.0577095     0.2722466   -0.1127992    1.0000000    0.9342209    0.2383239    0.2978862    1.0000000     0.9342209
## popover65        0.0733911    0.0030144     0.2014908   -0.1947994    0.9342209    1.0000000    0.2363786    0.4136327    0.9342209     1.0000000
## air_cond         0.2656498   -0.0858308     0.3210037   -0.1107561    0.2383239    0.2363786    1.0000000    0.5835344    0.2383239     0.2363786
## time_con         0.1961129    0.1328253     0.1699737   -0.1613095    0.2978862    0.4136327    0.5835344    1.0000000    0.2978862     0.4136327
## totpop1          0.2044237    0.0577095     0.2722466   -0.1127992    1.0000000    0.9342209    0.2383239    0.2978862    1.0000000     0.9342209
## popover65_1      0.0733911    0.0030144     0.2014908   -0.1947994    0.9342209    1.0000000    0.2363786    0.4136327    0.9342209     1.0000000
##         model_name cluster   lag
## 1         int_lag0       1 lag 0
## 2        time_lag0       1 lag 0
## 3 mod_add_all_lag0       1 lag 0
##                                                                                    mod
## 1                                                                          (Intercept)
## 2                                                  (Intercept) + time1 + time2 + time4
## 3 (Intercept) + air_cond + summer_mean + time1 + time2 + time4 + totpop1 + popover65_1
##           ref        Q df_q       Q_p       I2       AIC       BIC
## 1             198.6080  204 0.5933798 0.000000 -213.4489 -166.7234
## 2 (Intercept) 181.8878  192 0.6883580 0.000000 -202.5256 -115.7496
## 3 (Intercept) 217.9350  208 0.3042829 4.558691 -560.4659 -414.2791
##   model_name cluster  lag                                 mod         ref
## 1   int_ma15       1 MA15                         (Intercept)            
## 2  time_ma15       1 MA15 (Intercept) + time1 + time2 + time4 (Intercept)
##          Q df_q       Q_p I2       AIC       BIC
## 1 171.8391  204 0.9506366  0 -212.5059 -165.7803
## 2 160.7221  192 0.9513013  0 -199.0393 -112.2633
meta_reg_tab_cluster_9  %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/Cluster9/cluster_9_metareg.csv"))
meta_reg_tab_cluster_9
##          model_name cluster   lag
## 1          int_lag0       1 lag 0
## 2         time_lag0       1 lag 0
## 3  mod_add_all_lag0       1 lag 0
## 4    modac_int_lag0       1 lag 0
## 5        modac_lag0       1 lag 0
## 6         mod1_lag0       1 lag 0
## 7         mod2_lag0       1 lag 0
## 8         mod3_lag0       1 lag 0
## 9          int_ma15       1  MA15
## 10        time_ma15       1  MA15
## 11        moda_ma15       1  MA15
## 12        mod1_ma15       1  MA15
## 13        mod2_ma15       1  MA15
## 14        mod3_ma15       1  MA15
##                                                                                                   mod
## 1                                                                                         (Intercept)
## 2                                                                 (Intercept) + time1 + time2 + time4
## 3                (Intercept) + air_cond + summer_mean + time1 + time2 + time4 + totpop1 + popover65_1
## 4                 (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + air_cond:summer_mean
## 5  (Intercept) + time1 + time2 + time4 + air_cond + summer_mean + winter_mean + totpop1 + popover65_1
## 6                                     (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean
## 7                                                      (Intercept) + time1 + time2 + time4 + air_cond
## 8                                         (Intercept) + time1 + time2 + time4 + totpop1 + popover65_1
## 9                                                                                         (Intercept)
## 10                                                                (Intercept) + time1 + time2 + time4
## 11            (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean + totpop1 + popover65_1
## 12                                    (Intercept) + time1 + time2 + time4 + summer_mean + winter_mean
## 13                                                     (Intercept) + time1 + time2 + time4 + air_cond
## 14          (Intercept) + time1 + time2 + time4 + winter_mean + popover65_1 + winter_mean:popover65_1
##            ref        Q df_q       Q_p       I2       AIC        BIC
## 1              198.6080  204 0.5933798 0.000000 -213.4489 -166.72338
## 2  (Intercept) 181.8878  192 0.6883580 0.000000 -202.5256 -115.74957
## 3  (Intercept) 217.9350  208 0.3042829 4.558691 -560.4659 -414.27909
## 4  (Intercept) 285.3658  265 0.1862717 7.136738 -532.7325 -347.54338
## 5  (Intercept) 146.1928  172 0.9238431 0.000000 -196.2975  -42.77075
## 6  (Intercept) 157.4825  184 0.9221583 0.000000 -209.6957  -96.21943
## 7  (Intercept) 178.2531  188 0.6833469 0.000000 -197.7292  -97.60310
## 8  (Intercept) 167.5943  184 0.8015045 0.000000 -198.8747  -85.39838
## 9              171.8391  204 0.9506366 0.000000 -212.5059 -165.78034
## 10 (Intercept) 160.7221  192 0.9513013 0.000000 -199.0393 -112.26327
## 11 (Intercept) 147.5325  176 0.9420005 0.000000 -180.2066  -40.02997
## 12 (Intercept) 151.1997  184 0.9631286 0.000000 -192.5394  -79.06307
## 13 (Intercept) 158.4390  188 0.9426909 0.000000 -193.3001  -93.17392
## 14 (Intercept) 152.9035  180 0.9293812 0.000000 -182.8355  -56.00910
meta_reg_tab <-  bind_rows(meta_reg_tab_cluster1, meta_reg_tab_cluster2, meta_reg_tab_cluster_3, meta_reg_tab_cluster_4,  
          meta_reg_tab_cluster_6, meta_reg_tab_cluster_7, meta_reg_tab_cluster_9)

meta_reg_tab %>% write_csv(here("output/y4-final-manuscript-output/2020-03-03-Manu-meta-reg/2020-03-03-Meta-Compiled.csv"))
meta_reg_tab %>% filter(cluster ==3 & lag == "lag 0")
##  [1] model_name cluster    lag        mod        ref        Q         
##  [7] df_q       Q_p        I2         AIC        BIC       
## <0 rows> (or 0-length row.names)
print(cluster_load)
## function (c_cluster = 6, nsim = 1000, out_folder = "output/vodonos_out/2019-06-13_mvcl_4_time", 
##     time_periods = c(3, 4)) 
## {
##     require(here)
##     require(dplyr)
##     require(tidyr)
##     require(purrr)
##     require(MASS)
##     print(c_cluster)
##     cluster_f <- tibble(cluster = c_cluster, time_period = time_periods, 
##         cp_fname = here::here(out_folder, gsub("X", c_cluster, 
##             paste0("mv_clX/mv_clX/cp.clusterX.years", time_period, 
##                 ".Rdata"))), cp_out = map(cp_fname, .f = load_parms), 
##         mv_fname = here::here(out_folder, gsub("X", c_cluster, 
##             paste0("mv_clX/mv_clX/mv_clX_y", time_period, ".Rdata"))), 
##         mv_out = map(mv_fname, .f = load_parms), daily_arg = map(cp_out, 
##             ~pluck_all(., col = "cp_0", ns = nsim)), ma15_arg = map(cp_out, 
##             ~pluck_all(., col = "cp_05", ns = nsim)))
## }
## <bytecode: 0x00000000127c5778>
print(lag_calc)
## function (tmean_2lag = tmean_daily$doy_tmean0[[1]]) 
## {
##     lagged_tmean <- tmean_2lag %>% mutate(tmeanL1 = ifelse(is.na(lag(tmean0)), 
##         tmean0, lag(tmean0)), tmeanL2 = ifelse(is.na(lag(tmean0, 
##         2)), tmean0, lag(tmean0, 2)), tmeanL3 = ifelse(is.na(lag(tmean0, 
##         3)), tmean0, lag(tmean0, 3)), tmeanL4 = ifelse(is.na(lag(tmean0, 
##         4)), tmean0, lag(tmean0, 4)), tmean15 = (tmean0 + tmeanL1 + 
##         tmeanL2 + tmeanL3 + tmeanL4)/5) %>% select(-starts_with("tmeanL"))
##     return(lagged_tmean)
## }
print(calc_rr_se)
## function (tmean = seq(-25, 42, 1), coef_mat = pluck(cl4_y3, "cp_out", 
##     1, "cp_0", "coefficients"), vcov_mat = pluck(cl4_y3, "cp_out", 
##     1, "cp_0", "vcov"), center_val = 15.6, degree_val = 2, knot_vec = pluck(cl4_y3, 
##     "cp_out", 1, "knots_1"), bound_vec = c(-25, 42)) 
## {
##     require(dlnm)
##     require(splines)
##     bs_ce <- bs(center_val, degree = degree_val, knots = knot_vec, 
##         Bound = bound_vec, intercept = FALSE)
##     bs_var <- suppressWarnings(bs(tmean, degree = degree_val, 
##         knots = knot_vec, Bound = bound_vec, intercept = FALSE))
##     bs_vc <- scale(bs_var, center = bs_ce, scale = FALSE)
##     rr_log <- as.numeric(bs_vc %*% coef_mat)
##     rr_log_se <- sqrt(pmax(0, rowSums((bs_vc %*% vcov_mat) * 
##         bs_vc)))
##     RR_an <- 1 - exp(-rr_log)
##     RR_an_low <- 1 - exp(-(rr_log - 1.96 * rr_log_se))
##     RR_an_high <- 1 - exp(-(rr_log + 1.96 * rr_log_se))
##     tibble(rr_log = rr_log, rr_log_se = rr_log_se, RR = RR_an, 
##         RR_lci = RR_an_low, RR_uci = RR_an_high)
## }
print(ftab)
## function (m = time_lag0, mref = int_lag0, tests = c("LR", "wald")) 
## {
##     q <- qtest(m)
##     het <- data.frame(Q = q$Q[1], df = q$df[1], Q_p = q$pvalue[1], 
##         I2 = (q$Q[1] - q$df[1])/q$Q[1] * 100)
##     if (het$I2 < 0) {
##         het$I2 = 0
##     }
##     if (m$method != "reml") {
##         het <- cbind(het, data.frame(AIC = AIC(m), BIC = BIC(m)))
##     }
##     if ("LR" %in% tests & !is.null(mref) & m$method != "reml") {
##         lrstat <- -2 * (logLik(mref) - logLik(m))[1]
##         df <- attr(logLik(m), "df") - attr(logLik(mref), "df")
##         pvalue <- 1 - pchisq(lrstat, df)
##         lr <- data.frame(LR = lrstat, LR_df = df, LR_p = pvalue, 
##             stringsAsFactors = F)
##     }
##     if ("wald" %in% tests & !is.null(mref) & m$method != "reml") {
##         coef <- coef(m)[-grep("Int", names(coef(m)))]
##         vcov <- vcov(m)[-grep("Int", names(coef(m))), -grep("Int", 
##             names(coef(m)))]
##         waldstat <- coef %*% solve(vcov) %*% coef
##         df <- length(coef)
##         pvalue <- 1 - pchisq(waldstat, df)
##         wald <- data.frame(wald = waldstat, wald_df = df, wald_p = pvalue, 
##             stringsAsFactors = F)
##     }
##     if (class(m) == "mvmeta") {
##         labs = data.frame(mod = paste(rownames(m$coefficients), 
##             collapse = " + "), ref = paste(rownames(mref$coefficients), 
##             collapse = " + "))
##     }
##     else if (class(m) == "mixmeta") {
##         labs = data.frame(mod = paste(m$call$formula[3], collapse = ""), 
##             ref = paste(mref$call$formula[3], collapse = ""))
##     }
##     else {
##         stop("model class not supported")
##     }
##     if (m$method != "reml") {
##         if (!is.null(mref)) {
##             return(cbind(labs, het, lr, wald))
##         }
##         else {
##             return(cbind(labs, het, data.frame(LR = NA, LR_df = NA, 
##                 LR_p = NA, wald = NA, wald_df = NA, wald_p = NA, 
##                 stringsAsFactors = F)))
##         }
##     }
##     else {
##         if (!is.null(mref)) {
##             return(cbind(labs, het))
##         }
##         else {
##             return(cbind(labs, het))
##         }
##     }
## }
## <bytecode: 0x000000001db777d8>
devtools::session_info()
## - Session info ---------------------------------------------------------------
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       Windows 10 x64              
##  system   x86_64, mingw32             
##  ui       RTerm                       
##  language (EN)                        
##  collate  English_United States.1252  
##  ctype    English_United States.1252  
##  tz       America/Denver              
##  date     2020-03-04                  
## 
## - Packages -------------------------------------------------------------------
##  package     * version  date       lib source        
##  assertthat    0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
##  backports     1.1.5    2019-10-02 [1] CRAN (R 3.6.1)
##  broom         0.5.4    2020-01-27 [1] CRAN (R 3.6.2)
##  callr         3.4.2    2020-02-12 [1] CRAN (R 3.6.2)
##  cellranger    1.1.0    2016-07-27 [1] CRAN (R 3.6.0)
##  cli           2.0.1    2020-01-08 [1] CRAN (R 3.6.2)
##  colorspace    1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
##  crayon        1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
##  data.table  * 1.12.8   2019-12-09 [1] CRAN (R 3.6.2)
##  DBI           1.1.0    2019-12-15 [1] CRAN (R 3.6.2)
##  dbplyr        1.4.2    2019-06-17 [1] CRAN (R 3.6.2)
##  desc          1.2.0    2018-05-01 [1] CRAN (R 3.6.1)
##  devtools      2.2.2    2020-02-17 [1] CRAN (R 3.6.2)
##  digest        0.6.24   2020-02-12 [1] CRAN (R 3.6.2)
##  dlnm        * 2.3.9    2019-03-11 [1] CRAN (R 3.6.0)
##  dplyr       * 0.8.4    2020-01-31 [1] CRAN (R 3.6.2)
##  ellipsis      0.3.0    2019-09-20 [1] CRAN (R 3.6.2)
##  evaluate      0.14     2019-05-28 [1] CRAN (R 3.6.2)
##  fansi         0.4.1    2020-01-08 [1] CRAN (R 3.6.2)
##  forcats     * 0.4.0    2019-02-17 [1] CRAN (R 3.6.0)
##  fs            1.3.1    2019-05-06 [1] CRAN (R 3.6.0)
##  generics      0.0.2    2018-11-29 [1] CRAN (R 3.6.0)
##  ggplot2     * 3.2.1    2019-08-10 [1] CRAN (R 3.6.2)
##  glue          1.3.1    2019-03-12 [1] CRAN (R 3.6.0)
##  gridExtra   * 2.3      2017-09-09 [1] CRAN (R 3.6.0)
##  gtable        0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
##  haven         2.2.0    2019-11-08 [1] CRAN (R 3.6.2)
##  here        * 0.1      2017-05-28 [1] CRAN (R 3.6.0)
##  highr         0.8      2019-03-20 [1] CRAN (R 3.6.0)
##  hms           0.5.3    2020-01-08 [1] CRAN (R 3.6.2)
##  htmltools     0.4.0    2019-10-04 [1] CRAN (R 3.6.2)
##  httr          1.4.1    2019-08-05 [1] CRAN (R 3.6.2)
##  jsonlite      1.6.1    2020-02-02 [1] CRAN (R 3.6.2)
##  knitr         1.28     2020-02-06 [1] CRAN (R 3.6.2)
##  lattice       0.20-38  2018-11-04 [2] CRAN (R 3.6.2)
##  lazyeval      0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
##  lifecycle     0.1.0    2019-08-01 [1] CRAN (R 3.6.2)
##  lubridate     1.7.4    2018-04-11 [1] CRAN (R 3.6.0)
##  magrittr      1.5      2014-11-22 [1] CRAN (R 3.6.0)
##  MASS        * 7.3-51.5 2019-12-20 [1] CRAN (R 3.6.2)
##  Matrix        1.2-18   2019-11-27 [2] CRAN (R 3.6.2)
##  memoise       1.1.0    2017-04-21 [1] CRAN (R 3.6.1)
##  mgcv        * 1.8-31   2019-11-09 [1] CRAN (R 3.6.2)
##  mixmeta     * 1.0.7    2019-12-09 [1] CRAN (R 3.6.2)
##  modelr        0.1.5    2019-08-08 [1] CRAN (R 3.6.2)
##  munsell       0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
##  mvmeta      * 1.0.3    2019-12-10 [1] CRAN (R 3.6.2)
##  nlme        * 3.1-144  2020-02-06 [1] CRAN (R 3.6.2)
##  pillar        1.4.3    2019-12-20 [1] CRAN (R 3.6.2)
##  pkgbuild      1.0.6    2019-10-09 [1] CRAN (R 3.6.2)
##  pkgconfig     2.0.3    2019-09-22 [1] CRAN (R 3.6.2)
##  pkgload       1.0.2    2018-10-29 [1] CRAN (R 3.6.1)
##  prettyunits   1.1.1    2020-01-24 [1] CRAN (R 3.6.2)
##  processx      3.4.2    2020-02-09 [1] CRAN (R 3.6.2)
##  ps            1.3.2    2020-02-13 [1] CRAN (R 3.6.2)
##  purrr       * 0.3.3    2019-10-18 [1] CRAN (R 3.6.2)
##  R6            2.4.1    2019-11-12 [1] CRAN (R 3.6.2)
##  Rcpp          1.0.3    2019-11-08 [1] CRAN (R 3.6.2)
##  readr       * 1.3.1    2018-12-21 [1] CRAN (R 3.6.0)
##  readxl      * 1.3.1    2019-03-13 [1] CRAN (R 3.6.0)
##  remotes       2.1.1    2020-02-15 [1] CRAN (R 3.6.2)
##  reprex        0.3.0    2019-05-16 [1] CRAN (R 3.6.0)
##  rlang         0.4.4    2020-01-28 [1] CRAN (R 3.6.2)
##  rmarkdown     2.1      2020-01-20 [1] CRAN (R 3.6.2)
##  rprojroot     1.3-2    2018-01-03 [1] CRAN (R 3.6.0)
##  rstudioapi    0.11     2020-02-07 [1] CRAN (R 3.6.2)
##  rvest         0.3.5    2019-11-08 [1] CRAN (R 3.6.2)
##  scales      * 1.1.0    2019-11-18 [1] CRAN (R 3.6.2)
##  sessioninfo   1.1.1    2018-11-05 [1] CRAN (R 3.6.1)
##  stringi       1.4.3    2019-03-12 [1] CRAN (R 3.6.0)
##  stringr     * 1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
##  survival    * 3.1-8    2019-12-03 [1] CRAN (R 3.6.2)
##  testthat      2.3.1    2019-12-01 [1] CRAN (R 3.6.2)
##  tibble      * 2.1.3    2019-06-06 [1] CRAN (R 3.6.2)
##  tictoc      * 1.0      2014-06-17 [1] CRAN (R 3.6.0)
##  tidyr       * 1.0.2    2020-01-24 [1] CRAN (R 3.6.2)
##  tidyselect    1.0.0    2020-01-27 [1] CRAN (R 3.6.2)
##  tidyverse   * 1.3.0    2019-11-21 [1] CRAN (R 3.6.2)
##  tsModel     * 0.6      2013-06-24 [1] CRAN (R 3.6.0)
##  usethis       1.5.1    2019-07-04 [1] CRAN (R 3.6.1)
##  vctrs         0.2.3    2020-02-20 [1] CRAN (R 3.6.2)
##  withr         2.1.2    2018-03-15 [1] CRAN (R 3.6.0)
##  xfun          0.12     2020-01-13 [1] CRAN (R 3.6.2)
##  xml2          1.2.2    2019-08-09 [1] CRAN (R 3.6.2)
##  yaml          2.2.1    2020-02-01 [1] CRAN (R 3.6.2)
## 
## [1] C:/Users/LayC/R.Lib
## [2] C:/Program Files/R/R-3.6.2/library
## [3] C:/Users/LayC/Documents/R/win-library/3.6