Seasonal flow timeseries model
from paired flow and meterological/weather record
Seasonal mean flow
Headwaters of the Connecticut River Basin states
ALR July 2014
Conte Anadromous Fish Research Center
Fixed effects model summaries
# summaries of fixed models
for (i in season.names) {
print(paste("#####",i,"#####"))
print(summary(m.fixed.season[[i]]))
sort(vif(m.fixed.season[[i]]))
AIC(m.fixed.season[[i]])
}[1] "##### winter #####"
Call:
lm(formula = log(flow.mean) ~ log(drain_area_va) + log(non.zero(Agriculture)) +
log(non.zero(OffChannelWaterSqKM)) + log(non.zero(ReachSlopePCNT)) +
log(non.zero(Impervious)) + precip.e + precip.e.lag1, data = d.calib,
subset = season == "winter", na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-1.8802 -0.1965 0.0269 0.2421 2.1654
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.591345 0.093548 -6.32 4.8e-10
log(drain_area_va) 0.971027 0.021614 44.93 < 2e-16
log(non.zero(Agriculture)) -0.012747 0.004758 -2.68 0.00757
log(non.zero(OffChannelWaterSqKM)) 0.020969 0.006269 3.35 0.00087
log(non.zero(ReachSlopePCNT)) 0.075527 0.018521 4.08 5.1e-05
log(non.zero(Impervious)) -0.026342 0.008213 -3.21 0.00140
precip.e 0.003552 0.000164 21.64 < 2e-16
precip.e.lag1 0.001649 0.000180 9.16 < 2e-16
(Intercept) ***
log(drain_area_va) ***
log(non.zero(Agriculture)) **
log(non.zero(OffChannelWaterSqKM)) ***
log(non.zero(ReachSlopePCNT)) ***
log(non.zero(Impervious)) **
precip.e ***
precip.e.lag1 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.427 on 657 degrees of freedom
Multiple R-squared: 0.824, Adjusted R-squared: 0.822
F-statistic: 440 on 7 and 657 DF, p-value: <2e-16
[1] "##### spring #####"
Call:
lm(formula = log(flow.mean) ~ log(drain_area_va) + log(non.zero(PercentSandy)) +
log(non.zero(Forest)) + log(non.zero(Agriculture)) + log(non.zero(ReachSlopePCNT)) +
log(non.zero(precip.e)) + log(non.zero(precip.e.lag1)) +
log(non.zero(gdd)) + log(pet), data = d.calib, subset = season ==
"spring", na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-1.4318 -0.1673 0.0308 0.1936 1.5957
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.24571 0.53736 -11.62 < 2e-16 ***
log(drain_area_va) 0.95985 0.01490 64.42 < 2e-16 ***
log(non.zero(PercentSandy)) -0.00639 0.00200 -3.19 0.00147 **
log(non.zero(Forest)) 0.25214 0.04075 6.19 1.1e-09 ***
log(non.zero(Agriculture)) -0.01123 0.00299 -3.75 0.00019 ***
log(non.zero(ReachSlopePCNT)) 0.03012 0.01756 1.72 0.08680 .
log(non.zero(precip.e)) 0.66673 0.03923 17.00 < 2e-16 ***
log(non.zero(precip.e.lag1)) 0.11878 0.01866 6.37 3.6e-10 ***
log(non.zero(gdd)) 0.13374 0.03008 4.45 1.0e-05 ***
log(pet) -0.96215 0.24985 -3.85 0.00013 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.311 on 660 degrees of freedom
Multiple R-squared: 0.885, Adjusted R-squared: 0.884
F-statistic: 566 on 9 and 660 DF, p-value: <2e-16
[1] "##### summer #####"
Call:
lm(formula = log(flow.mean) ~ log(drain_area_va) + log(non.zero(DrainageClass)) +
log(non.zero(Forest)) + log(non.zero(Agriculture)) + log(non.zero(ReachSlopePCNT)) +
log(non.zero(precip.e)) + log(non.zero(precip.e.lag1)) +
log(non.zero(precip.e.lag2)) + log(non.zero(gdd)) + log(pet),
data = d.calib, subset = season == "summer", na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-1.873 -0.337 0.037 0.361 1.740
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -28.14305 1.66340 -16.92 < 2e-16 ***
log(drain_area_va) 0.94399 0.02671 35.34 < 2e-16 ***
log(non.zero(DrainageClass)) -0.55772 0.11058 -5.04 5.9e-07 ***
log(non.zero(Forest)) 0.24277 0.06676 3.64 0.00030 ***
log(non.zero(Agriculture)) -0.04225 0.00494 -8.55 < 2e-16 ***
log(non.zero(ReachSlopePCNT)) -0.10694 0.03075 -3.48 0.00054 ***
log(non.zero(precip.e)) 1.38931 0.07721 17.99 < 2e-16 ***
log(non.zero(precip.e.lag1)) 0.36951 0.06832 5.41 8.9e-08 ***
log(non.zero(precip.e.lag2)) 0.03773 0.03617 1.04 0.29735
log(non.zero(gdd)) 1.84067 0.17853 10.31 < 2e-16 ***
log(pet) -7.01483 0.75970 -9.23 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.534 on 647 degrees of freedom
Multiple R-squared: 0.778, Adjusted R-squared: 0.774
F-statistic: 226 on 10 and 647 DF, p-value: <2e-16
[1] "##### fall #####"
Call:
lm(formula = log(flow.mean) ~ log(drain_area_va) + log(non.zero(PercentSandy)) +
log(non.zero(Forest)) + log(non.zero(Agriculture)) + log(non.zero(OffChannelWaterSqKM)) +
log(non.zero(precip.e)) + log(non.zero(precip.e.lag1)) +
log(non.zero(gdd)) + log(pet), data = d.calib, subset = season ==
"fall", na.action = NULL)
Residuals:
Min 1Q Median 3Q Max
-2.7708 -0.2531 0.0362 0.3093 2.6051
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.07908 1.36619 -12.50 < 2e-16
log(drain_area_va) 0.93872 0.02786 33.70 < 2e-16
log(non.zero(PercentSandy)) 0.01220 0.00352 3.47 0.00056
log(non.zero(Forest)) 0.10346 0.04462 2.32 0.02073
log(non.zero(Agriculture)) -0.04544 0.00551 -8.25 9.6e-16
log(non.zero(OffChannelWaterSqKM)) -0.01669 0.00731 -2.28 0.02281
log(non.zero(precip.e)) 1.73096 0.06640 26.07 < 2e-16
log(non.zero(precip.e.lag1)) 0.72060 0.07244 9.95 < 2e-16
log(non.zero(gdd)) 0.17609 0.09065 1.94 0.05250
log(pet) -0.89457 0.41524 -2.15 0.03159
(Intercept) ***
log(drain_area_va) ***
log(non.zero(PercentSandy)) ***
log(non.zero(Forest)) *
log(non.zero(Agriculture)) ***
log(non.zero(OffChannelWaterSqKM)) *
log(non.zero(precip.e)) ***
log(non.zero(precip.e.lag1)) ***
log(non.zero(gdd)) .
log(pet) *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.515 on 629 degrees of freedom
Multiple R-squared: 0.809, Adjusted R-squared: 0.806
F-statistic: 295 on 9 and 629 DF, p-value: <2e-16Mixed (random year and site effects) model summaries
# summaries of fixed models
for (i in season.names) {
print(paste("#####",i,"#####"))
print(summary(m.mixed.season[[i]]))
AIC(m.mixed.season[[i]])
}[1] "##### winter #####"
Linear mixed model fit by REML ['lmerMod']
Formula:
log(flow.mean) ~ log(drain_area_va) + log(non.zero(Agriculture)) +
log(non.zero(OffChannelWaterSqKM)) + log(non.zero(ReachSlopePCNT)) +
(1 | date) + (1 | site_no) + log(non.zero(Impervious)) +
precip.e + precip.e.lag1
Data: d.calib
Subset: season == "winter"
REML criterion at convergence: 534.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.064 -0.468 0.061 0.537 2.991
Random effects:
Groups Name Variance Std.Dev.
site_no (Intercept) 0.0821 0.287
date (Intercept) 0.0670 0.259
Residual 0.0818 0.286
Number of obs: 665, groups: site_no, 77; date, 56
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.096590 0.140455 -0.69
log(drain_area_va) 0.929971 0.041508 22.40
log(non.zero(Agriculture)) -0.021456 0.010661 -2.01
log(non.zero(OffChannelWaterSqKM)) 0.039686 0.013529 2.93
log(non.zero(ReachSlopePCNT)) 0.082115 0.045221 1.82
log(non.zero(Impervious)) -0.031714 0.020120 -1.58
precip.e 0.002362 0.000184 12.85
precip.e.lag1 0.001406 0.000201 7.00
Correlation of Fixed Effects:
(Intr) lg(__) l(.(A) l(.(OC l(.(RS l(.(I) precp.
lg(drn_r_v) -0.530
lg(nn.z(A)) -0.006 -0.114
l(.(OCWSKM) 0.666 -0.192 -0.123
l(.(RSPCNT) 0.008 0.178 -0.158 0.273
lg(nn.z(I)) -0.336 0.347 -0.477 -0.263 0.454
precip.e -0.341 0.123 0.021 -0.106 0.063 -0.018
precip..lg1 -0.423 -0.011 -0.028 -0.042 0.002 0.080 -0.133
[1] "##### spring #####"
Linear mixed model fit by REML ['lmerMod']
Formula:
log(flow.mean) ~ log(drain_area_va) + log(non.zero(PercentSandy)) +
log(non.zero(Forest)) + log(non.zero(Agriculture)) + log(non.zero(ReachSlopePCNT)) +
log(non.zero(precip.e)) + log(non.zero(precip.e.lag1)) +
(1 | date) + (1 | site_no) + log(non.zero(gdd)) + log(pet)
Data: d.calib
Subset: season == "spring"
REML criterion at convergence: -152.9
Scaled residuals:
Min 1Q Median 3Q Max
-6.699 -0.532 0.032 0.520 3.756
Random effects:
Groups Name Variance Std.Dev.
site_no (Intercept) 0.0763 0.276
date (Intercept) 0.0139 0.118
Residual 0.0280 0.167
Number of obs: 670, groups: site_no, 78; date, 56
Fixed effects:
Estimate Std. Error t value
(Intercept) -5.65453 0.57895 -9.77
log(drain_area_va) 0.93740 0.03387 27.68
log(non.zero(PercentSandy)) -0.01055 0.00509 -2.07
log(non.zero(Forest)) 0.34255 0.08228 4.16
log(non.zero(Agriculture)) -0.01129 0.00807 -1.40
log(non.zero(ReachSlopePCNT)) -0.00553 0.04129 -0.13
log(non.zero(precip.e)) 0.66864 0.03733 17.91
log(non.zero(precip.e.lag1)) 0.09474 0.01483 6.39
log(non.zero(gdd)) 0.07368 0.02878 2.56
log(pet) -0.50894 0.27713 -1.84
Correlation of Fixed Effects:
(Intr) lg(__) l(.(PS l(.(F) l(.(A) l(.(RS l(.(.) l(.(.. l(.())
lg(drn_r_v) 0.006
lg(nn.(PS)) 0.017 0.011
lg(nn.z(F)) -0.594 -0.014 0.081
lg(nn.z(A)) -0.004 0.001 -0.157 0.006
l(.(RSPCNT) 0.354 0.062 0.189 -0.575 0.147
lg(nn.z(.)) -0.235 -0.019 0.023 -0.006 -0.019 -0.005
lg(n.(..1)) -0.198 0.047 0.043 -0.037 -0.032 0.083 0.284
lg(nn.zr()) -0.555 0.008 0.001 0.002 0.001 0.037 -0.109 -0.014
log(pet) 0.663 0.119 0.081 -0.033 -0.007 0.065 0.246 0.034 -0.596
[1] "##### summer #####"
Linear mixed model fit by REML ['lmerMod']
Formula:
log(flow.mean) ~ log(drain_area_va) + log(non.zero(DrainageClass)) +
log(non.zero(Forest)) + log(non.zero(Agriculture)) + log(non.zero(ReachSlopePCNT)) +
log(non.zero(precip.e)) + log(non.zero(precip.e.lag1)) +
(1 | date) + (1 | site_no) + log(non.zero(precip.e.lag2)) +
log(non.zero(gdd)) + log(pet)
Data: d.calib
Subset: season == "summer"
REML criterion at convergence: 816.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.857 -0.539 0.002 0.591 3.458
Random effects:
Groups Name Variance Std.Dev.
site_no (Intercept) 0.1148 0.339
date (Intercept) 0.0771 0.278
Residual 0.1386 0.372
Number of obs: 658, groups: site_no, 76; date, 56
Fixed effects:
Estimate Std. Error t value
(Intercept) -21.8336 2.2853 -9.55
log(drain_area_va) 0.9922 0.0514 19.30
log(non.zero(DrainageClass)) -0.8736 0.2388 -3.66
log(non.zero(Forest)) 0.1761 0.1136 1.55
log(non.zero(Agriculture)) -0.0402 0.0109 -3.68
log(non.zero(ReachSlopePCNT)) -0.0752 0.0580 -1.30
log(non.zero(precip.e)) 1.5432 0.0787 19.60
log(non.zero(precip.e.lag1)) 0.3121 0.0800 3.90
log(non.zero(precip.e.lag2)) 0.0886 0.0328 2.70
log(non.zero(gdd)) 1.2418 0.2462 5.04
log(pet) -3.4558 0.8792 -3.93
Correlation of Fixed Effects:
(Intr) lg(__) l(.(DC l(.(F) l(.(A) l(.(RS l(.(.) l(.(..1
lg(drn_r_v) -0.176
lg(nn.(DC)) -0.196 0.060
lg(nn.z(F)) -0.140 -0.054 0.067
lg(nn.z(A)) 0.021 -0.008 -0.111 -0.001
l(.(RSPCNT) -0.101 0.165 0.195 -0.607 0.149
lg(nn.z(.)) -0.180 -0.005 -0.038 0.056 0.028 -0.054
lg(n.(..1)) -0.271 -0.049 0.050 0.002 -0.019 -0.010 -0.017
lg(n.(..2)) 0.035 -0.008 0.028 -0.031 -0.024 0.056 -0.018 0.282
lg(nn.zr()) -0.908 0.217 0.055 -0.097 -0.014 0.240 0.031 0.064
log(pet) 0.730 -0.016 -0.026 0.041 -0.015 -0.102 0.151 0.013
l(.(..2 l(.())
lg(drn_r_v)
lg(nn.(DC))
lg(nn.z(F))
lg(nn.z(A))
l(.(RSPCNT)
lg(nn.z(.))
lg(n.(..1))
lg(n.(..2))
lg(nn.zr()) -0.197
log(pet) 0.071 -0.697
[1] "##### fall #####"
Linear mixed model fit by REML ['lmerMod']
Formula:
log(flow.mean) ~ log(drain_area_va) + log(non.zero(PercentSandy)) +
log(non.zero(Forest)) + log(non.zero(Agriculture)) + log(non.zero(OffChannelWaterSqKM)) +
log(non.zero(precip.e)) + log(non.zero(precip.e.lag1)) +
(1 | date) + (1 | site_no) + log(non.zero(gdd)) + log(pet)
Data: d.calib
Subset: season == "fall"
REML criterion at convergence: 838.8
Scaled residuals:
Min 1Q Median 3Q Max
-4.285 -0.512 0.051 0.606 3.868
Random effects:
Groups Name Variance Std.Dev.
site_no (Intercept) 0.2004 0.448
date (Intercept) 0.0305 0.175
Residual 0.1509 0.388
Number of obs: 639, groups: site_no, 71; date, 56
Fixed effects:
Estimate Std. Error t value
(Intercept) -15.28470 1.82552 -8.37
log(drain_area_va) 0.92524 0.07020 13.18
log(non.zero(PercentSandy)) 0.01361 0.00930 1.46
log(non.zero(Forest)) 0.22498 0.12098 1.86
log(non.zero(Agriculture)) -0.05859 0.01709 -3.43
log(non.zero(OffChannelWaterSqKM)) 0.00984 0.02000 0.49
log(non.zero(precip.e)) 1.66432 0.07469 22.28
log(non.zero(precip.e.lag1)) 0.78361 0.07751 10.11
log(non.zero(gdd)) 0.06416 0.12141 0.53
log(pet) -0.11829 0.64389 -0.18
Correlation of Fixed Effects:
(Intr) lg(__) l(.(PS l(.(F) l(.(A) l(.(OC l(.(.) l(.(.. l(.())
lg(drn_r_v) 0.000
lg(nn.(PS)) 0.068 0.002
lg(nn.z(F)) -0.290 0.009 0.215
lg(nn.z(A)) -0.004 0.160 -0.262 0.021
l(.(OCWSKM) -0.112 -0.221 -0.038 0.280 -0.368
lg(nn.z(.)) 0.048 0.009 0.042 -0.002 0.053 -0.091
lg(n.(..1)) -0.169 -0.009 0.042 0.053 0.035 -0.036 -0.107
lg(nn.zr()) -0.830 0.053 -0.065 0.010 -0.039 0.001 -0.234 -0.006
log(pet) 0.871 0.127 0.154 -0.027 0.070 -0.201 0.257 0.090 -0.764Goodness of fit statistics
season sample.n mean RMSE NSEff bias percent.bias pearsonR CV.error 1 winter 665 2.001 0.300 0.824 0 0 0.908 0.212 2 spring 670 2.608 0.218 0.885 0 0 0.941 0.118 3 summer 658 1.141 0.374 0.778 0 0 0.882 0.464 4 fall 639 1.470 0.361 0.809 0 0 0.899 0.348 season sample.n mean RMSE NSEff bias percent.bias pearsonR CV.error 1 winter 665 2.001 0.185 0.933 0 0 0.966 0.131 2 spring 670 2.608 0.108 0.972 0 0 0.986 0.059 3 summer 658 1.141 0.242 0.907 0 0 0.953 0.300 4 fall 639 1.470 0.254 0.906 0 0 0.952 0.244 season sample.n mean RMSE NSEff bias percent.bias pearsonR CV.error 1 winter 120 1.922 0.274 0.794 -0.126 -6.549 0.904 0.191 2 spring 122 2.279 0.277 0.807 0.077 3.388 0.903 0.169 3 summer 121 1.086 0.382 0.780 -0.066 -6.095 0.887 0.496 4 fall 118 1.370 0.398 0.799 0.135 9.860 0.907 0.401 season sample.n mean RMSE NSEff bias percent.bias pearsonR CV.error 1 winter 120 1.922 0.272 0.797 -0.147 -7.655 0.910 0.185 2 spring 122 2.279 0.253 0.839 0.040 1.771 0.917 0.157 3 summer 121 1.086 0.348 0.818 -0.018 -1.688 0.906 0.454 4 fall 118 1.370 0.395 0.802 0.077 5.646 0.901 0.406| season | sample.n | mean | RMSE | NSEff | bias | percent.bias | pearsonR | CV.error | |
|---|---|---|---|---|---|---|---|---|---|
| 1 | winter | 665.00 | 2.00 | 0.30 | 0.82 | -0.00 | -0.00 | 0.91 | 0.21 |
| 2 | spring | 670.00 | 2.61 | 0.22 | 0.88 | 0.00 | 0.00 | 0.94 | 0.12 |
| 3 | summer | 658.00 | 1.14 | 0.37 | 0.78 | -0.00 | -0.00 | 0.88 | 0.46 |
| 4 | fall | 639.00 | 1.47 | 0.36 | 0.81 | 0.00 | 0.00 | 0.90 | 0.35 |
| 5 | winter | 665.00 | 2.00 | 0.18 | 0.93 | 0.00 | 0.00 | 0.97 | 0.13 |
| 6 | spring | 670.00 | 2.61 | 0.11 | 0.97 | -0.00 | -0.00 | 0.99 | 0.06 |
| 7 | summer | 658.00 | 1.14 | 0.24 | 0.91 | 0.00 | 0.00 | 0.95 | 0.30 |
| 8 | fall | 639.00 | 1.47 | 0.25 | 0.91 | -0.00 | -0.00 | 0.95 | 0.24 |
| 9 | winter | 120.00 | 1.92 | 0.27 | 0.79 | -0.13 | -6.55 | 0.90 | 0.19 |
| 10 | spring | 122.00 | 2.28 | 0.28 | 0.81 | 0.08 | 3.39 | 0.90 | 0.17 |
| 11 | summer | 121.00 | 1.09 | 0.38 | 0.78 | -0.07 | -6.09 | 0.89 | 0.50 |
| 12 | fall | 118.00 | 1.37 | 0.40 | 0.80 | 0.14 | 9.86 | 0.91 | 0.40 |
| 13 | winter | 120.00 | 1.92 | 0.27 | 0.80 | -0.15 | -7.66 | 0.91 | 0.18 |
| 14 | spring | 122.00 | 2.28 | 0.25 | 0.84 | 0.04 | 1.77 | 0.92 | 0.16 |
| 15 | summer | 121.00 | 1.09 | 0.35 | 0.82 | -0.02 | -1.69 | 0.91 | 0.45 |
| 16 | fall | 118.00 | 1.37 | 0.40 | 0.80 | 0.08 | 5.65 | 0.90 | 0.41 |
Plots of observed vs predicted value correlation
Hydrograph of sample sites from calibration set
Hydrograph of validation sites
Trial of viewing validation by year
Need to work on creating some better plots 
par(mfrow=c(1,1))
plot(m.fixed.season[["winter"]], which=1)
par(mfrow=c(2,2))
plot(m.fixed.season[["winter"]], which=c(3:6), main=i)
par(mfrow=c(2,2))
for (i in season.names)
plot(m.fixed.season[[i]], which=1, main=capitalize(i),cex=.6)
for (i in season.names)
plot(m.fixed.season[[i]], which=2, main=capitalize(i),cex=.6)
for (i in season.names)
plot(m.fixed.season[[i]], which=c(3:6), main=i)
par(mfrow=c(1,1))
i[1] "fall"d.resid<-d.calib
for (i in season.names){
d.resid[d.resid$season==i,"residuals"]<-m.fixed.season[[i]]$residuals
d.resid[d.resid$season==i,"residuals.mixed"]<-residuals(m.mixed.season[[i]])
}
# d.resid<-melt(d.resid,measure.vars = "season")
d.resid$observed<-log(d.resid$flow.mean)
g<-ggplot(d.resid, aes(x=da_sqkm,y=residuals,fill=qannual)) + theme_bw() +
geom_point(cex=3,pch=23) +
geom_abline(intercept = 0, slope = 0,col="black",lty=2) +
scale_fill_continuous(limits=c(1, 20),
low="green3",high="darkorange",
na.value = "darkorange3",
name="Years of \r records") +
ggtitle("Drainage Area vs Residuals") +
facet_wrap(~season)
gWarning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
g + aes(x=observed) + ggtitle("Observed vs Residuals")Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
g + aes(y=residuals.mixed) + ggtitle("Mixed Model Drainage Are vs Residuals")Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
g + aes(x=observed) + ggtitle(" Mixed Model Observed vs Residuals")Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd
Warning: font width unknown for character 0xd