<!-- ################################### -->

<!-- #### MAR 536 Biol. Stats II -->

<!-- # Spring 2023 -->

<!-- # Lecture 13 -->

<!-- # 14 March 2023 -->

<!-- # Instructor: Gavin Fay, gfay@umassd.edu -->

<!-- ################################### -->

```{r}
library(tidyverse)
library(lattice)
```

### STREAMS EXAMPLE

```{r}
#library(nlme)
library(lme4)
Streams <- scan('../data/streams.dat',what=list(Stream=0,Density=0),n=3*18,skip=5)
streams.lm1 <- lm(Density ~ 1, data = Streams)
streams.lm2 <- lm(Density ~ factor(Stream) - 1, data = Streams)
streams.lme1 <- lmer(Density ~ 1 + (1 | Stream), data = Streams, REML=FALSE)
```

```{r}
#help function for lme()
?(lmer)
```

```{r}
#summary of results
summary(streams.lme1)
```

```{r}
#compare residuals among models
par(las=0)
boxplot(split(residuals(streams.lm1)/summary(streams.lm1)$sigma,Streams$Stream))
mtext(side=3,"lm, intercept only",line=2)
mtext(side=2,"Residuals",line=2)
abline(h=0)
boxplot(split(residuals(streams.lm2)/summary(streams.lm2)$sigma,Streams$Stream))
mtext(side=3,"lm, stream effect",line=2)
mtext(side=1,"Stream",line=2)
abline(h=0)
boxplot(split(residuals(streams.lme1)/summary(streams.lme1)$sigma,Streams$Stream))
mtext(side=3,"lmer, random stream effect",line=2)
abline(h=0)
```

```{r}
# plot the fitted values
# library(lattice)
# plot(streams.lme1,Density~fitted(.)|Stream, abline = c(0,1))
library(broom.mixed)
lme1_results <- streams.lme1 |>
   augment() |>
  janitor::clean_names()
lme1_results |>
  ggplot() +
  aes(x = fitted, y = density) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_wrap(~stream)
```

```{r}
# alternative using lme
library(nlme)
streams.lme1 <- lme(Density~1,random=~1|Stream,data=Streams,meth="ML")
# diagnostic plots of residuals, random effects
par(mfrow=c(2,3))
# QQ plot
qqnorm(streams.lme1$residuals/streams.lme1$sigma,ylab="Quantiles of Residuals")
qqline(streams.lme1$residuals/streams.lme1$sigma)

plot(residuals(streams.lme1)/streams.lme1$sigma,ylab="Standardized Residuals")
hist(residuals(streams.lme1)/streams.lme1$sigma,xlab="Standardized Residuals",main="")

#homogeneity of within group variance
boxplot(split(residuals(streams.lme1)/streams.lme1$sigma,Streams$Stream),ylab="Standardized Residual",xlab="Stream",csi=0.2)
abline(0,0,lwd=3)
#normality of the between-group residuals
re.sigma <- as.numeric(VarCorr(streams.lme1)[1,2])
print(streams.lme1$coefficients$random$Stream)
re<-streams.lme1$coefficients$random$Stream/re.sigma  ##/0.6184
qqnorm(re,ylab="Quantiles of random effects")
qqline(re)
hist(streams.lme1$coefficients$random$Stream/re.sigma,xlab="random effects",main="")
```

```{r}
### RIKZ Species Richness
RIKZ <- readRDS("../data/RIKZ.rds")

# no random effects
m1 <- lm(Richness ~ NAP, data = RIKZ)
# random intercept
m5 <- lmer(Richness ~ NAP + (1|Beach), data = RIKZ, REML = FALSE)
# random intercepts and slopes
m6 <- lmer(Richness ~ NAP + (1 + NAP|Beach), data = RIKZ, REML = FALSE)

AIC(m1,m5,m6)

xyplot(Richness~NAP|factor(Beach),data=RIKZ)
plot(m6,Richness~fitted(.)|factor(Beach),abline = c(0,1))

xyplot(resid(m1,type="pearson")~fitted(m1)|factor(RIKZ$Beach),abline=0,pch=16,ylab="stdized residuals")

library(nlme)
m5<-lme(Richness~NAP ,random= ~ 1 | factor(Beach), data = RIKZ, method="REML")
m6<-lme(Richness~NAP ,random= ~ NAP | factor(Beach), data = RIKZ, method="REML")

plot(m6,form=resid(.,type="p")~fitted(.)|factor(Beach),abline=0,pch=16)

plot(m5,form=resid(.,type="p")~fitted(.)|factor(Beach),abline=0,pch=16)


#      Output for model 6 with random slope and intercept, page 129
summary(m6)
anova(m6)
#I am not sure why these results are slightly different from those in the book
#It could be due to different R versions (1.6 in the book, and 2.7 now)
plot(m6)

#      Model selection (comparing models 5 and 6), page 130
anova(m5, m6)
```

```{r}
###### Weight-Length relationships ######
library(lme4)
WtLen <- scan('../data/wtlen.txt',what=list(Subject=0,Length=0,Weight=0),n=3*100,skip=1)
WtLen$LogLen <- log(WtLen$Length)
WtLen$LogWt <- log(WtLen$Weight)
# linear model, no individual parameters
wtlen.lm1 <- lm(LogWt~LogLen,data=WtLen)
# linear model, individual ln_a parameters
wtlen.lm2 <- lm(LogWt~factor(Subject)+LogLen-1,data=WtLen)
## linear mixed effects models
# mixed effects model, ln_a as a random effect by fish
wtlen.lme1 <- lmer(LogWt~LogLen + (1|Subject) ,data=WtLen, REML = FALSE)
# mixed effects model, b as a random effect by fish
wtlen.lme2 <- lmer(LogWt~LogLen + (-1+LogLen|Subject),data=WtLen,REML=FALSE)
# mixed effects model, random intercepts and slopes
#lmc <- lmeControl(niterEM = 5000, msMaxIter = 1000)
wtlen.lme3 <- lmer(LogWt~LogLen + (1+LogLen|Subject), 
                  data=WtLen,REML = FALSE)
```

```{r}
library(lattice)
xyplot(LogWt~LogLen|factor(Subject),data=WtLen)
xyplot(Weight~Length|factor(Subject),data=WtLen)


```