Агрохимия, 2022, № 3, стр. 81-93

ВВЕДЕНИЕ

Хотя задача анализа данных многолетних полевых опытов стала актуальной с момента закладки подобных опытов (как краткосрочных, так и длительных), более или менее ее полноценного общепринятого решения пока не было представлено. В работе [1] был сделан обзор и обоснование возможных методических подходов к решению некоторых типовых задач, а в настоящей работе эти разработки представлены в виде проекта рекомендаций.

Ограничимся решением следующих частных задач, которые считаются наиболее важными: 1) устойчивость по годам эффектов доз факторов опыта, 2) сходство и различие многолетних динамик изученного показателя, 3) наличие многолетнего тренда изученного показателя, 4) сравнение многолетних средних изученного показателя. Поскольку краткосрочные опыты (<10 лет и/или <2-х ротаций севооборота) имеют небольшую длину временно́го ряда, то для них 3-я и 4-я задачи не могут быть адекватно решены.

В качестве экспериментальных примеров используем здесь те же результаты опытов, что и в работе [1]: 3-летний опыт из работы [2] (табл. 1) и фрагмент опыта с формами калийных удобрений на Люберецком опытном поле (ЛОП) НИУИФ, предоставленный В.Г. Граковским (табл. 2). В качестве изученного показателя в обоих опытах взята урожайность культур, но это могут быть и другие ежегодные характеристики, относящиеся к культурам, почвам, экономическим показателям.

КРАТКОСРОЧНЫЕ ПОЛЕВЫЕ ОПЫТЫ

Устойчивость по годам эффектов доз факторов опыта. Для каждого года оценивается зависимость урожайности от факторов опыта дисперсионным (ДА) и регрессионным (РА) (если возможно) анализами. По результатам ДА определяется варьирование от года к году доли влияния и уровня значимости факторов и их взаимодействий: чем меньше варьирование, тем более устойчиво влияние изученных факторов. По результатам РА можно увидеть и оценить характер зависимости урожайности от факторов опыта – если ежегодные регрессионные зависимости имеют однотипный характер по структуре уравнений, знакам коэффициентов, то можно говорить об устойчивом влиянии доз факторов, хотя количественно величины коэффициентов могут меняться от года к году.

Для данных табл. 1 результаты ДА представлены в табл. 3. Из нее видно, что первые 2 года, как по влиянию изученных факторов, так и по их значимости, были относительно близкими между собой. А на 3-й год резко увеличилось влияние доз азотных удобрений и уменьшилось влияние фунгицида, несмотря на снижение средней урожайности (табл. 1).

РА дал следующие результаты:

Все зависимости однотипны по структуре и по знакам коэффициентов регрессии. Общий вывод: влияние факторов опыта, хотя и менялось количественно от года к году, качественно было устойчивым.

Сходство и различие многолетних динамик урожайности. Общее представление о сходстве и различии временны́х динамик для вариантов и делянок можно получить построением дендрограмм сходства и таблиц попарного сходства. Для расчета сходства можно использовать, например, многомерное евклидово расстояние между объектами (вариантами или делянками) k и l

и меру (степень) сходства c_kl = 100/(100 + d_kl), %, где d_kl – расстояние между рядами динамики объектов k и l, n – число лет, x_ki и x_li – величины урожайности для i-го года объектов k и l. Если изученный показатель – фракционный состав (для почв или растений), то расстояния между объектами можно рассчитать по формуле из работы [3], а если необходимо использовать совокупность разнородных показателей, то их предварительно нормируют (переводят в относительные единицы). Для облегчения расчетов имеются соответствующие программы для компьютера в различных статистических пакетах.

В нашем примере дендрограмма сходства временны́х динамик урожайности для делянок содержит 32 объекта (8 вариантов в 4-х повторениях) (рис. 1). Показано, что даже на уровне сходства 70% (довольно высоком) и выше образующиеся группы делянок, как правило, содержат делянки разных вариантов. Следовательно, временны́е динамики урожайности в данном случае слабо связаны с вариантами опыта и определяются другими факторами (возможно, различием в почвенных свойствах).

Полную таблицу сходств временны́х динамик для всех возможных пар делянок мы не приводим из-за ее громоздкости. Детальное рассмотрение этой таблицы показало, что, используя уровень сходства 80% и выше, получаются 3 группы делянок, сходные динамикой урожайности. К первой группе относятся делянки вариантов 1 (N0Ф0) и 2 (N0Ф1), за исключением делянок 1-2 (вариант 1, 2-е повторение) и 2-1 (вариант 2, 1-е повторение). Ко 2-й группе относятся делянки всех остальных вариантов, за исключением делянок 3‑4 и 5-4; отдельную группу составили делянки 8‑1 и 8-2. Эту картину можно увидеть и на дендрограмме сходства в несколько завуалированном виде.

Дендрограмма сходства временны́х динамик средних величин урожайностей в вариантах показана на рис. 2. На уровне сходства 70% образующиеся группы вариантов сходной динамики (кроме первой с вариантами N0Ф0 и N0Ф1) также неоднородны – в них изменяются дозы обоих факторов. Другими словами, при дозе азота 60 кг/га и выше дозы факторов не определяли сходство временны́х динамик урожайности вариантов.

Полная матрица попарных сходств временны́х динамик средних урожайностей для вариантов показана в табл. 4. Пары вариантов с уровнем сходства 70% и выше получились следующие: 1 и 2 (N0), 5, 6, 7 и примыкающие к ним 4 и 8. Динамика варианта 3 отличалась от других. Это рассмотрение детализирует дендрограмму.

Теперь рассмотрим этот же вопрос сравнения динамик с точки зрения сопряженных выборок. Для этого можно использовать t-критерий [4] и непараметрический критерий Вилкоксона [5]. Эти критерии рассматривают экспериментальный материал с противоположных предыдущему позиций – оцениваются не сходства динамик, а их различия.

Непараметрический критерий в данном примере оказался не применимым, т.к. так 3-х опыто-лет для его применения недостаточно – нужно как минимум 6 лет опыта. Полная таблица при использовании t-критерия для делянок тоже слишком громоздка, поэтому в табл. 5 показаны только значимо отличающиеся временны́е динамики. Из этой таблицы видно следующее. Варианты 1, 2 и 5 имеют делянки, значимо различающиеся по временно́й динамике (делянки 1 и 4 в варианте 1, делянки 6 и 7 в варианте 2, делянки 17 и 19 в варианте 5). Делянки вариантов 7 и 8 не различаются значимо динамикой урожайности; аналогичная ситуация для пар вариантов 1 и 2, 1 и 7, 2 и 3, 2 и 5, 2 и 6, 2 и 7, 4 и 7, 5 и 6.

Если сопряженно по t-критерию сравнивать динамики средних для вариантов, то получим значимое различие для вариантов 6 и 7 (0.01, средняя разница +2.7) и вариантов 3 и 8 (0.05, средняя разница +8.0); это не противоречит данным табл. 5 и рис. 2.

Сравнение методами сходства и различия могут давать и несовпадающие результаты (разные алгоритмы с разными предпосылками), в таком случае выводы должны быть более осторожными.

ДЛИТЕЛЬНЫЕ ПОЛЕВЫЕ ОПЫТЫ

Устойчивость по годам эффектов доз факторов опыта. Подходы здесь те же, что для краткосрочных опытов. Результаты для использованного примера представлены в табл. 6. Регрессионный анализ проводить в этом случае нет особого смысла из-за наличия только 2-х доз удобрений и малого числа вариантов. Так как размещение делянок по организованным повторениям не известно, то оценить вклад организованных повторений не было возможности, и разброс урожайностей в зависимости от делянок соответствовал случайной ошибке.

Данные табл. 6 показали, что влияние вариантов опыта имело тенденцию к изменению во времени: наблюдали рост доли влияния и числа значимых лет влияния от первых 3-х ротаций ко 2-й тройке ротаций. Имелась также четкая тенденция к изменению варианта с наибольшей урожайностью: в первые 3 ротации это могли быть варианты как NP, так и NP + K_x, в последующие годы – только вариант NP + K_x.

Сходство и различие многолетних динамик урожайности. Дендрограмма сходства делянок по многолетним динамикам урожайности для опыта ЛОП показана на рис. 3. В этом случае уровни сходства значительно меньше, чем для опыта из работы [2]. На уровне сходства 60% в каждом варианте лишь по 2 делянки имели сходные динамики урожайности (5 и 29 в контроле, 6 и 13 в варианте NP, 18 и 25 в варианте NP + K_x). Остальные делянки, хотя и группируются с вышеуказанными внутри вариантов (кроме делянки 32), но на еще более низких уровнях сходства. Тем не менее, в отличие от опыта из работы [2], в этом случае варианты более существенно влияли на динамику урожайности (группировка делянок по вариантам в основном выполнялась). Полная матрица парных сходств делянок по динамике урожайностей показана в табл. 7. В данном случае она полностью подтверждает результаты дендрограммы.

Для средних урожайностей в вариантах вообще не получено существенных сходств многолетних динамик (табл. 8). Все же можно отметить, что крайние варианты (контроль и NP + K_x) менее сходны между собой, чем соседние.

Теперь, как и ранее, рассмотрим сравнение многолетних динамик урожайностей с точки зрения сопряженных выборок. Так как в этом случае усредняются и сопоставляются разности урожайностей разных культур, то возникает необходимость перейти к более сопоставимым единицам, например, к зерновым [1]. Если использование t-критерия для сопряженных выборок хорошо известно, то использование непараметрического критерия Вилкоксона следует подробно изложить. Необходимость использовать 2 критерия (параметрический и непараметрический) обусловлена тем, что они имеют разные предпосылки для своего применения, а это повышает надежность результатов, если они совпадают по обоим критериям.

Использование критерия Вилкоксона для сопряженных выборок излагаем по работе [5]. Для конкретности возьмем сравнение делянок 5 и 12 (табл. 9). Сначала считаем разности для каждого года, затем эти разности ранжируем по абсолютной величине от самых маленьких до самых больших числами натурального ряда (1, 2, …). Далее считаем суммы рангов для положительных и отрицательных разностей раздельно. Меньшая сумма рангов представляет собой статистику Вилкоксона для сопряженных пар (выборок). В нашем примере меньшая сумма рангов равна 61. Далее используем таблицу критических величин критерия (табл. 10) [5]. При n = 23 (длина временно́го ряда) величина 61 меньше критической при уровне значимости 0.05 (74), но больше 56 (уровень значимости 0.01), следовательно делянки 5 и 12 различаются по временно́й динамике урожайности на уровне значимости 0.05, но не различаются при уровне значимости 0.01. Если какие-то разности в сопряженных рядах равны нулю, то соответствующие пары исключаются из рассмотрения, и соответственно уменьшается объем выборки.

Итоги сравнения многолетних динамик урожайностей культур для делянок и для вариантов сведены в табл. 11 и 12. Результаты оценки различий динамик для делянок не противоречат картине, полученной при изучении уровней сходства: дополнительно отмечено значительно больше неразличающихся значимо пар делянок. При сравнении динамик вариантов получена однозначная картина: все варианты значимо различаются по обоим критериям, и это соответствует низкому сходству между ними, полученному выше.

Наличие многолетнего тренда урожайности. Сама постановка такой задачи возможна только при изучении длительных многолетних опытов, где эти тренды могут проявиться. Однако и в этом случае задача достаточно неопределенна. В реальных временны́х рядах мы имеем дело с наложением нескольких источников варьирования: природных длиннопериодных циклов, смены культур в севооборотах, смены сортов культур, технологий обработки почвы, агроистощения или повышения обеспеченности растений питательными веществами, смены систем защиты растений, почвоутомления и т.п. На все это накладываются ежегодные изменения как природного, так и антропогенного характера [6].

При анализе длительных опытов обычно не стои́т задача выявить все детали имеющихся изменений, но определить основной тренд и, по возможности, наиболее выраженные циклы необходимо. Эти характеристики важны для сравнения между собой различных длительных опытов, различных вариантов опытов, а их сопоставление с природными процессами и антропогенными воздействиями будет способствовать лучшему пониманию причин получаемой в опытах динамики. Кроме того, как показано ниже, учет наличия временно́го тренда необходим для обоснованной оценки ошибки многолетних средних.

Для временны́х рядов урожайности показатели динамики могут сильно искажаться в зависимости от уровней первого и последнего члена ряда [7, c. 11]. Характерная картина динамики видна на данных ЛОП: динамика урожайностей делянок 2-х крайних вариантов (рис. 4) и динамика средних урожайностей для всех трех вариантов (рис. 5). В обоих случаях можно лишь качественно говорить, что в вариантах с удобрениями по годам с максимальной урожайностью наблюдается параболический тренд (повышение и понижение), а в годы с минимальной урожайностью и в  контроле – другая закономерность – чуть ли не снижение урожайности со временем.

С другой стороны, четко видна одинаковая цикличность урожайности как для делянок, так и для вариантов. Сроки между максимумами урожайности составили 2–6 лет, между минимумами урожайности – 2.5–5 лет, а в среднем длина цикла составила ≈4 года. В данном примере это совпало с длительностью ротации, но может быть случайным совпадением. Максимальным урожайностям в основном соответствовали пропашные культуры, но встречались также рожь озимая и люпин. Минимальным урожайностям в основном соответствовал люпин, но встречались также рожь озимая, гречиха и пропашные. Такое разнообразие свидетельствует, вероятно, о влиянии погодных условий разных лет.

Более длительные циклы в данном примере таким же графическим образом оценить невозможно из-за недостаточно большой длины временно́го ряда.

Сравнение многолетних средних для делянок и вариантов опыта. В этом вопросе, как и в предыдущем разделе, желаемые решения носят приближенный характер, т.к. реальные многолетние ряды урожайности могут не всегда соответствовать теоретическим допущениям.

Чтобы сравнивать многолетние средние, надо оценить их ошибки. Для конкретной делянки ошибка связана непосредственно с временны́м изменением урожайности, а для варианта опыта с делянками-повторениями к временно́й ошибке добавляется пространственная. Временну́ю ошибку можно оценить для стационарных рядов; приблизить реальный ряд к стационарности можно, убрав временно́й тренд. В нашей задаче речь может идти о том, чтобы убрать основную (более или менее гладкую) временну́ю тенденцию, а критерием успеха должна служить не значимость самого тренда (если он значим, то это еще лучше), а заметное снижение остаточной дисперсии по сравнению с исходной (без снятия тренда).

Квадратичные модели трендов для отдельных делянок представлены в табл. 13. Все модели оказались незначимы, что не удивительно из-за сильных изменений урожайности от года к году. В то же время стандартные отклонения остатков (s_ост) были меньше исходных (s₀), а средние из разниц между рассчитанными по моделям урожайностям и измеренными в опыте достоверно не отличались от нуля (табл. 14).

Если принять полученные ряды разниц стационарными, то для них имеются теоретические формулы для оценки ошибки среднего в ряду ([8], с. 136), которые можно было бы принять для наших целей. Эти оценки считают с использованием коэффициентов автокорреляции между членами ряда. Однако в нашем примере рассчитанные таким образом ошибки средних оказываются часто либо близкими к нулю, либо отрицательными, что не имеет физического смысла. Можно предположить, что ряды либо недостаточно длинные, либо они все же не являются стационарными. Поэтому в качестве ошибки средних многолетних для делянок предлагается брать ошибки средних из табл. 14, либо из остаточного стандартного отклонения ($\sqrt {{{s}_{{{\text{ост}}}}}} $) табл. 13. Эти параметры несколько различаются из-за ошибок пересчетов.

Еще одна возможность оценить ошибку непосредственно по ряду разниц – непараметрический подход, один его вариант – через квартили [5]:

А при наличии подходящей компьютерной программы (например, программа KVANT [9]) можно рассчитать величины x, соответствующие квантилям эмпирической функции распределения 0.159 и 0.841, тогда

Эти непараметрические подходы относятся к распределениям, которые не очень сильно отличаются от нормального.

Все полученные оценки стандартного отклонения сведены в табл. 15. Можно видеть, что в данном примере из оценок s по 4-м способам несколько выпадают оценки, рассчитанные через квартили. Поэтому в качестве окончательной оценки s берем среднее из 3-х других способов.

Для перехода от ошибок для делянок к ошибкам для вариантов опыта предполагаем, что делянки каждого варианта (повторности) принадлежат к одной генеральной совокупности. Тогда ([10], с. 105–106)

где Y_вар – средняя многолетняя урожайность для варианта опыта, N – сумма лет опыта по всем делянкам (повторностям) варианта, k – число делянок (повторностей) в варианте, n_i – число лет для i-й делянки, Y_срi – средняя многолетняя урожайность для i-й делянки, $S_{{{\text{вар}}}}^{2}$ – дисперсия урожайности для варианта, $s_{i}^{2}$ – дисперсия урожайности для i-й делянки. Степень свободы для дисперсии варианта – N – 1. Результаты расчетов статистических показателей для вариантов также показаны в табл. 15.

Что и как можно сравнивать в результате? При сравнении многолетних средних арифметических используем известную формулу для t-критерия в общем случае, когда дисперсии средних различаются:

где M_i – сравниваемые средние, m_i – ошибки этих средних. Сложность здесь состоит в расчете числа степеней свободы для t-критерия: где n_i – объемы сравниваемых выборок. Если объемы выборок одинаковы и равны n, то

Сравнивать можно средние многолетние показатели для делянок и вариантов, используя данные табл. 14 и 15. Для примера сравним средние многолетние урожайности для вариантов:

А – Контроль (NP)

уровень значимости получился немного >0.05.

Б – Аналогично для вариантов Контроль и NP + K_x получаем t = 3.665, u = 0.317, ν = 38.8, уровень значимости <0.001.

В – Варианты NP и NP + K_x: t = 1.628, u = = 0.452, ν = 43.6; уровень значимости >0.10.

Может представить интерес сравнение средних не для всего срока наблюдения, а для отдельных его периодов (например, в начале, в середине, в конце имеющихся наблюдений). Так как имеются временны́е тренды урожайности, то возникает вопрос, какую надо взять длительность периода, по которой будут рассчитываться среднее и его ошибка? При этом среднее считается по модели тренда, а ошибка среднего как s/$\sqrt n $, где n – длительность выбранного периода, s – из табл. 15 [7]. Заранее предугадать длительность периода затруднительно, поэтому можно провести следующую процедуру: беря длительность периода в 2, 3, 4 и так далее лет, рассчитываем для них средние (по тренду) и ошибки средних. Сравнивая средние (для делянок или вариантов) по t-критерию, находим такую длину периода, где различия средних максимальны и возможно значимы.

Автор с благодарностью примет замечания к этому проекту.