Lab 7
Kirk J Schmitz
11/11/2019
##Homework: Lab 7 This lab uses the dataset “RIKZ” to look at a number of samples
###Loading
library(readr)
library(wesanderson)
library(lme4)## Loading required package: Matrixlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepIt is important to load the appropriate packages before running code. The package wesanderson carries a number of color palletes inspired by the films of Wes Anderson.
###Creating Objects
#Creating an object for the "RIKZ" file
dat <- read_csv("ANT291WD/RIKZ data.csv")## Parsed with column specification:
## cols(
## Sample = col_double(),
## Richness = col_double(),
## Exposure = col_double(),
## NAP = col_double(),
## Beach = col_double()
## )dat$Beach <- as.factor(dat$Beach)
#Creating a color pallete for
pal <- wes_palette("Zissou1",9, type = "continuous")First we will look at one categorical variable, the “Beach” that the sample was collected from, and two continuous variables, the number of species observed
plot(as.factor(dat$Beach),dat$Richness, xlab="Beach", ylab="Richness", las=1, col=pal, main="Species Richness by Beach", type="l")
Species richness looks to be variable by beach. Beaches one and two have an average species richness greater than the other beaches, but beach five has the highest single observation of richness.
The next figure looks at the interaction of richness and the mean tidal level (NAP).
plot(dat$Richness, dat$NAP, col=dat$Beach, pch=16, xlab="Richness", ylab="NAP", las=1)
There appears to be a negative corrlation between NAP and species richness, with fewer species above tidal level and greater species richness the farther the sample is taken below mean tidal level.
###Constructing the model
I used the model from class on Wednesday to format my model today. I am looking at the Richness variable and how it relates to NAP. I changed the family from “binomial” to “Gaussian”, beceause number of species is a continuous variable (as opposed to the binomial remission/nonremission response variable we were looking at in class). I used the beach site as the random component.
model1 adds the variable exposure into the model. I get an error, and I am not sure what the problem is. I added the variable in, just as we added them in lab last wednesday. Possibly my assumption that the Richness variable was continuous is incorrect and that it is a categorical variable, thus the gaussian model may be innapproprite. Changing the model to a poisson model allows me to create the object.
model <- glmer(Richness~NAP+(1|Beach), data=dat,family="poisson", control=glmerControl(optimizer = "bobyqa"))
summary(model)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Richness ~ NAP + (1 | Beach)
## Data: dat
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 220.8 226.2 -107.4 214.8 42
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9648 -0.6155 -0.2243 0.2236 3.1869
##
## Random effects:
## Groups Name Variance Std.Dev.
## Beach (Intercept) 0.2249 0.4743
## Number of obs: 45, groups: Beach, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.66234 0.17373 9.569 < 2e-16 ***
## NAP -0.50389 0.07535 -6.687 2.28e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## NAP 0.013plot(model)
model1 <- glmer(Richness~NAP+Exposure+(1|Beach), data=dat, family ="poisson", control=glmerControl(optimizer = "bobyqa"))
summary(model1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Richness ~ NAP + Exposure + (1 | Beach)
## Data: dat
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 211.0 218.2 -101.5 203.0 41
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8186 -0.5430 -0.2574 0.1204 3.9435
##
## Random effects:
## Groups Name Variance Std.Dev.
## Beach (Intercept) 0.03485 0.1867
## Number of obs: 45, groups: Beach, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.25227 0.93145 6.712 1.91e-11 ***
## NAP -0.50570 0.07252 -6.974 3.09e-12 ***
## Exposure -0.44831 0.09279 -4.832 1.35e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) NAP
## NAP 0.071
## Exposure -0.995 -0.066plot(model1)
model2 is the same as model1, but examines whether NAP and Exposure interact with each other
model2 <- glmer(Richness~NAP*Exposure+(1|Beach), data=dat, family ="poisson", control=glmerControl(optimizer = "bobyqa"))
summary(model2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Richness ~ NAP * Exposure + (1 | Beach)
## Data: dat
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 210.2 219.2 -100.1 200.2 40
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7492 -0.5404 -0.2704 0.1349 3.8268
##
## Random effects:
## Groups Name Variance Std.Dev.
## Beach (Intercept) 0.03284 0.1812
## Number of obs: 45, groups: Beach, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.31778 0.91332 6.917 4.60e-12 ***
## NAP 0.50025 0.61181 0.818 0.4136
## Exposure -0.45498 0.09108 -4.996 5.86e-07 ***
## NAP:Exposure -0.10533 0.06365 -1.655 0.0979 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) NAP Exposr
## NAP 0.033
## Exposure -0.995 -0.037
## NAP:Exposur -0.034 -0.993 0.039plot(model2)