РЕГРЕССИЯ В ПАССИВНОМ ЭКСПЕРИМЕНТЕ:
Some uncertainties related to passive experiment designing is discussed.
Математическая статистика изначально базировалась на четких предпосылках, соответствующих весьма идеализированным ситуациям. Практическое экспериментирование привело со временем к размыванию жесткости предпосылок и откровениям типа: «Нормальность превратилась не более чем в частный случай» [1; с. 12] либо вовсе: «Нормальный закон как закон ошибок неверен» [15; с. 18). Поэтому экспериментатору приходится переходить с «гладкой дороги нереальных предпосылок, произвольных критериев и абстрактных результатов» ….«на каменистый путь реальных проблем» [12; с. 68].
Обсуждая возможные варианты обработки результатов измерений, В. Н. Тутубалин с соавт. [15] считают наиболее приемлемым метод наименьших квадратов, составляющий основу регрессионного анализа. Хотя основы матричной алгебры, по алгоритмам которой рассчитываются коэффициенты многофакторной регрессии, были заложены еще в XVIII веке, обработка результатов измерений долгое время сводилась к установлению парных зависимостей. Со времен И. Ньютона в течение 200 лет точные науки имели дело лишь с хорошо организованными системами, когда результаты исследований описывались легко интерпретируемыми функциональными связями, которым приписывалась роль абсолютных законов, а методология однофакторного эксперимента считалась единственно правильной [12].
Преимущества многофакторного эксперимента в условиях плохо организованных (диффузных) систем были показаны Р. Фишером в 1920-е годы [19], а в 1930-е годы А. К. Митропольским [10] предложены алгоритмы расчета многофакторных регрессий. Однако подвижников, которые бы отважились на подобные расчеты с помощью арифмометра, долгое время не находилось. Наступление 1950-60-х годов ознаменовалось в экологии «мистикой электронных вычислительных машин» [15; с. 201], «регрессионным бумом» [1; с. 7] и развитием планирования экспериментов [14]. Регрессионный анализ стали считать «методом века» [2].
Как известно, планируемый эксперимент может быть активным и пассивным, но в обоих случаях основное требование – воспроизводимость его результатов. При активном эксперименте объект исследования управляется путем задания действующим на него факторам определенных значений, когда все уровни некоторого фактора комбинируются со всеми уровнями остальных факторов. В лесоводстве исследователь имеет дело с пассивным экспериментом, т.е. «пассивно наблюдает за тем, как эксперимент ведет природа» [12; с. 161]. При этом все факторы не только не управляемы, но и тесно коррелированы с неучтенными факторами и между собой, что ведет к смещению оценок. Последнее обстоятельство порождает некоторые потенциально опасные ситуации.
Первая ситуация. Ошибка в модели или плохая воспроизводимость эксперимента не является случайной, а обусловлена действием скрытых (латентных) или дрейфующих факторов, не регистрируемых в ходе эксперимента. В этом случае возможно привлечение специальных методов планирования [8].
Вторая ситуация. Действующие факторы при эксперименте измеряются в искусственно суженных диапазонах, что приводит к занижению коэффициента регрессии. С увеличением диапазона варьирования фактора абсолютная величина коэффициента регрессии возрастает [3]. Кроме того, при сужении естественных диапазонов два фактора (в принципе ортогональных при условии, если эксперимент спланирован) становятся тесно коррелированными, и один из коэффициентов регрессии при таких факторах оказывается статистически не значимым. По-видимому, по этой причине у В. В. Кузьмичева [6] при расчете зависимости среднего диаметра древостоя (Dср) от трех массообразующих показателей: возраста (А), средней высоты (Нср) и среднего расстояния между деревьями (величины, обратной густоте N) – возраст А оказался коррелированным с Нср и отсеялся как статистически не значимый.
Третья ситуация проявляется при выборе структуры модели. Поскольку в реальном эксперименте и соответствующей регрессии число статистически значимых факторов ограничено (обычно не более 5-6), а в реальном насаждении можно измерить гораздо большее их количество и бóльшая их часть коррелирована по определению (объективно), то возникает проблема отбора наиболее информативных из них. Имеющиеся алгоритмы отбора (планы отсеивания) [12, 14] требуют специальной подготовки и программного обеспечения, совмещенного со стандартными программами обработки данных, и поэтому не нашли пока применения у лесоводов.
Метод наименьших квадратов предполагает наличие линейности исследуемой зависимости, тогда как в условиях лесной экосистемы более обычны нелинейные связи, среди которых линейная – частный случай. Поэтому прибегают к приемам линеаризации, в частности, логарифмическому преобразованию. Некорректная линеаризация аргумента при построении структуры модели создает четвертую ситуацию, чреватую невоспроизводимостью эксперимента.
Но еще более сложной представляется пятая ситуация, связанная с учетом одновременно действующих факторов - синергизмов. А. М. Мауринь с соавт. [9] путем ввода в трехфакторную линейную модель прогноза плодоношения тсуги канадской произведений и квадратов предикторов повысили коэффициент детерминации R2 с 0,476 до 0,790. Однако достаточно надежного математического аппарата для анализа подобных взаимовлияний пока не разработано.
Наконец, шестая ситуация, опасная в аспекте воспроизводимости, связана с тем, что с целью сокращения числа действующих факторов некоторые из них стремятся объединить в один, включающий в себя в сжатой форме информацию об основных воздействиях на объект. Например, уравнение Спурра, широко применяемое в лесной таксации для оценки объема ствола, включает в себя в качестве предиктора так называемый видовой цилиндр D2H, или произведение квадрата диаметра на высоту ствола, т.е. основные массообразующие показатели ствола. Слабая изменчивость полнодревесности и учет ее параметрами данной зависимости обусловил высокую адекватность уравнения при оценке объема и массы ствола: при одном и том же значении D2 масса ствола увеличивается прямо пропорционально его высоте. По молчаливому согласию ряд зарубежных и отечественных исследователей стали повсеместно применять этот агрегированный предиктор и при оценке массы кроны, хотя каких-либо убедительных предпосылок для такого переноса нет: при одном и том же D2 масса кроны дерева с увеличением его высоты не возрастает подобно массе ствола, а напротив, снижается. В совокупности древостоев разного возраста это объясняется сдвигом рангового положения дерева одного и того же диаметра с возрастом древостоя, а в совокупности древостоев разной производительности – бóльшим развитием ассимиляционного аппарата с ухудшением условий произрастания при одинаковой толщине деревьев.
В условиях реального эксперимента некорректность структуры уравнения с видовым цилиндром в качестве предиктора при оценке массы кроны проявляется в меньшей его детерминированности по сравнению с уравнением, имеющим в качестве предиктора только диаметр D. Расчет таких уравнений по совокупности 320 модельных деревьев в сосняках, взятых в диапазоне пяти классов возраста и пяти классов бонитета дал значения R2 для массы хвои 0,669 – по D2H и 0,758 – по D, для массы ветвей – соответственно 0,810 и 0,870. Подобная закономерность проявилась и при оценке массы крон по совокупности 104 деревьев, сплошь вырубленных в 45-летнем березняке: R2 для массы листвы составил соответственно 0,925 и 0,939, для массы ветвей – 0,899 и 0,924 [16].
Аналогичная ситуация имеет место и на уровне древостоя: если при оценке запаса стволовой древесины (М, м3/га) использование полноты (0,785 D2срN) в формуле запаса (М = 0,785 D2срНсрNF, где Dср и Нср –соответственно средние диаметр и высота; N – число стволов на 1 га и F – видовое число) логически обосновано, то при оценке фитомассы крон (полога) это далеко не очевидно. Преимущество густоты (N) перед полнотой (G = 0,785 D2срN) при оценке массы охвоенных побегов сосны на 1 га было статистически доказано С. Б. Байзаковым [4], но позднее к такому же выводу приходит В. В. Кузьмичев [6] и при оценке запаса стволовой древесины.
Агрегированный фактор – полнота оказывается несостоятельным в качестве предиктора при оценке фитомассы древесного полога, поскольку при одном и том же D2ср масса кроны с увеличением N не возрастает монотонно, как это имеет место при оценке запаса стволов, а изменяется по колоколообразной кривой, т.е. вначале увеличивается, а затем снижается. В результате при одном и том же значении G, но противоположных по величине составляющих Dср и N фитомасса крон в древостое различается в 2-3 раза [17, 18].
Изложенное имеет целью объяснить логику структурирования многофакторных моделей лесной фитомассы, разработанных на обширном экспериментальном материале для 13 лесообразующих пород Северной Евразии:
Pi/M = f ( A, Dср, Нср, N ) (1)
или Pi = f ( A, Dср, Нср, N, M), (2)
где Pi - фитомасса i-й фракции (стволы, ветви, хвоя, корни) древостоя, т/га, и где сумма площадей сечений G расчленена на составляющие Dср и N. Попытка включить в модель агрегированный показатель G всегда оказывалась безуспешной вследствие его статистической недостоверности [17, 18].
Несмотря на очевидную полезность расчленения полноты древостоя на составляющие, при расчете регрессионных моделей фитомассы продолжается использование агрегированного предиктора – полноты. Отрицая правомерность расчленения G на Dср и N при оценке фитомассы древостоев согласно (1) и (2) под тем предлогом, что “этот подход может дать существенные отклонения”, З.Я. Нагимов с соавт. [11; с. 142] рассчитывают модели фитомассы, вводя вместо двух факторов (Dср и N) один: либо G, либо Dср. Оба варианта дают у них примерно одинаковые R2, но густота в них выпадает в первом случае потому, что она не “работает”, будучи жестко вписанной в агрегированный фактор G, а во втором – потому, что так “мы склонны считать” (там же; с. 142).
Что же касается опасений З.Я. Нагимова с соавт. о возможных смещениях моделей (1)-(2), то здесь всё обстоит с точностью до наоборот: именно их “упрощенные” модели при экстраполяции на лесопокрытую площадь Невьянского лесхоза Свердловской области дали совершенно неприемлемые смещения (до 50 %) в основном вследствие неучета густоты древостоев [17].
З. Я. Нагимов с соавт. [11] игнорируют необходимость содержательного анализа на заключительном этапе пассивного эксперимента. Дело в том, что “оценки коэффициентов регрессии в реальных задачах так называемого пассивного (т.е. непланируемого) эксперимента всегда все же оказываются смещенными в силу того обстоятельства, что никогда нельзя включить в рассмотрение все независимые переменные, ответственные за изучаемое явление” [13; с. 281]. Поэтому прежде, чем отстаивать и рекомендовать ту или иную “упрощенную” регрессионную модель, вместо довода: “Мы склонны считать…”, неплохо бы оценить ее по минимуму упомянутых В. В. Налимовым смещений. Говоря словами Дж. Бернала, видимо “желание заработать почести исключительно для себя” [5; с. 689] двигает иногда исследователем, когда он увлекается необоснованными и сомнительными новациями, получая в итоге артефакт.
Таким образом, все возрастающие возможности вычислительной техники провоцируют исследователя на довольно рискованные пассивные эксперименты. Загружая в компьютер многофакторные матрицы данных, надо иметь в виду, что математическое моделирование, по В.В. Налимову [12] и Р. Мак-Лоуну [7], - это искусство применения логики и математики, а не просто умение запускать стандартную программу.
1. Адлер Ю.П., Горский В.Г. Предисловие к русскому изданию // Прикладной регрессионный анализ: В 2-х кн. Кн. 1 /Дрейпер Н., Смит Г.- М.: Финансы и статистика, 1986.- 366 с.
2. Адлер Ю.П., Горский В.Г. Предисловие к русскому изданию // Прикладной регрессионный анализ: В 2-х кн. Кн. 2 /Дрейпер Н, Смит Г.- М.: Финансы и статистика, 1987.- 352 с.
3. Адлер Ю.П., Маркова Е.В., Грановский Ю.В. Планирование эксперимента при поиске оптимальных условий.- М.: Наука, 1976.- 279 с.
4. Байзаков С. Б. Некоторые закономерности накопления древесной зелени в сосновых лесах Казахстана и перспективы ее промышленного использования: Автореф. дис… к.с.-х.н.- Алма-Ата: КазСХИ, 1969.- 28 с.
5. Бернал Дж. Наука в истории общества.- М.: Изд-во иностранной литературы, 1956.- 735 с.
6. Кузьмичев В.В. Закономерности роста древостоев.- Новосибирск: Наука, 1977.- 160 с.
7. Мак-Лоун Р.Р. Математическое моделирование – искусство применения математики // Математическое моделирование.- М.: Мир, 1979.- С. 9-20.
8. Маркова Е.В., Лисенков А.Н. Планирование эксперимента в условиях неоднородностей.- М.: Наука, 1973.- 219 с.
9. Мауринь А.М., Лиепа И.Я., Дрике А.Я., Поспелова Г.Е. Прогнозирование плодоношения древесных растений // Оптимизация использ. и воспроизводства лесов СССР.- М.: Наука, 1977.- С. 50-53.
10. Митропольский А.К. О множественных нелинейных корреляционных уравнениях // Изв. АН СССР, сер. матем.- 1939.- Т. 3.- С. 399-406.
11. Нагимов З.Я., Шевелина И.В., Анчугова Г.В., Сальникова И.С. Регрессионные модели переводных коэффициентов надземной фитомассы древостоев // Леса Урала и хоз-во в них. Вып. 21.- Екатеринбург: УГЛТУ, 2001.- С. 141-147.
12. Налимов В.В. Теория эксперимента.- М.: Наука, 1971.- 208 с.
13. Налимов В.В. Канатоходец.- М.: Прогресс, 1994.- 456 с.
14. Налимов В.В., Чернова Н.А. Статистические методы планирования экстремальных экспериментов.- М.: Наука, 1965.- 340 с.
15. Тутубалин В.Н., Барабашева Ю.М., Григорян А.А. и др. Математическое моделирование в экологии: Историко-методологический анализ.- М.: Языки русской культуры.- 1999.- 208 с.
16. Усольцев В.А. Рост и структура фитомассы древостоев. -Новосибирск: Наука, 1988. -253 с.
17. Усольцев В.А. Формирование банков данных о фитомассе лесов.- Екатеринбург: Изд-во УрО РАН, 1998. –541 с.
18. Усольцев В.А. Фитомасса лесов Северной Евразии: база данных и география.- Екатеринбург: Изд-во УрО РАН, 2001.- 708 с.
19. Fisher R.A. The influence of rainfall on the yield of wheat at Rothamsted // Phil. Trans. Roy. Soc.- 1924.- Vol. 213.- P. 89-142.
Работа поддержана РФФИ (гранты №№ 00-05-64532 и 01-04-96424).