ERG final
2024-09-26
Load ggplot2 library
library(ggplot2)A-wave
Read data
steps1to6 <- read.table(file = "ERG Steps 1-6.csv", header = TRUE, sep = ",")Subset data for a-wave
steps1to6_a <- steps1to6[steps1to6$Name == "a", ]Take average of left and right eye
steps1to6_a_average <- aggregate(uV_average ~ PD + log_intensity, data = steps1to6_a, FUN = mean)Re-arrange table to sort by PD
steps1to6_a_average$PD <- factor(steps1to6_a_average$PD, levels = c("Control ", "PD5", "PD7", "PD14", "PD21"))
steps1to6_a_average_sorted <- steps1to6_a_average[order(steps1to6_a_average$PD), ]Add SD column
steps1to6_a_average_sorted$SD <- c(10.06, 10.77, 23.08, 52.43, 68.63, 57.48, 16.29, 18.26, 12.59, 15.92, 20.49, 30.54, 9.26, 12.31, 26.61, 59.25, 70.02, 71.93, 12.37, 14.60, 20.37, 45.11, 50.91, 50.22, 9.02, 16.10, 23.92, 42.13, 39.56, 36.08)Plot a-wave with error bars
ggplot(data = steps1to6_a_average_sorted, aes(x = log_intensity, y = abs(uV_average), colour = PD)) +
geom_point(size = 2) +
geom_line() +
scale_x_continuous(breaks = seq(-2.5, 2, by = 0.5)) +
scale_y_continuous(breaks = seq(0, 300, by = 50)) +
labs(title = "A-wave ERG", x = "Flash Intensity (Log cd.s/m2)", y = "Average Amplitude (µV)") +
geom_errorbar(aes(x = log_intensity, ymin = abs(uV_average) - SD, ymax = abs(uV_average) + SD, width = 0.05)) +
theme_light() +
theme(legend.title = element_blank())
Plot a-wave without error bars
ggplot(data = steps1to6_a_average_sorted, aes(x = log_intensity, y = abs(uV_average), colour = PD)) +
geom_point(size = 2) +
geom_line() +
scale_x_continuous(breaks = seq(-2.5, 2, by = 0.5)) +
scale_y_continuous(breaks = seq(0, 300, by = 50)) +
labs(title = "A-wave ERG", x = "Flash Intensity (Log cd.s/m2)", y = "Average Amplitude (µV)") +
theme_light() +
theme(legend.title = element_blank())
Fitting a linear model
steps1to6_a_average_sorted_lm_model <- lm(uV_average ~ log_intensity * PD, data = steps1to6_a_average_sorted)
summary(steps1to6_a_average_sorted_lm_model)##
## Call:
## lm(formula = uV_average ~ log_intensity * PD, data = steps1to6_a_average_sorted)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.737 -7.674 2.853 10.141 27.936
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -152.173 9.367 -16.246 5.47e-13 ***
## log_intensity -66.208 6.959 -9.514 7.26e-09 ***
## PDPD5 98.034 13.246 7.401 3.80e-07 ***
## PDPD7 32.776 13.246 2.474 0.022429 *
## PDPD14 39.104 13.246 2.952 0.007881 **
## PDPD21 34.708 13.246 2.620 0.016395 *
## log_intensity:PDPD5 42.626 9.841 4.331 0.000324 ***
## log_intensity:PDPD7 14.596 9.841 1.483 0.153637
## log_intensity:PDPD14 19.661 9.841 1.998 0.059526 .
## log_intensity:PDPD21 16.455 9.841 1.672 0.110093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.29 on 20 degrees of freedom
## Multiple R-squared: 0.9371, Adjusted R-squared: 0.9088
## F-statistic: 33.13 on 9 and 20 DF, p-value: 4.291e-10library(emmeans)## Warning: package 'emmeans' was built under R version 4.3.3## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'em = emmeans(steps1to6_a_average_sorted_lm_model, "PD")## NOTE: Results may be misleading due to involvement in interactionspairs(em)## contrast estimate SE df t.ratio p.value
## Control - PD5 -84.39 12.9 20 -6.559 <.0001
## Control - PD7 -28.11 12.9 20 -2.184 0.2257
## Control - PD14 -32.81 12.9 20 -2.550 0.1189
## Control - PD21 -29.44 12.9 20 -2.288 0.1897
## PD5 - PD7 56.29 12.9 20 4.375 0.0024
## PD5 - PD14 51.58 12.9 20 4.009 0.0055
## PD5 - PD21 54.95 12.9 20 4.271 0.0031
## PD7 - PD14 -4.71 12.9 20 -0.366 0.9959
## PD7 - PD21 -1.34 12.9 20 -0.104 1.0000
## PD14 - PD21 3.37 12.9 20 0.262 0.9989
##
## P value adjustment: tukey method for comparing a family of 5 estimatesInterpretation:
The interaction term tests whether the slope between flash intensity and amplitude is significantly different for each PD group compared to the control.
PD5 has a significantly different rate of change (slope) compared to the control. For every unit increase in flash intensity, the amplitude in the PD5 group increases by 42.6 units less than the control group.
The other groups (PD7, PD14, PD21) show some differences, but they are not statistically significant. ________________________________________________________________________________
B-wave
Subset data for b-wave
steps1to6_b <- steps1to6[steps1to6$Name == "B", ]Take average of left and right eye
steps1to6_b_average <- aggregate(uV_average ~ PD + log_intensity, data = steps1to6_b, FUN = mean)Re-arrange table to sort by PD
steps1to6_b_average$PD <- factor(steps1to6_b_average$PD, levels = c("Control ", "PD5", "PD7", "PD14", "PD21"))
steps1to6_b_average_sorted <- steps1to6_b_average[order(steps1to6_b_average$PD), ]Add SD column
steps1to6_b_average_sorted$SD <- c(86.71, 85.11, 92.45, 124.15, 134.70, 120.61, 41.45, 42.02, 47.20, 56.51, 62.87, 73.80, 102.15, 97.66, 101.27, 131.22, 148.76, 146.31, 78.46, 74.58, 77.72, 99.71, 108.01, 107.43, 34.24, 32.97, 43.34, 54.45, 54.22, 48.73)Plot b-wave with error bars
ggplot(data = steps1to6_b_average_sorted, aes(x = log_intensity, y = abs(uV_average), colour = PD)) +
geom_point(size = 2) +
geom_line() +
scale_x_continuous(breaks = seq(-2.5, 2, by = 0.5)) +
scale_y_continuous(breaks = seq(0, 600, by = 50)) +
labs(title = "B-wave ERG", x = "Flash Intensity (Log cd.s/m2)", y = "Average Amplitude (µV)") +
geom_errorbar(aes(x = log_intensity, ymin = abs(uV_average) - SD, ymax = abs(uV_average) + SD, width = 0.05)) +
theme_light() +
theme(legend.title = element_blank())
Plot b-wave without error bars
ggplot(data = steps1to6_b_average_sorted, aes(x = log_intensity, y = abs(uV_average), colour = PD)) +
geom_point(size = 2) +
geom_line() +
scale_x_continuous(breaks = seq(-2.5, 2, by = 0.5)) +
scale_y_continuous(breaks = seq(0, 600, by = 50)) +
labs(title = "B-wave ERG", x = "Flash Intensity (Log cd.s/m2)", y = "Average Amplitude (µV)") +
theme_light() +
theme(legend.title = element_blank())
Fitting a linear model
steps1to6_b_average_sorted_lm_model <- lm(uV_average ~ log_intensity * PD, data = steps1to6_b_average_sorted)
summary(steps1to6_b_average_sorted_lm_model)##
## Call:
## lm(formula = uV_average ~ log_intensity * PD, data = steps1to6_b_average_sorted)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.52 -13.38 -1.78 13.40 32.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 306.798 8.394 36.551 < 2e-16 ***
## log_intensity 43.619 6.236 6.995 8.69e-07 ***
## PDPD5 -93.750 11.870 -7.898 1.42e-07 ***
## PDPD7 -5.440 11.870 -0.458 0.65169
## PDPD14 -23.453 11.870 -1.976 0.06215 .
## PDPD21 -36.126 11.870 -3.043 0.00642 **
## log_intensity:PDPD5 -19.566 8.819 -2.219 0.03825 *
## log_intensity:PDPD7 -4.369 8.819 -0.495 0.62568
## log_intensity:PDPD14 -17.268 8.819 -1.958 0.06433 .
## log_intensity:PDPD21 -14.210 8.819 -1.611 0.12279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.97 on 20 degrees of freedom
## Multiple R-squared: 0.9163, Adjusted R-squared: 0.8786
## F-statistic: 24.32 on 9 and 20 DF, p-value: 7.014e-09Interpretation:
The interaction term tests whether the slope between flash intensity and amplitude is significantly different for each PD group compared to the control.
PD5 has a significantly different rate of change (slope) compared to the control. For every unit increase in flash intensity, the amplitude in the PD5 group increases by 19.5 units less than the control group.
The other groups (PD7, PD14, PD21) show some differences, but they are not statistically significant.
C-Wave
Read data
step7 <- read.table(file = "ERG Step 7.csv", header = TRUE, sep = ",")Re-arrange table to sort by PD
unique(step7$PD) #check the way the control has been spelt with spaces## [1] "Control " "PD5" "PD7" "PD14" "PD21"step7$PD <- factor(step7$PD, levels = c("Control ", "PD5", "PD7", "PD14", "PD21"))Plot C-wave ERG
ggplot(step7, aes(x = PD, y = uV, fill = PD)) +
geom_boxplot() +
geom_jitter(width = 0.15) +
labs(title = "C-Wave ERG", x = "Mice", y = "C-Wave Amplitude (µV)") +
scale_y_continuous(breaks = seq(0, 400, by = 50)) +
theme_light() +
theme(legend.title = element_blank())
Testing for normality
shapiro.test(step7$uV)##
## Shapiro-Wilk normality test
##
## data: step7$uV
## W = 0.95117, p-value = 0.03529Interpretation:
Null Hypothesis: The data is normally distributed.
Since the p-value is below 0.05, the test indicates that the amplitude values do not follow a normal distribution.
Using Kruskal-Wallis test to check whether the median amplitude values differ between the groups.
kruskal_test <- kruskal.test(uV ~ PD, data = step7)
print(kruskal_test)##
## Kruskal-Wallis rank sum test
##
## data: uV by PD
## Kruskal-Wallis chi-squared = 4.905, df = 4, p-value = 0.2972Interpretation
The Kruskal-Wallis test assesses whether there are statistically significant differences in the median amplitude values across all groups (Control, PD5, PD7, PD14, PD21).
The null hypothesis of the Kruskal-Wallis test is that there is no significant difference in the distribution of amplitude values between the groups.
Given that the p-value is 0.2972 (which is greater than 0.05), we fail to reject the null hypothesis. This indicates that there is no significant evidence to suggest differences in the median amplitude values between any of the groups, including comparisons between the control and each PD group.