Supplemental Code for Verhoeven et al.
Mike Verhoeven
2019-12-30
This code will load in the MN CLP dataset and analyze and visualize the material presented in Verhoeven et al. 2019. Additonal dataset assembly and processing code are available upon request from the authors.
# header ------------------------------------------------------------------Tools
Load libraries:
#+warning=FALSE, message=FALSE
library(ggplot2)
library(lme4)## Loading required package: Matrix library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union library(data.table)## Warning: package 'data.table' was built under R version 3.6.2##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last library(numDeriv)
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3 library(mapdata)## Warning: package 'mapdata' was built under R version 3.6.2## Loading required package: maps library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine library(rlang)## Warning: package 'rlang' was built under R version 3.6.2##
## Attaching package: 'rlang'## The following object is masked from 'package:data.table':
##
## := library(devtools)## Warning: package 'devtools' was built under R version 3.6.2## Loading required package: usethis## Warning: package 'usethis' was built under R version 3.6.2 devtools::install_github('oswaldosantos/ggsn')## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools custom from http://cran.r-project.org/bin/windows/Rtools/.## Skipping install of 'ggsn' from a github remote, the SHA1 (1e208845) has not changed since last install.
## Use `force = TRUE` to force installation library(ggsn)## Loading required package: grid library(raster)## Loading required package: sp## Warning: package 'sp' was built under R version 3.6.2##
## Attaching package: 'raster'## The following object is masked from 'package:ggsn':
##
## scalebar## The following object is masked from 'package:data.table':
##
## shift## The following object is masked from 'package:dplyr':
##
## select library(ggthemes)## Warning: package 'ggthemes' was built under R version 3.6.2 library(mapproj)
library(maps)
library(GISTools) ## Loading required package: maptools## Warning: package 'maptools' was built under R version 3.6.2## Checking rgeos availability: TRUE## Loading required package: RColorBrewer## Loading required package: MASS##
## Attaching package: 'MASS'## The following objects are masked from 'package:raster':
##
## area, select## The following object is masked from 'package:dplyr':
##
## select## Loading required package: rgeos## rgeos version: 0.5-2, (SVN revision 621)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-1
## Polygon checking: TRUE##
## Attaching package: 'GISTools'## The following object is masked from 'package:maps':
##
## map.scale library(ggmap)## Warning: package 'ggmap' was built under R version 3.6.2## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details. #' Load custom functions:
# diagnostics
mv_diag <- function(model_to_diag, plots = c(1:4)) {
par(mfrow = c(2,2))
plot(model_to_diag, which = plots)
par(mfrow = c(1,1))
}
# ben bolkers CI calc fn (https://github.com/bbolker/asaglmm/blob/master/papers/bolker_chap.rmd)
easyPredCI <- function(model,newdata,alpha=0.05) {
## baseline prediction, on the linear predictor (log-odds) scale:
pred0 <- predict(model,re.form=NA,newdata=newdata)
## fixed-effects model matrix for new data
X <- model.matrix(formula(model,fixed.only=TRUE)[-2],
newdata)
beta <- fixef(model) ## fixed-effects coefficients
V <- vcov(model) ## variance-covariance matrix of beta
pred.se <- sqrt(diag(X %*% V %*% t(X))) ## std errors of predictions
## inverse-link (logistic) function: could also use plogis()
plogis <- model@resp$family$plogis
## construct 95% Normal CIs on the link scale and
## transform back to the response (probability) scale:
crit <- -qnorm(alpha/2)
plogis(cbind(lwr=pred0-crit*pred.se,
upr=pred0+crit*pred.se))
}
# load in data ------------------------------------------------------------Load Data
Load in clp dataset from the p_crispus dataset_prep script:
pcri.e <- fread("data/output_data/pcri.e_akDNRmv_combined.csv", drop = 1:2)
pcri.p <- fread("data/output_data/pcri.p_akDNRmv_combined.csv", drop = 1:2)Visualize the distributions of the responses we are interested in:
#foc: shown with my calcs(gray) and Adam's calcs (red)
ggplot(pcri.e, aes(foc_clp))+
geom_histogram(alpha = .5)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(pcri.p, aes(foc_clp))+
geom_histogram(alpha = .5)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#rake density
ggplot(pcri.e, aes(rd_pcri))+
geom_histogram(alpha = .5 )## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 38 rows containing non-finite values (stat_bin).
ggplot(pcri.p, aes(rd_pcri))+
geom_histogram(alpha = .5)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 65 rows containing non-finite values (stat_bin).
# cleaning data -----------------------------------------------------------Data Cleaning
Are our datasets complete for the parameters we want to estimate on? Note that the rake densities will be pres for a smaller proportion of the data (see descriptive stats for these values).
#check 6 num summaries
summary(pcri.e[,c("foc_clp","rd_pcri","nsamp","perc.cheme","cytrt","Snow","IceCover","PrevAUGSecchi")])## foc_clp rd_pcri nsamp perc.cheme
## Min. :0.0000 Min. :1.000 Min. : 8.0 Min. :0.0000
## 1st Qu.:0.2083 1st Qu.:1.087 1st Qu.:106.0 1st Qu.:0.0000
## Median :0.3431 Median :1.286 Median :126.0 Median :0.0000
## Mean :0.3760 Mean :1.562 Mean :149.4 Mean :0.2970
## 3rd Qu.:0.5209 3rd Qu.:1.884 3rd Qu.:165.0 3rd Qu.:0.3777
## Max. :1.0000 Max. :3.947 Max. :430.0 Max. :2.2932
## NA's :38
## cytrt Snow IceCover PrevAUGSecchi
## Min. :0.0000 Min. : 1.50 Min. :104.0 Min. :0.200
## 1st Qu.:0.0000 1st Qu.: 6.20 1st Qu.:122.0 1st Qu.:0.600
## Median :0.0000 Median :11.50 Median :137.0 Median :0.800
## Mean :0.5545 Mean :13.25 Mean :133.4 Mean :1.063
## 3rd Qu.:1.0000 3rd Qu.:19.70 3rd Qu.:143.0 3rd Qu.:1.400
## Max. :5.0000 Max. :39.20 Max. :162.0 Max. :3.500
## summary(pcri.p[,c("foc_clp","rd_pcri","nsamp","perc.lit.chem","cytrt","Snow","IceCover","PrevAUGSecchi")])## foc_clp rd_pcri nsamp perc.lit.chem
## Min. :0.00000 Min. :1.000 Min. : 5.0 Min. :0.0000
## 1st Qu.:0.02878 1st Qu.:1.000 1st Qu.: 91.0 1st Qu.:0.0000
## Median :0.17000 Median :1.359 Median :119.0 Median :0.0000
## Mean :0.23823 Mean :1.614 Mean :143.8 Mean :0.2589
## 3rd Qu.:0.39640 3rd Qu.:2.033 3rd Qu.:156.0 3rd Qu.:0.3219
## Max. :0.97479 Max. :3.358 Max. :653.0 Max. :2.2932
## NA's :65
## cytrt Snow IceCover PrevAUGSecchi
## Min. :0.0000 Min. : 1.1 Min. :103 Min. :0.200
## 1st Qu.:0.0000 1st Qu.: 6.2 1st Qu.:124 1st Qu.:0.500
## Median :0.0000 Median :12.6 Median :137 Median :0.800
## Mean :0.6345 Mean :14.2 Mean :134 Mean :1.181
## 3rd Qu.:1.0000 3rd Qu.:22.6 3rd Qu.:143 3rd Qu.:1.400
## Max. :8.0000 Max. :35.6 Max. :162 Max. :6.900
## #how many cases?
pcri.p[,unique(.(Lake, Year))] # 145## V1 V2
## 1: reshanau 2008
## 2: reshanau 2009
## 3: reshanau 2010
## 4: reshanau 2011
## 5: reshanau 2012
## ---
## 141: clearw 2013
## 142: clearw 2014
## 143: clearw 2015
## 144: stjames 2009
## 145: stjames 2011 pcri.e[,unique(.(Lake, Year))] # 101## V1 V2
## 1: belle 2010
## 2: blueberry 2009
## 3: clearm 2009
## 4: coal 2006
## 5: coal 2007
## ---
## 97: vails 2009
## 98: weaver 2006
## 99: weaver 2007
## 100: weaver 2008
## 101: weaver 2009 #which lakes?
pcri.p[,unique(Lake)] # 45## [1] "reshanau" "riley" "lucy" "minnewashta"
## [5] "susan" "hill" "lowermission" "keller"
## [9] "carlos" "hyland" "mitchell" "round"
## [13] "staring" "gleason" "medicine" "weaver"
## [17] "graystonka" "maxwelltonka" "halsteadtonka" "stubbstonka"
## [21] "christmas" "sarahwest" "saraheast" "rebecca"
## [25] "portage" "carrie" "belle" "clearm"
## [29] "crookneck" "cedar" "silver" "kohlman"
## [33] "crosbyupper" "cleary" "fish" "madison"
## [37] "julia" "rush" "pearl" "vails"
## [41] "coal" "blueberry" "stolaf" "clearw"
## [45] "stjames" pcri.e[,unique(Lake)] # 35## [1] "belle" "blueberry" "clearm" "coal" "como"
## [6] "crookneck" "fish" "gleason" "hill" "hyland"
## [11] "julia" "keller" "kohlman" "loeb" "lowermission"
## [16] "mccarron" "medicine" "mitchell" "pearl" "peltier"
## [21] "portage" "rebecca" "reshanau" "riley" "rush"
## [26] "saraheast" "sarahwest" "silver" "southcenter" "staring"
## [31] "stjames" "stolaf" "susan" "vails" "weaver" #Trt as factor:
str(pcri.e)## Classes 'data.table' and 'data.frame': 101 obs. of 33 variables:
## $ Lake.ID : int 47004900 80003400 47009500 77004600 77004600 77004600 77004600 62005500 49013300 49013300 ...
## $ Year : int 2010 2009 2009 2006 2007 2008 2009 2015 2006 2007 ...
## $ Lake : chr "belle" "blueberry" "clearm" "coal" ...
## $ Latitude.x : num 45 46.8 45.3 46.1 46.1 ...
## $ Trt : int 0 1 1 0 0 0 0 0 1 1 ...
## $ YrsTrt : int 0 2 4 0 0 0 0 0 1 2 ...
## $ PrevAUGSecchi : num 0.7 0.6 0.5 1.8 1.7 2.3 3 1.1 2.3 1.8 ...
## $ IceIn : chr "5-Dec" "16-Nov" "29-Nov" "29-Nov" ...
## $ IceOut : chr "30-Mar" "21-Apr" "13-Apr" "7-Apr" ...
## $ IceCover : int 116 157 136 130 141 157 141 138 135 142 ...
## $ Snow : num 24.3 21.3 13.7 15.4 7 14.2 20.1 4.2 1.5 2.9 ...
## $ Longitude : num -94.4 -95.1 -94.4 -94.8 -94.8 ...
## $ Ecoregion : chr "NCHF" "NCHF" "NCHF" "NLF" ...
## $ Lake.area..ha. : int 334 216 214 69 69 69 69 28 73 73 ...
## $ Littoral.area..ha. : num 155 216 178 29 29 29 29 28 53 53 ...
## $ Proportion.littoral: num 0.464 1 0.832 0.42 0.42 ...
## $ Max.depth..m. : num 7.6 4.6 5.5 7.6 7.6 7.6 7.6 4.6 6.7 6.7 ...
## $ Acres_Chem : num 0 190 60.5 0 0 0 0 0 60 60 ...
## $ Acres_Harvested : num NA NA NA NA NA NA NA 10 NA NA ...
## $ litt.ac : num 383 533.7 439.8 71.7 71.7 ...
## $ perc.lit.chem : num 0 0.356 0.138 0 0 ...
## $ Snow.sc : num 1.20042 0.86291 0.00788 0.19914 -0.74589 ...
## $ PrevAUGSecchi.sc : num -0.516 -0.616 -0.716 0.587 0.487 ...
## $ IceCover.sc : num -1.239 1.633 0.162 -0.258 0.512 ...
## $ nsamp : int 272 396 211 77 78 79 90 60 178 137 ...
## $ nocc_pcri : int 1 108 61 30 36 24 27 47 65 47 ...
## $ rd_pcri : num NA 1.61 1.43 1.33 1.9 ...
## $ datasource : chr "dustin" "newman" "newman" "johnson" ...
## $ foc_clp : num 0.00368 0.27273 0.2891 0.38961 0.46154 ...
## $ cytrt : int 0 0 0 0 0 0 0 0 0 0 ...
## $ yutrt : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Trte : int 0 1 0 0 0 0 0 0 0 1 ...
## $ perc.cheme : num 0 0.337 0 0 0 ...
## - attr(*, ".internal.selfref")=<externalptr> pcri.e[ , Trt := as.factor(Trt)]
pcri.e[ , Trte := as.factor(Trte)]
pcri.p[ , Trt := as.factor(Trt)]
#Lake as a factor:
pcri.e[ , Lake := as.factor(Lake)]
pcri.p[ , Lake := as.factor(Lake)]
#convert inches to centimeters in the Snow column:
pcri.e[,Snow:= Snow*2.54]
pcri.p[,Snow:= Snow*2.54]
#whats in each dataset?
names(pcri.e)## [1] "Lake.ID" "Year" "Lake"
## [4] "Latitude.x" "Trt" "YrsTrt"
## [7] "PrevAUGSecchi" "IceIn" "IceOut"
## [10] "IceCover" "Snow" "Longitude"
## [13] "Ecoregion" "Lake.area..ha." "Littoral.area..ha."
## [16] "Proportion.littoral" "Max.depth..m." "Acres_Chem"
## [19] "Acres_Harvested" "litt.ac" "perc.lit.chem"
## [22] "Snow.sc" "PrevAUGSecchi.sc" "IceCover.sc"
## [25] "nsamp" "nocc_pcri" "rd_pcri"
## [28] "datasource" "foc_clp" "cytrt"
## [31] "yutrt" "Trte" "perc.cheme" names(pcri.p)## [1] "Lake.ID" "Year" "Lake"
## [4] "Latitude.x" "Trt" "YrsTrt"
## [7] "PrevAUGSecchi" "IceIn" "IceOut"
## [10] "IceCover" "Snow" "Longitude"
## [13] "Ecoregion" "Lake.area..ha." "Littoral.area..ha."
## [16] "Proportion.littoral" "Max.depth..m." "Acres_Chem"
## [19] "Acres_Harvested" "litt.ac" "perc.lit.chem"
## [22] "Snow.sc" "PrevAUGSecchi.sc" "IceCover.sc"
## [25] "nsamp" "nocc_pcri" "rd_pcri"
## [28] "datasource" "foc_clp" "cytrt"
## [31] "yutrt"# distances to nearest NWS Station ----------------------------------------
nws_dist <- fread(file = "data/raw/SnowDepth/Lakes_NWS.txt")
lk_dist <- nws_dist[ , .(dowlknum, NEAR_DIST)]
e_dist <- merge(pcri.e, lk_dist, all.x = TRUE, by.x = "Lake.ID", by.y = "dowlknum")
p_dist <- merge(pcri.p, lk_dist, all.x = TRUE, by.x = "Lake.ID", by.y = "dowlknum")
all_dist <- rbind(e_dist, p_dist, fill = TRUE)
distances <- all_dist[, max(NEAR_DIST), .(Lake.ID) ]
summary(distances$V1)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1232 8294 11800 13103 17501 36945 5# descriptive stats -----------------------------------------------------------Summary statistics
# some simple summary stats first two data frames are tables 1 & 2 in the text.
data.summary <- list(
"Table_1_dat" = data.frame("data_type" = rep(c("freqoccurrence", "meanrakedens"), each = 3),
"group" = rep(c("total", "treated", "untreated"), 2),
"n lakes" = c(nrow(unique(rbind(pcri.p[,c("Lake","Year")],pcri.e[,c("Lake","Year")])[,1])), # total foc
nrow(unique(subset(rbind(pcri.p[,c("Lake","Year","Trt")],pcri.e[,c("Lake","Year","Trt")]), Trt == 1)[,1])), #treated foc nlakes
nrow(unique(subset(rbind(pcri.p[,c("Lake","Year","Trt")],pcri.e[,c("Lake","Year","Trt")]), Trt == 0)[,1])), #untrt foc nlakes
nrow(unique(rbind(subset(pcri.e, is.na(pcri.e$rd_pcri) == F)[,c("Lake","Trt")], subset(pcri.p, is.na(pcri.p$rd_pcri) == F)[,c("Lake","Trt")])[,1])), # total rake density nlakes
nrow(unique(subset(rbind(subset(pcri.e, is.na(pcri.e$rd_pcri) == F)[,c("Lake","Trt")], subset(pcri.p, is.na(pcri.p$rd_pcri) == F)[,c("Lake","Trt")]), Trt == 1))), # treated rake density nlakes
nrow(unique(subset(rbind(subset(pcri.e, is.na(pcri.e$rd_pcri) == F)[,c("Lake","Trt")], subset(pcri.p, is.na(pcri.p$rd_pcri) == F)[,c("Lake","Trt")]), Trt == 0))) # untrt rake dens nlakes
),
"n lake-years" = c(nrow(unique(rbind(pcri.p[,c("Lake","Year","Trt")],pcri.e[,c("Lake","Year","Trt")]))[,1]),# foc total n lakeyr
nrow(unique(subset(rbind(pcri.p[,c("Lake","Year","Trt")],pcri.e[,c("Lake","Year","Trt")]), Trt == 1))[,1]), # trt foc n lakeyr
nrow(unique(subset(rbind(pcri.p[,c("Lake","Year","Trt")],pcri.e[,c("Lake","Year","Trt")]), Trt == 0))[,1]), # untrt foc n lakeyr
nrow(unique(rbind(subset(pcri.e, is.na(pcri.e$rd_pcri) == F)[,c("Lake","Year")], subset(pcri.p, is.na(pcri.p$rd_pcri) == F)[,c("Lake","Year")]))),# rd total n lakeyr
nrow(unique(rbind(subset(pcri.e, is.na(pcri.e$rd_pcri) == F & pcri.e$Trt == 1)[,c("Lake","Year")], subset(pcri.p, is.na(pcri.p$rd_pcri) == F & pcri.p$Trt == 1)[,c("Lake","Year")]))),# rd trt n lakeyr
nrow(unique(rbind(subset(pcri.e, is.na(pcri.e$rd_pcri) == F & pcri.e$Trt == 0)[,c("Lake","Year")], subset(pcri.p, is.na(pcri.p$rd_pcri) == F & pcri.p$Trt == 0)[,c("Lake","Year")])))# rd untrt n lakeyr
),
"n late surveys" = c(nrow(unique(pcri.p[,c("Lake","Year")])),# foc total nsurveys
nrow(unique(pcri.p[Trt == 1,c("Lake","Year")])),# foc trt nsurveys
nrow(unique(pcri.p[Trt == 0,c("Lake","Year")])),# foc untrt nsurveys
nrow(pcri.p[is.na(pcri.p$rd_pcri) == F,c("Lake","Year")]),# rd total nsurveys
nrow(pcri.p[is.na(pcri.p$rd_pcri) == F & Trt == 1,c("Lake","Year")]),# rd trt nsurveys
nrow(pcri.p[is.na(pcri.p$rd_pcri) == F & Trt == 0,c("Lake","Year")])# rd untrt nsurveys
),
"n early surveys" = c(nrow(unique(pcri.e[,c("Lake","Year")])[,1]),# foc total nsurveys
nrow(unique(pcri.e[Trte == 1,c("Lake","Year")])[,1]),# foc trt nsurveys
nrow(unique(pcri.e[Trte == 0,c("Lake","Year")])[,1]),# foc untrt nsurveys
nrow(pcri.e[is.na(pcri.e$rd_pcri) == F ,c("Lake","Year")]),# rd total nsurveys
nrow(pcri.e[is.na(pcri.e$rd_pcri) == F & Trte == 1,c("Lake","Year")]),# rd trt nsurveys
nrow(pcri.e[is.na(pcri.e$rd_pcri) == F & Trte == 0,c("Lake","Year")]))# rd untrt nsurveys
),
"variables_summary" = data.frame(
"group" = c(rep(c("trt", "untrt"), each = 2)),
"season" = c(rep(c("early", "late"),2)),
"foc_mean" = c(mean(subset(pcri.e, Trte == 1)[,foc_clp]),
mean(subset(pcri.p, Trt == 1)[,foc_clp]),
mean(subset(pcri.e, Trte == 0)[,foc_clp]),
mean(subset(pcri.p, Trt == 0)[,foc_clp])),
"foc_sd" = c(sd(subset(pcri.e, Trte == 1)[,foc_clp]),
sd(subset(pcri.p, Trt == 1)[,foc_clp]),
sd(subset(pcri.e, Trte == 0)[,foc_clp]),
sd(subset(pcri.p, Trt == 0)[,foc_clp])),
"rd_mean" = c(mean(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
mean(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
mean(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
mean(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri])),
"rd_sd" = c(sd(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
sd(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
sd(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
sd(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri])),
"cytrt_mean" = c(mean(subset(pcri.e, Trte == 1)[,cytrt]),
mean(subset(pcri.p, Trt == 1)[,cytrt]),
NA,
NA),
"cytrt_sd" = c(sd(subset(pcri.e, Trte == 1)[,cytrt]),
sd(subset(pcri.p, Trt == 1)[,cytrt]),
NA,
NA),
"Secchi_mean_meter" = c(mean(subset(pcri.e, Trte == 1)[,PrevAUGSecchi]),
mean(subset(pcri.p, Trt == 1)[,PrevAUGSecchi]),
mean(subset(pcri.e, Trte == 0)[,PrevAUGSecchi]),
mean(subset(pcri.p, Trt == 0)[,PrevAUGSecchi])),
"Secchi_sd_meter" = c(sd(subset(pcri.e, Trte == 1)[,PrevAUGSecchi]),
sd(subset(pcri.p, Trt == 1)[,PrevAUGSecchi]),
sd(subset(pcri.e, Trte == 0)[,PrevAUGSecchi]),
sd(subset(pcri.p, Trt == 0)[,PrevAUGSecchi])),
"Snow_mean_in" = c(mean(subset(pcri.e, Trte == 1)[,Snow]),
mean(subset(pcri.p, Trt == 1)[,Snow]),
mean(subset(pcri.e, Trte == 0)[,Snow]),
mean(subset(pcri.p, Trt == 0)[,Snow])),
"Snow_sd_in" = c(sd(subset(pcri.e, Trte == 1)[,Snow]),
sd(subset(pcri.p, Trt == 1)[,Snow]),
sd(subset(pcri.e, Trte == 0)[,Snow]),
sd(subset(pcri.p, Trt == 0)[,Snow])),
"ice_mean_days" = c(mean(subset(pcri.e, Trte == 1)[,IceCover]),
mean(subset(pcri.p, Trt == 1)[,IceCover]),
mean(subset(pcri.e, Trte == 0)[,IceCover]),
mean(subset(pcri.p, Trt == 0)[,IceCover])),
"ice_sd_days" = c(sd(subset(pcri.e, Trte == 1)[,IceCover]),
sd(subset(pcri.p, Trt == 1)[,IceCover]),
sd(subset(pcri.e, Trte == 0)[,IceCover]),
sd(subset(pcri.p, Trt == 0)[,IceCover]))
),
"Descriptive_stats" = data.frame(
"Group" = rep(c("treated", "untreated"), each = 4),
"Metric" = rep(c("early rd", "peak rd", "early foc", "peak foc"), 2),
"Mean" = c(
mean(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
mean(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
mean(subset(pcri.e, Trte == 1)[,foc_clp]),
mean(subset(pcri.p, Trt == 1)[,foc_clp]),
mean(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
mean(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri]),
mean(subset(pcri.e, Trte == 0)[,foc_clp]),
mean(subset(pcri.p, Trt == 0)[,foc_clp])
),
"stdev" = c(
sd(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
sd(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
sd(subset(pcri.e, Trte == 1)[,foc_clp]),
sd(subset(pcri.p, Trt == 1)[,foc_clp]),
sd(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
sd(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri]),
sd(subset(pcri.e, Trte == 0)[,foc_clp]),
sd(subset(pcri.p, Trt == 0)[,foc_clp])
),
"Median" = c(
median(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
median(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
median(subset(pcri.e, Trte == 1)[,foc_clp]),
median(subset(pcri.p, Trt == 1)[,foc_clp]),
median(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
median(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri]),
median(subset(pcri.e, Trte == 0)[,foc_clp]),
median(subset(pcri.p, Trt == 0)[,foc_clp])
),
"Quartile1st" = c(
quantile(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri])[2],
quantile(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri])[2],
quantile(subset(pcri.e, Trte == 1)[,foc_clp])[2],
quantile(subset(pcri.p, Trt == 1)[,foc_clp])[2],
quantile(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri])[2],
quantile(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri])[2],
quantile(subset(pcri.e, Trte == 0)[,foc_clp])[2],
quantile(subset(pcri.p, Trt == 0)[,foc_clp])[2]
),
"Quartile3rd" = c(
quantile(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri])[4],
quantile(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri])[4],
quantile(subset(pcri.e, Trte == 1)[,foc_clp])[4],
quantile(subset(pcri.p, Trt == 1)[,foc_clp])[4],
quantile(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri])[4],
quantile(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri])[4],
quantile(subset(pcri.e, Trte == 0)[,foc_clp])[4],
quantile(subset(pcri.p, Trt == 0)[,foc_clp])[4]
),
"min" = c(
min(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
min(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
min(subset(pcri.e, Trte == 1)[,foc_clp]),
min(subset(pcri.p, Trt == 1)[,foc_clp]),
min(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
min(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri]),
min(subset(pcri.e, Trte == 0)[,foc_clp]),
min(subset(pcri.p, Trt == 0)[,foc_clp])
),
"max" = c(
max(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,rd_pcri]),
max(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,rd_pcri]),
max(subset(pcri.e, Trte == 1)[,foc_clp]),
max(subset(pcri.p, Trt == 1)[,foc_clp]),
max(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,rd_pcri]),
max(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,rd_pcri]),
max(subset(pcri.e, Trte == 0)[,foc_clp]),
max(subset(pcri.p, Trt == 0)[,foc_clp])
)
),
"nsamples_per_survey" = list(
"foc" = data.frame( "grp" = c("trt", "untrt"),
"late mean" = c(mean(subset(pcri.p, Trt == 1)[,nsamp]),
mean(subset(pcri.p, Trt == 0)[,nsamp])),
"late sd" = c(sd(subset(pcri.p, Trt == 1)[,nsamp]),
sd(subset(pcri.p, Trt == 0)[,nsamp])),
"early mean" = c(mean(subset(pcri.e, Trte == 1)[,nsamp]),
mean(subset(pcri.e, Trte == 0)[,nsamp])),
"early sd" = c(sd(subset(pcri.e, Trt == 1)[,nsamp]),
sd(subset(pcri.e, Trt == 0)[,nsamp]))
),
"rd"= data.frame( "grp" = c("trt", "untrt"),
"late mean" = c(mean(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,nocc_pcri]),
mean(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,nocc_pcri])),
"late sd" = c(sd(subset(pcri.p, is.na(rd_pcri) == F & Trt == 1)[,nocc_pcri]),
sd(subset(pcri.p, is.na(rd_pcri) == F & Trt == 0)[,nocc_pcri])),
"early mean" = c(mean(subset(pcri.e, is.na(rd_pcri) == F & Trte == 1)[,nocc_pcri]),
mean(subset(pcri.e, is.na(rd_pcri) == F & Trte == 0)[,nocc_pcri])),
"early sd" = c(sd(subset(pcri.e, is.na(rd_pcri) == F & Trt == 1)[,nocc_pcri]),
sd(subset(pcri.e, is.na(rd_pcri) == F & Trt == 0)[,nocc_pcri]))
))
)
data.summary## $Table_1_dat
## data_type group n.lakes n.lake.years n.late.surveys n.early.surveys
## 1 freqoccurrence total 50 177 145 101
## 2 freqoccurrence treated 22 73 63 48
## 3 freqoccurrence untreated 42 104 82 53
## 4 meanrakedens total 28 91 80 63
## 5 meanrakedens treated 16 49 40 35
## 6 meanrakedens untreated 20 42 40 28
##
## $variables_summary
## group season foc_mean foc_sd rd_mean rd_sd cytrt_mean cytrt_sd
## 1 trt early 0.2889259 0.2151302 1.250231 0.2650372 1.166667 1.243422
## 2 trt late 0.1315744 0.1336390 1.323346 0.5436439 1.460317 1.730424
## 3 untrt early 0.4548323 0.2453204 1.952101 0.8441731 NA NA
## 4 untrt late 0.3201809 0.2676298 1.904123 0.7056931 NA NA
## Secchi_mean_meter Secchi_sd_meter Snow_mean_in Snow_sd_in ice_mean_days
## 1 1.0020833 0.6446142 35.45946 20.73935 134.8125
## 2 0.8603175 0.5734839 32.63698 21.40893 134.0000
## 3 1.1188679 0.7465450 32.01886 22.15263 132.1132
## 4 1.4280488 1.2737884 38.68761 23.85164 134.0000
## ice_sd_days
## 1 14.94480
## 2 13.54562
## 3 14.17563
## 4 14.08440
##
## $Descriptive_stats
## Group Metric Mean stdev Median Quartile1st Quartile3rd
## 1 treated early rd 1.2502307 0.2650372 1.15000000 1.05777778 1.3333333
## 2 treated peak rd 1.3233459 0.5436439 1.00000000 1.00000000 1.4136364
## 3 treated early foc 0.2889259 0.2151302 0.26862170 0.11381579 0.4248610
## 4 treated peak foc 0.1315744 0.1336390 0.09090909 0.01222421 0.2072750
## 5 untreated early rd 1.9521012 0.8441731 1.88401328 1.24695122 2.4353448
## 6 untreated peak rd 1.9041232 0.7056931 1.78409091 1.33333333 2.4637097
## 7 untreated early foc 0.4548323 0.2453204 0.43462898 0.29931973 0.5714286
## 8 untreated peak foc 0.3201809 0.2676298 0.26836085 0.08012300 0.4609252
## min max
## 1 1.000000000 1.9315068
## 2 1.000000000 3.0000000
## 3 0.000000000 0.9516129
## 4 0.000000000 0.5232558
## 5 1.000000000 3.9473684
## 6 1.000000000 3.3577236
## 7 0.003676471 1.0000000
## 8 0.000000000 0.9747899
##
## $nsamples_per_survey
## $nsamples_per_survey$foc
## grp late.mean late.sd early.mean early.sd
## 1 trt 120.1270 48.76557 139.2292 55.00702
## 2 untrt 161.9878 139.44223 158.6038 122.42867
##
## $nsamples_per_survey$rd
## grp late.mean late.sd early.mean early.sd
## 1 trt 17.925 18.06170 41.08571 30.65920
## 2 untrt 43.250 38.83149 63.00000 34.86515Export Tables 1 & 2
# export data for Tables in manuscript
# write.csv(data.summary$variables_summary, file = "data/output_data/summary_stats_Tbl2.csv")
# write.csv(data.summary$Table_1_dat, file = "data/output_data/data_summary_Tbl1.csv")
# mixed eff models --------------------------------------------------------Mixed effects models
No interactions
Here we will start our modeling without considering any interactionsusing a glmme eval. First, we’ll need to add an observation level random effect to account for overdispersion (Harrsion, 2015). The following section uses centered & scaled predictor variables (done using function scale()) in order to converge during fitting. glmme for frequency of occurrence data:
# add obs level key to dataset
pcri.e[,obs := seq(nrow(pcri.e))]
pcri.p[,obs := seq(nrow(pcri.p))]Early Survey Freq
glme.foc.e <- glmer(foc_clp ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc + (1 | Lake) + (1 | obs), family = binomial(logit), data = pcri.e, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! summary(glme.foc.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake) + (1 | obs)
## Data: pcri.e
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 914.0 934.9 -449.0 898.0 93
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.36412 -0.12555 -0.00206 0.10055 1.18301
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6224 0.7889
## Lake (Intercept) 0.9775 0.9887
## Number of obs: 101, groups: obs, 101; Lake, 35
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.564845 0.216086 -2.614 0.00895 **
## Trte1 -0.219452 0.236166 -0.929 0.35277
## cytrt -0.529424 0.116909 -4.529 5.94e-06 ***
## PrevAUGSecchi.sc -0.464942 0.204855 -2.270 0.02323 *
## Snow.sc -0.508611 0.113872 -4.467 7.95e-06 ***
## IceCover.sc -0.008906 0.103509 -0.086 0.93143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.349
## cytrt -0.029 -0.441
## PrvAUGScch. 0.124 0.055 -0.267
## Snow.sc 0.061 -0.111 -0.097 0.015
## IceCover.sc 0.087 -0.049 -0.050 0.080 -0.305 plot(glme.foc.e) # not too bad
plot(glme.foc.e, resid(.)~ fitted(.) |Lake) # by group
plot(glme.foc.e, resid(.)~ Snow) # residuals versus Snow
plot(glme.foc.e, resid(.)~ PrevAUGSecchi) #secchi
plot(glme.foc.e, resid(.)~ IceCover) # icecover
plot(glme.foc.e, resid(.)~ cytrt) # cytrt - def more var at low cytrt. That is also where most of our data lie
Late Survey Freq
glme.foc.p <- glmer(foc_clp~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.p, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! summary(glme.foc.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 1201.1 1224.9 -592.6 1185.1 137
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38306 -0.20861 0.00594 0.09084 0.93707
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.794 1.339
## Lake (Intercept) 1.882 1.372
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.11214 0.27885 -3.988 6.65e-05 ***
## Trt1 -1.58941 0.36972 -4.299 1.72e-05 ***
## cytrt -0.18458 0.14372 -1.284 0.1990
## PrevAUGSecchi.sc -0.44663 0.18191 -2.455 0.0141 *
## Snow.sc 0.07494 0.13829 0.542 0.5879
## IceCover.sc -0.13352 0.14893 -0.897 0.3700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc
## Trt1 -0.407
## cytrt -0.021 -0.386
## PrvAUGScch. -0.135 0.193 -0.054
## Snow.sc -0.086 0.102 -0.084 -0.028
## IceCover.sc 0.077 -0.193 -0.027 -0.037 -0.347 plot(glme.foc.p) # this also loks o.k. but note the tails again are a bit wonky 
plot(glme.foc.p, resid(.)~ fitted(.) |Lake) # by group
plot(glme.foc.p, resid(.)~ Snow) # residuals versus Snow
plot(glme.foc.p, resid(.)~ PrevAUGSecchi) #secchi- here we can see that the secchi data have relatively poor coverage past ~3.5m
plot(glme.foc.p, resid(.)~ IceCover) # icecover
plot(glme.foc.p, resid(.)~ cytrt) # cytrt
Early Survey Rake Density:
glme.rd.e <- glmer(rd_pcri~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
(1 | Lake), family = Gamma(inverse), data = pcri.e)
summary(glme.rd.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( inverse )
## Formula: rd_pcri ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.e
##
## AIC BIC logLik deviance df.resid
## 41.8 58.9 -12.9 25.8 55
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.66066 -0.45020 0.02484 0.72711 1.82951
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.01357 0.1165
## Residual 0.03853 0.1963
## Number of obs: 63, groups: Lake, 19
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.64090 0.05345 11.991 < 2e-16 ***
## Trte1 0.14935 0.03965 3.767 0.000165 ***
## cytrt 0.03142 0.01989 1.579 0.114226
## PrevAUGSecchi.sc 0.05625 0.03573 1.574 0.115423
## Snow.sc 0.01794 0.02046 0.876 0.380789
## IceCover.sc 0.02648 0.01764 1.501 0.133369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.266
## cytrt -0.037 -0.456
## PrvAUGScch. 0.216 -0.044 -0.194
## Snow.sc 0.147 -0.045 -0.092 0.111
## IceCover.sc -0.065 -0.130 0.071 -0.110 -0.611 plot(glme.rd.e) # this looks o.k. but note the lack of data at large fitted vals 
Late Survey Rake Densities
glme.rd.p <- glmer(rd_pcri~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
(1 | Lake), family = Gamma(inverse), data = pcri.p)
summary(glme.rd.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( inverse )
## Formula: rd_pcri ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.p
##
## AIC BIC logLik deviance df.resid
## 116.2 135.3 -50.1 100.2 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8620 -0.5004 -0.2257 0.6633 2.3714
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.01745 0.1321
## Residual 0.08545 0.2923
## Number of obs: 80, groups: Lake, 28
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 5.649e-01 5.400e-02 10.461 < 2e-16 ***
## Trt1 2.612e-01 6.232e-02 4.191 2.78e-05 ***
## cytrt 2.881e-02 2.137e-02 1.349 0.1775
## PrevAUGSecchi.sc 2.072e-02 2.895e-02 0.716 0.4740
## Snow.sc 4.676e-02 2.404e-02 1.945 0.0518 .
## IceCover.sc -7.778e-06 2.551e-02 0.000 0.9998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc
## Trt1 -0.391
## cytrt -0.005 -0.330
## PrvAUGScch. -0.218 0.159 -0.098
## Snow.sc -0.097 0.316 -0.088 0.015
## IceCover.sc 0.090 -0.335 0.043 -0.031 -0.574 plot(glme.rd.p) # probably the best so far, but there is a weird lower bound on this one
# mixed models with interactions by trt ------------------------------------------------Interactions (trt X env. var.)
In order to evaluate whether or not chemical treatements can “harness” the natural fluctuations of CLP, we need to test for interactions between treatment and our environmental variables. Our approach is thus to test for interactions, but retain only significant (p = .05) interactions for estimation of fixed effects.
Late Surveys Freq
#peak foc
glme.foc.p <- glmer(foc_clp~ Trt + cytrt + Trt*PrevAUGSecchi.sc + Trt*Snow.sc + Trt*IceCover.sc +
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.p, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.01035 (tol = 0.001, component 1) summary(glme.foc.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + Trt * PrevAUGSecchi.sc + Trt * Snow.sc +
## Trt * IceCover.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 1196.0 1228.7 -587.0 1174.0 134
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38820 -0.19301 -0.00007 0.09729 0.88716
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.627 1.276
## Lake (Intercept) 1.737 1.318
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03954 0.26826 -3.875 0.000107 ***
## Trt1 -1.07594 0.39012 -2.758 0.005816 **
## cytrt -0.29859 0.14546 -2.053 0.040101 *
## PrevAUGSecchi.sc -0.54937 0.18678 -2.941 0.003269 **
## Snow.sc -0.15291 0.16653 -0.918 0.358504
## IceCover.sc -0.08506 0.18113 -0.470 0.638646
## Trt1:PrevAUGSecchi.sc 1.14252 0.48057 2.377 0.017434 *
## Trt1:Snow.sc 0.75621 0.28976 2.610 0.009059 **
## Trt1:IceCover.sc -0.11486 0.29774 -0.386 0.699666
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU Tr1:S.
## Trt1 -0.337
## cytrt -0.039 -0.431
## PrvAUGScch. -0.137 0.055 -0.026
## Snow.sc -0.080 -0.008 0.031 -0.088
## IceCover.sc 0.025 -0.080 0.001 -0.091 -0.262
## Trt1:PAUGS. 0.061 0.390 -0.154 -0.318 0.049 0.058
## Trt1:Snw.sc 0.025 0.208 -0.194 0.051 -0.595 0.133 0.084
## Trt1:IcCvr. 0.055 -0.059 -0.032 0.049 0.158 -0.607 0.032 -0.338
## convergence code: 0
## Model failed to converge with max|grad| = 0.01035 (tol = 0.001, component 1)Why the warnings? There are a few weird values in there, thus “fit with non-integer data”. These are likely the lakes where I used Adam’s calculated values.
pcri.p[ , nsamp*foc_clp] # note non-integer values## [1] 0.00 2.00 95.00 22.00 4.00 7.00 0.00 54.00 44.00 31.00
## [11] 35.00 18.00 28.00 50.00 18.00 60.86 35.00 37.00 39.00 26.00
## [21] 10.00 11.00 38.00 18.00 57.00 63.00 28.00 31.00 0.00 0.00
## [31] 2.00 116.00 100.00 11.00 0.00 0.00 122.00 85.00 3.00 13.00
## [41] 9.00 12.00 29.00 65.00 1.00 1.00 1.00 5.00 19.00 2.00
## [51] 0.00 3.00 0.00 25.00 45.00 35.00 2.00 0.00 1.00 39.00
## [61] 19.00 54.00 9.00 26.00 44.00 36.00 5.04 77.00 40.00 27.00
## [71] 8.00 36.00 51.00 35.00 24.00 2.00 114.00 126.00 5.00 33.00
## [81] 62.00 6.00 0.00 0.00 1.00 8.00 13.00 2.00 23.00 25.00
## [91] 37.00 22.00 0.00 0.00 0.00 62.78 34.00 69.00 101.00 1.00
## [101] 1.00 32.00 4.00 17.92 0.00 2.00 14.00 58.00 0.00 120.00
## [111] 94.00 99.00 47.00 24.00 20.00 28.00 264.00 308.00 0.00 0.00
## [121] 17.00 5.00 32.00 3.00 1.00 2.00 180.00 180.00 77.00 38.00
## [131] 50.00 35.00 36.00 72.00 21.00 6.00 5.00 8.00 3.00 7.00
## [141] 3.00 7.00 14.00 26.00 7.00Why non-convergence? Explore a few reasons (adapted from https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b615d8633c01d51.html):
# a singularity?
tt <- getME(glme.foc.p,"theta")
ll <- getME(glme.foc.p,"lower")
min(tt[ll==0]) # likely not because of a singularity## [1] 1.275533 # a gradient calculation error?
derivs1 <- glme.foc.p@optinfo$derivs
sc_grad1 <- with(derivs1,solve(Hessian,gradient))
max(abs(sc_grad1))## [1] 0.002362967 max(pmin(abs(sc_grad1),abs(derivs1$gradient)))## [1] 0.002362967 dd <- update(glme.foc.p,devFunOnly=TRUE)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! pars <- unlist(getME(glme.foc.p,c("theta","fixef")))
grad2 <- grad(dd,pars)
hess2 <- hessian(dd,pars)
sc_grad2 <- solve(hess2,grad2)
max(pmin(abs(sc_grad2),abs(grad2))) # not lack of gradient (usually set to require > 0.001)## [1] 0.00235688 # Restart from previous fit:
ss <- getME(glme.foc.p,c("theta","fixef"))
glme.foc.p2 <- update(glme.foc.p,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! # that worked
summary(glme.foc.p2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + Trt * PrevAUGSecchi.sc + Trt * Snow.sc +
## Trt * IceCover.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
## Control: glmerControl(optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 1196.0 1228.7 -587.0 1174.0 134
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38863 -0.19293 -0.00005 0.09743 0.88757
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.626 1.275
## Lake (Intercept) 1.739 1.319
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.03967 0.26834 -3.874 0.000107 ***
## Trt1 -1.07719 0.39005 -2.762 0.005750 **
## cytrt -0.29873 0.14542 -2.054 0.039955 *
## PrevAUGSecchi.sc -0.54976 0.18684 -2.942 0.003257 **
## Snow.sc -0.15325 0.16648 -0.921 0.357308
## IceCover.sc -0.08469 0.18108 -0.468 0.639978
## Trt1:PrevAUGSecchi.sc 1.14019 0.48053 2.373 0.017654 *
## Trt1:Snow.sc 0.75560 0.28966 2.609 0.009092 **
## Trt1:IceCover.sc -0.11417 0.29765 -0.384 0.701292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU Tr1:S.
## Trt1 -0.337
## cytrt -0.039 -0.430
## PrvAUGScch. -0.137 0.055 -0.026
## Snow.sc -0.079 -0.008 0.031 -0.088
## IceCover.sc 0.025 -0.080 0.001 -0.091 -0.262
## Trt1:PAUGS. 0.061 0.390 -0.154 -0.318 0.049 0.058
## Trt1:Snw.sc 0.025 0.208 -0.194 0.051 -0.595 0.133 0.084
## Trt1:IcCvr. 0.055 -0.059 -0.032 0.049 0.158 -0.607 0.032 -0.338 #then retain significant interactions (PAUGSecchi & Snow):
glme.foc.p <- glmer(foc_clp~ Trt + cytrt + PrevAUGSecchi.sc+ Snow.sc + IceCover.sc + Trt:PrevAUGSecchi.sc + Trt:Snow.sc+
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.p, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! #model.matrix(glme.foc.p)
summary(glme.foc.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## Trt:PrevAUGSecchi.sc + Trt:Snow.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 1194.1 1223.9 -587.1 1174.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38658 -0.19658 0.00133 0.09733 0.89071
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.629 1.276
## Lake (Intercept) 1.728 1.315
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0343 0.2676 -3.865 0.000111 ***
## Trt1 -1.0859 0.3897 -2.786 0.005329 **
## cytrt -0.3002 0.1453 -2.067 0.038746 *
## PrevAUGSecchi.sc -0.5457 0.1862 -2.930 0.003388 **
## Snow.sc -0.1430 0.1644 -0.870 0.384452
## IceCover.sc -0.1270 0.1439 -0.883 0.377422
## Trt1:PrevAUGSecchi.sc 1.1458 0.4801 2.387 0.016993 *
## Trt1:Snow.sc 0.7181 0.2729 2.632 0.008500 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU
## Trt1 -0.336
## cytrt -0.037 -0.434
## PrvAUGScch. -0.141 0.059 -0.024
## Snow.sc -0.089 0.000 0.036 -0.097
## IceCover.sc 0.075 -0.147 -0.023 -0.077 -0.211
## Trt1:PAUGS. 0.060 0.391 -0.153 -0.320 0.044 0.098
## Trt1:Snw.sc 0.045 0.202 -0.217 0.072 -0.583 -0.098 0.101So we can see that the secchi and snow slopes are different in Trt v Untrt lakes, but but now we want to ask if the values for the Trt interactions are significantly non-zero (vs. sig diff from UNtrt). To do this we’ll flip our intercept value, and estimate the slopes for treated lakes first:
str(pcri.p)## Classes 'data.table' and 'data.frame': 145 obs. of 32 variables:
## $ Lake.ID : int 2000900 2000900 2000900 2000900 2000900 2000900 2000900 10000200 10000200 10000200 ...
## $ Year : int 2008 2009 2010 2011 2012 2013 2014 2011 2012 2013 ...
## $ Lake : Factor w/ 45 levels "belle","blueberry",..: 32 32 32 32 32 32 32 33 33 33 ...
## $ Latitude.x : num 45.2 45.2 45.2 45.2 45.2 ...
## $ Trt : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 1 1 2 ...
## $ YrsTrt : int 4 5 0 1 2 0 1 0 0 1 ...
## $ PrevAUGSecchi : num 0.5 0.3 0.4 0.5 0.5 0.3 0.3 1 0.9 0.8 ...
## $ IceIn : chr "28-Nov" "22-Nov" "6-Dec" "24-Nov" ...
## $ IceOut : chr "22-Apr" "9-Apr" "31-Mar" "12-Apr" ...
## $ IceCover : int 147 139 116 140 110 143 152 137 106 147 ...
## $ Snow : num 27.18 32.77 57.4 78.99 5.33 ...
## $ Longitude : num -93.1 -93.1 -93.1 -93.1 -93.1 ...
## $ Ecoregion : chr "NCHF" "NCHF" "NCHF" "NCHF" ...
## $ Lake.area..ha. : num 136 136 136 136 136 136 136 120 120 120 ...
## $ Littoral.area..ha. : num 136 136 136 136 136 136 136 44 44 44 ...
## $ Proportion.littoral: num 1 1 1 1 1 ...
## $ Max.depth..m. : num 3 3 3 3 3 3 3 14.9 14.9 14.9 ...
## $ Acres_Chem : num 282 282 0 49 48 48 48 0 0 35 ...
## $ Acres_Harvested : num NA NA NA NA NA NA NA NA NA NA ...
## $ litt.ac : num 336 336 336 336 336 ...
## $ perc.lit.chem : num 0.839 0.839 0 0.146 0.143 ...
## $ Snow.sc : num -0.3296 -0.0821 1.0092 1.9654 -1.2972 ...
## $ PrevAUGSecchi.sc : num -0.716 -0.917 -0.817 -0.716 -0.716 ...
## $ IceCover.sc : num 0.933 0.372 -1.239 0.442 -1.659 ...
## $ nsamp : int 130 14 127 130 129 72 47 167 171 169 ...
## $ nocc_pcri : int 0 2 95 22 4 7 0 54 44 31 ...
## $ rd_pcri : num NA 1 2.74 1.05 1 ...
## $ datasource : chr "mccomas" "mccomas" "mccomas" "mccomas" ...
## $ foc_clp : num 0 0.143 0.748 0.169 0.031 ...
## $ cytrt : int 2 3 0 0 1 2 3 0 0 0 ...
## $ yutrt : int 0 0 1 0 0 0 0 0 0 0 ...
## $ obs : int 1 2 3 4 5 6 7 8 9 10 ...
## - attr(*, ".internal.selfref")=<externalptr> pcri.p$TrtN <- as.numeric(pcri.p$Trt == 0)
pcri.p[,.(Trt,TrtN)]## Trt TrtN
## 1: 1 0
## 2: 1 0
## 3: 0 1
## 4: 1 0
## 5: 1 0
## ---
## 141: 0 1
## 142: 0 1
## 143: 0 1
## 144: 0 1
## 145: 0 1 glme.foc.p.a <- glmer(foc_clp ~
TrtN + cytrt +
+ IceCover.sc +
TrtN*PrevAUGSecchi.sc + TrtN*Snow.sc+
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.p, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00984756 (tol = 0.001, component 1) # model.matrix(glme.foc.p.a)
# Restart from previous fit:
ss <- getME(glme.foc.p.a ,c("theta","fixef"))
glme.foc.p.a2 <- update(glme.foc.p.a,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! summary(glme.foc.p.a2) # ests for treated slopes as intercepts## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ TrtN + cytrt + +IceCover.sc + TrtN * PrevAUGSecchi.sc +
## TrtN * Snow.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
## Control: glmerControl(optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 1194.1 1223.9 -587.1 1174.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38657 -0.19658 0.00133 0.09733 0.89070
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.629 1.276
## Lake (Intercept) 1.728 1.315
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1203 0.3917 -5.414 6.18e-08 ***
## TrtN 1.0859 0.3897 2.787 0.00533 **
## cytrt -0.3002 0.1453 -2.067 0.03877 *
## IceCover.sc -0.1270 0.1439 -0.883 0.37744
## PrevAUGSecchi.sc 0.6001 0.4560 1.316 0.18818
## Snow.sc 0.5751 0.2218 2.593 0.00951 **
## TrtN:PrevAUGSecchi.sc -1.1458 0.4801 -2.387 0.01699 *
## TrtN:Snow.sc -0.7181 0.2729 -2.631 0.00850 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtN cytrt IcCvr. PAUGS. Snw.sc TN:PAU
## TrtN -0.765
## cytrt -0.457 0.434
## IceCover.sc -0.095 0.147 -0.023
## PrvAUGScch. 0.438 -0.436 -0.171 0.071
## Snow.sc 0.240 -0.248 -0.240 -0.277 0.172
## TrtN:PAUGS. -0.430 0.391 0.153 -0.098 -0.922 -0.157
## TrtN:Snw.sc -0.232 0.202 0.217 0.098 -0.136 -0.798 0.101 summary(glme.foc.p) # ests for untreated slopes as intercepts## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## Trt:PrevAUGSecchi.sc + Trt:Snow.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 1194.1 1223.9 -587.1 1174.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38658 -0.19658 0.00133 0.09733 0.89071
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.629 1.276
## Lake (Intercept) 1.728 1.315
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0343 0.2676 -3.865 0.000111 ***
## Trt1 -1.0859 0.3897 -2.786 0.005329 **
## cytrt -0.3002 0.1453 -2.067 0.038746 *
## PrevAUGSecchi.sc -0.5457 0.1862 -2.930 0.003388 **
## Snow.sc -0.1430 0.1644 -0.870 0.384452
## IceCover.sc -0.1270 0.1439 -0.883 0.377422
## Trt1:PrevAUGSecchi.sc 1.1458 0.4801 2.387 0.016993 *
## Trt1:Snow.sc 0.7181 0.2729 2.632 0.008500 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU
## Trt1 -0.336
## cytrt -0.037 -0.434
## PrvAUGScch. -0.141 0.059 -0.024
## Snow.sc -0.089 0.000 0.036 -0.097
## IceCover.sc 0.075 -0.147 -0.023 -0.077 -0.211
## Trt1:PAUGS. 0.060 0.391 -0.153 -0.320 0.044 0.098
## Trt1:Snw.sc 0.045 0.202 -0.217 0.072 -0.583 -0.098 0.101So secchi is non-significant in the treated lakes, whereas snow has a significantly positive effect on clp growth in the treated subset.From these models we can say that the influence on peak clp freq are: trt(-), cytrt (-), Prev Secchi (- in untrt lakes, non-sig in trt lakes), and Snow (non-sig in untrt lakes, but + in Trt lakes), ice cover (non-sig).
Early Surveys Freq
#early foc
glme.foc.e <- glmer(foc_clp~ Trte + cytrt + PrevAUGSecchi.sc*Trte + Snow.sc*Trte + IceCover.sc*Trte +
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.e, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! summary(glme.foc.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trte + cytrt + PrevAUGSecchi.sc * Trte + Snow.sc *
## Trte + IceCover.sc * Trte + (1 | Lake) + (1 | obs)
## Data: pcri.e
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 918.2 947.0 -448.1 896.2 90
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.30810 -0.13048 0.00094 0.11270 1.24141
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6120 0.7823
## Lake (Intercept) 0.9521 0.9757
## Number of obs: 101, groups: obs, 101; Lake, 35
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5424718 0.2149516 -2.524 0.01161 *
## Trte1 -0.3303027 0.2597987 -1.271 0.20359
## cytrt -0.4966936 0.1204203 -4.125 3.71e-05 ***
## PrevAUGSecchi.sc -0.3539969 0.2273198 -1.557 0.11941
## Snow.sc -0.4924664 0.1522249 -3.235 0.00122 **
## IceCover.sc 0.0649542 0.1445845 0.449 0.65325
## Trte1:PrevAUGSecchi.sc -0.3600133 0.3210735 -1.121 0.26217
## Trte1:Snow.sc 0.0004073 0.2287247 0.002 0.99858
## Trte1:IceCover.sc -0.1484899 0.2041444 -0.727 0.46700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU Tr1:S.
## Trte1 -0.340
## cytrt -0.011 -0.500
## PrvAUGScch. 0.125 -0.115 -0.128
## Snow.sc 0.083 -0.244 0.035 0.014
## IceCover.sc 0.119 -0.035 -0.026 0.031 -0.150
## Trt1:PAUGS. -0.045 0.393 -0.247 -0.443 -0.127 -0.002
## Trt1:Snw.sc -0.039 0.256 -0.155 -0.025 -0.651 0.147 0.212
## Trt1:IcCvr. -0.073 -0.025 0.006 0.022 0.082 -0.697 -0.008 -0.348 #no interactions are significant so include none:
glme.foc.e <- glmer(foc_clp~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.e, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! # model.matrix(glme.foc.e)
summary(glme.foc.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake) + (1 | obs)
## Data: pcri.e
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 914.0 934.9 -449.0 898.0 93
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.36412 -0.12555 -0.00206 0.10055 1.18301
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6224 0.7889
## Lake (Intercept) 0.9775 0.9887
## Number of obs: 101, groups: obs, 101; Lake, 35
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.564845 0.216086 -2.614 0.00895 **
## Trte1 -0.219452 0.236166 -0.929 0.35277
## cytrt -0.529424 0.116909 -4.529 5.94e-06 ***
## PrevAUGSecchi.sc -0.464942 0.204855 -2.270 0.02323 *
## Snow.sc -0.508611 0.113872 -4.467 7.95e-06 ***
## IceCover.sc -0.008906 0.103509 -0.086 0.93143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.349
## cytrt -0.029 -0.441
## PrvAUGScch. 0.124 0.055 -0.267
## Snow.sc 0.061 -0.111 -0.097 0.015
## IceCover.sc 0.087 -0.049 -0.050 0.080 -0.305Early freq affected by cumulative yrs of treatment (-), productivity (+), and snow cover (-), whereas ice had no effect.
Early Rake Density
glme.rd.e <- glmer(rd_pcri~ Trte + cytrt + PrevAUGSecchi.sc*Trte + Snow.sc*Trte + IceCover.sc*Trte +
(1 | Lake), family = Gamma(link = "log"), data = pcri.e)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0306323 (tol = 0.001, component 1) #restart at last fit:
ss <- getME(glme.rd.e,c("theta","fixef"))
glme.rd.e2 <- update(glme.rd.e,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
# worked again
summary(glme.rd.e2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trte + cytrt + PrevAUGSecchi.sc * Trte + Snow.sc *
## Trte + IceCover.sc * Trte + (1 | Lake)
## Data: pcri.e
## Control: glmerControl(optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 56.1 79.7 -17.1 34.1 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7827 -0.6004 0.0420 0.6780 1.8570
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03438 0.1854
## Residual 0.04193 0.2048
## Number of obs: 63, groups: Lake, 19
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.51905 0.08773 5.917 3.29e-09 ***
## Trte1 -0.20024 0.07655 -2.616 0.0089 **
## cytrt -0.03689 0.03075 -1.200 0.2303
## PrevAUGSecchi.sc -0.12413 0.07334 -1.693 0.0905 .
## Snow.sc -0.01873 0.05094 -0.368 0.7131
## IceCover.sc -0.11262 0.04596 -2.450 0.0143 *
## Trte1:PrevAUGSecchi.sc 0.02424 0.08723 0.278 0.7811
## Trte1:Snow.sc 0.00971 0.07128 0.136 0.8916
## Trte1:IceCover.sc 0.11180 0.06042 1.850 0.0643 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU Tr1:S.
## Trte1 -0.303
## cytrt -0.024 -0.454
## PrvAUGScch. 0.231 -0.119 -0.136
## Snow.sc 0.116 -0.253 0.106 -0.022
## IceCover.sc -0.040 0.017 -0.060 -0.074 -0.470
## Trt1:PAUGS. -0.094 0.374 -0.203 -0.379 -0.232 0.098
## Trt1:Snw.sc -0.058 0.345 -0.278 0.013 -0.702 0.339 0.318
## Trt1:IcCvr. 0.019 -0.125 0.082 0.100 0.317 -0.720 -0.136 -0.483 # retain no interactions
glme.rd.e <- glmer(rd_pcri~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc + (1 | Lake), family = Gamma(link = "log"), data = pcri.e)
summary(glme.rd.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.e
##
## AIC BIC logLik deviance df.resid
## 54.8 71.9 -19.4 38.8 55
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.69128 -0.60532 0.00492 0.58267 2.09493
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03594 0.1896
## Residual 0.04510 0.2124
## Number of obs: 63, groups: Lake, 19
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.52467 0.08930 5.875 4.22e-09 ***
## Trte1 -0.21914 0.07106 -3.084 0.00204 **
## cytrt -0.03077 0.03057 -1.007 0.31412
## PrevAUGSecchi.sc -0.13701 0.06903 -1.985 0.04715 *
## Snow.sc -0.01116 0.03788 -0.295 0.76833
## IceCover.sc -0.05131 0.03328 -1.542 0.12309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.297
## cytrt -0.057 -0.379
## PrvAUGScch. 0.227 -0.018 -0.213
## Snow.sc 0.106 -0.007 -0.132 -0.019
## IceCover.sc -0.032 -0.133 -0.010 -0.005 -0.504Only Trt last year (-) and productivity (+) are affecting rake densities
Late Rake Density
# peak rake density
glme.rd.p <- glmer(rd_pcri~ Trt + cytrt + PrevAUGSecchi.sc*Trt + Snow.sc*Trt + IceCover.sc*Trt +
(1 | Lake), family = Gamma(link = "log"), data = pcri.p)
# #restart at last fit:
# ss <- getME(glme.rd.p,c("theta","fixef"))
# glme.rd.p2 <- update(glme.rd.p,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
# # worked again
summary(glme.rd.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trt + cytrt + PrevAUGSecchi.sc * Trt + Snow.sc * Trt +
## IceCover.sc * Trt + (1 | Lake)
## Data: pcri.p
##
## AIC BIC logLik deviance df.resid
## 123.6 149.8 -50.8 101.6 69
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8505 -0.5382 -0.2189 0.6545 2.3310
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.04117 0.2029
## Residual 0.08504 0.2916
## Number of obs: 80, groups: Lake, 28
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.631120 0.083805 7.531 5.04e-14 ***
## Trt1 -0.414851 0.101878 -4.072 4.66e-05 ***
## cytrt -0.040711 0.033490 -1.216 0.224
## PrevAUGSecchi.sc -0.052605 0.047452 -1.109 0.268
## Snow.sc -0.062227 0.053340 -1.167 0.243
## IceCover.sc -0.000906 0.059984 -0.015 0.988
## Trt1:PrevAUGSecchi.sc 0.054079 0.116131 0.466 0.641
## Trt1:Snow.sc -0.055736 0.084470 -0.660 0.509
## Trt1:IceCover.sc 0.030229 0.095003 0.318 0.750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU Tr1:S.
## Trt1 -0.408
## cytrt -0.068 -0.365
## PrvAUGScch. -0.216 0.109 -0.022
## Snow.sc -0.218 0.211 0.018 0.095
## IceCover.sc 0.126 -0.200 -0.022 -0.159 -0.519
## Trt1:PAUGS. 0.135 0.182 -0.191 -0.287 -0.067 0.144
## Trt1:Snw.sc 0.137 0.153 -0.130 -0.107 -0.591 0.266 0.069
## Trt1:IcCvr. -0.028 -0.064 -0.022 0.148 0.292 -0.577 0.096 -0.556 #retain no interactions:
glme.rd.p <- glmer(rd_pcri~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
(1 | Lake), family = Gamma(link = "log"), data = pcri.p)
summary(glme.rd.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.p
##
## AIC BIC logLik deviance df.resid
## 118.3 137.4 -51.2 102.3 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8474 -0.5613 -0.2213 0.8268 2.4093
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03687 0.1920
## Residual 0.08600 0.2933
## Number of obs: 80, groups: Lake, 28
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.634152 0.080237 7.903 2.71e-15 ***
## Trt1 -0.414004 0.099660 -4.154 3.26e-05 ***
## cytrt -0.040986 0.032754 -1.251 0.2108
## PrevAUGSecchi.sc -0.048128 0.043548 -1.105 0.2691
## Snow.sc -0.082609 0.043393 -1.904 0.0569 .
## IceCover.sc 0.001175 0.047325 0.025 0.9802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc
## Trt1 -0.491
## cytrt -0.023 -0.346
## PrvAUGScch. -0.195 0.202 -0.078
## Snow.sc -0.178 0.395 -0.094 0.054
## IceCover.sc 0.127 -0.337 -0.003 -0.026 -0.596Only Trt reduced rd, but snow was marginally significant (-)
Final models:
summary(glme.foc.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## Trt:PrevAUGSecchi.sc + Trt:Snow.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 1194.1 1223.9 -587.1 1174.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38658 -0.19658 0.00133 0.09733 0.89071
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.629 1.276
## Lake (Intercept) 1.728 1.315
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0343 0.2676 -3.865 0.000111 ***
## Trt1 -1.0859 0.3897 -2.786 0.005329 **
## cytrt -0.3002 0.1453 -2.067 0.038746 *
## PrevAUGSecchi.sc -0.5457 0.1862 -2.930 0.003388 **
## Snow.sc -0.1430 0.1644 -0.870 0.384452
## IceCover.sc -0.1270 0.1439 -0.883 0.377422
## Trt1:PrevAUGSecchi.sc 1.1458 0.4801 2.387 0.016993 *
## Trt1:Snow.sc 0.7181 0.2729 2.632 0.008500 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU
## Trt1 -0.336
## cytrt -0.037 -0.434
## PrvAUGScch. -0.141 0.059 -0.024
## Snow.sc -0.089 0.000 0.036 -0.097
## IceCover.sc 0.075 -0.147 -0.023 -0.077 -0.211
## Trt1:PAUGS. 0.060 0.391 -0.153 -0.320 0.044 0.098
## Trt1:Snw.sc 0.045 0.202 -0.217 0.072 -0.583 -0.098 0.101 summary(glme.foc.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake) + (1 | obs)
## Data: pcri.e
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 914.0 934.9 -449.0 898.0 93
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.36412 -0.12555 -0.00206 0.10055 1.18301
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6224 0.7889
## Lake (Intercept) 0.9775 0.9887
## Number of obs: 101, groups: obs, 101; Lake, 35
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.564845 0.216086 -2.614 0.00895 **
## Trte1 -0.219452 0.236166 -0.929 0.35277
## cytrt -0.529424 0.116909 -4.529 5.94e-06 ***
## PrevAUGSecchi.sc -0.464942 0.204855 -2.270 0.02323 *
## Snow.sc -0.508611 0.113872 -4.467 7.95e-06 ***
## IceCover.sc -0.008906 0.103509 -0.086 0.93143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.349
## cytrt -0.029 -0.441
## PrvAUGScch. 0.124 0.055 -0.267
## Snow.sc 0.061 -0.111 -0.097 0.015
## IceCover.sc 0.087 -0.049 -0.050 0.080 -0.305 summary(glme.rd.e)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.e
##
## AIC BIC logLik deviance df.resid
## 54.8 71.9 -19.4 38.8 55
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.69128 -0.60532 0.00492 0.58267 2.09493
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03594 0.1896
## Residual 0.04510 0.2124
## Number of obs: 63, groups: Lake, 19
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.52467 0.08930 5.875 4.22e-09 ***
## Trte1 -0.21914 0.07106 -3.084 0.00204 **
## cytrt -0.03077 0.03057 -1.007 0.31412
## PrevAUGSecchi.sc -0.13701 0.06903 -1.985 0.04715 *
## Snow.sc -0.01116 0.03788 -0.295 0.76833
## IceCover.sc -0.05131 0.03328 -1.542 0.12309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.297
## cytrt -0.057 -0.379
## PrvAUGScch. 0.227 -0.018 -0.213
## Snow.sc 0.106 -0.007 -0.132 -0.019
## IceCover.sc -0.032 -0.133 -0.010 -0.005 -0.504 summary(glme.rd.p)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.p
##
## AIC BIC logLik deviance df.resid
## 118.3 137.4 -51.2 102.3 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8474 -0.5613 -0.2213 0.8268 2.4093
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03687 0.1920
## Residual 0.08600 0.2933
## Number of obs: 80, groups: Lake, 28
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.634152 0.080237 7.903 2.71e-15 ***
## Trt1 -0.414004 0.099660 -4.154 3.26e-05 ***
## cytrt -0.040986 0.032754 -1.251 0.2108
## PrevAUGSecchi.sc -0.048128 0.043548 -1.105 0.2691
## Snow.sc -0.082609 0.043393 -1.904 0.0569 .
## IceCover.sc 0.001175 0.047325 0.025 0.9802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc
## Trt1 -0.491
## cytrt -0.023 -0.346
## PrvAUGScch. -0.195 0.202 -0.078
## Snow.sc -0.178 0.395 -0.094 0.054
## IceCover.sc 0.127 -0.337 -0.003 -0.026 -0.596# unscaled models for use in visualization ---------------------------------------Here we will try to re-construct the same models with unscaled data. This will make sure that our visualizations are on the actual scales of the data.
Models on unscaled predictors
Late & Early Survey Freq
# peak foc
glme.foc.pu <- glmer(foc_clp~ Trt + cytrt + PrevAUGSecchi+ Snow + IceCover + Trt:PrevAUGSecchi + Trt:Snow+
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.p, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0165531 (tol = 0.001, component 1) # restart at last fit:
ss <- getME(glme.foc.pu,c("theta","fixef"))
glme.foc.pu2 <- update(glme.foc.pu,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm! summary(glme.foc.pu2)# worked again## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## foc_clp ~ Trt + cytrt + PrevAUGSecchi + Snow + IceCover + Trt:PrevAUGSecchi +
## Trt:Snow + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
## Control: glmerControl(optCtrl = list(maxfun = 20000))
##
## AIC BIC logLik deviance df.resid
## 1194.1 1223.9 -587.1 1174.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38656 -0.19657 0.00133 0.09733 0.89071
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.629 1.276
## Lake (Intercept) 1.728 1.315
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.039169 1.333049 0.780 0.43566
## Trt1 -3.582182 0.688451 -5.203 1.96e-07 ***
## cytrt -0.300219 0.145266 -2.067 0.03876 *
## PrevAUGSecchi -0.547247 0.186771 -2.930 0.00339 **
## Snow -0.006336 0.007284 -0.870 0.38440
## IceCover -0.008899 0.010082 -0.883 0.37741
## Trt1:PrevAUGSecchi 1.148978 0.481414 2.387 0.01700 *
## Trt1:Snow 0.031807 0.012087 2.631 0.00850 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PrAUGS Snow IceCvr T1:PAU
## Trt1 -0.050
## cytrt 0.013 0.017
## PrevAUGScch -0.103 0.261 -0.024
## Snow 0.023 0.317 0.036 -0.097
## IceCover -0.943 -0.107 -0.023 -0.077 -0.211
## Trt1:PrAUGS -0.041 -0.689 -0.153 -0.320 0.044 0.098
## Trt1:Snow 0.206 -0.579 -0.217 0.072 -0.583 -0.098 0.101 # early foc
glme.foc.eu <- glmer(foc_clp~ Trte + cytrt + PrevAUGSecchi + Snow + IceCover +
(1 | Lake)+ (1 | obs), family = binomial(logit), data = pcri.e, weights = nsamp)## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables? summary(glme.foc.eu)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trte + cytrt + PrevAUGSecchi + Snow + IceCover + (1 |
## Lake) + (1 | obs)
## Data: pcri.e
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 914.0 934.9 -449.0 898.0 93
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.36413 -0.12555 -0.00206 0.10056 1.18302
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6224 0.7889
## Lake (Intercept) 0.9776 0.9887
## Number of obs: 101, groups: obs, 101; Lake, 35
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8644954 0.9806268 0.882 0.3780
## Trte1 -0.2194573 0.2361646 -0.929 0.3528
## cytrt -0.5294258 0.1169077 -4.529 5.94e-06 ***
## PrevAUGSecchi -0.4662558 0.2054352 -2.270 0.0232 *
## Snow -0.0225266 0.0050436 -4.466 7.96e-06 ***
## IceCover -0.0006233 0.0072507 -0.086 0.9315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PrAUGS Snow
## Trte1 -0.022
## cytrt 0.128 -0.441
## PrevAUGScch -0.309 0.055 -0.267
## Snow 0.133 -0.111 -0.097 0.015
## IceCover -0.936 -0.049 -0.050 0.080 -0.305
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?Late & Early Rake Densities
# early rd
glme.rd.eu <- glmer(rd_pcri~ Trte + cytrt + PrevAUGSecchi + Snow + IceCover + (1 | Lake), family = Gamma(inverse), data = pcri.e)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00455251 (tol = 0.001, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables? # the scaling issue was induced when we unscaled the predictor vars. Because they lie on radically dif scales, it doesn't like them. Try restart at last fit:
ss <- getME(glme.rd.eu,c("theta","fixef"))
glme.rd.eu2 <- update(glme.rd.eu,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))# no worky## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00246188 (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables? #peak rd?
glme.rd.pu <- glmer(rd_pcri~ Trt + cytrt + PrevAUGSecchi + Snow + IceCover +
(1 | Lake), family = Gamma(inverse), data = pcri.p)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00248663 (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables? # Again, scaling issue induced when we unscaled the predictors. Try restart at last fit:
ss <- getME(glme.rd.pu,c("theta","fixef"))
glme.rd.pu2 <- update(glme.rd.pu,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))# no worky## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00160221 (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?We have actual scale models for frequency, but not rake density. This is likley because the gamma models do not do well with the variation among the predictor variables.
# visualizations ------------------------------------Visualizing results
# secchi ----------------------------------------------------Secchi
# Early Surveys
#Prev Secchi
length(pcri.e$cytrt)## [1] 101 newdat1<-data.frame(cytrt =rep(0,length(pcri.e$cytrt)),
Trte = factor(x = rep(1, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.e$Snow),length(pcri.e$Snow)),
PrevAUGSecchi = seq(min(pcri.e$PrevAUGSecchi), max(pcri.e$PrevAUGSecchi), length.out = length(pcri.e$cytrt)),
IceCover = rep(mean(pcri.e$IceCover),length(pcri.e$cytrt)),
Lake = pcri.e$Lake,
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
newdat1<-data.frame(cytrt =rep(0,length(pcri.e$cytrt)),
Trte = factor(x = rep(0, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.e$Snow),length(pcri.e$Snow)),
PrevAUGSecchi = seq(min(pcri.e$PrevAUGSecchi), max(pcri.e$PrevAUGSecchi), length.out = length(pcri.e$cytrt)),
IceCover = rep(mean(pcri.e$IceCover),length(pcri.e$cytrt)),
Lake = pcri.e$Lake,
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
predtype <- list("Untreated", "Treated")
pred <- rbindlist(list("Untrt" = UnTrt.pred,"Trt" = Trt.pred), idcol = TRUE)
secchi_early <- ggplot(pred, aes(PrevAUGSecchi,fitted))+
geom_line(size = 2, aes(linetype = .id))+
geom_ribbon(aes(ymin=low, ymax=up, group = .id), color = NA , fill = "black", alpha = .15)+
ylab("Frequency of Occurrence")+
xlab("Previous August Secchi (m)") +
scale_linetype_manual("Previous-year", values = c("Trt" = "solid", "Untrt" = "dashed" ), labels = c("Treated","Untreated"))+
theme(text = element_text(size=20), legend.position = c(.60, .85), legend.key.width = unit(1.25, "in"), legend.background = element_blank())+
ggtitle("Early Survey Period")+
theme(plot.title = element_text(hjust = 0.5))
secchi_early 
# mean change in foc per meter of secchi
lm(fitted~PrevAUGSecchi, data = rbind(UnTrt.pred,Trt.pred))##
## Call:
## lm(formula = fitted ~ PrevAUGSecchi, data = rbind(UnTrt.pred,
## Trt.pred))
##
## Coefficients:
## (Intercept) PrevAUGSecchi
## 0.45923 -0.09249 #' ON average, the early season freq declines 9.2% for each added meter ofSecchi clarity.
# LATE
#note that we cut off our predicitons at 3m Secchi because beyond that we don't
# have enough data to inform about trend/pattern
#Prev Secchi
length(pcri.p$cytrt)## [1] 145 #trtd
newdat1<-data.frame(cytrt =rep(0,length(pcri.p$cytrt)),
Trt = factor(x = rep(1, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.p$Snow),length(pcri.p$Snow)),
PrevAUGSecchi = seq(min(pcri.p$PrevAUGSecchi), 3, length.out = length(pcri.p$cytrt)),
IceCover = rep(mean(pcri.p$IceCover),length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
#untrtd
newdat1<-data.frame(cytrt =rep(0,length(pcri.p$cytrt)),
Trt = factor(x = rep(0, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.p$Snow),length(pcri.p$Snow)),
PrevAUGSecchi = seq(min(pcri.p$PrevAUGSecchi), 3, length.out = length(pcri.p$cytrt)),
IceCover = rep(mean(pcri.p$IceCover),length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
predtype <- list("Untreated", "Treated")
pred <- rbindlist(list("Untrt" = UnTrt.pred,"Trt" = Trt.pred), idcol = TRUE)
secchi_late <- ggplot(pred, aes(PrevAUGSecchi,fitted))+
geom_line(size = 2, aes(linetype = .id))+
geom_ribbon(aes(ymin=low, ymax=up, group = .id), color = NA , fill = "black", alpha = .15)+
ylab("Frequency of Occurrence")+
xlab("Previous August Secchi (m)") +
scale_linetype_manual("Within-year", values = c("Trt" = "solid", "Untrt" = "dashed" ), labels = c("Treated","Untreated"))+
theme(text = element_text(size=20), legend.position = c(.60, .85), legend.key.width = unit(1.25, "in"), legend.background = element_blank())+
ggtitle("Late Survey Period")+
theme(plot.title = element_text(hjust = 0.5))
secchi_late
# secchi_late+
# geom_point(aes(x = pcri.p$PrevAUGSecchi, y = pcri.p$foc_clp, shape = pcri.p$Trt ), size = 5, alpha = .4)+
# geom_line(aes(x = pcri.p$PrevAUGSecchi, y=pcri.p$foc_clp ,group = pcri.p$Lake), alpha = .3)
# mean change in foc per meter of secchi
lm(fitted~PrevAUGSecchi, data = rbind(UnTrt.pred,Trt.pred))##
## Call:
## lm(formula = fitted ~ PrevAUGSecchi, data = rbind(UnTrt.pred,
## Trt.pred))
##
## Coefficients:
## (Intercept) PrevAUGSecchi
## 0.20586 -0.01113Again, a 9.2% decline in frequency for each additonal meter of Secchi
# secchi_late+
# geom_point(aes(x = pcri.p$PrevAUGSecchi, y = pcri.p$foc_clp, shape = pcri.p$Trt ), size = 5, alpha = .4)+
# geom_line(aes(x = pcri.p$PrevAUGSecchi, y=pcri.p$foc_clp ,group = pcri.p$Lake), alpha = .3)
# mean change in foc per meter of secchi
lm(fitted~PrevAUGSecchi, data = UnTrt.pred)##
## Call:
## lm(formula = fitted ~ PrevAUGSecchi, data = UnTrt.pred)
##
## Coefficients:
## (Intercept) PrevAUGSecchi
## 0.38119 -0.09408Again, a 9.4% decline in frequency for each additonal meter of Secchi clarity.
# snow ------------------------------------------------------Snow
# EARLY Surveys
length(pcri.e$cytrt)## [1] 101 newdat1<-data.frame(cytrt =rep(0,length(pcri.e$cytrt)),
Trte = factor(x = rep(1, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = seq(min(pcri.e$Snow), max(pcri.e$Snow), length.out = length(pcri.e$cytrt)),
PrevAUGSecchi = rep(mean(pcri.e$PrevAUGSecchi),length(pcri.e$cytrt)),
IceCover = rep(mean(pcri.e$IceCover),length(pcri.e$cytrt)),
Lake = pcri.e$Lake,
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
newdat1<-data.frame(cytrt =rep(0,length(pcri.e$cytrt)),
Trte = factor(x = rep(0, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = seq(min(pcri.e$Snow), max(pcri.e$Snow), length.out = length(pcri.e$cytrt)),
PrevAUGSecchi = rep(mean(pcri.e$PrevAUGSecchi),length(pcri.e$cytrt)),
IceCover = rep(mean(pcri.e$IceCover),length(pcri.e$cytrt)),
Lake = pcri.e$Lake,
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
head(Trt.pred)## cytrt Trte Snow PrevAUGSecchi IceCover Lake obs fitted up
## 1 0 1 3.81000 1.063366 133.396 belle 1 0.4950654 0.6409608
## 2 0 1 4.76758 1.063366 133.396 blueberry 2 0.4896740 0.6347863
## 3 0 1 5.72516 1.063366 133.396 clearm 3 0.4842849 0.6285921
## 4 0 1 6.68274 1.063366 133.396 coal 4 0.4788995 0.6223812
## 5 0 1 7.64032 1.063366 133.396 coal 5 0.4735190 0.6161564
## 6 0 1 8.59790 1.063366 133.396 coal 6 0.4681447 0.6099208
## low
## 1 0.3500056
## 2 0.3462813
## 3 0.3425517
## 4 0.3388168
## 5 0.3350768
## 6 0.3313318 tail(Trt.pred)## cytrt Trte Snow PrevAUGSecchi IceCover Lake obs fitted up
## 96 0 1 94.78010 1.063366 133.396 susan 96 0.1121478 0.2122954
## 97 0 1 95.73768 1.063366 133.396 vails 97 0.1100179 0.2098915
## 98 0 1 96.69526 1.063366 133.396 weaver 98 0.1079235 0.2075160
## 99 0 1 97.65284 1.063366 133.396 weaver 99 0.1058642 0.2051683
## 100 0 1 98.61042 1.063366 133.396 weaver 100 0.1038396 0.2028480
## 101 0 1 99.56800 1.063366 133.396 weaver 101 0.1018494 0.2005547
## low
## 96 0.05589143
## 97 0.05439579
## 98 0.05293541
## 99 0.05150972
## 100 0.05011814
## 101 0.04876009 snow_early <- ggplot(Trt.pred, aes(Snow,fitted))+
geom_line( size = 2)+
geom_ribbon(aes(ymin=low, ymax=up), alpha=0.2)+
geom_line(data = UnTrt.pred,aes(Snow,fitted), size = 2, linetype = 2)+
geom_ribbon(data = UnTrt.pred,aes(ymin=low, ymax=up), alpha=0.2)+
theme(text = element_text(size=20), legend.position = "none")+
ylab("Frequency of Occurrence")+
xlab("Mean Winter Snow Cover (cm)")
snow_early
# show data and average value line
snow_early+
geom_line(aes(x = pcri.e$Snow, y=pcri.e$foc_clp ,group = pcri.e$Lake), alpha = .3)+
geom_point(aes(x = pcri.e$Snow, y = pcri.e$foc_clp, shape = pcri.e$Trte ), size = 5, alpha = .4)+
geom_smooth( method = "lm", color = "red")
# mean change in foc per centimeter of snow
bound <- bind_rows(UnTrt.pred,Trt.pred)
lm(fitted~Snow, data = bound)##
## Call:
## lm(formula = fitted ~ Snow, data = bound)
##
## Coefficients:
## (Intercept) Snow
## 0.512700 -0.004345 lm(fitted~Snow, data = bound)$'coefficients'[2]*2.54 #convert to inches## Snow
## -0.01103517 #' This comes out to approximately a 1.1% change in foc per addt'l inch of snow
#LATE Surveys
length(pcri.p$cytrt)## [1] 145 newdat1<-data.frame(cytrt =rep(0,length(pcri.p$cytrt)),
Trt = factor(x = rep(1, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = seq(min(pcri.p$Snow), max(pcri.p$Snow), length.out = length(pcri.p$cytrt)),
PrevAUGSecchi = rep(mean(pcri.p$PrevAUGSecchi),length(pcri.p$cytrt)),
IceCover = rep(mean(pcri.p$IceCover),length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
newdat1<-data.frame(cytrt =rep(0,length(pcri.p$cytrt)),
Trt = factor(x = rep(0, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = seq(min(pcri.p$Snow), max(pcri.p$Snow), length.out = length(pcri.p$cytrt)),
PrevAUGSecchi = rep(mean(pcri.p$PrevAUGSecchi),length(pcri.p$cytrt)),
IceCover = rep(mean(pcri.p$IceCover),length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
head(Trt.pred)## cytrt Trt Snow PrevAUGSecchi IceCover Lake obs fitted up
## 1 0 1 2.794000 1.181379 134 reshanau 1 0.04957300 0.1090345
## 2 0 1 3.402542 1.181379 134 reshanau 2 0.05030841 0.1099688
## 3 0 1 4.011083 1.181379 134 reshanau 3 0.05105415 0.1109221
## 4 0 1 4.619625 1.181379 134 reshanau 4 0.05181034 0.1118949
## 5 0 1 5.228167 1.181379 134 reshanau 5 0.05257711 0.1128879
## 6 0 1 5.836708 1.181379 134 reshanau 6 0.05335458 0.1139016
## low
## 1 0.02174709
## 2 0.02220745
## 3 0.02267463
## 4 0.02314860
## 5 0.02362933
## 6 0.02411674 tail(Trt.pred)## cytrt Trt Snow PrevAUGSecchi IceCover Lake obs fitted up
## 140 0 1 87.38129 1.181379 134 clearw 140 0.3102490 0.6460522
## 141 0 1 87.98983 1.181379 134 clearw 141 0.3135756 0.6518601
## 142 0 1 88.59838 1.181379 134 clearw 142 0.3169215 0.6576291
## 143 0 1 89.20692 1.181379 134 clearw 143 0.3202865 0.6633577
## 144 0 1 89.81546 1.181379 134 stjames 144 0.3236703 0.6690442
## 145 0 1 90.42400 1.181379 134 stjames 145 0.3270726 0.6746872
## low
## 140 0.09978262
## 141 0.10027825
## 142 0.10077367
## 143 0.10126893
## 144 0.10176403
## 145 0.10225902 snow_late <- ggplot(Trt.pred, aes(Snow,fitted))+
geom_line(size = 2)+
geom_ribbon(aes(ymin=low, ymax=up), alpha=0.2)+
geom_line(data = UnTrt.pred,aes(Snow,fitted), size = 2, linetype = 2)+
geom_ribbon(data = UnTrt.pred,aes(ymin=low, ymax=up), alpha=0.2)+
theme(text = element_text(size=20), legend.position = "none")+
ylab("Frequency of Occurrence")+
xlab("Mean Winter Snow Cover (cm)")
snow_late
# uncomment these lines to add data in background:
snow_late+
geom_point(aes(x = pcri.p$Snow, y = pcri.p$foc_clp, shape = pcri.p$Trt ), size = 5, alpha = .4)+
geom_line(aes(x = pcri.p$Snow, y=pcri.p$foc_clp ,group = pcri.p$Lake), alpha = .1)
# mean change in foc per inch of snow
bound <- bind_rows(Trt.pred)
lm(fitted~Snow, data = bound)##
## Call:
## lm(formula = fitted ~ Snow, data = bound)
##
## Coefficients:
## (Intercept) Snow
## 0.009798 0.003114 lm(fitted~Snow, data = bound)$'coefficients'[2]*2.54 #convert to inches## Snow
## 0.00790883 # summary(trt.snow.pred)Approximately a 0.7% increase in foc per addt’l inch of snow in the treated lakes
# ice -------------------------------------------------------Ice
# Early Surveys
length(pcri.e$cytrt)## [1] 101 newdat1<-data.frame(cytrt =rep(0,length(pcri.e$cytrt)),
Trte = factor(x = rep(1, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.e$Snow),length(pcri.e$Snow)),
PrevAUGSecchi = rep(mean(pcri.e$PrevAUGSecchi),length(pcri.e$cytrt)),
IceCover = seq(min(pcri.e$IceCover), max(pcri.e$IceCover), length.out = length(pcri.e$cytrt)),
Lake = pcri.e$Lake,
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
newdat1<-data.frame(cytrt =rep(0,length(pcri.e$cytrt)),
Trte = factor(x = rep(0, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.e$Snow),length(pcri.e$Snow)),
PrevAUGSecchi = rep(mean(pcri.e$PrevAUGSecchi),length(pcri.e$cytrt)),
IceCover = seq(min(pcri.e$IceCover), max(pcri.e$IceCover), length.out = length(pcri.e$cytrt)),
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
ice_early <- ggplot(Trt.pred, aes(IceCover,fitted))+
geom_line( size = 2)+
geom_ribbon(aes(ymin=low, ymax=up), alpha=0.2)+
geom_line(data = UnTrt.pred,aes(IceCover, fitted), size = 2, linetype = 2)+
geom_ribbon(data = UnTrt.pred,aes(ymin=low, ymax=up), alpha=0.2)+
theme(text = element_text(size=20), legend.position = "none")+
ylab("Frequency of Occurrence")+
xlab("Winter Ice Cover (days)") +
scale_color_manual(breaks = c(0,1), values = c("blue", "red"))
ice_early
# plots data in background:
ice_early+
geom_point(aes(x = pcri.e$IceCover, y = pcri.e$foc_clp, shape = pcri.e$Trte ), size = 5, alpha = .4)+
geom_line(aes(x = pcri.e$IceCover, y=pcri.e$foc_clp ,group = pcri.e$Lake), alpha = .3)
# Late Surveys
length(pcri.p$cytrt)## [1] 145 newdat1<-data.frame(cytrt =rep(0,length(pcri.p$cytrt)),
Trt = factor(x = rep(1, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.p$Snow),length(pcri.p$Snow)),
PrevAUGSecchi = rep(mean(pcri.p$PrevAUGSecchi),length(pcri.p$cytrt)),
IceCover = seq(min(pcri.p$IceCover), max(pcri.p$IceCover), length.out = length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
newdat1<-data.frame(cytrt =rep(0,length(pcri.p$cytrt)),
Trt = factor(x = rep(0, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.p$Snow),length(pcri.p$Snow)),
PrevAUGSecchi = rep(mean(pcri.p$PrevAUGSecchi),length(pcri.p$cytrt)),
IceCover = seq(min(pcri.p$IceCover), max(pcri.p$IceCover), length.out = length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
ice_late <- ggplot(Trt.pred, aes(IceCover,fitted))+
geom_line(size = 2)+
geom_ribbon(aes(ymin=low, ymax=up), alpha=0.2)+
geom_line(data = UnTrt.pred,aes(IceCover, fitted),size = 2, linetype = 2 )+
geom_ribbon(data = UnTrt.pred,aes(ymin=low, ymax=up), alpha=0.2)+
theme(text = element_text(size=20), legend.position = "none")+
ylab("Frequency of Occurrence")+
xlab("Winter Ice Cover (days)")
ice_late
ice_late+
geom_point(aes(x = pcri.p$IceCover, y = pcri.p$foc_clp, shape = pcri.p$Trt ), size = 5, alpha = .4)+
geom_line(aes(x = pcri.p$IceCover, y=pcri.p$foc_clp ,group = pcri.p$Lake), alpha = .3)+
scale_color_manual(breaks = c(0,1), values = c("blue", "red"))
these are much nicer, but they do not explain how the effect that is observed is after accounting for other drivers (background data are “raw”)
# cytrt -------------------------------------------------------------------Consecutive Years of Treatment
# EARLY Surveys
length(pcri.e$cytrt)## [1] 101 newdat1<-data.frame(cytrt = seq(min(pcri.e$cytrt),
max(pcri.e$cytrt), length.out = length(pcri.e$cytrt)),
Trte = factor(x = rep(1, length(pcri.e$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.e$Snow),length(pcri.e$cytrt)),
PrevAUGSecchi = rep(mean(pcri.e$PrevAUGSecchi),length(pcri.e$cytrt)),
IceCover = rep(mean(pcri.e$IceCover),length(pcri.e$cytrt)),
Lake = pcri.e$Lake,
obs = pcri.e$obs)
#fitted values
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1) # see fn documentation in load fns section
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1The line we’ll show is smoothed, but here are the declines for each added, consecutive year of treatment:
Trt.pred[Trt.pred$cytrt == 0,"fitted"] - Trt.pred[Trt.pred$cytrt == 1,"fitted"] # This means going from a single year of treating to two years (1 addt'l consecutive year)## [1] 0.1059016 Trt.pred[Trt.pred$cytrt == 1,"fitted"] - Trt.pred[Trt.pred$cytrt == 2,"fitted"]## [1] 0.07974466 Trt.pred[Trt.pred$cytrt == 2,"fitted"] - Trt.pred[Trt.pred$cytrt == 3,"fitted"]## [1] 0.05516925 Trt.pred[Trt.pred$cytrt == 3,"fitted"] - Trt.pred[Trt.pred$cytrt == 4,"fitted"]## [1] 0.03596682 Trt.pred[Trt.pred$cytrt == 4,"fitted"] - Trt.pred[Trt.pred$cytrt == 5,"fitted"]## [1] 0.02254865If we strike a line through the model we find mean change in freq per additonal year of treatment.
bound <- bind_rows(Trt.pred)
lm(fitted~cytrt, data = bound)##
## Call:
## lm(formula = fitted ~ cytrt, data = bound)
##
## Coefficients:
## (Intercept) cytrt
## 0.28591 -0.05813 # note that I have shifted the axis labels so that the first year of trt is not zero
cytrt_early <- ggplot(Trt.pred, aes(cytrt,fitted))+
geom_line(size = 1.5)+
geom_ribbon(aes(ymin=low, ymax=up), alpha=0.4)+
theme(text = element_text(size=20), legend.position = "none")+
ylab("Frequency of Occurrence")+
xlab("Consecutive Years Treated") +
scale_x_continuous( breaks = seq(0,5, by = 1),
labels = seq(1,6, by = 1),
limits = c(0,5))
cytrt_early 
# plot data in backgroud
cytrt_early +
geom_point(x = pcri.e$cytrt, y = pcri.e$foc_clp)+
geom_smooth(method = "lm")
# LATE Surveys
length(pcri.p$cytrt)## [1] 145 newdat1<-data.frame(cytrt = seq(min(pcri.p$cytrt),
max(pcri.p$cytrt), length.out = length(pcri.p$cytrt)),
Trt = factor(x = rep(1, length(pcri.p$cytrt)), levels = c(0,1)),
Snow = rep(mean(pcri.p$Snow),length(pcri.p$cytrt)),
PrevAUGSecchi = rep(mean(pcri.p$PrevAUGSecchi),length(pcri.p$cytrt)),
IceCover = rep(mean(pcri.p$IceCover),length(pcri.p$cytrt)),
Lake = pcri.p$Lake,
obs = pcri.p$obs)
#fitted values
phat1<-predict(glme.foc.pu2, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu2,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
Trt.pred <- newdat1
# declines for each added, consecutive year of treatment:
Trt.pred[Trt.pred$cytrt == 0,"fitted"] - Trt.pred[Trt.pred$cytrt == 1,"fitted"] # This means going from a single year of treating to two years (1 addt'l consecutive year)## [1] 0.0258115 Trt.pred[Trt.pred$cytrt == 1,"fitted"] - Trt.pred[Trt.pred$cytrt == 2,"fitted"]## [1] 0.020102 Trt.pred[Trt.pred$cytrt == 2,"fitted"] - Trt.pred[Trt.pred$cytrt == 3,"fitted"]## [1] 0.01546596 Trt.pred[Trt.pred$cytrt == 3,"fitted"] - Trt.pred[Trt.pred$cytrt == 4,"fitted"]## [1] 0.01178798 Trt.pred[Trt.pred$cytrt == 4,"fitted"] - Trt.pred[Trt.pred$cytrt == 5,"fitted"]## [1] 0.008920582 Trt.pred[Trt.pred$cytrt == 5,"fitted"] - Trt.pred[Trt.pred$cytrt == 6,"fitted"]## [1] 0.006714163 Trt.pred[Trt.pred$cytrt == 6,"fitted"] - Trt.pred[Trt.pred$cytrt == 7,"fitted"]## [1] 0.005032884 Trt.pred[Trt.pred$cytrt == 7,"fitted"] - Trt.pred[Trt.pred$cytrt == 8,"fitted"]## [1] 0.003761074 # avg change per year modeled
bound <- bind_rows(Trt.pred)
lm(fitted~cytrt, data = bound)##
## Call:
## lm(formula = fitted ~ cytrt, data = bound)
##
## Coefficients:
## (Intercept) cytrt
## 0.08908 -0.01143 cytrt_late <- ggplot(Trt.pred, aes(cytrt,fitted))+
geom_line( size = 1.5)+
geom_ribbon(aes(ymin=low, ymax=up), alpha=0.4)+
theme(text = element_text(size=20), legend.position = "none")+
ylab("Frequency of Occurrence")+
xlab("Consecutive Years Treated") +
scale_x_continuous( breaks = seq(0,8, by = 1),
labels = seq(1,9, by = 1),
limits = c(0,9))
cytrt_late
# mean change in foc per yr trt
cytrt_late+
geom_point(x = pcri.p$cytrt, y = pcri.p$foc_clp, alpha = .5)+
geom_smooth(method = "lm")
# treatment effects -------------------------------------------------------Treatments
#Trte switchover for EARLY freq:
length(pcri.e$cytrt)## [1] 101 newdat1<-data.frame(Trte = factor(x = c(rep(0,50),rep(1,51)), levels = c(0,1)),
IceCover = rep(mean(pcri.e$IceCover),length(pcri.e$IceCover)),
Snow = rep(mean(pcri.e$Snow),length(pcri.e$Snow)),
cytrt = rep(0,length(pcri.e$cytrt)),
PrevAUGSecchi = rep(mean(pcri.e$PrevAUGSecchi),length(pcri.e$PrevAUGSecchi)))
phat1<-predict(glme.foc.eu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.eu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
modvals = data.frame(Trte= c(0,1), modfoc=c(max(UnTrt.pred$fitted),min(UnTrt.pred$fitted)), lwr= c(.75, 1.75), upr = c(1.25, 2.25) )
# ggplot(pcri.e, aes(x = pcri.e$foc_clp))+
# geom_density(aes(group = Trte, lty = Trte))+
# geom_segment(data = modvals,aes(y = 0, yend = 2, x = modfoc, xend = modfoc, lty = as.factor(Trte)))+
# theme(text = element_text(size=20), legend.position = "none")+
# ylab("Density")+xlab("Frequency of Occurrence")+
# scale_fill_manual( values = c("blue", "red"))+
# theme_bw()
#
trt_foc_early <- ggplot(pcri.e, aes(x=foc_clp)) +
geom_density( aes(linetype = Trte), alpha=0.25, position="identity", lwd = 1.5)+
scale_linetype_manual("Previous-year", values = c("0" = "dotted", "1" = "solid"), labels = c("Untreated", "Treated"))+
geom_vline(xintercept = modvals$modfoc, size = 1.5, lty = c("0" = "dotted", "1" = "solid"))+
theme(legend.position = c(.75, .75), legend.background = element_blank())+
ylab("Kernel Density")+xlab("Frequency of Occurrence")+
geom_text( aes(.42,.75, label = "a"), size = 7)+
geom_text( aes(.30,.75, label = "a"), size = 7)+
theme(text = element_text(size=20))+
ggtitle("Early Survey Period")+
theme(plot.title = element_text(hjust = 0.5))
trt_foc_early
#Trt switchover for Late freq:
length(pcri.p$cytrt)## [1] 145 newdat1<-data.frame(Trt = factor(x = c(rep(0,72),rep(1,73)), levels = c(0,1)),
IceCover = rep(mean(pcri.p$IceCover),length(pcri.p$IceCover)),
Snow = rep(mean(pcri.p$Snow),length(pcri.p$Snow)),
cytrt = rep(0,length(pcri.p$cytrt)),
PrevAUGSecchi = rep(mean(pcri.p$PrevAUGSecchi),length(pcri.p$PrevAUGSecchi)))
phat1<-predict(glme.foc.pu, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-plogis(phat1)
cpred1.CI <- easyPredCI(glme.foc.pu,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
modvals = data.frame(Trt= c(0,1), modfoc=c(max(UnTrt.pred$fitted),min(UnTrt.pred$fitted)), lwr= c(.75, 1.75), upr = c(1.25, 2.25) )
# ggplot(pcri.p, aes(x = pcri.p$foc_clp))+
# geom_density(aes(group = Trt, lty = Trt))+
# theme(text = element_text(size=20), legend.position = "none")+
# ylab("Density")+xlab("Frequency of Occurrence")+
# geom_segment(data = modvals,aes(y = 0, yend = 3, x = modfoc, xend = modfoc, lty = as.factor(Trt)))+
# theme_bw()
trt_foc_late <- ggplot(pcri.p, aes(x=foc_clp)) +
geom_density(aes(linetype = Trt), alpha=0.25, position="identity", size = 1.5)+
theme(text = element_text(size=20))+
scale_linetype_manual("Within-year", values = c("0" = "dotted", "1" = "solid"), labels = c("Untreated","Treated"))+
geom_vline(xintercept = modvals$modfoc, size = 1.5, lty = c("0" = "dotted", "1" = "solid"))+
ylab("Kernel Density")+xlab("Frequency of Occurrence")+
geom_text( aes(.07,.75, label = "a"), size = 7)+
geom_text( aes(.30,.75, label = "b"), size = 7)+
theme(legend.position = c(.75, .75), legend.background = element_blank())+
ggtitle("Late Survey Period")+
theme(plot.title = element_text(hjust = 0.5))
trt_foc_late
# Trt switchover for early rd:
length(pcri.e$cytrt)## [1] 101 newdat1<-data.frame(Trte = factor(x = c(rep(0,50),rep(1,51)), levels = c(0,1)),
IceCover.sc = rep(mean(pcri.e$IceCover.sc),length(pcri.e$IceCover.sc)),
Snow.sc = rep(mean(pcri.e$Snow.sc),length(pcri.e$Snow.sc)),
cytrt = rep(0,length(pcri.e$cytrt)),
PrevAUGSecchi.sc = rep(mean(pcri.e$PrevAUGSecchi.sc),length(pcri.e$PrevAUGSecchi.sc)))
phat1<-predict(glme.rd.e, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-exp(phat1)
cpred1.CI <- easyPredCI(glme.rd.e,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
modvals = data.frame(Trte= c(0,1), modfoc=c(max(UnTrt.pred$fitted),min(UnTrt.pred$fitted)), lwr= c(.75, 1.75), upr = c(1.25, 2.25) )
trt_rd_early <- ggplot(pcri.e, aes(x=rd_pcri)) +
geom_density(aes(linetype=Trte), alpha=0.25, position="identity", size = 1.5)+
xlim(c(0,4))+
theme(text = element_text(size=20))+
ylab("Kernel Density")+xlab("Mean Rake Density")+
scale_linetype_manual("Previous-year", values = c("0" = "dotted", "1" = "solid"), labels = c("Untreated", "Treated"))+
geom_vline(xintercept = modvals$modfoc, size = 1.5, lty = c("0" = "dotted", "1" = "solid"))+
geom_text( aes(1.25,.75, label = "a"), size = 7)+
geom_text( aes(1.86,.75, label = "b"), size = 7)+
theme(legend.position = c(.75, .75), legend.background = element_blank())
trt_rd_early## Warning: Removed 38 rows containing non-finite values (stat_density).
#trt switchover for peak rd
length(pcri.p$cytrt)## [1] 145 newdat1<-data.frame(Trt = factor(x = c(rep(0,72),rep(1,73)), levels = c(0,1)),
IceCover.sc = rep(mean(pcri.p$IceCover.sc),length(pcri.p$IceCover.sc)),
Snow.sc = rep(mean(pcri.p$Snow.sc),length(pcri.p$Snow.sc)),
cytrt = rep(0,length(pcri.p$cytrt)),
PrevAUGSecchi.sc = rep(mean(pcri.p$PrevAUGSecchi.sc),length(pcri.p$PrevAUGSecchi.sc)))
phat1<-predict(glme.rd.p, newdat=newdat1, re.form=NA, type="link")
newdat1$fitted<-exp(phat1)
cpred1.CI <- easyPredCI(glme.rd.p,newdata = newdat1)
newdat1$up<-cpred1.CI[,2]
newdat1$low<-cpred1.CI[,1]
UnTrt.pred <- newdat1
modvals = data.frame(Trt= c(0,1), modfoc=c(max(UnTrt.pred$fitted),min(UnTrt.pred$fitted)), lwr= c(.75, 1.75), upr = c(1.25, 2.25) )
trt_rd_late <- ggplot(pcri.p, aes(x=rd_pcri)) +
geom_density(aes( linetype = Trt), alpha=0.25, position="identity", size = 1.5)+
xlim(c(0,4))+
theme(text = element_text(size=20))+
ylab("Kernel Density")+xlab("Mean Rake Density")+
scale_linetype_manual("Within-year", values = c("0" = "dotted", "1" = "solid"), labels = c("Untreated", "Treated"))+
geom_vline(xintercept = modvals$modfoc, size = 1.5, lty = c("0" = "dotted", "1" = "solid"))+
geom_text( aes(1.1,.75, label = "a"), size = 7)+
geom_text( aes(2.05,.75, label = "b"), size = 7)+
theme(legend.position = c(.75, .75), legend.background = element_blank())
trt_rd_late## Warning: Removed 65 rows containing non-finite values (stat_density).
# area treated declines -----------------------------------------------Area Managed
Evidence for decreased intervention thru time?
# does acreage decline with increased cum yrs trt (results in hectares)
ggplot(subset(pcri.p, Trt == 1), aes(jitter(cytrt), Acres_Chem*0.404686))+
geom_point(size = 4, alpha = .2)+
geom_smooth(method = "lm", color = "black")+
theme_bw()
area_trt <- summary(lmer(Acres_Chem ~ cytrt + ( 1 | Lake ), data = subset(pcri.p, Trt == 1)))
summary(area_trt)## Linear mixed model fit by REML ['lmerMod']
## Formula: Acres_Chem ~ cytrt + (1 | Lake)
## Data: subset(pcri.p, Trt == 1)
##
## REML criterion at convergence: 711
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6916 -0.4242 -0.1450 0.4189 3.0122
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 5526 74.34
## Residual 3395 58.26
## Number of obs: 63, groups: Lake, 21
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 111.417 19.113 5.829
## cytrt -13.684 5.512 -2.482
##
## Correlation of Fixed Effects:
## (Intr)
## cytrt -0.301#convert to ha
area_trt$coefficients*0.404686## Estimate Std. Error t value
## (Intercept) 45.088915 7.734835 2.359049
## cytrt -5.537669 2.230710 -1.004620so the acreage treated declines by 5.5 hectares +- 1.96*2.23 for each additonal year of treatment and avearge lakes treated 45.1 ha in first treatment
# treating based on spring populations? -----------------------------------Bias in Treated Late Surveys
Do you get more treatment in years with high early season foc?
ggplot(subset(pcri.e, Trte == 1), aes(foc_clp,perc.lit.chem))+
geom_point()+
geom_smooth(method = "lm")
ggplot(subset(pcri.e, Trte == 1), aes(Snow,perc.lit.chem))+
geom_point()+
geom_smooth(method = "lm")
summary(lm(data = subset(pcri.e, Trte == 1), perc.lit.chem ~ foc_clp))##
## Call:
## lm(formula = perc.lit.chem ~ foc_clp, data = subset(pcri.e, Trte ==
## 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6535 -0.2833 -0.1632 0.1403 1.8107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2613 0.1237 2.112 0.0401 *
## foc_clp 0.8497 0.3447 2.465 0.0175 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5084 on 46 degrees of freedom
## Multiple R-squared: 0.1167, Adjusted R-squared: 0.09747
## F-statistic: 6.076 on 1 and 46 DF, p-value: 0.0175summary(lm(data = subset(pcri.e, Trte == 1), perc.lit.chem ~ Snow))##
## Call:
## lm(formula = perc.lit.chem ~ Snow, data = subset(pcri.e, Trte ==
## 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6306 -0.3272 -0.1766 0.3011 1.9110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.667385 0.153466 4.349 7.52e-05 ***
## Snow -0.004528 0.003746 -1.209 0.233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5326 on 46 degrees of freedom
## Multiple R-squared: 0.03079, Adjusted R-squared: 0.009718
## F-statistic: 1.461 on 1 and 46 DF, p-value: 0.2329plot(pcri.p$perc.lit.chem ~ pcri.p$foc_clp)
plot(pcri.p$foc_clp~ pcri.p$perc.lit.chem)
ggplot(pcri.p, aes(cytrt, Acres_Chem))+
geom_point()+
geom_path(aes(group = Lake))
# export centered scaled models for reporting -----------------------------Review and export modeled estimates
Reported parameter estimates should be on a centered scaled model, thus allowing intuitive inference of sign.
summary(glme.foc.p)# June_FOC_split_secchi(notsigdifzero)&snow(nonzero) (centered scaled)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## Trt:PrevAUGSecchi.sc + Trt:Snow.sc + (1 | Lake) + (1 | obs)
## Data: pcri.p
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 1194.1 1223.9 -587.1 1174.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38658 -0.19658 0.00133 0.09733 0.89071
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 1.629 1.276
## Lake (Intercept) 1.728 1.315
## Number of obs: 145, groups: obs, 145; Lake, 45
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0343 0.2676 -3.865 0.000111 ***
## Trt1 -1.0859 0.3897 -2.786 0.005329 **
## cytrt -0.3002 0.1453 -2.067 0.038746 *
## PrevAUGSecchi.sc -0.5457 0.1862 -2.930 0.003388 **
## Snow.sc -0.1430 0.1644 -0.870 0.384452
## IceCover.sc -0.1270 0.1439 -0.883 0.377422
## Trt1:PrevAUGSecchi.sc 1.1458 0.4801 2.387 0.016993 *
## Trt1:Snow.sc 0.7181 0.2729 2.632 0.008500 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc IcCvr. T1:PAU
## Trt1 -0.336
## cytrt -0.037 -0.434
## PrvAUGScch. -0.141 0.059 -0.024
## Snow.sc -0.089 0.000 0.036 -0.097
## IceCover.sc 0.075 -0.147 -0.023 -0.077 -0.211
## Trt1:PAUGS. 0.060 0.391 -0.153 -0.320 0.044 0.098
## Trt1:Snw.sc 0.045 0.202 -0.217 0.072 -0.583 -0.098 0.101summary(glme.foc.e)# April_FOC_split_none (centered scaled)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: foc_clp ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake) + (1 | obs)
## Data: pcri.e
## Weights: nsamp
##
## AIC BIC logLik deviance df.resid
## 914.0 934.9 -449.0 898.0 93
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.36412 -0.12555 -0.00206 0.10055 1.18301
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6224 0.7889
## Lake (Intercept) 0.9775 0.9887
## Number of obs: 101, groups: obs, 101; Lake, 35
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.564845 0.216086 -2.614 0.00895 **
## Trte1 -0.219452 0.236166 -0.929 0.35277
## cytrt -0.529424 0.116909 -4.529 5.94e-06 ***
## PrevAUGSecchi.sc -0.464942 0.204855 -2.270 0.02323 *
## Snow.sc -0.508611 0.113872 -4.467 7.95e-06 ***
## IceCover.sc -0.008906 0.103509 -0.086 0.93143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.349
## cytrt -0.029 -0.441
## PrvAUGScch. 0.124 0.055 -0.267
## Snow.sc 0.061 -0.111 -0.097 0.015
## IceCover.sc 0.087 -0.049 -0.050 0.080 -0.305summary(glme.rd.e)# April_DENS_split_none (centered scaled)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trte + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.e
##
## AIC BIC logLik deviance df.resid
## 54.8 71.9 -19.4 38.8 55
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.69128 -0.60532 0.00492 0.58267 2.09493
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03594 0.1896
## Residual 0.04510 0.2124
## Number of obs: 63, groups: Lake, 19
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.52467 0.08930 5.875 4.22e-09 ***
## Trte1 -0.21914 0.07106 -3.084 0.00204 **
## cytrt -0.03077 0.03057 -1.007 0.31412
## PrevAUGSecchi.sc -0.13701 0.06903 -1.985 0.04715 *
## Snow.sc -0.01116 0.03788 -0.295 0.76833
## IceCover.sc -0.05131 0.03328 -1.542 0.12309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trte1 cytrt PAUGS. Snw.sc
## Trte1 -0.297
## cytrt -0.057 -0.379
## PrvAUGScch. 0.227 -0.018 -0.213
## Snow.sc 0.106 -0.007 -0.132 -0.019
## IceCover.sc -0.032 -0.133 -0.010 -0.005 -0.504summary(glme.rd.p)# June_DENS_split_none (centered scaled)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: rd_pcri ~ Trt + cytrt + PrevAUGSecchi.sc + Snow.sc + IceCover.sc +
## (1 | Lake)
## Data: pcri.p
##
## AIC BIC logLik deviance df.resid
## 118.3 137.4 -51.2 102.3 72
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8474 -0.5613 -0.2213 0.8268 2.4093
##
## Random effects:
## Groups Name Variance Std.Dev.
## Lake (Intercept) 0.03687 0.1920
## Residual 0.08600 0.2933
## Number of obs: 80, groups: Lake, 28
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.634152 0.080237 7.903 2.71e-15 ***
## Trt1 -0.414004 0.099660 -4.154 3.26e-05 ***
## cytrt -0.040986 0.032754 -1.251 0.2108
## PrevAUGSecchi.sc -0.048128 0.043548 -1.105 0.2691
## Snow.sc -0.082609 0.043393 -1.904 0.0569 .
## IceCover.sc 0.001175 0.047325 0.025 0.9802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trt1 cytrt PAUGS. Snw.sc
## Trt1 -0.491
## cytrt -0.023 -0.346
## PrvAUGScch. -0.195 0.202 -0.078
## Snow.sc -0.178 0.395 -0.094 0.054
## IceCover.sc 0.127 -0.337 -0.003 -0.026 -0.596summary(glme.foc.p)[["coefficients"]]## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0342917 0.2676094 -3.8649297 0.0001111212
## Trt1 -1.0859093 0.3897078 -2.7864701 0.0053285527
## cytrt -0.3002459 0.1452655 -2.0668768 0.0387457600
## PrevAUGSecchi.sc -0.5457330 0.1862460 -2.9301730 0.0033877333
## Snow.sc -0.1430248 0.1644487 -0.8697225 0.3844520754
## IceCover.sc -0.1270340 0.1439223 -0.8826567 0.3774217828
## Trt1:PrevAUGSecchi.sc 1.1458199 0.4800532 2.3868602 0.0169929537
## Trt1:Snow.sc 0.7181161 0.2728871 2.6315502 0.0084996312Build and export Table 3
early_mod_foc <- summary(glme.foc.e)[["coefficients"]]
peak_mod_foc <- summary(glme.foc.p)[["coefficients"]]
early_mod_rd <- summary(glme.rd.e)[["coefficients"]]
peak_mod_rd <- summary(glme.rd.p)[["coefficients"]]
final_model_coefs <- cbind(c(rep("early_foc",6),rep("peak_foc",8),rep("early_rd",6),rep("peak_rd",6)), rbind(early_mod_foc,peak_mod_foc,early_mod_rd,peak_mod_rd))
# write.csv(final_model_coefs, file = "data/output_data/final_model_coefs.csv")
# mapping data ------------------------------------------------------------Map of Study
states <- map_data("state")
mn_df <- subset(states, region == "minnesota")
mn_lakes <- st_read("G:/My Drive/Documents/UMN/Grad School/Larkin Lab/R_projects/curlyleaf_mgmt/data/spatial_data/MN_lakesonly_trimmed_strytrk.shp")## Reading layer `MN_lakesonly_trimmed_strytrk' from data source `G:\My Drive\Documents\UMN\Grad School\Larkin Lab\R_projects\curlyleaf_mgmt\data\spatial_data\MN_lakesonly_trimmed_strytrk.shp' using driver `ESRI Shapefile'
## Simple feature collection with 11715 features and 31 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 191375.6 ymin: 4817046 xmax: 913706 ymax: 5463414
## epsg (SRID): 26915
## proj4string: +proj=utm +zone=15 +datum=NAD83 +units=m +no_defsstudy_lakes <- st_read("G:/My Drive/Documents/UMN/Grad School/Larkin Lab/R_projects/curlyleaf_mgmt/data/spatial_data/lake_loc_fromRN.shp")## Reading layer `lake_loc_fromRN' from data source `G:\My Drive\Documents\UMN\Grad School\Larkin Lab\R_projects\curlyleaf_mgmt\data\spatial_data\lake_loc_fromRN.shp' using driver `ESRI Shapefile'
## Simple feature collection with 60 features and 10 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -96.251 ymin: 43.903 xmax: -92.826 ymax: 47.002
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defsclp_lakes <- st_read("G:/My Drive/Documents/UMN/Grad School/Larkin Lab/R_projects/curlyleaf_mgmt/data/spatial_data/clppoints.shp")## Reading layer `clppoints' from data source `G:\My Drive\Documents\UMN\Grad School\Larkin Lab\R_projects\curlyleaf_mgmt\data\spatial_data\clppoints.shp' using driver `ESRI Shapefile'
## Simple feature collection with 759 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 217230.6 ymin: 4819588 xmax: 747105.7 ymax: 5427245
## epsg (SRID): NA
## proj4string: NAclp_lakes <- clp_lakes %>% st_set_crs(2027)
study_map <- ggplot(study_lakes)+
geom_sf(data = clp_lakes, size = 3 ,shape = 1, alpha = 0.5)+
geom_sf(size = 8, shape = "+", color= "black")+
geom_polygon(data = mn_df,aes(x = long, y = lat), color = "black", alpha = .05)+
scale_shape_discrete(solid = FALSE)+
theme(text = element_text(size=20), legend.position = )+
ylab("Longitude")+
xlab("Latitude")
study_map
#add scalebar & North arrow
study_map+
north(study_lakes, location = "topright", anchor = c("x"=-90, "y"= 49), scale = 0.2)
#(https://stackoverflow.com/questions/39067838/parsimonious-way-to-add-north-arrow-and-scale-bar-to-ggmap):
scalebar = function(x,y,w,n,d, units="Degrees"){
# x,y = lower left coordinate of bar
# w = width of bar
# n = number of divisions on bar
# d = distance along each division
bar = data.frame(
xmin = seq(0.0, n*d, by=d) + x,
xmax = seq(0.0, n*d, by=d) + x + d,
ymin = y,
ymax = y+w,
z = rep(c(1,0),n)[1:(n+1)],
fill.col = rep(c("black","white"),n)[1:(n+1)])
labs = data.frame(
xlab = c(seq(0.0, (n+1)*d, by=d) + x, x),
ylab = c(rep(y-w*1.5, n+2), y-3*w),
text = c(as.character(seq(0.0, ((n+1)*d), by=d)), units)
)
list(bar, labs)
}
sb = scalebar( - 92.5, 49.1, .1, 2, .5 )
# Plot map
study_map+
north(study_lakes, location = "topright", anchor = c("x"=-90, "y"= 49), scale = 0.2)+
geom_rect(data=sb[[1]], aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, fill=z), inherit.aes=F,
show.legend = F, color = "black", fill = sb[[1]]$fill.col) +
geom_text(data=sb[[2]], aes(x=xlab, y=ylab, label=text), inherit.aes=F, show.legend = F) 
# add inset:
study_map <- ggplot(study_lakes)+
geom_sf(data = clp_lakes, size = 3 ,shape = 1, alpha = 0.5)+
geom_sf(size = 8, shape = "+", color= "black")+
geom_polygon(data = mn_df,aes(x = long, y = lat), color = "black", alpha = .05)+
scale_shape_discrete(solid = FALSE)+
theme(text = element_text(size=20), legend.position = )+
ylab("Longitude")+
xlab("Latitude")
study_map
#inset
foo <- map_data("state")
g2 <- ggplotGrob(
ggplot() +
geom_polygon(data = foo,
aes(x = long, y = lat, group = group),
fill = NA, color = "black", alpha = 0.5) +
geom_polygon(data = mn_df,aes(x = long, y = lat), color = "black", alpha = .95)+
theme(panel.background = element_rect(fill = NULL))+
coord_map("polyconic")+
theme_bw()+
theme_inset()
)
g3 <- study_map +
annotation_custom(grob = g2, xmin = -92.5, xmax = -89.1,
ymin = 44.4, ymax = 46.1)
sb = scalebar( - 91.7, 46.5, .1, 2, .5 )
study_map <- ggplot(study_lakes)+
geom_sf(data = clp_lakes, size = 3 ,shape = 1, alpha = 0.5)+
geom_sf(size = 8, shape = "+", color= "black")+
geom_polygon(data = mn_df,aes(x = long, y = lat), color = "black", alpha = .0)+
scale_shape_discrete(solid = FALSE)+
annotation_custom(grob = g2, xmin = -92.5, xmax = -89.1,
ymin = 44.4, ymax = 46.1)+
north(study_lakes, location = "topright", anchor = c("x"=-89.3, "y"= 46.7), scale = 0.2)+
geom_rect(data=sb[[1]], aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, fill=z), inherit.aes=F,
show.legend = F, color = "black", fill = sb[[1]]$fill.col) +
geom_text(data=sb[[2]], aes(x=xlab, y=ylab, label=text), inherit.aes=F, show.legend = F)+
theme_bw()+
theme(text = element_text(size=20), legend.position = )+
ylab("Longitude")+
xlab("Latitude")
# figure construction -----------------------------------------------------
#Figure 1
tiff("Fig1.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw")
plot(study_map) # Make plot
dev.off()## png
## 2# Figure 2
fig2 <- grid.arrange(trt_foc_early,trt_foc_late,trt_rd_early,trt_rd_late, cytrt_early, cytrt_late )## Warning: Removed 38 rows containing non-finite values (stat_density).## Warning: Removed 65 rows containing non-finite values (stat_density).
tiff("Fig2.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw")
plot(fig2) # Make plot
dev.off()## png
## 2#Figure 3
fig3 <- grid.arrange( secchi_early, secchi_late, snow_early, snow_late, ice_early, ice_late )
tiff("Fig3.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw")
plot(fig3) # Make plot
dev.off()## png
## 2# footer ------------------------------------------------------------------Session Information
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggmap_3.0.0 GISTools_0.7-4 rgeos_0.5-2
## [4] MASS_7.3-51.4 RColorBrewer_1.1-2 maptools_0.9-9
## [7] mapproj_1.2.6 ggthemes_4.2.0 raster_3.0-7
## [10] sp_1.3-2 ggsn_0.5.3 devtools_2.2.1
## [13] usethis_1.5.1 rlang_0.4.2 gridExtra_2.3
## [16] mapdata_2.3.0 maps_3.3.0 sf_0.8-0
## [19] numDeriv_2016.8-1.1 data.table_1.12.8 dplyr_0.8.3
## [22] lme4_1.1-21 Matrix_1.2-17 ggplot2_3.2.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-140 bitops_1.0-6 fs_1.3.1 httr_1.4.1
## [5] rprojroot_1.3-2 tools_3.6.1 backports_1.1.5 R6_2.4.1
## [9] KernSmooth_2.23-15 DBI_1.1.0 lazyeval_0.2.2 colorspace_1.4-1
## [13] withr_2.1.2 tidyselect_0.2.5 prettyunits_1.0.2 processx_3.4.1
## [17] curl_4.3 compiler_3.6.1 cli_2.0.0 desc_1.2.0
## [21] labeling_0.3 scales_1.1.0 classInt_0.4-2 callr_3.4.0
## [25] stringr_1.4.0 digest_0.6.23 foreign_0.8-71 minqa_1.2.4
## [29] rmarkdown_2.0 jpeg_0.1-8.1 pkgconfig_2.0.3 htmltools_0.4.0
## [33] sessioninfo_1.1.1 highr_0.8 farver_2.0.1 magrittr_1.5
## [37] Rcpp_1.0.3 munsell_0.5.0 fansi_0.4.0 lifecycle_0.1.0
## [41] stringi_1.4.3 yaml_2.2.0 pkgbuild_1.0.6 plyr_1.8.5
## [45] crayon_1.3.4 lattice_0.20-38 splines_3.6.1 zeallot_0.1.0
## [49] knitr_1.26 ps_1.3.0 pillar_1.4.3 boot_1.3-22
## [53] rjson_0.2.20 codetools_0.2-16 pkgload_1.0.2 glue_1.3.1
## [57] evaluate_0.14 remotes_2.1.0 png_0.1-7 vctrs_0.2.1
## [61] nloptr_1.2.1 testthat_2.3.1 RgoogleMaps_1.4.5 gtable_0.3.0
## [65] purrr_0.3.3 tidyr_1.0.0 assertthat_0.2.1 xfun_0.11
## [69] e1071_1.7-3 class_7.3-15 tibble_2.1.3 memoise_1.1.0
## [73] units_0.6-5 ellipsis_0.3.0Cite packages:
citation("ggplot2")##
## To cite ggplot2 in publications, please use:
##
## H. Wickham. ggplot2: Elegant Graphics for Data Analysis.
## Springer-Verlag New York, 2016.
##
## A BibTeX entry for LaTeX users is
##
## @Book{,
## author = {Hadley Wickham},
## title = {ggplot2: Elegant Graphics for Data Analysis},
## publisher = {Springer-Verlag New York},
## year = {2016},
## isbn = {978-3-319-24277-4},
## url = {https://ggplot2.tidyverse.org},
## }citation("lme4")##
## To cite lme4 in publications use:
##
## Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).
## Fitting Linear Mixed-Effects Models Using lme4. Journal of
## Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {Fitting Linear Mixed-Effects Models Using {lme4}},
## author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
## journal = {Journal of Statistical Software},
## year = {2015},
## volume = {67},
## number = {1},
## pages = {1--48},
## doi = {10.18637/jss.v067.i01},
## }citation("dplyr")##
## To cite package 'dplyr' in publications use:
##
## Hadley Wickham, Romain François, Lionel Henry and Kirill Müller
## (2019). dplyr: A Grammar of Data Manipulation. R package version
## 0.8.3. https://CRAN.R-project.org/package=dplyr
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {dplyr: A Grammar of Data Manipulation},
## author = {Hadley Wickham and Romain François and Lionel Henry and Kirill Müller},
## year = {2019},
## note = {R package version 0.8.3},
## url = {https://CRAN.R-project.org/package=dplyr},
## }citation("data.table")##
## To cite package 'data.table' in publications use:
##
## Matt Dowle and Arun Srinivasan (2019). data.table: Extension of
## `data.frame`. R package version 1.12.8.
## https://CRAN.R-project.org/package=data.table
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {data.table: Extension of `data.frame`},
## author = {Matt Dowle and Arun Srinivasan},
## year = {2019},
## note = {R package version 1.12.8},
## url = {https://CRAN.R-project.org/package=data.table},
## }citation("numDeriv")##
## To cite package 'numDeriv' in publications use:
##
## Paul Gilbert and Ravi Varadhan (2019). numDeriv: Accurate Numerical
## Derivatives. R package version 2016.8-1.1.
## https://CRAN.R-project.org/package=numDeriv
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {numDeriv: Accurate Numerical Derivatives},
## author = {Paul Gilbert and Ravi Varadhan},
## year = {2019},
## note = {R package version 2016.8-1.1},
## url = {https://CRAN.R-project.org/package=numDeriv},
## }
##
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.citation("sf")##
## To cite package sf in publications, please use:
##
## Pebesma, E., 2018. Simple Features for R: Standardized Support for
## Spatial Vector Data. The R Journal 10 (1), 439-446,
## https://doi.org/10.32614/RJ-2018-009
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## author = {Edzer Pebesma},
## title = {{Simple Features for R: Standardized Support for Spatial Vector Data}},
## year = {2018},
## journal = {{The R Journal}},
## doi = {10.32614/RJ-2018-009},
## url = {https://doi.org/10.32614/RJ-2018-009},
## pages = {439--446},
## volume = {10},
## number = {1},
## }citation("mapdata")##
## To cite package 'mapdata' in publications use:
##
## Original S code by Richard A. Becker and Allan R. Wilks. R version by
## Ray Brownrigg. (2018). mapdata: Extra Map Databases. R package
## version 2.3.0. https://CRAN.R-project.org/package=mapdata
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {mapdata: Extra Map Databases},
## author = {Original S code by Richard A. Becker and Allan R. Wilks. R version by Ray Brownrigg.},
## year = {2018},
## note = {R package version 2.3.0},
## url = {https://CRAN.R-project.org/package=mapdata},
## }
##
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.citation("ggmap")##
## To cite ggmap in publications, please use:
##
## D. Kahle and H. Wickham. ggmap: Spatial Visualization with ggplot2.
## The R Journal, 5(1), 144-161. URL
## http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## author = {David Kahle and Hadley Wickham},
## title = {ggmap: Spatial Visualization with ggplot2},
## journal = {The R Journal},
## year = {2013},
## volume = {5},
## number = {1},
## pages = {144--161},
## url = {https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf},
## }