Neural Cluster 10 and Comparison of Quantification Methods
2025-06-17
library(dplyr) ## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)#This dataframe has un-separated hemispheres
MomentsU <- read.csv("C:/Users/gibbsj1/OneDrive - College of Charleston/Jesi CofC/Thesis/Neurobiological Method/MomentsUnseparated.csv")
#create unique section names
MomentsU$section <- interaction(MomentsU$round, MomentsU$subject, MomentsU$section, sep = "_")Percent Area Averaged by Subject
subject_avg <- MomentsU %>%
group_by(subject, hot) %>%
summarise(mean_X_area = mean(X.area, na.rm = TRUE), .groups = "drop")
group_means <- subject_avg %>%
group_by(hot) %>%
summarise(group_mean = mean(mean_X_area), .groups = "drop")
ggplot(subject_avg, aes(x = hot, y = mean_X_area)) +
geom_jitter(width = 0.1, size = 1, alpha = 0.8) +
geom_point(data = group_means, aes(x = hot, y = group_mean),
shape = 17, size = 3, color = "black") + # triangle group means
labs(title = "Percent Area by Treatment",
x = "Treatment",
y = "Average Percent Area") +
theme_minimal()## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
Nested Analysis % Area
library(lmerTest)## Warning: package 'lmerTest' was built under R version 4.3.3## 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':
##
## stepmodel <- lmer(X.area ~ hot +
(1 | round) +
(1 | subject) +
(1 | subject:section),
data = MomentsU)
nullmodel <- lmer(X.area ~ 1 +
(1 | round) +
(1 | subject) +
(1 | subject:section),
data = MomentsU)
# Basic variance components
anova(model,nullmodel)## refitting model(s) with ML (instead of REML)## Data: MomentsU
## Models:
## nullmodel: X.area ~ 1 + (1 | round) + (1 | subject) + (1 | subject:section)
## model: X.area ~ hot + (1 | round) + (1 | subject) + (1 | subject:section)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 5 403.23 413.53 -196.61 393.23
## model 6 401.22 413.58 -194.61 389.22 4.011 1 0.0452 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Model including treatment is significantly better than null model and AIC is better. Likelihood ratio test is significant. Treatment has a significant effect on percent area stained.
model <- lmer(X.area ~ hot +
(1 | round) +
(1 | subject) +
(1 | section),
data = MomentsU)
nullmodel <- lmer(X.area ~ 1 +
(1 | round) +
(1 | subject) +
(1 | section),
data = MomentsU)
# Basic variance components
anova(model,nullmodel)## refitting model(s) with ML (instead of REML)## Data: MomentsU
## Models:
## nullmodel: X.area ~ 1 + (1 | round) + (1 | subject) + (1 | section)
## model: X.area ~ hot + (1 | round) + (1 | subject) + (1 | section)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodel 5 403.23 413.53 -196.61 393.23
## model 6 401.22 413.58 -194.61 389.22 4.011 1 0.0452 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#section and subject:section test the same thing as long as section names are unique
Testing sex effects:
modelsex <- lmer(X.area ~ hot + sex +
(1 | round) +
(1 | subject) +
(1 | section),
data = MomentsU)
nullmodelsex <- lmer(X.area ~ 1 +
(1 | round) +
(1 | subject) +
(1 | section),
data = MomentsU)
anova(modelsex,nullmodelsex)## refitting model(s) with ML (instead of REML)## Data: MomentsU
## Models:
## nullmodelsex: X.area ~ 1 + (1 | round) + (1 | subject) + (1 | section)
## modelsex: X.area ~ hot + sex + (1 | round) + (1 | subject) + (1 | section)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodelsex 5 403.23 413.53 -196.61 393.23
## modelsex 7 402.61 417.03 -194.30 388.61 4.6208 2 0.09922 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Sex accounts for some variation in the dataset, but is not significantly different from the null model with only random variation
# Basic variance components
anova(modelsex,model)## refitting model(s) with ML (instead of REML)## Data: MomentsU
## Models:
## model: X.area ~ hot + (1 | round) + (1 | subject) + (1 | section)
## modelsex: X.area ~ hot + sex + (1 | round) + (1 | subject) + (1 | section)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## model 6 401.22 413.58 -194.61 389.22
## modelsex 7 402.61 417.03 -194.30 388.61 0.6097 1 0.4349#Sex does not account for more variance than the treatment-only model
Checking Residuals
# Raw residuals
resid_raw <- resid(model)
# Fitted values
fitted_vals <- fitted(model)
# looking for random cloud of points - funnels or patterns indicate heteroscedastity
plot(fitted_vals, resid_raw,
xlab = "Fitted values",
ylab = "Residuals",
main = "Residuals vs Fitted")
abline(h = 0, lty = 2)
#Normal Q-Q plot
#Residuals should lie on the diagonal line if normally distributed.
qqnorm(resid_raw)
qqline(resid_raw, col = "red")
#Histogram of Residuals
#To visualize overall distribution.
hist(resid_raw, breaks = 20, main = "Histogram of Residuals",
xlab = "Residuals")
library(performance)
check_model(model)
#Residuals look pretty good. Going to try a transformation to see if it is better.
Arcsine Square Root Transformation
# Step 1: Convert to proportion (if your values are 0–100)
MomentsU$X.area_prop <- MomentsU$X.area / 100
# Step 2: Apply arcsine square root transformation
MomentsU$X.area_trans <- asin(sqrt(MomentsU$X.area_prop))model_trans <- lmer(X.area_trans ~ hot + (1 | round) + (1 | subject) + (1 | subject:section),
data = MomentsU)#Checking transformation
# Raw residuals
resid_raw <- resid(model_trans)
# Fitted values
fitted_vals <- fitted(model_trans)
# looking for random cloud of points - funnels or patterns indicate heteroscedastity
plot(fitted_vals, resid_raw,
xlab = "Fitted values",
ylab = "Residuals",
main = "Residuals vs Fitted")
abline(h = 0, lty = 2)
#Normal Q-Q plot
#Residuals should lie on the diagonal line if normally distributed.
qqnorm(resid_raw)
qqline(resid_raw, col = "red")
#Histogram of Residuals
#To visualize overall distribution.
hist(resid_raw, breaks = 20, main = "Histogram of Residuals",
xlab = "Residuals")
check_model(model_trans)
Influential observations is much improved. Otherwise looks the same. Influential observations are data points that disproportionately affect model estimates (e.g., high leverage or large residuals).The fact that these look better after transformation implies that the transformation reduced heteroscedasticity or outlier influence, which often happens with bounded data like proportions.
nullmodeltrans <- lmer(X.area_trans ~ 1 +
(1 | round) +
(1 | subject) +
(1 | section),
data = MomentsU)
anova(model_trans, nullmodeltrans)## refitting model(s) with ML (instead of REML)## Data: MomentsU
## Models:
## nullmodeltrans: X.area_trans ~ 1 + (1 | round) + (1 | subject) + (1 | section)
## model_trans: X.area_trans ~ hot + (1 | round) + (1 | subject) + (1 | subject:section)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodeltrans 5 -117.79 -107.49 63.896 -127.79
## model_trans 6 -119.91 -107.54 65.953 -131.91 4.1138 1 0.04253 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: X.area ~ hot + (1 | round) + (1 | subject) + (1 | section)
## Data: MomentsU
##
## REML criterion at convergence: 384
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0263 -0.5830 -0.2660 0.5456 2.1934
##
## Random effects:
## Groups Name Variance Std.Dev.
## section (Intercept) 6.2238 2.4947
## subject (Intercept) 4.4124 2.1006
## round (Intercept) 0.9812 0.9906
## Residual 40.6114 6.3727
## Number of obs: 58, groups: section, 33; subject, 12; round, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.998 1.661 1.283 15.654 0.0201 *
## hoty 4.590 2.356 4.785 1.949 0.1115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.580summary(model_trans)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## X.area_trans ~ hot + (1 | round) + (1 | subject) + (1 | subject:section)
## Data: MomentsU
##
## REML criterion at convergence: -119.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3448 -0.5815 -0.2556 0.5759 2.0793
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject:section (Intercept) 0.0009057 0.03010
## subject (Intercept) 0.0004788 0.02188
## round (Intercept) 0.0001342 0.01158
## Residual 0.0050237 0.07088
## Number of obs: 58, groups: subject:section, 33; subject, 12; round, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.53165 0.01849 1.15585 28.754 0.0136 *
## hoty 0.05202 0.02609 4.47671 1.994 0.1093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.570#The effect of hot is positive and marginally significant in both models.
#With data that are bounded percentages (0–100), the transformed model is statistically safer and likely preferable.
Use transformation for data.
Hemispheric Differences
#This df has absolute difference values
Moments <- read.csv("C:/Users/gibbsj1/OneDrive - College of Charleston/Jesi CofC/Thesis/Neurobiological Method/Moments.csv")
#create unique section names
Moments$section <- interaction(Moments$round, Moments$subject, Moments$section, sep = "_")
#Arcsine Square Root Transform
# Step 1: Convert to proportion (if your values are 0–100)
Moments$X.areadiff_prop <- Moments$X.areadiff / 100
# Step 2: Apply arcsine square root transformation
Moments$X.areadiff_trans <- asin(sqrt(Moments$X.areadiff_prop))group_means <- Moments %>%
group_by(hot) %>%
summarise(group_mean = mean(X.areadiff), .groups = "drop")
group_means## # A tibble: 2 × 2
## hot group_mean
## <chr> <dbl>
## 1 n 5.31
## 2 y 9.35subject_avg <- Moments %>%
group_by(subject, hot) %>%
summarise(mean_X_areadiff = mean(X.areadiff, na.rm = TRUE), .groups = "drop")
group_means <- subject_avg %>%
group_by(hot) %>%
summarise(group_mean = mean(mean_X_areadiff), .groups = "drop")
ggplot(subject_avg, aes(x = hot, y = mean_X_areadiff)) +
geom_jitter(width = 0.1, size = 1, alpha = 0.8) +
geom_point(data = group_means, aes(x = hot, y = group_mean),
shape = 17, size = 3, color = "black") + # triangle group means
labs(title = "Hemisphere Asymmetry by Treatment",
x = "Treatment",
y = "Average Percent Area Absolute Difference") +
theme_minimal()
Nested Anlaysis for Difference Scores (%area)
modeldiff <- lmer(X.areadiff ~ hot +
(1 | subject),
data = Moments)
nullmodeldiff <- lmer(X.areadiff ~ 1 +
(1 | subject),
data = Moments)
# Basic variance components
anova(modeldiff,nullmodeldiff)## refitting model(s) with ML (instead of REML)## Data: Moments
## Models:
## nullmodeldiff: X.areadiff ~ 1 + (1 | subject)
## modeldiff: X.areadiff ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodeldiff 3 161.99 165.64 -77.993 155.99
## modeldiff 4 161.54 166.41 -76.769 153.54 2.4476 1 0.1177modeldifftrans <- lmer(X.areadiff_trans ~ hot +
(1 | subject),
data = Moments)
nullmodeldifftrans <- lmer(X.areadiff_trans ~ 1 +
(1 | subject),
data = Moments)
# Basic variance components
anova(modeldifftrans,nullmodeldifftrans)## refitting model(s) with ML (instead of REML)## Data: Moments
## Models:
## nullmodeldifftrans: X.areadiff_trans ~ 1 + (1 | subject)
## modeldifftrans: X.areadiff_trans ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## nullmodeldifftrans 3 -32.298 -28.642 19.149 -38.298
## modeldifftrans 4 -33.735 -28.860 20.868 -41.735 3.4372 1 0.06374 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1check_model(modeldiff)
check_model(modeldifftrans)
#Multiple improvements. Definitely use transformation with %area data
Comparing quantification methods
meth<-MethodComparison <- read.csv("C:/Users/gibbsj1/OneDrive - College of Charleston/Jesi CofC/Thesis/Neurobiological Method/MethodComparison.csv")
#create unique section names
meth$section <- interaction(meth$round, meth$subject, meth$section, sep = "_")
#Arcsine Square Root Transform
# Step 1: Convert to proportion (if your values are 0–100)
meth$X.area_prop <- meth$X.area / 100
# Step 2: Apply arcsine square root transformation
meth$X.area_trans <- asin(sqrt(meth$X.area_prop))hotmodel <- lmer(X.area_trans ~ hot +
(1 | round) +
(1 | subject) +
(1 | section),
data = meth)## boundary (singular) fit: see help('isSingular')modelmethod <- lmer(X.area_trans ~ hot * method +
(1 | round) +
(1 | subject) +
(1 | section),
data = meth)## boundary (singular) fit: see help('isSingular')modelmethodrandom <- lmer(X.area_trans ~ hot +
(1 | round) +
(1 | subject) +
(1 | section) +
(1| method),
data = meth)## boundary (singular) fit: see help('isSingular')anova(modelmethod,hotmodel)## refitting model(s) with ML (instead of REML)## Data: meth
## Models:
## hotmodel: X.area_trans ~ hot + (1 | round) + (1 | subject) + (1 | section)
## modelmethod: X.area_trans ~ hot * method + (1 | round) + (1 | subject) + (1 | section)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## hotmodel 6 -151.40 -134.88 81.702 -163.40
## modelmethod 8 -155.46 -133.44 85.732 -171.46 8.0611 2 0.01776 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(hotmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: X.area_trans ~ hot + (1 | round) + (1 | subject) + (1 | section)
## Data: meth
##
## REML criterion at convergence: -152.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3825 -0.5542 -0.0075 0.6564 3.5653
##
## Random effects:
## Groups Name Variance Std.Dev.
## section (Intercept) 0.00000 0.00000
## subject (Intercept) 0.00221 0.04702
## round (Intercept) 0.00000 0.00000
## Residual 0.01328 0.11526
## Number of obs: 116, groups: section, 33; subject, 12; round, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.50157 0.02416 6.95786 20.761 1.62e-07 ***
## hoty 0.12322 0.03590 8.65931 3.433 0.00792 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.673
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')summary(modelmethod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: X.area_trans ~ hot * method + (1 | round) + (1 | subject) + (1 |
## section)
## Data: meth
##
## REML criterion at convergence: -149.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7395 -0.5876 -0.0374 0.4886 3.3536
##
## Random effects:
## Groups Name Variance Std.Dev.
## section (Intercept) 0.000000 0.00000
## subject (Intercept) 0.002324 0.04821
## round (Intercept) 0.000000 0.00000
## Residual 0.012526 0.11192
## Number of obs: 116, groups: section, 33; subject, 12; round, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.47337 0.02759 11.77264 17.157 1.08e-09 ***
## hoty 0.18430 0.04201 16.00549 4.387 0.000459 ***
## methodmoments 0.05609 0.02602 101.12131 2.156 0.033476 *
## hoty:methodmoments -0.12161 0.04324 101.12131 -2.812 0.005913 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hoty mthdmm
## hoty -0.657
## methodmmnts -0.472 0.310
## hty:mthdmmn 0.284 -0.515 -0.602
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')summary(modelmethodrandom)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: X.area_trans ~ hot + (1 | round) + (1 | subject) + (1 | section) +
## (1 | method)
## Data: meth
##
## REML criterion at convergence: -152.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3825 -0.5542 -0.0075 0.6564 3.5653
##
## Random effects:
## Groups Name Variance Std.Dev.
## section (Intercept) 2.062e-12 1.436e-06
## subject (Intercept) 2.210e-03 4.701e-02
## method (Intercept) 0.000e+00 0.000e+00
## round (Intercept) 1.724e-12 1.313e-06
## Residual 1.328e-02 1.153e-01
## Number of obs: 116, groups: section, 33; subject, 12; method, 2; round, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.50157 0.02416 6.95791 20.761 1.62e-07 ***
## hoty 0.12322 0.03590 8.65939 3.433 0.00792 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.673
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')library(ggplot2)
ggplot(meth, aes(x = method, y = X.area, group = interaction(subject, section))) +
geom_point(aes(color = method), size = 3) +
geom_line(aes(group = interaction(subject, section)), color = "gray") +
labs(title = "Percent Area Stained by Quantification Method",
y = "Percent Area Stained",
x = "Quantification Method") +
theme_minimal()
df_wide<- MethodComparisonWide <- read.csv("C:/Users/gibbsj1/OneDrive - College of Charleston/Jesi CofC/Thesis/Neurobiological Method/MethodComparisonWide.csv")
#create unique section names
df_wide$section <- interaction(df_wide$round, df_wide$subject, df_wide$section, sep = "_")
head(df_wide)## location hot sex round subject section hemi method moments manual avg
## 1 C10 n f 1 CC 1_CC_4 a moments 43.818 32.268 38.0430
## 2 C10 n f 1 CC 1_CC_4 b moments 24.950 26.512 25.7310
## 3 C10 n f 1 CC 1_CC_5 a moments 31.235 10.167 20.7010
## 4 C10 n f 1 CC 1_CC_5 b moments 32.220 15.747 23.9835
## 5 C10 n m 2 Z 2_Z_3 a moments 43.762 30.980 37.3710
## 6 C10 n m 2 Z 2_Z_4 a moments 27.381 11.799 19.5900
## diff
## 1 11.550
## 2 -1.562
## 3 21.068
## 4 16.473
## 5 12.782
## 6 15.582library(ggplot2)
ggplot(df_wide, aes(x = avg, y = diff)) +
geom_point(size = 2, alpha = 0.7) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_hline(yintercept = mean(df_wide$diff), color = "blue") +
geom_hline(yintercept = mean(df_wide$diff) + 1.96 * sd(df_wide$diff), linetype = "dotted", color = "red") +
geom_hline(yintercept = mean(df_wide$diff) - 1.96 * sd(df_wide$diff), linetype = "dotted", color = "red") +
labs(title = "Bland-Altman Plot: Moments vs Manual",
x = "Mean of Methods",
y = "Difference (Moments - Manual)") +
theme_minimal()
#No clear bias of either method (one doesn’t consistently over/underestimate) - blue line
#No clear agreement between methods (red lines) and a few outliers
ggplot(df_wide, aes(x = manual, y = moments)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
labs(title = "Scatterplot: Moments vs Manual",
x = "Manual Percent Area",
y = "Moments Percent Area") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
#No good linear relationship between methods. Test:
# Pearson: assumes linear relationship
cor.test(df_wide$manual, df_wide$moments, method = "pearson")##
## Pearson's product-moment correlation
##
## data: df_wide$manual and df_wide$moments
## t = 2.2963, df = 56, p-value = 0.02542
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.03793473 0.51279596
## sample estimates:
## cor
## 0.2933561# Spearman: for monotonic but not necessarily linear
cor.test(df_wide$manual, df_wide$moments, method = "spearman")## Warning in cor.test.default(df_wide$manual, df_wide$moments, method =
## "spearman"): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: df_wide$manual and df_wide$moments
## S = 26015, p-value = 0.1327
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1997631#Weak but existent positive correlation. Can’t use moments in place of manual.
#Compare method means:
library(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'# Get EMMs for each combination of hot and method
emm <- emmeans(modelmethod, ~ hot * method)library(ggplot2)
# Convert emmeans result to a data frame
emm_df <- as.data.frame(emm)
# Plot
ggplot(emm_df, aes(x = hot, y = emmean, color = method, group = method)) +
geom_point(position = position_dodge(width = 0.2), size = 3) +
geom_line(position = position_dodge(width = 0.2), size = 1) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
width = 0.1, position = position_dodge(width = 0.2)) +
labs(title = "Interaction of Hot Treatment and Quantification Method",
x = "Treatment (Hot vs Not)",
y = "Estimated Marginal Mean (Transformed Area)",
color = "Method") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#not great agreement between methods. Manual method shows stronger effect of treatment but could be biased.
#Continue to use both methods and include as a fixed effect in the model?
#Could also unbias the manual method by blinding myself to treatment or by having an additional person quantify the same areas and checking inter-rater reliability or averaging the results.