chapter1
Chapter 1 Analysis
Reading in data
## Load packages
library(dplyr) # data wrangling##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 3.6.2library(lme4) # glmer## Warning: package 'lme4' was built under R version 3.6.2## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 3.6.2library(MuMIn) # r.squaredGLMM## Warning: package 'MuMIn' was built under R version 3.6.2library(emmeans) #pairwise## Warning: package 'emmeans' was built under R version 3.6.2# Read the data
# all lidar only has primary unit
all<- read.csv("/Users/NonisA/Documents/Feb 18th - finalizing stats/Chapter.1/data/species.age1.sample.events.csv")
all$cover.dom<-as.character(all$cover.dom)
all$cover.dom[is.na(all$cover.dom)] <- "none"
all$cover.dom<-as.factor(all$cover.dom)
all$sub.dom<-as.factor(all$sub.dom)
all$level<-as.factor(all$level)
all$stream<-as.factor(all$stream)
all$location<-as.factor(all$location)
all$sample<-as.factor(all$sample)
#### ordering of the categorical factors; how do i choose to order them? does order matter?
all$sub.dom <- factor(all$sub.dom, levels = c("C", "B", "G", "R", "F"))
all$cover.dom <- factor(all$cover.dom, levels = c("none", "l.wood", "s.wood", "OV", "UCB", "B"))
all$stream <- factor(all$stream, levels = c("rainbow.creek", "head.water.creek", "elk.creek", "view.creek", "headache.creek"))
all$habitat <- factor(all$habitat, levels = c("POOL", "G", "RI", "C", "STEP"))
# set up labels for plots;
all$plot<- paste(all$stream, all$location, sep="-")
all$plot<-as.factor(all$plot)
str(all)## 'data.frame': 1136 obs. of 24 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ stream : Factor w/ 5 levels "rainbow.creek",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ location : Factor w/ 67 levels "1","10","10A",..: 12 14 18 20 28 30 51 59 64 3 ...
## $ level : Factor w/ 3 levels "P","S","T": 1 1 3 1 1 1 1 1 1 3 ...
## $ sample : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ total.fish : int 2 1 1 1 1 1 2 3 1 1 ...
## $ density.m2 : num 0.0275 0.0245 0.0581 0.085 0.0608 ...
## $ fish.p.a : int 1 1 1 1 1 1 1 1 1 1 ...
## $ species : Factor w/ 2 levels "DV","RB": 2 2 2 2 2 2 2 2 2 2 ...
## $ habitat : Factor w/ 5 levels "POOL","G","RI",..: 2 3 1 1 1 1 2 1 1 1 ...
## $ gradient : num 0.5 1 0 0 0 0 1 0 0 0 ...
## $ bfw : num 8 8 8.2 4.2 4.7 6.6 2.9 4.4 6.9 4.1 ...
## $ sub.dom : Factor w/ 5 levels "C","B","G","R",..: 3 1 2 2 4 2 1 1 2 1 ...
## $ cover.dom : Factor w/ 6 levels "none","l.wood",..: 4 4 4 4 6 4 4 4 2 2 ...
## $ pool.area : num NA NA 4.41 11.76 12.6 ...
## $ wetted.width.m: num 3.3 2.1 2.1 3.7 3.6 6.5 1.9 3.2 3.7 1 ...
## $ mean.depth.m : num 0.263 0.143 0.09 0.367 0.36 ...
## $ total.wood : int 0 0 0 0 0 2 0 0 1 1 ...
## $ small : int 0 0 0 0 0 0 0 0 0 0 ...
## $ medium : int 0 0 0 0 0 0 0 0 0 0 ...
## $ large : int 0 0 0 0 0 2 0 0 1 1 ...
## $ wood.m2 : num 0 0 0 0 0 ...
## $ wood.p.a : int 0 0 0 0 0 1 0 0 1 1 ...
## $ plot : Factor w/ 142 levels "elk.creek-1",..: 98 99 101 102 106 107 108 113 116 19 ...head(all)## X stream location level sample total.fish density.m2 fish.p.a species
## 1 1 rainbow.creek 12 P 1 2 0.02747253 1 RB
## 2 2 rainbow.creek 13 P 1 1 0.02450980 1 RB
## 3 3 rainbow.creek 14A T 1 1 0.05807201 1 RB
## 4 4 rainbow.creek 15 P 1 1 0.08503401 1 RB
## 5 5 rainbow.creek 19 P 1 1 0.06079027 1 RB
## 6 6 rainbow.creek 2 P 1 1 0.04329004 1 RB
## habitat gradient bfw sub.dom cover.dom pool.area wetted.width.m mean.depth.m
## 1 G 0.5 8.0 G OV NA 3.3 0.2633
## 2 RI 1.0 8.0 C OV NA 2.1 0.1433
## 3 POOL 0.0 8.2 B OV 4.41 2.1 0.0900
## 4 POOL 0.0 4.2 B OV 11.76 3.7 0.3667
## 5 POOL 0.0 4.7 R B 12.60 3.6 0.3600
## 6 POOL 0.0 6.6 B OV 23.10 6.5 0.3400
## total.wood small medium large wood.m2 wood.p.a plot
## 1 0 0 0 0 0.00000000 0 rainbow.creek-12
## 2 0 0 0 0 0.00000000 0 rainbow.creek-13
## 3 0 0 0 0 0.00000000 0 rainbow.creek-14A
## 4 0 0 0 0 0.00000000 0 rainbow.creek-15
## 5 0 0 0 0 0.00000000 0 rainbow.creek-19
## 6 2 0 0 2 0.08658009 1 rainbow.creek-2RAINBOW TROUT ANALYSIS 250 observations (total units 4 sample events) Age 1+
Gradient and bf –> both significant
no zeros within this analysis
rb<-all%>%
filter(species == "RB")%>%
filter(total.fish != "0")
rb$log.density<-log(rb$density.m2)
#explore the data a bit: general shape per stream and sample
ggplot(data = rb) +
geom_density(aes(x = log.density), col = "darkblue", fill = "lightblue") +
theme_bw()+
facet_wrap(~stream)
all.lm<- lmer(log.density ~ gradient+ bfw # significant
+ (1|stream/plot), data=rb)
summary(all.lm) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ gradient + bfw + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 468.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5681 -0.4165 0.0349 0.4086 2.2509
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.35272 0.5939
## stream (Intercept) 0.07122 0.2669
## Residual 0.19461 0.4411
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -2.03540 0.25394 -8.015
## gradient -0.25068 0.04807 -5.215
## bfw -0.15296 0.02719 -5.626
##
## Correlation of Fixed Effects:
## (Intr) gradnt
## gradient -0.160
## bfw -0.820 0.022#pseudo R^2
r.squaredGLMM(all.lm) #MuMIN package## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.3623537 0.7993864#100*(null_model$deviance - all.lm$deviance)/null_model$deviance
# model explains 76% of variation in density
# RMSE
rb$log.density.pred <- predict(all.lm, rb, type = "response")
sqrt(mean((rb$log.density - rb$log.density.pred)^2))## [1] 0.3716258# Intercept
#exp(X)
# X variables
#(exp(X)-1)*100
# or XX% decrease in density per one unit increase in gradient
#Check the Model Assumptions
# Linearity
### Residuals & heteroskedasticity
par(mfrow=c(2,2), cex=1.0, mai=c(1.0,1.0,0.6,0.6))
plot(all.lm)
par(mfrow=c(1,1), cex=1.0, mai=c(1.0,1.0,1.0,1.0))
# Normality
hist(residuals(all.lm))
# QQ Plot
dres <- resid(all.lm)
qqnorm(dres, main="Deviance residuals", las=1); qqline(dres)
## Partial F-test to test significance of terms
# bfw
all.null<- lmer(log.density ~ 1 + bfw
+ (1|stream/plot), data=rb)
summary(all.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + bfw + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 488
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5454 -0.4205 0.0266 0.3922 2.3152
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.49485 0.7035
## stream (Intercept) 0.02986 0.1728
## Residual 0.19535 0.4420
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -2.18144 0.24821 -8.789
## bfw -0.15624 0.02885 -5.416
##
## Correlation of Fixed Effects:
## (Intr)
## bfw -0.888anova(all.lm, all.null)## refitting model(s) with ML (instead of REML)## Data: rb
## Models:
## all.null: log.density ~ 1 + bfw + (1 | stream/plot)
## all.lm: log.density ~ gradient + bfw + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## all.null 5 489.89 507.50 -239.94 479.89
## all.lm 6 468.83 489.95 -228.41 456.83 23.063 1 1.568e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#gradient
all.null<- lmer(log.density ~ 1 + gradient
+ (1|stream/plot), data=rb)
summary(all.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + gradient + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 489
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5436 -0.3857 0.0052 0.4539 2.1774
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.4686 0.6845
## stream (Intercept) 0.2535 0.5035
## Residual 0.1957 0.4424
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.2158 0.2439 -13.183
## gradient -0.2510 0.0543 -4.622
##
## Correlation of Fixed Effects:
## (Intr)
## gradient -0.164anova(all.lm, all.null)## refitting model(s) with ML (instead of REML)## Data: rb
## Models:
## all.null: log.density ~ 1 + gradient + (1 | stream/plot)
## all.lm: log.density ~ gradient + bfw + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## all.null 5 493.89 511.50 -241.94 483.89
## all.lm 6 468.83 489.95 -228.41 456.83 27.064 1 1.968e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wood investigation –> no significance
log.density & total.fish (CPUE) tried
rb$wood.p.a <- factor(rb$wood.p.a, levels = c("0", "1"))
wood<-lmer(log.density ~ wood.p.a + (1|stream/plot), data=rb)
summary(wood)## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ wood.p.a + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 503.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5162 -0.3799 0.0183 0.4361 2.2641
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.5947 0.7712
## stream (Intercept) 0.1975 0.4444
## Residual 0.1960 0.4427
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.2679 0.2349 -13.915
## wood.p.a1 -0.2843 0.1797 -1.583
##
## Correlation of Fixed Effects:
## (Intr)
## wood.p.a1 -0.354anova(wood)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## wood.p.a 1 0.49087 0.49087 2.5047wood.null<- lmer(log.density ~ 1
+ (1|stream/plot), data=rb)
summary(wood.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 504.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5336 -0.3569 0.0288 0.4766 2.2426
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.6062 0.7786
## stream (Intercept) 0.2090 0.4572
## Residual 0.1959 0.4426
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.3998 0.2252 -15.1anova(wood, wood.null)## refitting model(s) with ML (instead of REML)## Data: rb
## Models:
## wood.null: log.density ~ 1 + (1 | stream/plot)
## wood: log.density ~ wood.p.a + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## wood.null 4 510.91 525.00 -251.46 502.91
## wood 5 510.36 527.97 -250.18 500.36 2.5461 1 0.1106(plots.lsm <- emmeans (wood, pairwise~wood.p.a,
adjust="bonferroni", side="two-sided"))## $emmeans
## wood.p.a emmean SE df lower.CL upper.CL
## 0 -3.27 0.235 5.16 -3.87 -2.67
## 1 -3.55 0.241 5.56 -4.15 -2.95
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 0 - 1 0.284 0.181 84.7 1.573 0.1194
##
## Degrees-of-freedom method: kenward-rogerNOTE: wood should be an influential factor; how to make this work?
Cover investigation –> not significant
density and CPUE
cover<-lmer(log.density ~ cover.dom + (1|stream/plot), data=rb)
summary(cover)## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ cover.dom + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 499.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5856 -0.4415 0.0319 0.4553 2.2130
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.5909 0.7687
## stream (Intercept) 0.2865 0.5353
## Residual 0.1953 0.4419
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.61524 0.64322 -5.621
## cover.doml.wood 0.53570 0.61755 0.867
## cover.doms.wood 0.36076 0.64478 0.560
## cover.domOV 0.05156 0.61040 0.084
## cover.domUCB 0.66615 0.72312 0.921
## cover.domB -0.03884 0.65700 -0.059
##
## Correlation of Fixed Effects:
## (Intr) cvr.dml. cvr.dms. cvr.OV cv.UCB
## covr.dml.wd -0.886
## covr.dms.wd -0.842 0.873
## cover.domOV -0.906 0.932 0.891
## cover.dmUCB -0.733 0.767 0.722 0.764
## cover.domB -0.849 0.872 0.835 0.896 0.706anova(cover)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## cover.dom 5 1.2661 0.25322 1.2967cover.null<- lmer(log.density ~ 1
+ (1|stream/plot), data=rb)
summary(cover.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 504.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5336 -0.3569 0.0288 0.4766 2.2426
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.6062 0.7786
## stream (Intercept) 0.2090 0.4572
## Residual 0.1959 0.4426
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.3998 0.2252 -15.1anova(cover, cover.null)## refitting model(s) with ML (instead of REML)## Data: rb
## Models:
## cover.null: log.density ~ 1 + (1 | stream/plot)
## cover: log.density ~ cover.dom + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## cover.null 4 510.91 525.00 -251.46 502.91
## cover 9 514.67 546.36 -248.34 496.67 6.2409 5 0.2835(plots.lsm <- emmeans (cover, pairwise~cover.dom,
adjust="bonferroni", side="two-sided"))## $emmeans
## cover.dom emmean SE df lower.CL upper.CL
## none -3.62 0.645 64.52 -4.90 -2.33
## l.wood -3.08 0.302 7.27 -3.79 -2.37
## s.wood -3.25 0.363 14.10 -4.03 -2.48
## OV -3.56 0.274 5.00 -4.27 -2.86
## UCB -2.95 0.512 36.21 -3.99 -1.91
## B -3.65 0.358 13.61 -4.42 -2.88
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## none - l.wood -0.5357 0.618 81.8 -0.866 1.0000
## none - s.wood -0.3608 0.646 80.5 -0.559 1.0000
## none - OV -0.0516 0.612 82.4 -0.084 1.0000
## none - UCB -0.6661 0.725 79.1 -0.919 1.0000
## none - B 0.0388 0.660 81.7 0.059 1.0000
## l.wood - s.wood 0.1749 0.320 77.1 0.547 1.0000
## l.wood - OV 0.4841 0.227 83.9 2.129 0.5426
## l.wood - UCB -0.1305 0.473 77.4 -0.276 1.0000
## l.wood - B 0.5745 0.326 79.3 1.760 1.0000
## s.wood - OV 0.3092 0.297 76.4 1.042 1.0000
## s.wood - UCB -0.3054 0.521 76.5 -0.586 1.0000
## s.wood - B 0.3996 0.376 75.7 1.063 1.0000
## OV - UCB -0.6146 0.479 78.7 -1.284 1.0000
## OV - B 0.0904 0.293 76.0 0.308 1.0000
## UCB - B 0.7050 0.542 78.3 1.300 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 15 testssubstrate investigation –> significant!! Pairwise.. no pairs significantly different?
sub<-lmer(log.density ~ sub.dom + (1|stream/plot), data=rb)
summary(sub)## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ sub.dom + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 495.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5561 -0.3389 0.0398 0.4617 2.2641
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.5632 0.7505
## stream (Intercept) 0.2266 0.4761
## Residual 0.1953 0.4420
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.2992 0.2458 -13.424
## sub.domB -0.4499 0.2225 -2.023
## sub.domG -0.2068 0.2795 -0.740
## sub.domR 0.4441 0.3909 1.136
## sub.domF 0.9144 0.5792 1.579
##
## Correlation of Fixed Effects:
## (Intr) sb.dmB sb.dmG sb.dmR
## sub.domB -0.292
## sub.domG -0.194 0.235
## sub.domR -0.161 0.211 0.181
## sub.domF -0.076 0.070 0.088 0.057anova(sub)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## sub.dom 4 1.9085 0.47713 2.4425sub.null<- lmer(log.density ~ 1
+ (1|stream/plot), data=rb)
summary(sub.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 504.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5336 -0.3569 0.0288 0.4766 2.2426
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.6062 0.7786
## stream (Intercept) 0.2090 0.4572
## Residual 0.1959 0.4426
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.3998 0.2252 -15.1anova(sub, sub.null)## refitting model(s) with ML (instead of REML)## Data: rb
## Models:
## sub.null: log.density ~ 1 + (1 | stream/plot)
## sub: log.density ~ sub.dom + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## sub.null 4 510.91 525.00 -251.46 502.91
## sub 8 509.26 537.43 -246.63 493.26 9.6494 4 0.04677 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(plots.lsm <- emmeans (sub, pairwise~sub.dom,
adjust="bonferroni", side="two-sided"))## $emmeans
## sub.dom emmean SE df lower.CL upper.CL
## C -3.30 0.247 4.86 -3.94 -2.66
## B -3.75 0.280 8.06 -4.39 -3.10
## G -3.51 0.336 15.59 -4.22 -2.79
## R -2.86 0.429 32.99 -3.73 -1.98
## F -2.38 0.614 65.93 -3.61 -1.16
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## C - B 0.450 0.225 81.3 1.997 0.4916
## C - G 0.207 0.281 84.3 0.736 1.0000
## C - R -0.444 0.395 78.1 -1.125 1.0000
## C - F -0.914 0.580 76.9 -1.576 1.0000
## B - G -0.243 0.316 84.1 -0.770 1.0000
## B - R -0.894 0.409 77.4 -2.186 0.3181
## B - F -1.364 0.609 78.2 -2.240 0.2792
## G - R -0.651 0.439 78.6 -1.484 1.0000
## G - F -1.121 0.622 78.0 -1.804 0.7512
## R - F -0.470 0.683 77.5 -0.689 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 10 tests# Intercept
#exp(X)
# X variables
#(exp(X)-1)*100
# or XX% decrease in density per one unit increase in gradient
#Check the Model Assumptions
# Linearity
### Residuals & heteroskedasticity
par(mfrow=c(2,2), cex=1.0, mai=c(1.0,1.0,0.6,0.6))
plot(sub)
par(mfrow=c(1,1), cex=1.0, mai=c(1.0,1.0,1.0,1.0))
# Normality
hist(residuals(sub))
# QQ Plot
dres <- resid(sub)
qqnorm(dres, main="Deviance residuals", las=1); qqline(dres)
Field Habitat Selection
habitat<-lmer(log.density ~ habitat + (1|stream/plot), data=rb)
summary(habitat)## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ habitat + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 471.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5906 -0.3213 -0.0039 0.4688 2.2477
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.3799 0.6163
## stream (Intercept) 0.1679 0.4098
## Residual 0.1970 0.4438
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -2.9184 0.2142 -13.624
## habitatG -0.7386 0.1806 -4.089
## habitatRI -0.9935 0.1980 -5.017
## habitatC -1.3623 0.2983 -4.566
##
## Correlation of Fixed Effects:
## (Intr) habttG hbttRI
## habitatG -0.315
## habitatRI -0.275 0.342
## habitatC -0.152 0.178 0.170anova(habitat)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## habitat 3 8.5858 2.862 14.531habitat.null<- lmer(log.density ~ 1
+ (1|stream/plot), data=rb)
summary(habitat.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + (1 | stream/plot)
## Data: rb
##
## REML criterion at convergence: 504.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5336 -0.3569 0.0288 0.4766 2.2426
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.6062 0.7786
## stream (Intercept) 0.2090 0.4572
## Residual 0.1959 0.4426
## Number of obs: 250, groups: plot:stream, 89; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.3998 0.2252 -15.1anova(habitat, habitat.null)## refitting model(s) with ML (instead of REML)## Data: rb
## Models:
## habitat.null: log.density ~ 1 + (1 | stream/plot)
## habitat: log.density ~ habitat + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## habitat.null 4 510.91 525.00 -251.46 502.91
## habitat 7 480.46 505.11 -233.23 466.46 36.454 3 6.004e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(plots.lsm <- emmeans (habitat, pairwise~habitat,
adjust="bonferroni", side="two-sided"))## $emmeans
## habitat emmean SE df lower.CL upper.CL
## POOL -2.92 0.215 5.27 -3.46 -2.37
## G -3.66 0.233 7.34 -4.20 -3.11
## RI -3.91 0.249 9.39 -4.47 -3.35
## C -4.28 0.341 27.44 -4.98 -3.58
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## POOL - G 0.739 0.181 79.7 4.072 0.0007
## POOL - RI 0.994 0.199 83.0 4.999 <.0001
## POOL - C 1.362 0.299 78.2 4.558 0.0001
## G - RI 0.255 0.218 82.3 1.169 1.0000
## G - C 0.624 0.322 78.8 1.938 0.3370
## RI - C 0.369 0.330 79.4 1.117 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 6 tests# Intercept
#exp(X)
# X variables
#(exp(X)-1)*100
# or XX% decrease in density per one unit increase in gradient
#Check the Model Assumptions
# Linearity
### Residuals & heteroskedasticity
par(mfrow=c(2,2), cex=1.0, mai=c(1.0,1.0,0.6,0.6))
plot(habitat)
par(mfrow=c(1,1), cex=1.0, mai=c(1.0,1.0,1.0,1.0))
# Normality
hist(residuals(habitat))
# QQ Plot
dres <- resid(habitat)
qqnorm(dres, main="Deviance residuals", las=1); qqline(dres)
########################## ## Field Pool area –> significant #########################
pool.area.data<-rb%>%
filter(habitat == "POOL")
foo<-pool.area.data%>%
distinct(plot, .keep_all = TRUE) # 42 distinct pools (all levels)
pool.area<-lmer(total.fish ~ pool.area + (1|stream/plot), data=pool.area.data)
summary(pool.area)## Linear mixed model fit by REML ['lmerMod']
## Formula: total.fish ~ pool.area + (1 | stream/plot)
## Data: pool.area.data
##
## REML criterion at convergence: 460.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5177 -0.4065 -0.1390 0.1709 4.1074
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.5652 0.7518
## stream (Intercept) 0.1157 0.3402
## Residual 2.0811 1.4426
## Number of obs: 120, groups: plot:stream, 42; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.187275 0.300543 3.950
## pool.area 0.037764 0.008781 4.301
##
## Correlation of Fixed Effects:
## (Intr)
## pool.area -0.586anova(pool.area)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## pool.area 1 38.489 38.489 18.494pool.null<- lmer(total.fish ~ 1
+ (1|stream/plot), data=pool.area.data)
summary(pool.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: total.fish ~ 1 + (1 | stream/plot)
## Data: pool.area.data
##
## REML criterion at convergence: 468.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1694 -0.3584 -0.1810 0.1465 4.4465
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 1.1133 1.0551
## stream (Intercept) 0.1622 0.4028
## Residual 2.0862 1.4444
## Number of obs: 120, groups: plot:stream, 42; stream, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.9025 0.2898 6.564anova(pool.area, pool.null)## refitting model(s) with ML (instead of REML)## Data: pool.area.data
## Models:
## pool.null: total.fish ~ 1 + (1 | stream/plot)
## pool.area: total.fish ~ pool.area + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## pool.null 4 475.35 486.50 -233.68 467.35
## pool.area 5 462.09 476.03 -226.05 452.09 15.258 1 9.377e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Intercept
#exp(X)
# X variables
#(exp(X)-1)*100
# or XX% decrease in density per one unit increase in gradient
#Check the Model Assumptions
# some of the assumptions look funky
# Linearity
### Residuals & heteroskedasticity
par(mfrow=c(2,2), cex=1.0, mai=c(1.0,1.0,0.6,0.6))
plot(pool.area)
par(mfrow=c(1,1), cex=1.0, mai=c(1.0,1.0,1.0,1.0))
# Normality
hist(residuals(pool.area))
# QQ Plot
dres <- resid(pool.area)
qqnorm(dres, main="Deviance residuals", las=1); qqline(dres)
Dolly Varden Same steps and analysis done with dolly varden age 1+ 57 observations (total units over 4 sample)
dv<-all %>%
filter(species == "DV")%>%
filter(total.fish != "0")
dv$log.density<-log(dv$density.m2+0.01)
#explore the data a bit: general shape per stream and sample
ggplot(data = dv) +
geom_density(aes(x = log.density), col = "darkblue", fill = "lightblue") +
theme_bw()+
facet_wrap(~stream)
Results: Gradient and BFW significant Habitat (pool result) Pool area (assumptions also look funky)
Substrate
sub<-lmer(log.density ~ sub.dom + (1|stream/plot), data=dv)
summary(sub)## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ sub.dom + (1 | stream/plot)
## Data: dv
##
## REML criterion at convergence: 77
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15446 -0.29826 -0.02394 0.33907 2.45853
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.20478 0.4525
## stream (Intercept) 0.14672 0.3830
## Residual 0.06956 0.2637
## Number of obs: 57, groups: plot:stream, 36; stream, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.2026 0.2190 -14.624
## sub.domB -0.5493 0.2200 -2.497
## sub.domG -0.2603 0.2793 -0.932
##
## Correlation of Fixed Effects:
## (Intr) sb.dmB
## sub.domB -0.205
## sub.domG -0.182 0.159anova(sub)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## sub.dom 2 0.45398 0.22699 3.2635sub.null<- lmer(log.density ~ 1
+ (1|stream/plot), data=dv)
summary(sub.null) ## Linear mixed model fit by REML ['lmerMod']
## Formula: log.density ~ 1 + (1 | stream/plot)
## Data: dv
##
## REML criterion at convergence: 80.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15226 -0.31996 -0.03135 0.32237 2.53439
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot:stream (Intercept) 0.25252 0.5025
## stream (Intercept) 0.10319 0.3212
## Residual 0.06739 0.2596
## Number of obs: 57, groups: plot:stream, 36; stream, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.3290 0.1878 -17.73anova(sub, sub.null)## refitting model(s) with ML (instead of REML)## Data: dv
## Models:
## sub.null: log.density ~ 1 + (1 | stream/plot)
## sub: log.density ~ sub.dom + (1 | stream/plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## sub.null 4 87.287 95.459 -39.643 79.287
## sub 6 85.508 97.766 -36.754 73.508 5.7793 2 0.0556 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1(plots.lsm <- emmeans (sub, pairwise~sub.dom,
adjust="bonferroni", side="two-sided"))## $emmeans
## sub.dom emmean SE df lower.CL upper.CL
## C -3.20 0.219 3.28 -3.87 -2.54
## B -3.75 0.280 7.62 -4.40 -3.10
## G -3.46 0.330 11.85 -4.18 -2.74
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## C - B 0.549 0.223 28.7 2.468 0.0594
## C - G 0.260 0.289 34.2 0.902 1.0000
## B - G -0.289 0.337 33.1 -0.858 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 3 tests