# Invertebrate data from a simple before/after analysis.
# Two streams were measured at a small hydro electric project. At each stream
# multiple samples were taken on the steam, and the number of invertebrates
# was counted.
# After 5 years, the hydro electric plant started up, and an additional 5 years
# of data were collected.
# Is there evidence of a change in abundance after the plant starts?
# Note that this is NOT a BACI design, and is a VERY poor substitute for such.
# The key problem is that changes may have occured over time unrelated
# to the impact, so a "change" may simple be natural.
# Also, in this analysis we have ignored autocorrelation which can be
# a serious problem in long-term studies.
library(ggfortify)
## Loading required package: ggplot2
library(ggplot2)
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(plyr)
source("../../schwarz.functions.r")
# Read in the data; declare Period as a factor; create factor for Year
# sink('invert-R-001.txt', split=TRUE)
##***part001b;
invert <- read.csv('invert.csv', header=TRUE, as.is=TRUE, strip.white=TRUE)
# Declare Period as an ordered factor so that sorts in proper order
invert$Period <- factor(invert$Period, levels=c("Before","After"), ordered=TRUE)
# We also want Year to be categorical variable so make a new variable
invert$YearF <- factor(invert$Year)
# Stream should also be defined as a factor
invert$StreamF <- factor(invert$Stream)
# It turns out that missing values cause subtle problems in some R
# functions that deal with the fitted object. Sigh... Usually, R
# is sensible enough, but not all package authors did a good job.
cat("Dimensions of full dataset :",dim(invert), "\n")
## Dimensions of full dataset : 100 7
invert <- invert[complete.cases(invert),]
cat("Dimensions of reduced dataset after removing missing values: ", dim(invert), "\n")
## Dimensions of reduced dataset after removing missing values: 90 7
invert[1:12,]
## Year Period Stream quadrat count YearF StreamF
## 1 1 Before 1 1 58 1 1
## 2 1 Before 1 2 62 1 1
## 3 1 Before 1 3 48 1 1
## 4 1 Before 1 4 72 1 1
## 5 1 Before 1 5 41 1 1
## 6 1 Before 2 1 76 1 2
## 7 1 Before 2 2 73 1 2
## 9 1 Before 2 4 69 1 2
## 10 1 Before 2 5 66 1 2
## 11 2 Before 1 1 61 2 1
## 12 2 Before 1 2 48 2 1
## 13 2 Before 1 3 60 2 1
##***part001e;
# sink()
str(invert) # check the structure
## 'data.frame': 90 obs. of 7 variables:
## $ Year : int 1 1 1 1 1 1 1 1 1 2 ...
## $ Period : Ord.factor w/ 2 levels "Before"<"After": 1 1 1 1 1 1 1 1 1 1 ...
## $ Stream : int 1 1 1 1 1 2 2 2 2 1 ...
## $ quadrat: int 1 2 3 4 5 1 2 4 5 1 ...
## $ count : int 58 62 48 72 41 76 73 69 66 61 ...
## $ YearF : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 2 ...
## $ StreamF: Factor w/ 2 levels "1","2": 1 1 1 1 1 2 2 2 2 1 ...
# look at the std dev to mean plot
# Compute the averages etc
# sink('invert-R-020.txt', split=TRUE)
##***part020b;
mean.invert <- plyr::ddply(invert, c("Year","YearF","Stream","StreamF","Period"), function(x){
# compute sample size, missing values, mean, and sd
n <- nrow(x)
n.miss <- sum(is.na(x$count))
count.mean <- mean(x$count, na.rm=TRUE)
count.sd <- sd (x$count, na.rm=TRUE)
res <- c(n=n, n.miss=n.miss, count.mean=count.mean, count.sd=count.sd) # put results together
res
})
mean.invert
## Year YearF Stream StreamF Period n n.miss count.mean count.sd
## 1 1 1 1 1 Before 5 0 56.20 12.091319
## 2 1 1 2 2 Before 4 0 71.00 4.396969
## 3 2 2 1 1 Before 5 0 56.20 5.167204
## 4 2 2 2 2 Before 5 0 60.00 6.204837
## 5 3 3 1 1 Before 4 0 49.75 7.889867
## 6 3 3 2 2 Before 5 0 56.60 7.733046
## 7 4 4 1 1 Before 3 0 50.00 8.888194
## 8 4 4 2 2 Before 5 0 59.20 12.111978
## 9 5 5 1 1 Before 4 0 41.25 3.774917
## 10 5 5 2 2 Before 4 0 52.25 12.867919
## 11 6 6 1 1 After 5 0 68.80 9.782638
## 12 6 6 2 2 After 5 0 70.60 8.734987
## 13 7 7 1 1 After 5 0 52.00 7.314369
## 14 7 7 2 2 After 5 0 61.20 5.718391
## 15 8 8 1 1 After 5 0 58.80 6.760178
## 16 8 8 2 2 After 3 0 73.00 14.000000
## 17 9 9 1 1 After 5 0 67.80 4.919350
## 18 9 9 2 2 After 5 0 80.60 8.203658
## 19 10 10 1 1 After 4 0 56.25 6.396614
## 20 10 10 2 2 After 4 0 64.25 7.182154
##***part020e;
# sink()
##***part060b;
# Plot the log(std) vs. log(mean)
# Because there is no relationship between the mean and the sd,
# this plot shows that indicates that no transformation is needed.
sd.plot <- ggplot(data=mean.invert, aes(x=log(count.mean), y=log(count.sd),
shape=StreamF, color=StreamF))+
ggtitle("SD vs. Mean")+
geom_point(size=4)
sd.plot
##***part060e;
# ggsave(plot=sd.plot, file='invert-R-060-varmeanplot.png', height=4, width=6, units="in", dpi=300)
# Plot the mean count over time by stream
##***part070b;
plot.mean <- ggplot(data=invert, aes(x=Year, y=count,
shape=StreamF, group=StreamF, color=StreamF))+
ggtitle("Preliminary plot of data and means")+
geom_point(size=2, position=position_dodge(w=0.2))+
geom_vline(xintercept=.5+max(invert$Year[as.character(invert$Period)=="Before"],na.rm=TRUE))+
stat_summary(fun="mean", geom="line") # , position=position_dodge(w=0.2))
plot.mean
##***part070e;
# ggsave(plot=plot.mean, file='invert-R-070-trendplot.png', height=4, width=6, units='in', dpi=300)
############################################################################################
# Analysis of Stream 1 alone to illustrate how to analyze one stream's data
##***part99b;
# Subset the data
stream1 <- invert[ invert$Stream ==1,]
head(stream1)
## Year Period Stream quadrat count YearF StreamF
## 1 1 Before 1 1 58 1 1
## 2 1 Before 1 2 62 1 1
## 3 1 Before 1 3 48 1 1
## 4 1 Before 1 4 72 1 1
## 5 1 Before 1 5 41 1 1
## 11 2 Before 1 1 61 2 1
# The analysis of the mean(Count)) for each year
mean.invert1 <- mean.invert[ mean.invert$Stream ==1, ]
mean.invert1
## Year YearF Stream StreamF Period n n.miss count.mean count.sd
## 1 1 1 1 1 Before 5 0 56.20 12.091319
## 3 2 2 1 1 Before 5 0 56.20 5.167204
## 5 3 3 1 1 Before 4 0 49.75 7.889867
## 7 4 4 1 1 Before 3 0 50.00 8.888194
## 9 5 5 1 1 Before 4 0 41.25 3.774917
## 11 6 6 1 1 After 5 0 68.80 9.782638
## 13 7 7 1 1 After 5 0 52.00 7.314369
## 15 8 8 1 1 After 5 0 58.80 6.760178
## 17 9 9 1 1 After 5 0 67.80 4.919350
## 19 10 10 1 1 After 4 0 56.25 6.396614
##***part99e;
# Test for a period effect using the mean count data
# sink('invert-R-100.txt', split=TRUE)
##***part100b;
result100 <- lm( count.mean ~ Period, data=mean.invert1)
anova(result100)
## Analysis of Variance Table
##
## Response: count.mean
## Df Sum Sq Mean Sq F value Pr(>F)
## Period 1 252.51 252.506 5.5146 0.0468 *
## Residuals 8 366.31 45.789
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##***part100e;
# sink()
# Estimate the period effect
# Do NOT use the summary output because it gives odd results
# with ordered factors. Compare the two following models summary output.
# They should give identical estimates but are not even though the p-value is the same
summary(result100) #with ordered factor
##
## Call:
## lm(formula = count.mean ~ Period, data = mean.invert1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.430 -3.842 -0.805 5.520 8.070
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 55.705 2.140 26.032 5.09e-09 ***
## Period.L 7.106 3.026 2.348 0.0468 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.767 on 8 degrees of freedom
## Multiple R-squared: 0.408, Adjusted R-squared: 0.3341
## F-statistic: 5.515 on 1 and 8 DF, p-value: 0.0468
summary(lm(count.mean ~ factor(as.character(Period)), data=mean.invert1))
##
## Call:
## lm(formula = count.mean ~ factor(as.character(Period)), data = mean.invert1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.430 -3.842 -0.805 5.520 8.070
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.730 3.026 20.068 3.97e-08 ***
## factor(as.character(Period))Before -10.050 4.280 -2.348 0.0468 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.767 on 8 degrees of freedom
## Multiple R-squared: 0.408, Adjusted R-squared: 0.3341
## F-statistic: 5.515 on 1 and 8 DF, p-value: 0.0468
# It is better to formally estimate the difference using the
# emmeans package and use the pairs to estimate the difference
# sink('invert-R-105.txt', split=TRUE)
##***part105b;
# Estimate the marginal means and the differences.
result100.emmo <- emmeans::emmeans(result100, ~Period)
cat("\nEstimated marginal means for each Period \n\n")
##
## Estimated marginal means for each Period
summary(result100.emmo)
## Period emmean SE df lower.CL upper.CL
## Before 50.7 3.03 8 43.7 57.7
## After 60.7 3.03 8 53.8 67.7
##
## Confidence level used: 0.95
cat("\n\nEstimated difference between means in Before and After periods\n")
##
##
## Estimated difference between means in Before and After periods
summary(pairs(result100.emmo), infer=TRUE)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Before - After -10.1 4.28 8 -19.9 -0.181 -2.348 0.0468
##
## Confidence level used: 0.95
##***part105e;
# sink()
# Residual and other plots
##***part106b;
plot.diag100 <- ggplot2::autoplot(result100)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
## Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot.diag100
##***part106e;
# ggsave.ggmultiplot(plot.diag100, # see https://github.com/sinhrks/ggfortify/issues/98 for bug in autoplot
# file='invert-R-106.png', height=4, width=6, units="in", dpi=300)
#################################################################################
# Analysis of the individual quadrat data for Stream 1.
#
# sink('invert-R-200-type3.txt', split=TRUE)
##***part200b;
# Don't forget that you need a Year as a factor.
# We created YearF as a factor when we read in the data
# We don't need to transform counts because they are largish.
result200 <- lmerTest::lmer(count ~ Period + (1|YearF), data=stream1)
anova(result200, ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 309.89 309.89 1 7.9842 5.2373 0.05145 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##***part200e;
# sink()
# Again, don't trust the estimates from the summary table.
summary(result200)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: count ~ Period + (1 | YearF)
## Data: stream1
##
## REML criterion at convergence: 314.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79514 -0.81619 0.02489 0.67490 2.23492
##
## Random effects:
## Groups Name Variance Std.Dev.
## YearF (Intercept) 33.80 5.814
## Residual 59.17 7.692
## Number of obs: 45, groups: YearF, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 55.810 2.173 8.210 25.68 3.91e-09 ***
## Period.L 7.036 3.073 8.210 2.29 0.0505 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Period.L -0.022
# sink('invert-R-200-vc.txt', split=TRUE)
##***part200-vcb;
# Extract the variance components.
vc <- VarCorr(result200)
vc
## Groups Name Std.Dev.
## YearF (Intercept) 5.8137
## Residual 7.6922
##***part200-vce;
# sink()
# emmeans after a lmer() fit
# sink('invert-R-200LSM-Period.txt', split=TRUE)
##***part200LSM-Periodb;
# Estimate the marginal means.
result200.emmo <- emmeans::emmeans(result200, ~Period)
cat("\n\nEstimated marginal means\n")
##
##
## Estimated marginal means
summary(result200.emmo)
## Period emmean SE df lower.CL upper.CL
## Before 50.8 3.11 8.31 43.7 58.0
## After 60.8 3.04 7.66 53.7 67.8
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part200LSM-Periode;
# sink()
# sink("invert-R-200baci.txt", split=TRUE)
##***part200bacib;
# Estimate the period effect along with a se
cat("\n\nEstimated difference between means in Before and After periods\n")
##
##
## Estimated difference between means in Before and After periods
summary(pairs(result200.emmo), infer=TRUE)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Before - After -9.95 4.35 7.98 -20 0.0795 -2.289 0.0514
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part200bacie;
# sink()
# Diagnostic plots for mixed models are not as standardized as for
# ordinary linear models so there are several than can be produced.
##***part200diagnosticb;
plot.diag200 <- sf.autoplot.lmer(result200)
## Loading required package: broom.mixed
## Loading required package: grid
## Loading required package: gridExtra
## Loading required package: lattice
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
plot(plot.diag200)
##***part200diagnostice;
# ggsave.ggmultiplot(plot.diag200, # see https://github.com/sinhrks/ggfortify/issues/98 for bug in autoplot
# file='invert-R-200-diagnostic.png', height=4, width=6, units="in", dpi=300)
#--------------------------------------------------- Analyze all streams together --------------
# The analysis of the means values
# sink('invert-R-300-type3.txt', split=TRUE)
##***part300b;
# Model for all streams together
# Note that Stream, Period, and Year all need to be factors
# These were defined as such when the data were read in.
result300 <- lmerTest::lmer(count.mean~ Period + (1|Stream) + (1|YearF), data=mean.invert)
anova(result300, ddf='Kenward-Roger')
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 52.033 52.033 1 8 5.7309 0.04359 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##***part300e;
# sink()
# Again be careful with the summary table because Period is an ordered factor
summary(result300)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: count.mean ~ Period + (1 | Stream) + (1 | YearF)
## Data: mean.invert
##
## REML criterion at convergence: 118.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0409 -0.4726 -0.1592 0.4327 1.3376
##
## Random effects:
## Groups Name Variance Std.Dev.
## YearF (Intercept) 39.828 6.311
## Stream (Intercept) 41.089 6.410
## Residual 9.079 3.013
## Number of obs: 20, groups: YearF, 10; Stream, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60.287 4.998 1.407 12.062 0.0217 *
## Period.L 7.131 2.979 8.000 2.394 0.0436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Period.L 0.000
# sink('invert-R-300-vc.txt', split=TRUE)
##***part300-vcb;
# Extract the variance components
vc <-VarCorr(result300)
vc
## Groups Name Std.Dev.
## YearF (Intercept) 6.3110
## Stream (Intercept) 6.4101
## Residual 3.0132
##***part300-vce;
# sink()
# Estimate the marginal means
# sink('invert-R-300LSM-Period.txt', split=TRUE)
##***part300LSM-Periodb;
# Estimate the marginal means
result300.emmo <- emmeans::emmeans(result300, ~Period)
cat("\n\n Estimated marginal means for each Period \n\n")
##
##
## Estimated marginal means for each Period
summary(result300.emmo)
## Period emmean SE df lower.CL upper.CL
## Before 55.2 5.42 1.92 30.9 79.5
## After 65.3 5.42 1.92 41.0 89.6
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part300LSM-Periode;
# sink()
# sink("invert-R-300baci.txt", split=TRUE)
##***part300bacib;
# Estimate the Period effect
cat("\n\nEstimated difference between means in Before and After periods\n")
##
##
## Estimated difference between means in Before and After periods
pairs(result300.emmo)
## contrast estimate SE df t.ratio p.value
## Before - After -10.1 4.21 8 -2.394 0.0436
##
## Degrees-of-freedom method: kenward-roger
summary(pairs(result300.emmo), infer=TRUE)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Before - After -10.1 4.21 8 -19.8 -0.37 -2.394 0.0436
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part300bacie;
# # sink()
# Check the lmer residuals
##***part300diagnosticb;
plot.diag300 <- sf.autoplot.lmer(result300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
plot(plot.diag300)
##***part300diagnostice;
# ggsave.ggmultiplot(plot.diag300, # see https://github.com/sinhrks/ggfortify/issues/98 for bug in autoplot
# file='invert-R-300-diagnostic.png', height=4, width=6, units="in", dpi=300)
#######################################################################################
# Now for the analysis of the entire dataset (whew)
# # sink('invert-R-400-type3.txt', split=TRUE)
##***part400b;
# Year, Period, Stream must all be defined as factors.
# This was done when the data were read in.
result400 <- lmerTest::lmer(count ~ Period + (1|StreamF)+(1|YearF)+(1|StreamF:YearF), data=invert)
## boundary (singular) fit: see help('isSingular')
anova(result400, ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Period 376.31 376.31 1 7.9986 5.6712 0.04445 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##***part400e;
# sink()
# Again be very careful about using the results
# from the summary() because Period is defined as an ordered factor
summary(result400)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: count ~ Period + (1 | StreamF) + (1 | YearF) + (1 | StreamF:YearF)
## Data: invert
##
## REML criterion at convergence: 644.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0397 -0.7550 0.1514 0.7071 2.2441
##
## Random effects:
## Groups Name Variance Std.Dev.
## StreamF:YearF (Intercept) 0.00 0.000
## YearF (Intercept) 36.65 6.054
## StreamF (Intercept) 38.31 6.190
## Residual 66.35 8.146
## Number of obs: 90, groups: StreamF:YearF, 20; YearF, 10; StreamF, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60.223 4.854 1.393 12.406 0.0215 *
## Period.L 7.077 2.970 8.081 2.383 0.0440 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Period.L -0.001
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# sink('invert-R-400-vc.txt', split=TRUE)
##***part400-vcb;
# Extract the variance components
vc <- VarCorr(result400)
vc
## Groups Name Std.Dev.
## StreamF:YearF (Intercept) 0.0000
## YearF (Intercept) 6.0537
## StreamF (Intercept) 6.1897
## Residual 8.1459
##***part400-vce;
# sink()
# sink('invert-R-400LSM-Period.txt', split=TRUE) # Note these have the wrong standard error
##***part400LSM-Periodb;
# Estimate the marginal means
result400.emmo <- emmeans::emmeans(result400, ~Period)
cat("\n\n Estimated marginal means for each Period \n\n")
##
##
## Estimated marginal means for each Period
summary(result400.emmo)
## Period emmean SE df lower.CL upper.CL
## Before 55.2 5.29 1.93 31.7 78.8
## After 65.2 5.29 1.92 41.6 88.9
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part400LSM-Periode;
# sink()
# sink("invert-R-400baci.txt", split=TRUE) # note these have the wrong df
##***part400bacib;
# Estimate the Period effect
cat("\n\nEstimated difference between means in Before and After periods\n")
##
##
## Estimated difference between means in Before and After periods
summary(pairs(result400.emmo), infer=TRUE)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Before - After -10 4.2 8 -19.7 -0.317 -2.381 0.0445
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part400bacie;
# sink()
# diagnostic plots for the lmer fit
##***part400diagnosticb;
plot.diag400 <- sf.autoplot.lmer(result400)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
plot(plot.diag400)
##***part400diagnostice;
# ggsave.ggmultiplot(plot.diag400, # see https://github.com/sinhrks/ggfortify/issues/98 for bug in autoplot
# file='invert-R-400-diagnostic.png', height=4, width=6, units="in", dpi=300)