# BACI design with multiple controls; 2 factor; interaction;
# A BACI design was used to assess the impact
# of cooling water discharge on the density of
# shore crabs.
# The beach near the outlet of the cooling water
# was sampled using several quadrats
# before and after the plant started operation.
# Two control beaches at other locations
# in the inlet were also sampled.
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)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
#source("http://www.stat.sfu.ca/~cschwarz/Stat-Ecology-Datasets/schwarz.functions.r")
source("../../schwarz.functions.r")
source("../baci-power.r")
# Read in the actual data
# sink("baci-crabs-mod-R-001.txt", split=TRUE)
##***part001b;
crabs <- read.csv("baci-crabs-mod.csv", header=TRUE, as.is=TRUE, strip.white=TRUE)
crabs$trt <- interaction(crabs$SiteClass, crabs$Site, crabs$Period)
crabs$Site <- factor(crabs$Site)
crabs$SiteClass <- factor(crabs$SiteClass)
crabs$Period <- factor(crabs$Period, levels=c("Before","After"), ordered=TRUE)
crabs$trt <- factor(crabs$trt)
crabs$logDensity<- log(crabs$Density)
cat("Listing of part of the raw data \n")
## Listing of part of the raw data
crabs[1:10,]
## SiteClass Site Period Quadrat Density trt logDensity
## 1 Impact I1 Before Q1 23 Impact.I1.Before 3.135494
## 2 Impact I1 Before Q2 30 Impact.I1.Before 3.401197
## 3 Impact I1 Before Q3 27 Impact.I1.Before 3.295837
## 4 Impact I1 Before Q4 25 Impact.I1.Before 3.218876
## 5 Impact I1 After Q5 17 Impact.I1.After 2.833213
## 6 Impact I1 After Q6 19 Impact.I1.After 2.944439
## 7 Impact I1 After Q7 23 Impact.I1.After 3.135494
## 8 Impact I1 After Q8 25 Impact.I1.After 3.218876
## 9 Impact I1 After Q9 17 Impact.I1.After 2.833213
## 10 Control C1 Before Q10 32 Control.C1.Before 3.465736
##***part001e;
# sink()
str(crabs)
## 'data.frame': 30 obs. of 7 variables:
## $ SiteClass : Factor w/ 2 levels "Control","Impact": 2 2 2 2 2 2 2 2 2 1 ...
## $ Site : Factor w/ 3 levels "C1","C2","I1": 3 3 3 3 3 3 3 3 3 1 ...
## $ Period : Ord.factor w/ 2 levels "Before"<"After": 1 1 1 1 2 2 2 2 2 1 ...
## $ Quadrat : chr "Q1" "Q2" "Q3" "Q4" ...
## $ Density : int 23 30 27 25 17 19 23 25 17 32 ...
## $ trt : Factor w/ 6 levels "Control.C1.After",..: 6 6 6 6 3 3 3 3 3 4 ...
## $ logDensity: num 3.14 3.4 3.3 3.22 2.83 ...
# Preliminary plot
# Get side-by-side dot plots
##***part010b;
prelimplot <- ggplot(data=crabs, aes(x=trt, y=logDensity))+
ggtitle("Preliminary plot to look for outliers etc")+
geom_point(position=position_jitter(w=0.1))+
geom_boxplot(alpha=0.1)
prelimplot
##***part010e;
#ggsave(plot=prelimplot, file="baci-crabs-mod-R-010.png",
# h=4, w=6, units="in", dpi=300)
# Get some simple summary statistics
# sink('baci-crabs-mod-R-020.txt', split=TRUE)
##***part020b;
report <- plyr::ddply(crabs, c("Period","SiteClass","Site"), function(x){
res <- sf.simple.summary(x, "logDensity", crd=TRUE)
return(res)
})
report
## Period SiteClass Site n nmiss mean sd se lcl
## 1 Before Control C1 6 0 3.412531 0.11191792 0.04569030 3.295081
## 2 Before Control C2 6 0 3.634335 0.08851389 0.03613565 3.541445
## 3 Before Impact I1 4 0 3.262851 0.11310962 0.05655481 3.082868
## 4 After Control C1 4 0 3.309636 0.10870652 0.05435326 3.136660
## 5 After Control C2 5 0 3.436145 0.10386916 0.04645170 3.307174
## 6 After Impact I1 5 0 2.993047 0.17659714 0.07897664 2.773773
## ucl
## 1 3.529982
## 2 3.727224
## 3 3.442834
## 4 3.482613
## 5 3.565115
## 6 3.212321
##***part020e;
# sink()
# Draw a profile plot
##***part030b;
profileplot <- ggplot(data=report, aes(x=Period, y=mean,
group=Site, color=SiteClass, shape=Site))+
ggtitle("Profile plot of crab density")+
ylab("Density with mean and 95% ci")+
geom_point(position=position_dodge(w=0.1))+
geom_line(position=position_dodge(w=0.1))+
geom_errorbar(aes(ymax=ucl, ymin=lcl), width=0.1, position=position_dodge(w=0.1))+
geom_point(data=crabs, aes(y=logDensity), position=position_dodge(w=0.2))
profileplot
##***part030e;
#ggsave(plot=profileplot, file="baci-crabs-mod-R-030.png",
# h=4, w=6, units="in", dpi=300)
# There are several ways in which this BACI design can be analyzed.
#################################################################################
#################################################################################
# Do a t-test on the differences of the averages for each site
# sink('baci-crabs-mod-R-100.txt', split=TRUE)
##***part100b;
# The averages are available in the report dataframe
crabs.avg <- report[,c("Site","SiteClass","Period","mean")]
crabs.avg
## Site SiteClass Period mean
## 1 C1 Control Before 3.412531
## 2 C2 Control Before 3.634335
## 3 I1 Impact Before 3.262851
## 4 C1 Control After 3.309636
## 5 C2 Control After 3.436145
## 6 I1 Impact After 2.993047
##***part100e;
# )
# sink('baci-crabs-mod-R-101.txt', split=TRUE)
##***part101b;
# Reshape the file to get the Before and After measurements on the same line
#
crabs.site.wide <- tidyr::pivot_wider(crabs.avg,
id_cols=c("Site","SiteClass"),
names_from="Period",
values_from="mean"
)
crabs.site.wide$diff <- crabs.site.wide$After - crabs.site.wide$Before
crabs.site.wide
## # A tibble: 3 × 5
## Site SiteClass Before After diff
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 C1 Control 3.41 3.31 -0.103
## 2 C2 Control 3.63 3.44 -0.198
## 3 I1 Impact 3.26 2.99 -0.270
##***part101e;
# sink()
# do the two sample t-test not assuming equal variances
# Unfortunately, you need at least 2 sites in EACH SiteClass, so this does not work here
# sink('baci-crabs-mod-R-104.txt', split=TRUE)
##***part104b;
result <- try(t.test(diff ~ SiteClass, data=crabs.site.wide),silent=TRUE)
if(class(result)=="try-error")
{cat("Unable to do unequal variance t-test because of small sample size\n")} else
{
result$diff.in.means <- sum(result$estimate*c(1,-1))
names(result$diff.in.means)<- "diff.in.means"
print(result)
cat("Estimated difference in the means on log scale: ",result$diff.in.means,
"( SE ", result$stderr, ") \n")
}
## Unable to do unequal variance t-test because of small sample size
##***part104e;
# sink()
# do the two sample t-test assuming equal variances
# This only requires at least 2 sites in at least one of the SiteClasses
# sink('baci-crabs-mod-R-105.txt', split=TRUE)
##***part105b;
result <- t.test(diff ~ SiteClass, data=crabs.site.wide, var.equal=TRUE)
result$diff.in.means <- sum(result$estimate*c(1,-1))
names(result$diff.in.means)<- "diff.in.means"
result
##
## Two Sample t-test
##
## data: diff by SiteClass
## t = 1.4451, df = 1, p-value = 0.3854
## alternative hypothesis: true difference in means between group Control and group Impact is not equal to 0
## 95 percent confidence interval:
## -0.9293537 1.1678764
## sample estimates:
## mean in group Control mean in group Impact
## -0.1505426 -0.2698039
cat("Estimated difference in the means on log scale: ",result$diff.in.means,
"( SE ", result$stderr, ") \n")
## Estimated difference in the means on log scale: 0.1192614 ( SE 0.0825278 )
##***part105e;
# sink()
#################################################################################
#################################################################################
# Do a Mixed effect linear model on the averages for each site
crabs.avg <- report[,c("Site","SiteClass","Period","mean")]
crabs.avg
## Site SiteClass Period mean
## 1 C1 Control Before 3.412531
## 2 C2 Control Before 3.634335
## 3 I1 Impact Before 3.262851
## 4 C1 Control After 3.309636
## 5 C2 Control After 3.436145
## 6 I1 Impact After 2.993047
# sink('baci-crabs-mod-R-200-type3.txt', split=TRUE)
##***part200b;
result.avg.lmer <- lmerTest::lmer(mean ~ SiteClass+Period+SiteClass:Period +
(1|Site), data=crabs.avg)
##***part200e;
##***part201b;
anova(result.avg.lmer, ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## SiteClass 0.010233 0.010233 1 1 4.5075 0.2802
## Period 0.058897 0.058897 1 1 25.9427 0.1234
## SiteClass:Period 0.004741 0.004741 1 1 2.0883 0.3854
##***part201e;
# sink()
# sink('baci-crabs-mod-R-200-vc.txt', split=TRUE)
##***part200vcb;
# Extract the variance components
vc.avg <- VarCorr(result.avg.lmer)
vc.avg
## Groups Name Std.Dev.
## Site (Intercept) 0.118448
## Residual 0.047647
##***part200vce;
# sink()
# Extract the BACI effect
# You could get this from the summary table looking at the
# interaction term, but it is more robust to get this from the
# emmeans
summary(result.avg.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ SiteClass + Period + SiteClass:Period + (1 | Site)
## Data: crabs.avg
##
## REML criterion at convergence: -1.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6368 -0.2724 0.0000 0.2724 0.6368
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.01403 0.11845
## Residual 0.00227 0.04765
## Number of obs: 6, groups: Site, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.44816 0.08708 1.00000 39.599 0.0161 *
## SiteClassImpact -0.32021 0.15082 1.00000 -2.123 0.2802
## Period.L -0.10645 0.03369 1.00000 -3.160 0.1951
## SiteClassImpact:Period.L -0.08433 0.05836 1.00000 -1.445 0.3854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StClsI Perd.L
## StClssImpct -0.577
## Period.L 0.000 0.000
## StClssI:P.L 0.000 0.000 -0.577
# emmeans after a lm() fit
# Note that there is a emmeans() function in both the emmeans and lmerTest package
# so we must specify which one we want
# sink('baci-crabs-mod-R-s00LSM-SiteClass.txt', split=TRUE)
##***parts200LSM-SiteClassb;
result.avg.lmer.emmo.S <- emmeans::emmeans(result.avg.lmer, ~SiteClass)
## NOTE: Results may be misleading due to involvement in interactions
cat("\nEstimated marginal means for SiteClass \n\n")
##
## Estimated marginal means for SiteClass
summary(result.avg.lmer.emmo.S)
## SiteClass emmean SE df lower.CL upper.CL
## Control 3.45 0.0871 1 2.34 4.55
## Impact 3.13 0.1230 1 1.56 4.69
##
## Results are averaged over the levels of: Period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***parts200LSM-SiteClasse;
# sink()
# sink('baci-crabs-mod-R-200LSM-Period.txt', split=TRUE)
##***part200LSM-Periodb;
result.avg.lmer.emmo.P <- emmeans::emmeans(result.avg.lmer, ~Period)
## NOTE: Results may be misleading due to involvement in interactions
cat("\nEstimated marginal means for Period \n\n")
##
## Estimated marginal means for Period
summary(result.avg.lmer.emmo.P)
## Period emmean SE df lower.CL upper.CL
## Before 3.39 0.0782 1.15 2.66 4.13
## After 3.18 0.0782 1.15 2.45 3.92
##
## Results are averaged over the levels of: SiteClass
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part200LSM-Periode;
# sink()
# sink('baci-crabs-mod-R-200LSM-int.txt', split=TRUE)
##***part200LSM-intb;
result.avg.lmer.emmo.SP <- emmeans::emmeans(result.avg.lmer, ~SiteClass:Period)
cat("\nEstimated marginal means for SiteClass:Period \n\n")
##
## Estimated marginal means for SiteClass:Period
summary(result.avg.lmer.emmo.SP)
## SiteClass Period emmean SE df lower.CL upper.CL
## Control Before 3.52 0.0903 1.15 2.67 4.37
## Impact Before 3.26 0.1280 1.15 2.06 4.46
## Control After 3.37 0.0903 1.15 2.52 4.22
## Impact After 2.99 0.1280 1.15 1.79 4.19
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part200LSM-inte;
# sink()
# You could look at the entry in the summary table from the model fit, but
# this is dangerous as these entries depend on the contrast matrix.
# It is far safer to the contrast function applied to an emmeans object
temp <- summary(result.avg.lmer)$coefficients # get all the coefficients
temp[grepl("SiteClass",rownames(temp)) & grepl("Period", rownames(temp)),]
## Estimate Std. Error df t value Pr(>|t|)
## -0.08433052 0.05835596 1.00000005 -1.44510551 0.38536539
# sink("baci-crabs-mod-R-200baci.txt", split=TRUE)
##***part200bacib;
# Estimate the BACI contrast along with a se
summary(contrast(result.avg.lmer.emmo.SP, list(baci=c(1,-1,-1,1))), infer=TRUE)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## baci -0.119 0.0825 1 -1.17 0.929 -1.445 0.3854
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##***part200bacie;
# sink()
# Check the residuals etc
##***part200diagnosticb;
diagplot <- sf.autoplot.lmer(result.avg.lmer)
## Loading required package: broom.mixed
## Loading required package: grid
## Loading required package: gridExtra
## Loading required package: lattice
## Warning in indices[which(stats::complete.cases(original))] <- seq_len(nrow(x)):
## number of items to replace is not a multiple of replacement length
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : Chernobyl! trL>n 6
## Warning in sqrt(sum.squares/one.delta): NaNs produced
plot(diagplot)
##***part200diagnostice;
#ggsave.ggmultiplot(plotdiag, # see https://github.com/sinhrks/ggfortify/issues/98 for bug in autoplot
# file='baci-crabs-mod-R-200-diagnostic.png',
# h=4, w=6, unit="in", dpi=300)
#################################################################################
#################################################################################
# Do a Mixed effect linear model on the individual values
# Results will differ slightly from above analyeses on the averages because
# the design is not balanced.
# sink('baci-crabs-mod-R-300-type3.txt', split=TRUE)
##***part300b;
result.all.lmer <- lmerTest::lmer(log(Density) ~ SiteClass+Period+SiteClass:Period +
(1|Site) + (1|Site:Period), data=crabs)
## boundary (singular) fit: see help('isSingular')
##***part300e;
##***part301b;
anova(result.all.lmer, ddf="Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## SiteClass 0.058088 0.058088 1 1.0204 4.0892 0.2885
## Period 0.275358 0.275358 1 1.2061 19.3844 0.1096
## SiteClass:Period 0.021130 0.021130 1 1.2061 1.4875 0.4095
##***part301e;
# sink()
summary(result.all.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Density) ~ SiteClass + Period + SiteClass:Period + (1 | Site) +
## (1 | Site:Period)
## Data: crabs
##
## REML criterion at convergence: -25.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52883 -0.69753 0.05302 0.63592 1.89477
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site:Period (Intercept) 0.00000 0.0000
## Site (Intercept) 0.01507 0.1227
## Residual 0.01421 0.1192
## Number of obs: 30, groups: Site:Period, 6; Site, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.44706 0.09069 1.00298 38.007 0.01657 *
## SiteClassImpact -0.31911 0.15776 1.02041 -2.023 0.28840
## Period.L -0.10801 0.03721 25.01251 -2.902 0.00762 **
## SiteClassImpact:Period.L -0.08277 0.06768 25.00379 -1.223 0.23277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StClsI Perd.L
## StClssImpct -0.575
## Period.L 0.042 -0.024
## StClssI:P.L -0.023 -0.010 -0.550
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# sink('baci-crabs-mod-R-300-vc.txt', split=TRUE)
##***part300vcb;
# Extract the variance components
vc <- VarCorr(result.all.lmer)
vc
## Groups Name Std.Dev.
## Site:Period (Intercept) 0.00000
## Site (Intercept) 0.12274
## Residual 0.11919
##***part300vce;
# sink()
# emmeans after a lm() fit
# sink('baci-crabs-mod-R-s300LSM-SiteClass.txt', split=TRUE)
##***parts300LSM-SiteClassb;
result.all.lmer.emmo.S <- emmeans::emmeans(result.all.lmer, ~SiteClass)
## NOTE: Results may be misleading due to involvement in interactions
cat("\nEstimated marginal means for SiteClass\n\n")
##
## Estimated marginal means for SiteClass
summary(result.all.lmer.emmo.S)
## SiteClass emmean SE df lower.CL upper.CL
## Control 3.45 0.0908 1.00 2.30 4.59
## Impact 3.13 0.1290 1.03 1.59 4.66
##
## Results are averaged over the levels of: Period
## Degrees-of-freedom method: kenward-roger
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##***parts300LSM-SiteClasse;
# sink()
# sink('baci-crabs-mod-R-300LSM-Period.txt', split=TRUE)
##***part300LSM-Periodb;
result.all.lmer.emmo.P <- emmeans::emmeans(result.all.lmer, ~ Period)
## NOTE: Results may be misleading due to involvement in interactions
cat("\nEstimated marginal means for Period \n\n")
##
## Estimated marginal means for Period
summary(result.all.lmer.emmo.P)
## Period emmean SE df lower.CL upper.CL
## Before 3.39 0.0827 1.22 2.70 4.08
## After 3.18 0.0823 1.19 2.46 3.90
##
## Results are averaged over the levels of: SiteClass
## Degrees-of-freedom method: kenward-roger
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##***part300LSM-Periode;
# sink()
# sink('baci-crabs-mod-R-300LSM-int.txt', split=TRUE)
##***part300LSM-intb;
result.all.lmer.emmo.SP <- emmeans::emmeans(result.all.lmer, ~SiteClass:Period)
cat("\nEstimated marginal means for SiteClass:Period \n\n")
##
## Estimated marginal means for SiteClass:Period
summary(result.all.lmer.emmo.SP)
## SiteClass Period emmean SE df lower.CL upper.CL
## Control Before 3.52 0.0934 1.12 2.60 4.45
## Impact Before 3.26 0.1360 1.28 2.21 4.32
## Control After 3.37 0.0958 1.22 2.57 4.17
## Impact After 2.99 0.1340 1.18 1.80 4.19
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##***part300LSM-inte;
# sink()
# You could look at the entry in the summary table from the model fit, but
# this is dangerous as these entries depend on the contrast matrix.
# It is far safer to the contrast function applied to an emmeans object
temp <- summary(result.all.lmer)$coefficients # get all the coefficients
temp[grepl("SiteClass",rownames(temp)) & grepl("Period", rownames(temp)),]
## Estimate Std. Error df t value Pr(>|t|)
## -0.08276953 0.06768356 25.00379042 -1.22288972 0.23277289
# sink("baci-crabs-mod-R-300baci.txt", split=TRUE)
##***part300bacib;
# Estimate the BACI contrast along with a se
summary(contrast(result.all.lmer.emmo.SP, list(baci=c(1,-1,-1,1))), infer=TRUE)
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## baci -0.117 0.096 1.21 -0.939 0.705 -1.220 0.4095
##
## Degrees-of-freedom method: kenward-roger
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##***part300bacie;
# sink()
# Check the residuals etc
##***part300diagnosticb;
diagplot <- sf.autoplot.lmer(result.all.lmer)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
plot(diagplot)
##***part300diagnostice;
#ggsave.ggmultiplot(plotdiag, # see https://github.com/sinhrks/ggfortify/issues/98 for bug in autoplot
# file='baci-crabs-mod-R-300-diagnostic.png',
# h=4, w=6, units="in", dpi=300)
###############################################################################################################
###############################################################################################################
###############################################################################################################
###############################################################################################################
# Power analysis for BACI design with multiple sites, but one year before/after
# Example of crab - modifieds.
# Density is about 35 and a BACI difference on the log scale of .2 is important
# We use mu_IA=.2, mu_IB=0, mu_CB=0, mu_CA=0
#
# We don't have to worry about the sdSite or the sdYear because these "cancel"
# because every site is measured every year.
# The residual standard deviation is obtained from the above it.
#
# We examine various scenarios for the number of sites and number o subsamples
# An illustration of how to run for single case
baci.power(n_IA=5, n_IB=5, n_CA=5, n_CB=5,
ns_I=1, ns_C=2,
ny_B=1, ny_A=1,
mu_IA=30, mu_IB=35, mu_CA=35, mu_CB=35,
sdYear=0, sdSite=3.831, sdSiteYear=1.296, sdResid=3.30)
## alpha sdSite sdYear sdSiteYear sdResid n_IA n_IB n_CA n_CB ns_I ns_C ny_B
## 1 0.05 3.831 0 1.296 3.3 5 5 5 5 1 2 1
## ny_A mu_IA mu_IB mu_CA mu_CB baci baci.se dfdenom ncp Fcrit
## 1 1 30 35 35 35 -5 3.401889 1 2.160229 161.4476
## power Tcrit os.power1 os.power2
## 1 0.09574448 6.313752 0.1858008 0.003952477
# get the variance components
vc <- as.data.frame(VarCorr(result.all.lmer))
sdResid <- vc[ vc$grp=="Residual", "sdcor"]
sdSite <- vc[ vc$grp=="Site", "sdcor"]
sdSiteYear<- vc[ vc$grp=="Site:Period", "sdcor"]
##***part500b;
cat("Estimated variance components are \nsdSiteYear ", round(sdSiteYear,2),
";\n sdSite ", round(sdSite,2),
"; \nsdResid ",round(sdResid,2),"\n")
## Estimated variance components are
## sdSiteYear 0 ;
## sdSite 0.12 ;
## sdResid 0.12
cat("\n")
scenarios <- expand.grid(n_quad=seq(5,40,5),
ns_I = c(1,2),
ns_C = seq(2,10,2),
baci_effect=0.2,
sdSiteYear=sdSiteYear,
sdSite =sdSite,
sdResid=sdResid)
cat("Some scenarios\n")
## Some scenarios
head(scenarios)
## n_quad ns_I ns_C baci_effect sdSiteYear sdSite sdResid
## 1 5 1 2 0.2 0 0.1227443 0.1191853
## 2 10 1 2 0.2 0 0.1227443 0.1191853
## 3 15 1 2 0.2 0 0.1227443 0.1191853
## 4 20 1 2 0.2 0 0.1227443 0.1191853
## 5 25 1 2 0.2 0 0.1227443 0.1191853
## 6 30 1 2 0.2 0 0.1227443 0.1191853
power <- plyr::adply(scenarios,1,function(x){
#browser()
power <- baci.power(
n_IA=x$n_quad, n_IB=x$n_quad, n_CA=x$n_quad, n_CB=x$n_quad,
ns_I=x$ns_I, ns_C=x$ns_C,
ny_B=1, ny_A=1,
mu_IA=x$baci_effect, mu_IB=0, mu_CA=0, mu_CB=0,
sdYear=0, sdSite=x$sdSite, sdSiteYear=x$sdSiteYear, sdResid=x$sdResid)
power
})
cat("\nSome of the power computations\n")
##
## Some of the power computations
head(power)
## n_quad baci_effect alpha sdSite sdYear sdSiteYear sdResid n_IA n_IB n_CA
## 1 5 0.2 0.05 0.1227443 0 0 0.1191853 5 5 5
## 2 10 0.2 0.05 0.1227443 0 0 0.1191853 10 10 10
## 3 15 0.2 0.05 0.1227443 0 0 0.1191853 15 15 15
## 4 20 0.2 0.05 0.1227443 0 0 0.1191853 20 20 20
## 5 25 0.2 0.05 0.1227443 0 0 0.1191853 25 25 25
## 6 30 0.2 0.05 0.1227443 0 0 0.1191853 30 30 30
## n_CB ns_I ns_C ny_B ny_A mu_IA mu_IB mu_CA mu_CB baci baci.se dfdenom
## 1 5 1 2 1 1 0.2 0 0 0 0.2 0.09232057 1
## 2 10 1 2 1 1 0.2 0 0 0 0.2 0.06528050 1
## 3 15 1 2 1 1 0.2 0 0 0 0.2 0.05330131 1
## 4 20 1 2 1 1 0.2 0 0 0 0.2 0.04616029 1
## 5 25 1 2 1 1 0.2 0 0 0 0.2 0.04128701 1
## 6 30 1 2 1 1 0.2 0 0 0 0.2 0.03768972 1
## ncp Fcrit power Tcrit os.power1 os.power2
## 1 4.693135 161.4476 0.1356417 6.313752 0.2659881 6.776132e-04
## 2 9.386270 161.4476 0.1899989 6.313752 0.3682909 3.843098e-05
## 3 14.079405 161.4476 0.2315491 6.313752 0.4427867 2.629404e-06
## 4 18.772539 161.4476 0.2661000 6.313752 0.5020942 1.961622e-07
## 5 23.465674 161.4476 0.2961041 6.313752 0.5514234 1.539306e-08
## 6 28.158809 161.4476 0.3228404 6.313752 0.5935277 1.248884e-09
cat("\n")
head(power[,c("alpha","ns_I","ns_C","n_IA","n_IB","n_CA","n_CB","baci","power")])
## alpha ns_I ns_C n_IA n_IB n_CA n_CB baci power
## 1 0.05 1 2 5 5 5 5 0.2 0.1356417
## 2 0.05 1 2 10 10 10 10 0.2 0.1899989
## 3 0.05 1 2 15 15 15 15 0.2 0.2315491
## 4 0.05 1 2 20 20 20 20 0.2 0.2661000
## 5 0.05 1 2 25 25 25 25 0.2 0.2961041
## 6 0.05 1 2 30 30 30 30 0.2 0.3228404
##***part500e;
# sink("baci-crabs-mod-power-R-500.txt", split=TRUE) # get a smaller report
power[,c("alpha","ns_I","ns_C","n_IA","n_IB","n_CA","n_CB","baci","power")]
## alpha ns_I ns_C n_IA n_IB n_CA n_CB baci power
## 1 0.05 1 2 5 5 5 5 0.2 0.1356417
## 2 0.05 1 2 10 10 10 10 0.2 0.1899989
## 3 0.05 1 2 15 15 15 15 0.2 0.2315491
## 4 0.05 1 2 20 20 20 20 0.2 0.2661000
## 5 0.05 1 2 25 25 25 25 0.2 0.2961041
## 6 0.05 1 2 30 30 30 30 0.2 0.3228404
## 7 0.05 1 2 35 35 35 35 0.2 0.3470739
## 8 0.05 1 2 40 40 40 40 0.2 0.3693062
## 9 0.05 2 2 5 5 5 5 0.2 0.3259687
## 10 0.05 2 2 10 10 10 10 0.2 0.5217704
## 11 0.05 2 2 15 15 15 15 0.2 0.6606929
## 12 0.05 2 2 20 20 20 20 0.2 0.7592594
## 13 0.05 2 2 25 25 25 25 0.2 0.8291930
## 14 0.05 2 2 30 30 30 30 0.2 0.8788113
## 15 0.05 2 2 35 35 35 35 0.2 0.9140158
## 16 0.05 2 2 40 40 40 40 0.2 0.9389936
## 17 0.05 1 4 5 5 5 5 0.2 0.3765797
## 18 0.05 1 4 10 10 10 10 0.2 0.6186940
## 19 0.05 1 4 15 15 15 15 0.2 0.7750979
## 20 0.05 1 4 20 20 20 20 0.2 0.8702541
## 21 0.05 1 4 25 25 25 25 0.2 0.9262601
## 22 0.05 1 4 30 30 30 30 0.2 0.9585400
## 23 0.05 1 4 35 35 35 35 0.2 0.9768786
## 24 0.05 1 4 40 40 40 40 0.2 0.9871878
## 25 0.05 2 4 5 5 5 5 0.2 0.6364791
## 26 0.05 2 4 10 10 10 10 0.2 0.8921661
## 27 0.05 2 4 15 15 15 15 0.2 0.9713214
## 28 0.05 2 4 20 20 20 20 0.2 0.9928246
## 29 0.05 2 4 25 25 25 25 0.2 0.9982735
## 30 0.05 2 4 30 30 30 30 0.2 0.9995957
## 31 0.05 2 4 35 35 35 35 0.2 0.9999072
## 32 0.05 2 4 40 40 40 40 0.2 0.9999790
## 33 0.05 1 6 5 5 5 5 0.2 0.5090098
## 34 0.05 1 6 10 10 10 10 0.2 0.7913120
## 35 0.05 1 6 15 15 15 15 0.2 0.9199517
## 36 0.05 1 6 20 20 20 20 0.2 0.9711913
## 37 0.05 1 6 25 25 25 25 0.2 0.9900739
## 38 0.05 1 6 30 30 30 30 0.2 0.9966867
## 39 0.05 1 6 35 35 35 35 0.2 0.9989204
## 40 0.05 1 6 40 40 40 40 0.2 0.9996548
## 41 0.05 2 6 5 5 5 5 0.2 0.7723324
## 42 0.05 2 6 10 10 10 10 0.2 0.9661684
## 43 0.05 2 6 15 15 15 15 0.2 0.9958936
## 44 0.05 2 6 20 20 20 20 0.2 0.9995553
## 45 0.05 2 6 25 25 25 25 0.2 0.9999552
## 46 0.05 2 6 30 30 30 30 0.2 0.9999957
## 47 0.05 2 6 35 35 35 35 0.2 0.9999996
## 48 0.05 2 6 40 40 40 40 0.2 1.0000000
## 49 0.05 1 8 5 5 5 5 0.2 0.5767271
## 50 0.05 1 8 10 10 10 10 0.2 0.8568460
## 51 0.05 1 8 15 15 15 15 0.2 0.9580912
## 52 0.05 1 8 20 20 20 20 0.2 0.9887922
## 53 0.05 1 8 25 25 25 25 0.2 0.9971843
## 54 0.05 1 8 30 30 30 30 0.2 0.9993244
## 55 0.05 1 8 35 35 35 35 0.2 0.9998435
## 56 0.05 1 8 40 40 40 40 0.2 0.9999648
## 57 0.05 2 8 5 5 5 5 0.2 0.8352975
## 58 0.05 2 8 10 10 10 10 0.2 0.9847948
## 59 0.05 2 8 15 15 15 15 0.2 0.9989229
## 60 0.05 2 8 20 20 20 20 0.2 0.9999344
## 61 0.05 2 8 25 25 25 25 0.2 0.9999964
## 62 0.05 2 8 30 30 30 30 0.2 0.9999998
## 63 0.05 2 8 35 35 35 35 0.2 1.0000000
## 64 0.05 2 8 40 40 40 40 0.2 1.0000000
## 65 0.05 1 10 5 5 5 5 0.2 0.6162147
## 66 0.05 1 10 10 10 10 10 0.2 0.8882422
## 67 0.05 1 10 15 15 15 15 0.2 0.9725569
## 68 0.05 1 10 20 20 20 20 0.2 0.9939468
## 69 0.05 1 10 25 25 25 25 0.2 0.9987609
## 70 0.05 1 10 30 30 30 30 0.2 0.9997600
## 71 0.05 1 10 35 35 35 35 0.2 0.9999555
## 72 0.05 1 10 40 40 40 40 0.2 0.9999920
## 73 0.05 2 10 5 5 5 5 0.2 0.8692245
## 74 0.05 2 10 10 10 10 10 0.2 0.9913696
## 75 0.05 2 10 15 15 15 15 0.2 0.9995810
## 76 0.05 2 10 20 20 20 20 0.2 0.9999829
## 77 0.05 2 10 25 25 25 25 0.2 0.9999994
## 78 0.05 2 10 30 30 30 30 0.2 1.0000000
## 79 0.05 2 10 35 35 35 35 0.2 1.0000000
## 80 0.05 2 10 40 40 40 40 0.2 1.0000000
#sink()
##***part510b;
power.plot <- ggplot(data=power, aes(x=n_quad, y=power, color=as.factor(ns_C)))+
ggtitle("Estimated power",
subtitle=paste("alpha: ", power$alpha[1]))+
geom_point()+
geom_line()+
ylim(0,1)+
geom_hline(yintercept=0.8, color="blue")+
xlab("Number of quadrats at each site:year")+
scale_color_discrete(name="# Site\nControl")+
facet_wrap(~ns_I, ncol=2, labeller=label_both)
power.plot
##***part510e;
##***part520b;
# you could also have different sample sizes in impact an control areas.
results <- baci.power(n_IA=20, n_IB=20, n_CA=5, n_CB=5,
ns_I=1, ns_C=4,
ny_B=1, ny_A=1,
mu_IA=0.2, mu_IB=0, mu_CA=0, mu_CB=0,
sdYear=0, sdSite=sdSite, sdSiteYear=sdSiteYear, sdResid=sdResid)
results
## alpha sdSite sdYear sdSiteYear sdResid n_IA n_IB n_CA n_CB ns_I ns_C
## 1 0.05 0.1227443 0 0 0.1191853 20 20 5 5 1 4
## ny_B ny_A mu_IA mu_IB mu_CA mu_CB baci baci.se dfdenom ncp Fcrit
## 1 1 1 0.2 0 0 0 0.2 0.05330131 3 14.0794 10.12796
## power Tcrit os.power1 os.power2
## 1 0.7061564 2.353363 0.876833 5.764102e-07
##***part520e;