Bat data models

Set up stuff

rm(list=ls())
pacman::p_load(plyr,dplyr,ggplot2,lme4,car,multcomp, tidyr, reshape2)
setwd("C:/Users/Katy/Desktop/thesisr")

Data organization

# Read in data, add column calculating Calls Per Night
data.allparks <- read.csv("./MixedModel_allparks.csv")
data.allparks$Year <- as.factor(data.allparks$Year)
data.allparks$CallsPN <- round(data.allparks$Calls/data.allparks$Nights, 3)
str(data.allparks)

## 'data.frame':    2460 obs. of  7 variables:
##  $ Park   : Factor w/ 8 levels "APIS","GRPO",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Site   : Factor w/ 185 levels "APIS001A","APIS002A",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ Year   : Factor w/ 2 levels "2016","2017": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Nights : int  12 12 12 12 12 12 9 9 9 9 ...
##  $ Species: Factor w/ 9 levels "EPFU","LABO",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ Calls  : int  130 92 95 1496 2311 1 8 2 7 47 ...
##  $ CallsPN: num  10.83 7.67 7.92 124.67 192.58 ...

head(data.allparks)

##   Park     Site Year Nights Species Calls CallsPN
## 1 APIS APIS001A 2017     12    EPFU   130  10.833
## 2 APIS APIS001A 2017     12    LABO    92   7.667
## 3 APIS APIS001A 2017     12    LACI    95   7.917
## 4 APIS APIS001A 2017     12    LANO  1496 124.667
## 5 APIS APIS001A 2017     12    MYLU  2311 192.583
## 6 APIS APIS001A 2017     12    MYSE     1   0.083

# Temporarily subset the data by year. Rename columns to include year information.
# and remove the actual "Year" column. Then merge all data back to one wide
# data frame, with rows for each combination of Park, Site, and Species. Add a new
# column showing difference between Calls Per Night in 2017 vs 2016. 
data.allparks.2016 <- data.allparks %>% 
                      subset(Year == "2016") %>%
                      rename(Nights.2016 = Nights, Calls.2016 = Calls, 
                             CallsPN.2016 = CallsPN) %>%
                      subset(select = -c(Year))

data.allparks.2017 <- data.allparks %>% 
                      subset(Year == "2017") %>%
                      rename(Nights.2017 = Nights, Calls.2017 = Calls, 
                           CallsPN.2017 = CallsPN) %>%
                      subset(select = -c(Year))

data.allparks.wide <- merge(data.allparks.2016, data.allparks.2017, 
                            by = c("Park","Site","Species")) %>%
                      mutate(Diff = CallsPN.2017 - CallsPN.2016)
str(data.allparks.wide)

## 'data.frame':    1230 obs. of  10 variables:
##  $ Park        : Factor w/ 8 levels "APIS","GRPO",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Site        : Factor w/ 185 levels "APIS001A","APIS002A",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ Species     : Factor w/ 9 levels "EPFU","LABO",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ Nights.2016 : int  7 7 7 7 7 7 10 10 10 10 ...
##  $ Calls.2016  : int  17 91 35 144 2380 16 4 1 4 33 ...
##  $ CallsPN.2016: num  2.43 13 5 20.57 340 ...
##  $ Nights.2017 : int  12 12 12 12 12 12 9 9 9 9 ...
##  $ Calls.2017  : int  130 92 95 1496 2311 1 8 2 7 47 ...
##  $ CallsPN.2017: num  10.83 7.67 7.92 124.67 192.58 ...
##  $ Diff        : num  8.4 -5.33 2.92 104.1 -147.42 ...

head(data.allparks.wide)

##   Park     Site Species Nights.2016 Calls.2016 CallsPN.2016 Nights.2017
## 1 APIS APIS001A    EPFU           7         17        2.429          12
## 2 APIS APIS001A    LABO           7         91       13.000          12
## 3 APIS APIS001A    LACI           7         35        5.000          12
## 4 APIS APIS001A    LANO           7        144       20.571          12
## 5 APIS APIS001A    MYLU           7       2380      340.000          12
## 6 APIS APIS001A    MYSE           7         16        2.286          12
##   Calls.2017 CallsPN.2017     Diff
## 1        130       10.833    8.404
## 2         92        7.667   -5.333
## 3         95        7.917    2.917
## 4       1496      124.667  104.096
## 5       2311      192.583 -147.417
## 6          1        0.083   -2.203

Plot 1 - histogram of counts

gg.diff1 <- ggplot(data.allparks.wide, aes(data.allparks.wide$Diff)) + 
            geom_histogram(binwidth = 20, colour="black", fill="gray") + 
            theme_bw()
gg.diff1

Plot 2 - histogram of densities, density estimate, and probability density function

gg.diff2 <- ggplot(data.allparks.wide, aes(x=Diff)) + 
            geom_histogram(aes(y=..density..), binwidth = 20, 
                           colour="black", fill="gray") + 
            geom_density(alpha=0.4, color= "coral2", size=.7, fill= "coral2") + 
            theme_bw() 

mean(data.allparks.wide$Diff)

## [1] -6.78478

sd(data.allparks.wide$Diff)

## [1] 60.77969

gg.diff3 <- gg.diff2 + 
            stat_function(fun = dnorm, colour="blue", size=1.2,
                    args = list(mean = -6.78478, sd = 60.77969))
gg.diff3

# normal distribution seems to fit okay

Test significance of “Park” variable

# Using a linear model
mod.null <- lm(data.allparks.wide$Diff ~ 1 , data= data.allparks.wide)
mod.park <- lm(data.allparks.wide$Diff ~ Park, data= data.allparks.wide)
anova(mod.null, mod.park)

## Analysis of Variance Table
## 
## Model 1: data.allparks.wide$Diff ~ 1
## Model 2: data.allparks.wide$Diff ~ Park
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
## 1   1229 4540136                                
## 2   1222 4469108  7     71028 2.7745 0.007275 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Using a linear mixed model
mod.null2 <- lmer(data.allparks.wide$Diff ~ 1 + (1|Site), data= data.allparks.wide)
mod.park2 <- lmer(data.allparks.wide$Diff ~ Park + (1|Site), data= data.allparks.wide)
anova(mod.null2, mod.park2)

## refitting model(s) with ML (instead of REML)

## Data: data.allparks.wide
## Models:
## mod.null2: data.allparks.wide$Diff ~ 1 + (1 | Site)
## mod.park2: data.allparks.wide$Diff ~ Park + (1 | Site)
##           Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mod.null2  3 13597 13612 -6795.6    13591                           
## mod.park2 10 13594 13645 -6786.9    13574 17.368      7    0.01517 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# I wasn't sure whether to used the regular linear vs. mixed model for this test? Ran both and both gave a significant p-value. Therefore, conclude that there are significant differences by park, so I will run separate models for each park.

Subset data by parks

park.list <- list("APIS", "GRPO", "INDU", "ISRO","MISS", "SACN", "SLBE", "VOYA")
for (i in park.list) {
  temp <- subset(data.allparks.wide, Park == i)
  temp <- droplevels(temp)
  assign(paste("data",i,sep="."), temp)
        }

Run model and hypothesis tests for each of the 8 parks

Apostle Islands

# This model yields a "singular fit" and lots of warnings when I get the confidence intervals. So still need to figure out what the problem is and how to resolve it.
model.APIS <- lmer(Diff ~ 0 + Species + (1|Site), data=data.APIS)

## singular fit

summary(model.APIS)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.APIS
## 
## REML criterion at convergence: 1582
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.1232 -0.0631  0.0088  0.3139  2.6408 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)    0      0.00   
##  Residual             4856     69.69   
## Number of obs: 144, groups:  Site, 24
## 
## Fixed effects:
##             Estimate Std. Error t value
## SpeciesEPFU   0.4025    14.2250   0.028
## SpeciesLABO  -2.4514    14.2250  -0.172
## SpeciesLACI -24.2901    14.2250  -1.708
## SpeciesLANO   4.1382    14.2250   0.291
## SpeciesMYLU -35.0799    14.2250  -2.466
## SpeciesMYSE  -1.2812    14.2250  -0.090
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU
## SpeciesLABO 0.000                             
## SpeciesLACI 0.000  0.000                      
## SpeciesLANO 0.000  0.000  0.000               
## SpeciesMYLU 0.000  0.000  0.000  0.000        
## SpeciesMYSE 0.000  0.000  0.000  0.000  0.000 
## convergence code: 0
## singular fit

car::Anova(model.APIS)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)
## Species 9.1206  6     0.1669

hyptest.APIS <- glht(model.APIS, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0")))
summary(hyptest.APIS)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.APIS)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)  
## SpeciesEPFU == 0   0.4025    14.2250   0.028   1.0000  
## SpeciesLABO == 0  -2.4514    14.2250  -0.172   1.0000  
## SpeciesLACI == 0 -24.2901    14.2250  -1.708   0.4235  
## SpeciesLANO == 0   4.1382    14.2250   0.291   0.9999  
## SpeciesMYLU == 0 -35.0799    14.2250  -2.466   0.0792 .
## SpeciesMYSE == 0  -1.2812    14.2250  -0.090   1.0000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.APIS, oldNames=FALSE)

## Computing profile confidence intervals ...

## Warning in zeta(shiftpar, start = opt[seqpar1][-w]): slightly lower
## deviances (diff=-2.27374e-13) detected

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep

## Warning in FUN(X[[i]], ...): non-monotonic profile for sd_(Intercept)|Site

## Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
## sd_(Intercept)|Site: falling back to linear interpolation

##                         2.5 %    97.5 %
## sd_(Intercept)|Site   0.00000 20.430679
## sigma                61.04023 76.927177
## SpeciesEPFU         -27.07394 27.878856
## SpeciesLABO         -29.92777 25.025023
## SpeciesLACI         -51.76652  3.186273
## SpeciesLANO         -23.33823 31.614565
## SpeciesMYLU         -62.55627 -7.603477
## SpeciesMYSE         -28.75756 26.195231

Grand Portage

model.GRPO <- lmer(Diff ~ 0 + Species + (1|Site), data=data.GRPO)
summary(model.GRPO)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.GRPO
## 
## REML criterion at convergence: 1087.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.1459 -0.1493 -0.0292  0.2432  2.6616 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)  565.7   23.78   
##  Residual             3656.7   60.47   
## Number of obs: 102, groups:  Site, 17
## 
## Fixed effects:
##             Estimate Std. Error t value
## SpeciesEPFU   0.3791    15.7599   0.024
## SpeciesLABO -16.7073    15.7599  -1.060
## SpeciesLACI   1.6259    15.7599   0.103
## SpeciesLANO  12.9774    15.7599   0.823
## SpeciesMYLU -91.8349    15.7599  -5.827
## SpeciesMYSE -26.4316    15.7599  -1.677
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU
## SpeciesLABO 0.134                             
## SpeciesLACI 0.134  0.134                      
## SpeciesLANO 0.134  0.134  0.134               
## SpeciesMYLU 0.134  0.134  0.134  0.134        
## SpeciesMYSE 0.134  0.134  0.134  0.134  0.134

car::Anova(model.GRPO)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##         Chisq Df Pr(>Chisq)    
## Species 39.18  6    6.6e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hyptest.GRPO <- glht(model.GRPO, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0")))
summary(hyptest.GRPO)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.GRPO)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## SpeciesEPFU == 0   0.3791    15.7599   0.024    1.000    
## SpeciesLABO == 0 -16.7073    15.7599  -1.060    0.864    
## SpeciesLACI == 0   1.6259    15.7599   0.103    1.000    
## SpeciesLANO == 0  12.9774    15.7599   0.823    0.955    
## SpeciesMYLU == 0 -91.8349    15.7599  -5.827   <1e-04 ***
## SpeciesMYSE == 0 -26.4316    15.7599  -1.677    0.436    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.GRPO, oldNames=FALSE)

## Computing profile confidence intervals ...

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

##                          2.5 %    97.5 %
## sd_(Intercept)|Site    0.00000  42.31964
## sigma                 50.84040  68.72344
## SpeciesEPFU          -29.89961  30.65773
## SpeciesLABO          -46.98596  13.57137
## SpeciesLACI          -28.65273  31.90461
## SpeciesLANO          -17.30126  43.25608
## SpeciesMYLU         -122.11361 -61.55627
## SpeciesMYSE          -56.71031   3.84702

Indiana Dunes

model.INDU <- lmer(Diff ~ 0 + Species + (1|Site), data=data.INDU)
summary(model.INDU)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.INDU
## 
## REML criterion at convergence: 2370.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9357 -0.0992  0.0156  0.1085  7.9895 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)   67.39   8.209  
##  Residual             1873.34  43.282  
## Number of obs: 234, groups:  Site, 26
## 
## Fixed effects:
##              Estimate Std. Error t value
## SpeciesEPFU  57.19138    8.63967   6.620
## SpeciesLABO -23.56919    8.63967  -2.728
## SpeciesLACI   6.66581    8.63967   0.772
## SpeciesLANO  14.80400    8.63967   1.713
## SpeciesMYLU  -2.11131    8.63967  -0.244
## SpeciesMYSE  -0.01527    8.63967  -0.002
## SpeciesMYSO  -0.12923    8.63967  -0.015
## SpeciesNYHU  -3.81546    8.63967  -0.442
## SpeciesPESU  -1.07900    8.63967  -0.125
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU SpMYSE SpMYSO SpNYHU
## SpeciesLABO 0.035                                                  
## SpeciesLACI 0.035  0.035                                           
## SpeciesLANO 0.035  0.035  0.035                                    
## SpeciesMYLU 0.035  0.035  0.035  0.035                             
## SpeciesMYSE 0.035  0.035  0.035  0.035  0.035                      
## SpeciesMYSO 0.035  0.035  0.035  0.035  0.035  0.035               
## SpeciesNYHU 0.035  0.035  0.035  0.035  0.035  0.035  0.035        
## SpeciesPESU 0.035  0.035  0.035  0.035  0.035  0.035  0.035  0.035

car::Anova(model.INDU)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)    
## Species 56.178  9   7.26e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hyptest.INDU <- glht(model.INDU, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0", "SpeciesMYSO = 0", "SpeciesNYHU = 0", "SpeciesPESU = 0")))
summary(hyptest.INDU)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.INDU)
## 
## Linear Hypotheses:
##                   Estimate Std. Error z value Pr(>|z|)    
## SpeciesEPFU == 0  57.19138    8.63967   6.620   <1e-05 ***
## SpeciesLABO == 0 -23.56919    8.63967  -2.728   0.0558 .  
## SpeciesLACI == 0   6.66581    8.63967   0.772   0.9945    
## SpeciesLANO == 0  14.80400    8.63967   1.713   0.5565    
## SpeciesMYLU == 0  -2.11131    8.63967  -0.244   1.0000    
## SpeciesMYSE == 0  -0.01527    8.63967  -0.002   1.0000    
## SpeciesMYSO == 0  -0.12923    8.63967  -0.015   1.0000    
## SpeciesNYHU == 0  -3.81546    8.63967  -0.442   0.9999    
## SpeciesPESU == 0  -1.07900    8.63967  -0.125   1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.INDU, oldNames=FALSE)

## Computing profile confidence intervals ...

##                          2.5 %    97.5 %
## sd_(Intercept)|Site   0.000000 16.844888
## sigma                38.668172 46.871222
## SpeciesEPFU          40.517729 73.865040
## SpeciesLABO         -40.242848 -6.895536
## SpeciesLACI         -10.007848 23.339464
## SpeciesLANO          -1.869656 31.477656
## SpeciesMYLU         -18.784964 14.562348
## SpeciesMYSE         -16.688925 16.658387
## SpeciesMYSO         -16.802887 16.544425
## SpeciesNYHU         -20.489117 12.858194
## SpeciesPESU         -17.752656 15.594656

Isle Royale

model.ISRO <- lmer(Diff ~ 0 + Species + (1|Site), data=data.ISRO)
summary(model.ISRO)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.ISRO
## 
## REML criterion at convergence: 1477.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.5007 -0.0222  0.0062  0.0856  1.4926 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)   69.62   8.344  
##  Residual             2211.23  47.024  
## Number of obs: 144, groups:  Site, 24
## 
## Fixed effects:
##             Estimate Std. Error t value
## SpeciesEPFU   -1.342      9.749  -0.138
## SpeciesLABO   -5.914      9.749  -0.607
## SpeciesLACI   -1.212      9.749  -0.124
## SpeciesLANO   -2.188      9.749  -0.224
## SpeciesMYLU  -72.283      9.749  -7.415
## SpeciesMYSE   -3.560      9.749  -0.365
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU
## SpeciesLABO 0.031                             
## SpeciesLACI 0.031  0.031                      
## SpeciesLANO 0.031  0.031  0.031               
## SpeciesMYLU 0.031  0.031  0.031  0.031        
## SpeciesMYSE 0.031  0.031  0.031  0.031  0.031

car::Anova(model.ISRO)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)    
## Species 55.162  6  4.299e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hyptest.ISRO <- glht(model.ISRO, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0")))
summary(hyptest.ISRO)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.ISRO)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## SpeciesEPFU == 0   -1.342      9.749  -0.138    1.000    
## SpeciesLABO == 0   -5.914      9.749  -0.607    0.991    
## SpeciesLACI == 0   -1.212      9.749  -0.124    1.000    
## SpeciesLANO == 0   -2.188      9.749  -0.224    1.000    
## SpeciesMYLU == 0  -72.283      9.749  -7.415   <1e-08 ***
## SpeciesMYSE == 0   -3.560      9.749  -0.365    0.999    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.ISRO, oldNames=FALSE)

## Computing profile confidence intervals ...

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

##                         2.5 %    97.5 %
## sd_(Intercept)|Site   0.00000  20.92750
## sigma                40.77121  52.45007
## SpeciesEPFU         -20.17314  17.48823
## SpeciesLABO         -24.74498  12.91639
## SpeciesLACI         -20.04231  17.61906
## SpeciesLANO         -21.01889  16.64248
## SpeciesMYLU         -91.11323 -53.45186
## SpeciesMYSE         -22.39069  15.27069

Mississippi River

model.MISS <- lmer(Diff ~ 0 + Species + (1|Site), data=data.MISS)
summary(model.MISS)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.MISS
## 
## REML criterion at convergence: 1740.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.7940 -0.1851  0.0102  0.2840  3.2979 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)  245.4   15.67   
##  Residual             2328.2   48.25   
## Number of obs: 168, groups:  Site, 24
## 
## Fixed effects:
##             Estimate Std. Error t value
## SpeciesEPFU  12.5875    10.3555   1.216
## SpeciesLABO  -4.9254    10.3555  -0.476
## SpeciesLACI   4.7345    10.3555   0.457
## SpeciesLANO  -6.3638    10.3555  -0.615
## SpeciesMYLU -44.6007    10.3555  -4.307
## SpeciesMYSE  -0.2728    10.3555  -0.026
## SpeciesPESU  -2.0285    10.3555  -0.196
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU SpMYSE
## SpeciesLABO 0.095                                    
## SpeciesLACI 0.095  0.095                             
## SpeciesLANO 0.095  0.095  0.095                      
## SpeciesMYLU 0.095  0.095  0.095  0.095               
## SpeciesMYSE 0.095  0.095  0.095  0.095  0.095        
## SpeciesPESU 0.095  0.095  0.095  0.095  0.095  0.095

car::Anova(model.MISS)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)   
## Species 22.036  7   0.002504 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hyptest.MISS <- glht(model.MISS, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0", "SpeciesPESU = 0")))
summary(hyptest.MISS)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.MISS)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## SpeciesEPFU == 0  12.5875    10.3555   1.216 0.825845    
## SpeciesLABO == 0  -4.9254    10.3555  -0.476 0.999067    
## SpeciesLACI == 0   4.7345    10.3555   0.457 0.999278    
## SpeciesLANO == 0  -6.3638    10.3555  -0.615 0.995296    
## SpeciesMYLU == 0 -44.6007    10.3555  -4.307 0.000116 ***
## SpeciesMYSE == 0  -0.2728    10.3555  -0.026 1.000000    
## SpeciesPESU == 0  -2.0285    10.3555  -0.196 0.999998    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.MISS, oldNames=FALSE)

## Computing profile confidence intervals ...

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

##                          2.5 %    97.5 %
## sd_(Intercept)|Site   0.000000  26.82173
## sigma                42.264240  53.26436
## SpeciesEPFU          -7.402145  32.57723
## SpeciesLABO         -24.915062  15.06431
## SpeciesLACI         -15.255228  24.72415
## SpeciesLANO         -26.353478  13.62590
## SpeciesMYLU         -64.590395 -24.61102
## SpeciesMYSE         -20.262520  19.71685
## SpeciesPESU         -22.018187  17.96119

Saint Croix

model.SACN <- lmer(Diff ~ 0 + Species + (1|Site), data=data.SACN)
summary(model.SACN)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.SACN
## 
## REML criterion at convergence: 1354.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7251 -0.1015  0.0454  0.2406  5.8159 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)  226.1   15.04   
##  Residual             4145.8   64.39   
## Number of obs: 126, groups:  Site, 18
## 
## Fixed effects:
##              Estimate Std. Error t value
## SpeciesEPFU -15.86394   15.58479  -1.018
## SpeciesLABO -11.86317   15.58479  -0.761
## SpeciesLACI  -3.93667   15.58479  -0.253
## SpeciesLANO   0.83044   15.58479   0.053
## SpeciesMYLU -41.06939   15.58479  -2.635
## SpeciesMYSE  -5.64994   15.58479  -0.363
## SpeciesPESU  -0.04722   15.58479  -0.003
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU SpMYSE
## SpeciesLABO 0.052                                    
## SpeciesLACI 0.052  0.052                             
## SpeciesLANO 0.052  0.052  0.052                      
## SpeciesMYLU 0.052  0.052  0.052  0.052               
## SpeciesMYSE 0.052  0.052  0.052  0.052  0.052        
## SpeciesPESU 0.052  0.052  0.052  0.052  0.052  0.052

car::Anova(model.SACN)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)
## Species 8.2038  7      0.315

hyptest.SACN <- glht(model.SACN, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0", "SpeciesPESU = 0")))
summary(hyptest.SACN)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.SACN)
## 
## Linear Hypotheses:
##                   Estimate Std. Error z value Pr(>|z|)  
## SpeciesEPFU == 0 -15.86394   15.58479  -1.018   0.9236  
## SpeciesLABO == 0 -11.86317   15.58479  -0.761   0.9838  
## SpeciesLACI == 0  -3.93667   15.58479  -0.253   1.0000  
## SpeciesLANO == 0   0.83044   15.58479   0.053   1.0000  
## SpeciesMYLU == 0 -41.06939   15.58479  -2.635   0.0573 .
## SpeciesMYSE == 0  -5.64994   15.58479  -0.363   0.9999  
## SpeciesPESU == 0  -0.04722   15.58479  -0.003   1.0000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.SACN, oldNames=FALSE)

## Computing profile confidence intervals ...

##                         2.5 %    97.5 %
## sd_(Intercept)|Site   0.00000  32.51508
## sigma                55.07339  71.94799
## SpeciesEPFU         -45.78052  14.05263
## SpeciesLABO         -41.77974  18.05341
## SpeciesLACI         -33.85324  25.97991
## SpeciesLANO         -29.08613  30.74702
## SpeciesMYLU         -70.98596 -11.15281
## SpeciesMYSE         -35.56652  24.26663
## SpeciesPESU         -29.96380  29.86935

Sleeping Bear Dunes

model.SLBE <- lmer(Diff ~ 0 + Species + (1|Site), data=data.SLBE)
summary(model.SLBE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.SLBE
## 
## REML criterion at convergence: 2109.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7327 -0.1780 -0.0113  0.2552  6.3065 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)   55.17   7.428  
##  Residual             2176.45  46.652  
## Number of obs: 204, groups:  Site, 34
## 
## Fixed effects:
##              Estimate Std. Error t value
## SpeciesEPFU   5.07206    8.10160   0.626
## SpeciesLABO -18.99468    8.10160  -2.345
## SpeciesLACI   6.86500    8.10160   0.847
## SpeciesLANO   8.03294    8.10160   0.992
## SpeciesMYLU -14.02665    8.10160  -1.731
## SpeciesMYSE  -0.08391    8.10160  -0.010
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU
## SpeciesLABO 0.025                             
## SpeciesLACI 0.025  0.025                      
## SpeciesLANO 0.025  0.025  0.025               
## SpeciesMYLU 0.025  0.025  0.025  0.025        
## SpeciesMYSE 0.025  0.025  0.025  0.025  0.025

car::Anova(model.SLBE)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)  
## Species 10.797  6    0.09486 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hyptest.SLBE <- glht(model.SLBE, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0")))
summary(hyptest.SLBE)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.SLBE)
## 
## Linear Hypotheses:
##                   Estimate Std. Error z value Pr(>|z|)
## SpeciesEPFU == 0   5.07206    8.10160   0.626    0.989
## SpeciesLABO == 0 -18.99468    8.10160  -2.345    0.109
## SpeciesLACI == 0   6.86500    8.10160   0.847    0.952
## SpeciesLANO == 0   8.03294    8.10160   0.992    0.902
## SpeciesMYLU == 0 -14.02665    8.10160  -1.731    0.407
## SpeciesMYSE == 0  -0.08391    8.10160  -0.010    1.000
## (Adjusted p values reported -- single-step method)

confint(model.SLBE, oldNames=FALSE)

## Computing profile confidence intervals ...

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

##                          2.5 %    97.5 %
## sd_(Intercept)|Site   0.000000 18.291135
## sigma                41.477071 51.243103
## SpeciesEPFU         -10.645686 20.789804
## SpeciesLABO         -34.712421 -3.276932
## SpeciesLACI          -8.852745 22.582745
## SpeciesLANO          -7.684804 23.750686
## SpeciesMYLU         -29.744392  1.691098
## SpeciesMYSE         -15.801657 15.633833

Voyageurs

model.VOYA <- lmer(Diff ~ 0 + Species + (1|Site), data=data.VOYA)
summary(model.VOYA)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Diff ~ 0 + Species + (1 | Site)
##    Data: data.VOYA
## 
## REML criterion at convergence: 1215.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.5840 -0.1220  0.0136  0.0766  5.8464 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept)    2.961  1.721  
##  Residual             7364.439 85.816  
## Number of obs: 108, groups:  Site, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## SpeciesEPFU  -0.6048    20.2312  -0.030
## SpeciesLABO  -4.8387    20.2312  -0.239
## SpeciesLACI  14.3252    20.2312   0.708
## SpeciesLANO  -3.2594    20.2312  -0.161
## SpeciesMYLU -66.6106    20.2312  -3.292
## SpeciesMYSE  -8.2444    20.2312  -0.408
## 
## Correlation of Fixed Effects:
##             SpEPFU SpLABO SpLACI SpLANO SpMYLU
## SpeciesLABO 0.000                             
## SpeciesLACI 0.000  0.000                      
## SpeciesLANO 0.000  0.000  0.000               
## SpeciesMYLU 0.000  0.000  0.000  0.000        
## SpeciesMYSE 0.000  0.000  0.000  0.000  0.000

car::Anova(model.VOYA)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Diff
##          Chisq Df Pr(>Chisq)  
## Species 11.592  6    0.07172 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hyptest.VOYA <- glht(model.VOYA, linfct= (Species = c("SpeciesEPFU = 0", "SpeciesLABO = 0", "SpeciesLACI = 0", "SpeciesLANO = 0", "SpeciesMYLU = 0", "SpeciesMYSE = 0")))
summary(hyptest.VOYA)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = Diff ~ 0 + Species + (1 | Site), data = data.VOYA)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)   
## SpeciesEPFU == 0  -0.6048    20.2312  -0.030  1.00000   
## SpeciesLABO == 0  -4.8387    20.2312  -0.239  0.99995   
## SpeciesLACI == 0  14.3252    20.2312   0.708  0.97998   
## SpeciesLANO == 0  -3.2594    20.2312  -0.161  1.00000   
## SpeciesMYLU == 0 -66.6106    20.2312  -3.292  0.00594 **
## SpeciesMYSE == 0  -8.2444    20.2312  -0.408  0.99900   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model.VOYA, oldNames=FALSE)

## Computing profile confidence intervals ...

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

## Warning in optwrap(optimizer, par = start, fn = function(x)
## dd(mkpar(npar1, : convergence code 3 from bobyqa: bobyqa -- a trust region
## step failed to reduce q

##                          2.5 %    97.5 %
## sd_(Intercept)|Site    0.00000  35.85976
## sigma                 72.55307  95.90790
## SpeciesEPFU          -39.48513  38.27557
## SpeciesLABO          -43.71907  34.04163
## SpeciesLACI          -24.55513  53.20557
## SpeciesLANO          -42.13980  35.62091
## SpeciesMYLU         -105.49091 -27.73020
## SpeciesMYSE          -47.12474  30.63596

Bat data models

Katy Goodwin

February 5, 2019

Set up stuff

Data organization

Plot 1 - histogram of counts

Plot 2 - histogram of densities, density estimate, and probability density function

Test significance of “Park” variable

Subset data by parks

Run model and hypothesis tests for each of the 8 parks

Apostle Islands

Grand Portage

Indiana Dunes

Isle Royale

Mississippi River

Saint Croix

Sleeping Bear Dunes

Voyageurs