Федосеева Е. В., Терехова В. А., Цесаренко О. В., Гладкова М. М. Обработка результатов токсикологических исследований в статистической программе R // Принципы экологии. 2015. № 3. С. 12–26. DOI: 10.15393/j1.art.2015.4381 # Блок 1. ИНИЦИАЛИЗАЦИЯ ДАННЫХ # Создаем вектор с данными об уровне флуоресценции хлорофилла fluor <- c(0.0114, 0.0109, 0.0109, 0.0076, 0.0075, 0.0079, 0.0092, 0.0094,0.0094, 0.0067, 0.0067, 0.0068, 0.0075, 0.0071, 0.0075, 0.0057, 0.0060, 0.0056, 0.0121, 0.0119, 0.0121, 0.0117, 0.0113, 0.0115, 0.0102, 0.0105, 0.0107, 0.0086, 0.0085, 0.0084, 0.0100, 0.0095, 0.0095, 0.0095, 0.0104, 0.0105) # Превращаем вектор в матрицу fluorescence <- matrix(fluor, ncol=4, dimnames = list(1:9,c("контроль", "исходная", "1к50", "1к100"))) # Создаем вектор с данными о численности клеток водорослей quant <- c(880000, 930000, 630000, 530000, 620000, 720000, 270000, 230000, 380000, 520000, 310000, 250000, 1050000, 1180000, 610000, 830000, 970000, 950000, 940000, 1030000, 810000, 1090000, 1040000, 840000) # Превращаем вектор в матрицу quantity <- matrix(quant, ncol=4, dimnames = list(1:6,c("контроль", "исходная", "1к50", "1к100" ))) # Создаем два вектора с группировочными факторами и изменяем их уровни trt<-as.factor( rep(c("контроль", "исходная", "1к50", "1к100"), c(9,9,9,9))) trt.fluor=relevel(trt, ref="контроль") # для данных об уровне флуоресценции # Создаем фактор и изменяем его уровни trt= as.factor( rep(c("контроль", "исходная", "1к50", "1к100"), c(6,6,6,6))) trt.quant=relevel(trt, ref="контроль") # для данных о численности клеток # ------------------------------------------ # Блок 2. ОЦЕНКА ВЫБОРОЧНЫХ ПАРАМЕТРОВ # Загружаем пакеты, необходимые для расчетов library(car); library(pastecs); library(multcomp); library(ggplot2) # Если пакет отсутствует в среде, используйте install.packages(<имя пакета>) # Собственная функция расчета выборочных характеристик по столбцам матрицы StatMat = function(x) { SEmpg = function(y) sd(y)/sqrt(length(y)) # Функция ошибки среднего # Рассчитываем среднее значение (mean), дисперсию (var) и стандартное отклонение (sd) data.frame(mean=apply (x, 2, mean),var=apply (x, 2, var), sd=apply (x, 2, sd), SEm=apply (x, 2, SEmpg)) } # Анализируем данные об уровне флуоресценции хлорофилла stat.fluor = StatMat(fluorescence) (stat.desc(fluorescence)) # Анализируем данные о численности клеток водорослей stat.quant = StatMat(quantity) (stat.desc(quantity)) # ------------------------------------------ # Блок 3. ОЦЕНКА ТОКСИЧНОСТИ # Создаем матрицу процентов прироста/ингибирования уровня флуоресценции и численности клеток gentoxicity = t(data.frame(100*stat.fluor[,1]/stat.fluor[1,1],100*stat.quant[,1]/stat.quant[1,1])) dimnames(gentoxicity) = list(c("уровень флуоресценции","численность клеток"), c("контроль", "исходная", "1/50", "1/100")) # Рисуем столбчатую диаграмму barplot(gentoxicity, beside = TRUE, col = c("gray30", "gray80"), ylab = "% от контроля", xlab = "Наименование проб", family="serif", font=1, cex.lab=1.5, cex.axis=1.5, ylim = c(0, 160), legend.text = rownames(gentoxicity)) # Добавляем линию острой токсичности lines(rep("50", 12), lwd = 4, lty = 5, col="black") # ------------------------------------------ # Блок 4. ДИСПЕРСИОННЫЙ АНАЛИЗ # Проверка предположения о нормальности распределения # Тест Шапиро-Уилка для данных по уровню флуоресценции хлорофилла shapiro.test(fluorescence[,1]); shapiro.test(fluorescence[,2]) shapiro.test(fluorescence[,3]); shapiro.test(fluorescence[,4]) # Тест Шапиро-Уилка для данных по численности клеток водорослей shapiro.test(quantity[,1]); shapiro.test(quantity[,2]) shapiro.test(quantity[,3]); shapiro.test(quantity[,4]) # Проверяем гомогенность дисперсий в группах по тесту Левина leveneTest(fluor, trt.fluor, center = mean) leveneTest(quant, trt.quant, center = mean) # Выполняем дисперсионный анализ для данных по уровню флуоресценции Mod.fluor <- aov(fluor~trt.fluor); summary(Mod.fluor) # То же для данных по численности клеток водорослей Mod.quant <- aov(quant~trt.quant); summary(Mod.quant) # ------------------------------------------ # Блок 5. МНОЖЕСТВЕННЫЕ СРАВНЕНИЯ # Проверка достоверно значимой разности HSD Тьюки (пакет car) TukeyHSD(Mod.fluor); TukeyHSD(Mod.quant) # Рисуем графики с доверительными интервалами Тьюки # Чтобы увеличить левое поле подписи, корректируем графические параметры par(mar = c(5.1, 12.1, 4.1, 2.1)) plot(TukeyHSD(Mod.fluor), family="serif", font=1, cex.axis=1.5, las = 1) par(mar = c(5.1, 12.1, 4.1, 2.1)) plot(TukeyHSD(Mod.quant), family="serif", font=1, cex.axis=1.5, las = 1) # Используем для множественных сравнений пакет multcomp # Проверку осуществляем с использованием матриц контрастов Тьюки и Даннета summary(glht(Mod.fluor, linfct = mcp(trt.fluor = "Tukey"))) summary(glht(Mod.fluor, linfct = mcp(trt.fluor = "Dunnett"))) summary(glht(Mod.quant, linfct = mcp(trt.quant = "Tukey"))) summary(glht(Mod.quant, linfct = mcp(trt.quant = "Dunnett"))) # Рисуем графики с доверительными интервалами Даннета par(mar = c(5.1, 12.1, 4.1, 2.1)) plot(confint(glht(Mod.fluor, linfct = mcp(trt.fluor = "Dunnett")) , level=0.95), family="serif", font=1, cex.axis=1.5) par(mar = c(5.1, 12.1, 4.1, 2.1)) plot(confint(glht(Mod.quant, linfct = mcp(trt.quant = "Dunnett")) , level=0.95), family="serif", font=1, cex.axis=1.5) # ------------------------------------------ # Блок 6. ОЦЕНКА КОРРЕЛЯЦИИ НАБЛЮДАЕМЫХ ПЕРЕМЕННЫХ # Создаем вектор со средними значениями уровня флуоресценции для каждой повторности проб cor.fluor<-c(0.01107, 0.00767, 0.00933, 0.00673, 0.00737, 0.00577, 0.01203, 0.01150, 0.00741, 0.00850, 0.00967, 0.01013) # Создаем вектор со средними значениями численности клеток для каждой повторности проб cor.quant<-c(905000, 580000, 670000, 250000, 450000, 280000, 1115000, 720000, 960000, 985000, 950000, 940000) # Рассчитываем коэффициент корреляции Пирсона cor.test(cor.fluor, cor.quant) # Рассчитываем коэффициент корреляции Спирмена cor.test(cor.fluor, cor.quant, method = "spearman") # Рисуем график зависимости ggplot(,aes(x=cor.quant, y=cor.fluor)) + xlab("Численность клеток, тыс.", family="serif", font=1) + ylab("уровень флуоресценции", family="serif", font=1) + geom_point(shape=16) + # Точки наблюдений "ромбиками" geom_smooth(method=lm,se=FALSE, size = 2 ) + # Добавить линейную регрессию geom_smooth(colour = "red", linetype = 2, size = 2) + # Добавить кривую сглаживания theme_bw()