library(emon)
library(tidyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(glue)
library(ggplot2)
options(jupyter.plot_scale=1)
ff = read.csv(
file="forest_fuel_pilot.csv",
header=T
)
head(ff)
## PROJECT.CODE BA CI Plot.Num Point Transect Litter Canopy.cover LLB CBH x
## 1 DRYL_062421 B C 1 1 1 0.75 NA NA 7.5
## 2 DRYL_062421 B C 1 2 1 1.00 NA NA 7.5
## 3 DRYL_062421 B C 1 3 1 0.25 17.510 30 28.8 7.5
## 4 DRYL_062421 B C 1 4 1 1.00 NA NA 7.5
## 5 DRYL_062421 B C 1 5 1 0.25 NA NA 7.5
## 6 DRYL_062421 B C 1 6 1 0.75 42.199 NA 7.5
## y
## 1 3
## 2 6
## 3 9
## 4 12
## 5 15
## 6 18
# Format CBH measurements into separate columns and calculate mean for that microplot
cbh_cols = c('CBH1', 'CBH2', 'CBH3', 'CBH4', 'CBH5', 'CBH6')
ff2 = separate(ff, CBH, cbh_cols, sep=",", convert=TRUE)
## Warning: Expected 6 pieces. Missing pieces filled with `NA` in 360 rows [1, 2, 3, 4, 5,
## 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
ff2$CBH_mean = rowMeans(select(ff2, all_of(cbh_cols)), na.rm=TRUE)
head(ff2)
## PROJECT.CODE BA CI Plot.Num Point Transect Litter Canopy.cover LLB CBH1 CBH2
## 1 DRYL_062421 B C 1 1 1 0.75 NA NA NA NA
## 2 DRYL_062421 B C 1 2 1 1.00 NA NA NA NA
## 3 DRYL_062421 B C 1 3 1 0.25 17.510 30 28.8 NA
## 4 DRYL_062421 B C 1 4 1 1.00 NA NA NA NA
## 5 DRYL_062421 B C 1 5 1 0.25 NA NA NA NA
## 6 DRYL_062421 B C 1 6 1 0.75 42.199 NA NA NA
## CBH3 CBH4 CBH5 CBH6 x y CBH_mean
## 1 NA NA NA NA 7.5 3 NaN
## 2 NA NA NA NA 7.5 6 NaN
## 3 NA NA NA NA 7.5 9 28.8
## 4 NA NA NA NA 7.5 12 NaN
## 5 NA NA NA NA 7.5 15 NaN
## 6 NA NA NA NA 7.5 18 NaN
p <- ggplot(
ff2,
aes(x=Litter, fill=CI, color=CI)
) + geom_histogram(alpha=0.5, position="dodge", binwidth=0.25) +
facet_grid(Plot.Num ~ .) + geom_density(alpha=.2)
p
## Warning: Removed 120 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 120 rows containing non-finite values (`stat_density()`).
p <- ggplot(
ff2,
aes(x=Canopy.cover, fill=CI, color=CI)
) + geom_histogram(alpha=0.5, position="dodge", binwidth=5) +
facet_grid(Plot.Num ~ .) + geom_density(alpha=.2)
p
## Warning: Removed 240 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 240 rows containing non-finite values (`stat_density()`).
p <- ggplot(
ff2,
aes(x=LLB, fill=CI, color=CI)
) + geom_histogram(alpha=0.5, position="dodge", binwidth=5) +
facet_grid(Plot.Num ~ .) + geom_density(alpha=.2)
p
## Warning: Removed 298 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 298 rows containing non-finite values (`stat_density()`).
p <- ggplot(
ff2,
aes(x=CBH1, fill=CI, color=CI)
) + geom_histogram(alpha=0.5, position="dodge", binwidth=5) +
facet_grid(Plot.Num ~ .)
p
## Warning: Removed 295 rows containing non-finite values (`stat_bin()`).
p <- ggplot(
ff2,
aes(x=CBH_mean, fill=CI, color=CI)
) + geom_histogram(alpha=0.5, position="dodge", binwidth=5) +
facet_grid(Plot.Num ~ .)
p
## Warning: Removed 295 rows containing non-finite values (`stat_bin()`).
litter_df = setNames(
aggregate(
ff2$Litter,
list(ff2$BA, ff2$CI, ff2$Plot.Num, ff2$Transect),
mean,
na.rm=TRUE,
na.action=NULL
),
c('BA', 'CI', 'Plot', 'Transect', 'mean')
)
litter_df
## BA CI Plot Transect mean
## 1 B C 1 1 0.7833333
## 2 B I 1 1 1.3166667
## 3 B C 2 1 1.2000000
## 4 B I 2 1 0.8250000
## 5 B C 1 2 1.0000000
## 6 B I 1 2 0.8250000
## 7 B C 2 2 0.7000000
## 8 B I 2 2 0.7750000
## 9 B C 1 3 1.2250000
## 10 B I 1 3 0.9750000
## 11 B C 2 3 1.2750000
## 12 B I 2 3 1.1500000
litter = litter_df$mean
print(litter)
## [1] 0.7833333 1.3166667 1.2000000 0.8250000 1.0000000 0.8250000 0.7000000
## [8] 0.7750000 1.2250000 0.9750000 1.2750000 1.1500000
cc_df = setNames(
aggregate(
ff2$Canopy.cover,
list(ff2$Plot.Num, ff2$Transect, ff2$BA, ff2$CI),
mean,
na.rm=TRUE,
na.action=NULL
),
c('Plot', 'Transect', 'BA', 'CI', 'mean')
)
cc_df
## Plot Transect BA CI mean
## 1 1 1 B C 38.0627
## 2 2 1 B C 50.5257
## 3 1 2 B C 45.7746
## 4 2 2 B C 48.4897
## 5 1 3 B C 52.5877
## 6 2 3 B C 55.4446
## 7 1 1 B I 50.0808
## 8 2 1 B I 56.0394
## 9 1 2 B I 49.9197
## 10 2 2 B I 50.3985
## 11 1 3 B I 48.5185
## 12 2 3 B I 45.4680
cc = cc_df$mean
print(cc)
## [1] 38.0627 50.5257 45.7746 48.4897 52.5877 55.4446 50.0808 56.0394 49.9197
## [10] 50.3985 48.5185 45.4680
cbh_df = setNames(
aggregate(
ff2$CBH_mean,
list(ff2$Plot.Num, ff2$Transect, ff2$BA, ff2$CI),
mean,
na.rm=TRUE,
na.action=NULL
),
c('Plot', 'Transect', 'BA', 'CI', 'mean')
)
cbh_df
## Plot Transect BA CI mean
## 1 1 1 B C 22.44000
## 2 2 1 B C 16.60000
## 3 1 2 B C 21.30000
## 4 2 2 B C 17.68000
## 5 1 3 B C 24.48000
## 6 2 3 B C 10.46000
## 7 1 1 B I 19.80000
## 8 2 1 B I 16.40000
## 9 1 2 B I 25.68000
## 10 2 2 B I 21.60000
## 11 1 3 B I 17.18571
## 12 2 3 B I 19.44000
cbh = cbh_df$mean
cbh
## [1] 22.44000 16.60000 21.30000 17.68000 24.48000 10.46000 19.80000 16.40000
## [9] 25.68000 21.60000 17.18571 19.44000
llb_df = setNames(
aggregate(
ff2$LLB,
list(ff2$Plot.Num, ff2$Transect, ff2$BA, ff2$CI),
mean,
na.rm=TRUE,
na.action=NULL
),
c('Plot', 'Transect', 'BA', 'CI', 'mean')
)
llb_df
## Plot Transect BA CI mean
## 1 1 1 B C 21.45000
## 2 2 1 B C 11.33333
## 3 1 2 B C 27.83333
## 4 2 2 B C 13.20000
## 5 1 3 B C 20.06250
## 6 2 3 B C 11.30000
## 7 1 1 B I 14.83333
## 8 2 1 B I 13.45833
## 9 1 2 B I 12.00000
## 10 2 2 B I 13.62500
## 11 1 3 B I 15.78571
## 12 2 3 B I 11.60000
llb = llb_df$mean
llb
## [1] 21.45000 11.33333 27.83333 13.20000 20.06250 11.30000 14.83333 13.45833
## [9] 12.00000 13.62500 15.78571 11.60000
par(mfrow=c(2,4))
hist(litter)
hist(cc)
hist(llb)
hist(cbh)
qqnorm(litter, main='Original', pch=1, frame=FALSE); qqline(litter)
qqnorm(cc, main='Original', pch=1, frame=FALSE); qqline(cc)
qqnorm(llb, main='Original', pch=1, frame=FALSE); qqline(llb)
qqnorm(cbh, main='Original', pch=1, frame=FALSE); qqline(cbh)
par(mfrow=c(2,4))
hist(log(litter))
hist(log(cc))
hist(log(llb))
hist(log(cbh))
qqnorm((log(litter)), main='Log Transformed', pch=1, frame=FALSE); qqline(log(litter))
qqnorm((log(cc)), main='Log Transformed', pch=1, frame=FALSE); qqline(log(cc))
qqnorm((log(llb)), main='Log Transformed', pch=1, frame=FALSE); qqline(log(llb))
qqnorm((log(cbh)), main='Log Transformed', pch=1, frame=FALSE); qqline(log(cbh))
par(mfrow=c(2,4))
hist(sqrt(litter))
hist(sqrt(cc))
hist(sqrt(llb))
hist(sqrt(cbh))
qqnorm((sqrt(litter)), main='Sqrt Transformed', pch=1, frame=FALSE); qqline(sqrt(litter))
qqnorm((sqrt(cc)), main='Sqrt Transformed',pch=1, frame=FALSE); qqline(sqrt(cc))
qqnorm((sqrt(llb)), main='Sqrt Transformed',pch=1, frame=FALSE); qqline(sqrt(llb))
qqnorm((sqrt(cbh)), main='Sqrt Transformed',pch=1, frame=FALSE); qqline(sqrt(cbh))
litter_ttest = t.test(
x=litter_df[which (litter_df$CI == 'C'), ]$mean,
y=litter_df[which (litter_df$CI == 'I'), ]$mean,
alternative = 'two.sided',
)
litter_ttest
##
## Welch Two Sample t-test
##
## data: litter_df[which(litter_df$CI == "C"), ]$mean and litter_df[which(litter_df$CI == "I"), ]$mean
## t = 0.39711, df = 9.8505, p-value = 0.6998
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2439619 0.3495175
## sample estimates:
## mean of x mean of y
## 1.0305556 0.9777778
cc_ttest = t.test(
x=cc_df[which (cc_df$CI == 'C'), ]$mean,
y=cc_df[which (cc_df$CI == 'I'), ]$mean,
alternative = 'two.sided',
)
cc_ttest
##
## Welch Two Sample t-test
##
## data: cc_df[which(cc_df$CI == "C"), ]$mean and cc_df[which(cc_df$CI == "I"), ]$mean
## t = -0.55666, df = 7.9036, p-value = 0.5932
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -8.190629 5.010662
## sample estimates:
## mean of x mean of y
## 48.48083 50.07082
## Note: not always based on 30 observations per group
cbh_ttest = t.test(
x=cbh_df[which (cbh_df$CI == 'C'), ]$mean,
y=cbh_df[which (cbh_df$CI == 'I'), ]$mean,
alternative = 'two.sided',
)
cbh_ttest
##
## Welch Two Sample t-test
##
## data: cbh_df[which(cbh_df$CI == "C"), ]$mean and cbh_df[which(cbh_df$CI == "I"), ]$mean
## t = -0.48159, df = 8.6861, p-value = 0.642
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.816087 4.434182
## sample estimates:
## mean of x mean of y
## 18.82667 20.01762
## Note: not always based on 30 observations per group
llb_ttest = t.test(
x=llb_df[which (llb_df$CI == 'C'), ]$mean,
y=llb_df[which (llb_df$CI == 'I'), ]$mean,
alternative = 'two.sided',
)
llb_ttest
##
## Welch Two Sample t-test
##
## data: llb_df[which(llb_df$CI == "C"), ]$mean and llb_df[which(llb_df$CI == "I"), ]$mean
## t = 1.4165, df = 5.5727, p-value = 0.21
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.024445 10.983373
## sample estimates:
## mean of x mean of y
## 17.52986 13.55040
ssize = seq(2, 50, 3)
parst.litter = rep(0,2)
parsc.litter = parst.litter
parst.litter[1] = mean(litter)
parst.litter[2] = sd(litter)
parsc.litter = parst.litter
parst.cc = rep(0,2)
parsc.cc = parst.cc
parst.cc[1] = mean(cc)
parst.cc[2] = sd(cc)
parsc.cc = parst.cc
parst.llb = rep(0,2)
parsc.llb = parst.llb
parst.llb[1] = mean(llb)
parst.llb[2] = sd(llb)
parsc.llb = parst.llb
parst.cbh = rep(0,2)
parsc.cbh = parst.cbh
parst.cbh[1] = mean(cbh)
parst.cbh[2] = sd(cbh)
parsc.cbh = parst.cbh
power.litter = rep(0, length(ssize))
power.cc = rep(0, length(ssize))
power.llb = rep(0, length(ssize))
power.cbh = rep(0, length(ssize))
pct_change = 10
power.litter = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.litter, parsc=parsc.litter,
distribution="Normal", test="P", nsims=1000)$power
power.cc = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cc, parsc=parsc.cc, distribution="Normal",
test="P", nsims=1000)$power
power.llb = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.llb, parsc=parsc.llb,
distribution="Normal", test="P", nsims=1000)$power
power.cbh = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cbh, parsc=parsc.cbh, distribution="Normal",
test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.litter, ylim=c(0,1), ylab="Power", xlab="Number of Transects", type="l")
lines(ssize, power.cc, lty=2, col=2)
lines(ssize, power.llb, lty=3, col=3)
lines(ssize, power.cbh, lty=6, col=4)
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Litter Depth Power",
"Canopy Cover Power",
"Lowest Live Branch Power",
"CBH Power",
"Transects needed for 80% power"
),
lty=c(1,2, 3, 6, 4),
col=c(1,2, 3, 4, 'grey'))
title(glue("BACI power plots - {pct_change}% change"))
pct_change = 20
power.litter = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.litter, parsc=parsc.litter,
distribution="Normal", test="P", nsims=1000)$power
power.cc = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cc, parsc=parsc.cc, distribution="Normal",
test="P", nsims=1000)$power
power.llb = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.llb, parsc=parsc.llb,
distribution="Normal", test="P", nsims=1000)$power
power.cbh = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cbh, parsc=parsc.cbh, distribution="Normal",
test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.litter, ylim=c(0,1), ylab="Power", xlab="Number of Transects", type="l")
lines(ssize, power.cc, lty=2, col=2)
lines(ssize, power.llb, lty=3, col=3)
lines(ssize, power.cbh, lty=6, col=4)
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Litter Depth Power",
"Canopy Cover Power",
"Lowest Live Branch Power",
"CBH Power",
"Transects needed for 80% power"
),
lty=c(1,2, 3, 6, 4),
col=c(1,2, 3, 4, 'grey'))
title(glue("BACI power plots - {pct_change}% change"))
pct_change = 25
power.litter = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.litter, parsc=parsc.litter,
distribution="Normal", test="P", nsims=1000)$power
power.cc = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cc, parsc=parsc.cc, distribution="Normal",
test="P", nsims=1000)$power
power.llb = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.llb, parsc=parsc.llb,
distribution="Normal", test="P", nsims=1000)$power
power.cbh = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cbh, parsc=parsc.cbh, distribution="Normal",
test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.litter, ylim=c(0,1), ylab="Power", xlab="Number of Transects", type="l")
lines(ssize, power.cc, lty=2, col=2)
lines(ssize, power.llb, lty=3, col=3)
lines(ssize, power.cbh, lty=6, col=4)
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Litter Depth Power",
"Canopy Cover Power",
"Lowest Live Branch Power",
"CBH Power",
"Transects needed for 80% power"
),
lty=c(1,2, 3, 6, 4),
col=c(1,2, 3, 4, 'grey'))
title(glue("BACI power plots - {pct_change}% change"))
pct_change = 30
power.litter = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.litter, parsc=parsc.litter,
distribution="Normal", test="P", alpha=0.05, nsims=1000)$power
power.cc = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cc, parsc=parsc.cc, distribution="Normal",
test="P", alpha=0.05, nsims=1000)$power
power.llb = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.llb, parsc=parsc.llb,
distribution="Normal", test="P", alpha=0.05, nsims=1000)$power
power.cbh = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cbh, parsc=parsc.cbh, distribution="Normal",
test="P", alpha=0.05, nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.litter, ylim=c(0,1), ylab="Power", xlab="Number of Transects", type="l")
lines(ssize, power.cc, lty=2, col=2)
lines(ssize, power.llb, lty=3, col=3)
lines(ssize, power.cbh, lty=6, col=4)
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Litter Depth Power",
"Canopy Cover Power",
"Lowest Live Branch Power",
"CBH Power",
"Transects needed for 80% power"
),
lty=c(1,2, 3, 6, 4),
col=c(1,2, 3, 4, 'grey'))
title(glue("BACI power plots - {pct_change}% change"))
pct_change = 30
ssize = seq(2, 50, 3)
litter2_c = litter_df[litter_df$CI == 'C', ]$mean
litter2_i = litter_df[litter_df$CI == 'I', ]$mean
parst.litter2 = rep(0,2)
parsc.litter2 = parst.litter2
parst.litter2[1] = mean(litter2_i)
parst.litter2[2] = sd(litter2_i)
parsc.litter2[1] = mean(litter2_c)
parsc.litter2[2] = sd(litter2_c)
parst.litter
## [1] 1.0041667 0.2212082
parsc.litter2
## [1] 1.0305556 0.2439642
parst.litter2
## [1] 0.9777778 0.2155527
power.litter2 = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.litter2, parsc=parsc.litter2,
distribution="Normal", test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.litter, ylim=c(0,1), ylab="Power", xlab="Number of Transects", type="l")
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Litter Depth Power",
"Transects needed for 80% power"
),
lty=c(1,4),
col=c(1,'grey'))
title(glue("BACI power plots - {pct_change}% change"))
pct_change = 50
power.litter = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.litter, parsc=parsc.litter,
distribution="Normal", test="P", nsims=1000)$power
power.cc = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cc, parsc=parsc.cc, distribution="Normal",
test="P", nsims=1000)$power
power.llb = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.llb, parsc=parsc.llb,
distribution="Normal", test="P", nsims=1000)$power
power.cbh = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.cbh, parsc=parsc.cbh, distribution="Normal",
test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.litter, ylim=c(0,1), ylab="Power", xlab="Number of Transects", type="l")
lines(ssize, power.cc, lty=2, col=2)
lines(ssize, power.llb, lty=3, col=3)
lines(ssize, power.cbh, lty=6, col=4)
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Litter Depth Power",
"Canopy Cover Power",
"Lowest Live Branch Power",
"CBH Power",
"Transects needed for 80% power"
),
lty=c(1,2, 3, 6, 4),
col=c(1,2, 3, 4, 'grey'))
title(glue("BACI power plots - {pct_change}% change"))
plot = read.csv(
file="forest_fuel_plot_data.csv",
header=T
)
plot
## PROJECT.CODE BA CI Plot.Num Dead_Decline_Trees Dom_Species Lon
## 1 DRYL_062421 B C 1 10 Lodgepole pine -106.2189
## 2 DRYL_062421 B C 2 58 Lodgepole pine -106.2212
## 3 DRYL_062421 B I 1 29 Lodgepole pine -106.2499
## 4 DRYL_062421 B I 2 18 Lodgepole pine -106.2200
## Lat
## 1 38.56509
## 2 38.56418
## 3 38.56484
## 4 38.56454
num_trees = plot$Dead_Decline_Trees
num_trees
## [1] 10 58 29 18
ssize = seq(2, 20, 2)
parst.tree = mean(num_trees)
parsc.tree = parst.tree
power.tree = rep(0, length(ssize))
pct_change = 50
power.tree = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.tree, parsc=parsc.tree, distribution="Poisson",
test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.tree, ylim=c(0,1), ylab="Power", xlab="Number of Plots", type="l")
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Number Trees Power",
"Transects needed for 80% power"
),
lty=c(1,4),
col=c(1,'grey'))
title(glue("BACI power plots - {pct_change}% change"))
pct_change = 90
power.tree = power.BACI(change= pct_change, change.type="M", nt=ssize, nc=ssize,
parst=parst.tree, parsc=parsc.tree, distribution="Poisson",
test="P", nsims=1000)$power
par(mfrow=c(1,1))
plot(ssize, power.tree, ylim=c(0,1), ylab="Power", xlab="Number of Plots", type="l")
abline(h=0.8, col="grey", lty=4)
legend("bottomright", legend=c(
"Number Trees Power",
"Transects needed for 80% power"
),
lty=c(1,4),
col=c(1,'grey'))
title(glue("BACI power plots - {pct_change}% change"))
power.groups(
change=50, # % change
change.type="M",
n1=12,
n2=12,
pars1=c(80, 8), # mean before, standard deviation before
pars2=(35), # mean after
test='P',
distribution='Normal',
nsims=1000
)
## [1] 0.942