낙엽활엽수의 수고-흉고직경 조건부 분포 분석을 위한 수종별 분위수 회귀식 구축^†

본 연구에서는 수종별 수고–흉고직경 관계의 조건부 분포를 분석하기 위해 국가산림자원조사(National Forest Inventory, NFI) 자료를 활용하였다(산림청, 2022a; 2022b; 2022c). NFI 자료는 전국 산림을 대상으로 표준화된 조사체계에 따라 구축된 대규모 자료로서, 다양한 지역과 입지 및 임분조건에 분포하는 수목의 구조적 특성을 수종별로 비교하는 데 활용될 수 있다(강진택, 2019; Sharma et al., 2019). 다만 NFI 자료는 도시숲이나 가로수를 대상으로 수집된 자료가 아니므로, 본 연구에서는 분석 결과를 도시수목의 생육상태를 직접 판정하는 기준이 아니라 산림수목의 수고–흉고직경 관계에 관한 탐색적 비교자료로 한정하였다.

분석대상 기간은 제5차, 제6차 및 제7차 국가산림자원조사에 해당하는 2006년부터 2020년까지로 설정하였다. 연도별 NFI 원자료의 임목조사표에서 집락번호, 표본점번호, 조사차기, 조사연도, 개체목 번호, 학명번호, 수종명, 흉고직경, 실측수고, 개체목 구분코드 및 수령을 추출하였다.

NFI 임목조사표에는 현장에서 기록된 실측수고와 수고곡선식으로 산정된 추정수고가 함께 수록되어 있다. 추정수고는 흉고직경을 포함한 수고곡선식으로 산정되므로, 이를 다시 흉고직경의 설명변수로 사용할 경우 두 변수 사이에 순환논리가 발생할 수 있다. 이에 본 연구에서는 추정수고를 분석에서 제외하고, 현장에서 측정된 실측수고와 흉고직경이 함께 기록된 관측값만을 사용하였다. 실측수고는 모든 조사 개체목에 기록된 변수가 아니므로, 본 연구의 분석대상은 NFI 전체 조사목이 아니라 실측수고와 흉고직경이 모두 기록된 개체목으로 한정하였다.

수종 후보군은 2025년도 도시숲경관 사업계획에 제시된 온대중부 적용 가능 수종을 기준으로 설정하였다. 해당 계획은 국내 도시숲 적용 가능 수종을 온대중부, 온대남부 및 난온대 유형으로 구분하여 제시하고 있다. 본 연구에서는 이 가운데 온대중부 수종군에 포함된 낙엽활엽교목 및 낙엽활엽아교목 52종을 1차 후보수종으로 선정하였다(표 1 참조). 침엽교목은 낙엽활엽수와 수관구조, 줄기 형태 및 수고–흉고직경 관계가 다르게 나타날 수 있으므로, 생육형이 비교적 유사한 수종군 내에서 결과를 비교하기 위해 후보수종에서 제외하였다.

분석자료를 구축하기 위해 표 1에 제시된 52개 1차 후보수종과 NFI 원자료의 수종명을 대조한 후, 후보수종에 해당하는 개체목 자료만을 추출하였다. 이에 따라 NFI 원자료에 존재하더라도 표 1의 후보수종 목록에 포함되지 않은 수종은 분석자료 구축 단계에서 제외하였다. 추출된 자료 가운데 실측수고 또는 흉고직경이 결측이거나 0 이하인 관측값을 제외하였으며 수종별 실측수고와 흉고직경의 최솟값과 최댓값을 검토하여 입력 오류 또는 생물학적으로 해석하기 어려운 극단값의 존재 여부를 확인하였다.

이후 수종별 분위회귀모형의 추정에 필요한 관측값과 표본점 단위 검증에 필요한 집단 수를 확보하기 위해 유효 관측값이 100개 이상이고 고유 표본점이 30개 이상인 수종을 주 분석대상으로 선정하였다. 유효 관측값 100개 기준은 절대적인 통계 임계값이라기보다, q10 및 q90과 같은 꼬리 분위수를 추정할 때 지나치게 작은 표본으로 인해 회귀선이 불안정하게 추정되는 문제를 줄이기 위한 기준으로 적용하였다. 또한 고유 표본점 30개 기준은 특정 수종의 관측값이 소수의 표본점에 편중되는 상황을 완화하고, 표본점 단위 5-fold 교차검증과 군집 부트스트랩을 수행할 수 있을 만큼의 집단 수를 확보하기 위해 설정하였다. 이러한 자료 선별 절차를 거쳐 유효 관측값 수와 고유 표본점 수 기준을 모두 충족한 수종을 최종 분석대상으로 선정하였다.

모든 자료처리와 통계분석은 Python 3.9.6 환경에서 수행하였다. 국가산림자원조사 원자료의 불러오기, 결측자료 처리, 변수 변환, 수종별 자료 분할, 반복 표본점 및 반복 개체 후보 확인, 기초통계량 산정과 결과표 작성에는 pandas 2.3.1과 NumPy 2.0.2를 활용하였다.

수종별 q10 및 q90 분위회귀모형은 statsmodels 0.14.6의 QuantReg를 이용하여 추정하였다. 수종별 선형, 이차 및 cubic B-spline 후보모형의 예측성능을 비교하기 위한 표본점 단위 5-fold 교차검증에는 scikit-learn 1.6.1의 GroupKFold를 사용하였으며, cubic B-spline 기저함수의 생성에는 동일 라이브러리의 SplineTransformer를 활용하였다. 표본점 단위 군집 부트스트랩의 난수 생성과 복원추출에는 NumPy의 default_rng를 사용하였으며, 분석의 재현성을 확보하기 위해 고정된 난수 시드를 적용하였다.

Cubic B-spline 후보모형에서는 SplineTransformer로 생성한 기저변수를 각 훈련 fold에서 표준화하고, 동일한 변환을 해당 검증 fold에 적용하였다. 이 절차는 기저변수의 척도 차이로 인한 수치적 불안정성을 줄이면서 검증자료의 정보가 훈련과정에 유입되는 것을 방지하기 위해 적용하였다.

본 연구는 실측수고가 주어졌을 때 흉고직경 조건부 분포의 특정 위치를 추정하기 위해 분위회귀분석을 적용하였다. 조건부 평균을 추정하는 최소제곱회귀와 달리 분위회귀는 q10, q50, q90과 같이 관심 있는 조건부 분위함수를 직접 추정하므로, 분포가 비대칭이거나 수고 구간에 따라 산포가 달라지는 자료에서 하위와 상위 관계를 구분할 수 있다(Koenker, 2017; Zhang et al., 2020).

본 연구에서는 흉고직경을 종속변수로, 실측수고를 설명변수로 설정하였다. 회귀계수의 해석 편의성과 비선형모형의 수치적 안정성을 고려하여 실측수고는 원자료에 기록된 cm 단위 값을 m 단위로 환산하였으며, 흉고직경은 원자료의 cm 단위를 유지하였다. 본 연구에서 추정한 모형은 동일한 개체목의 시간에 따른 변화를 분석하는 생장모형이 아니라, 2006–2020년 여러 조사차수를 통합한 자료에서 시간적 변화가 아닌 개체목 간 실측수고와 흉고직경의 구조적 관계를 분석하는 모형이다. 따라서 모형의 계수와 예측값은 개별 수목의 미래 생장을 의미하는 것이 아니라, 통합된 조사자료에서 동일한 실측수고 조건에 나타나는 흉고직경의 조건부 분포로 해석하였다.

흉고직경 조건부 분포의 하위 및 상위 위치를 함께 파악하기 위해 q10과 q90을 분석 분위로 설정하였다. q10은 동일한 실측수고 조건에서 흉고직경의 조건부 10% 분위값으로, 관측값의 약 10%가 그 이하에 위치한다. q90은 조건부 90% 분위값으로, 이론적으로 관측값의 약 90%가 그 이하에 위치한다. 따라서 q10 이하와 q90 초과에는 각각 약 10%의 관측값이 위치하며, q10과 q90 사이에는 조건부 분포의 중앙 약 80%가 포함된다.

이러한 분위 설정은 직경생장 예측 및 산림 생육범위 추정에 분위회귀식이 활용되어 온 선행연구를 고려해 반영한 것이다. 선행연구에서는 복수의 분위회귀식을 기반으로 직경생장 예측이 가능함을 제시하였으며, 99%와 같은 극단 분위수는 대규모 자료에서 이상치의 영향을 받을 수 있어 90% 분위수가 보다 안정적인 상위 경계 추정에 적합할 수 있다고 언급하였다(Bohora and Cao, 2014; Raptis et al., 2022). 또한 스페인 대규모 산림 인벤토리 자료를 활용한 연구에서는 생육범위의 경계를 설정하기 위해 수관폭의 90% 분위수와 추정수고의 10% 분위수를 기준으로 개방 생육목을 선별하고, 이들의 수관폭 관계를 분석하였다(Condés and Sterba, 2005). 이러한 선행연구를 바탕으로 본 연구에서는 극단 분위수에 비해 일부 관측값의 영향을 상대적으로 줄이면서 조건부 분포의 하위와 상위 위치를 함께 비교할 수 있도록 10%와 90% 분위수를 적용하였다.

수고–흉고직경 관계는 수종과 개체의 크기 범위에 따라 직선 또는 곡선 형태로 나타날 수 있으므로, 선형모형, 평균중심화 이차모형 및 cubic B-spline 모형을 후보모형으로 설정하였다. 선형모형은 실측수고가 서로 다른 개체들 사이에서 흉고직경의 조건부 분위값 차이가 일정한 비율로 나타난다고 가정하는 가장 단순한 모형이다. 회귀계수의 의미가 명확하고 수종 간 비교가 용이하지만, 수고 구간에 따라 흉고직경 분위값의 차이가 달라지는 관계를 충분히 표현하지 못할 수 있다.

이차모형은 실측수고와 그 제곱항을 사용하여 전체 수고 범위에 하나의 완만한 곡선을 적용하는 모형이다. 이 모형은 수고가 큰 개체에서 흉고직경 조건부 분위값의 차이가 점차 커지거나 작아지는 관계를 표현하기 위해 사용하였다. 실측수고와 그 제곱항은 동일한 변수에서 생성되므로 구조적 공선성이 발생할 수 있다. 이를 완화하기 위해 각 개체목의 실측수고에서 해당 수종의 평균 실측수고를 뺀 평균중심화 변수와 그 제곱항을 사용하였다. 이때 1차항은 수종별 평균 수고 부근에서 수고가 서로 다른 개체들 사이에 나타나는 흉고직경 분위값의 차이를, 제곱항은 그 관계가 휘어지는 방향과 정도를 나타낸다.

Cubic B-spline 모형은 실측수고의 범위를 여러 구간으로 나누어 각 구간에 3차 다항식 형태의 곡선을 적용한 후, 구간 사이를 부드럽게 연결하는 방법이다. 이 모형은 낮은 수고 구간과 높은 수고 구간에 속한 개체들에서 흉고직경의 조건부 분위값 차이가 서로 다른 복잡한 비선형 관계를 검토하기 위해 사용하였다. 본 연구에서는 실측수고 범위에 균등하게 설정된 4개의 절점과 3차 B-spline을 적용하고, 표본점 단위 교차검증 예측성능을 선형모형 및 이차모형과 비교하였다. 즉, B-spline은 모든 수종의 기본모형이 아니라 선형 또는 이차모형보다 복잡한 함수형태가 실제로 필요한지를 확인하기 위한 비교 후보모형으로 사용하였다.

후보모형의 비교에는 표본점번호를 그룹으로 하는 5-fold 교차검증을 적용하였다. 교차검증의 fold 수는 모형 추정에 충분한 훈련자료를 확보하는 동시에 각 검증자료에도 일정 수 이상의 표본점이 포함되도록 하기 위해 5개로 설정하였다. 본 연구에서는 수종별 고유 표본점이 30개 이상인 경우를 분석대상으로 선정하였으므로, 5-fold 적용 시 각 반복에서 전체 표본점의 약 80%를 훈련자료로 활용하면서 나머지 약 20%를 독립적인 검증자료로 사용할 수 있다. 또한 최소 기준인 30개 표본점을 기준으로 할 경우 각 fold에 약 6개 표본점이 포함되어, 10-fold와 같이 fold 수를 증가시켰을 때 개별 검증 fold에 포함되는 표본점 수가 지나치게 감소하는 문제를 완화할 수 있다.

따라서 본 연구에서는 전체 표본점은 서로 겹치지 않는 5개 집단으로 나누고, 이 가운데 4개 집단으로 모형을 추정한 후 나머지 1개 집단으로 예측성능을 평가하였다. 검증 집단을 차례로 변경하여 동일한 과정을 5회 반복함으로써 모든 표본점이 한 번씩 검증자료로 사용되도록 하였다. 동일한 표본점에서 조사된 수목은 입지 및 임분환경 등 일부 생육조건을 공유할 수 있으므로, 같은 표본점의 관측값이 훈련자료와 검증자료에 동시에 포함될 경우 자료의 군집구조로 인해 예측성능이 과대평가될 가능성이 있다. 따라서 동일한 표본점에 속한 모든 개체목 관측값을 하나의 fold에 배정하여 훈련자료와 검증자료 간 표본점의 중복을 방지하였다. 이를 통해 모형 추정에 사용되지 않은 표본점의 개체목을 대상으로 각 후보모형이 목표 분위값을 얼마나 잘 추정하는지를 비교하였다.

후보모형의 예측오차는 분위회귀에 적합한 손실함수인 pinball loss를 이용하여 평가하였다. Pinball loss는 목표 분위에 따라 과대예측과 과소예측에 서로 다른 가중치를 부여한다. q10에서는 하위 분위값을 지나치게 높게 예측하는 경우에 상대적으로 큰 손실을 부여하고, q90에서는 상위 분위값을 지나치게 낮게 예측하는 경우에 상대적으로 큰 손실을 부여한다. 따라서 pinball loss가 작을수록 해당 모형이 검증자료에서 목표로 설정한 조건부 분위의 위치를 더 정확하게 추정한 것으로 해석하였다.

각 후보모형에서 q10과 q90의 교차검증 pinball loss를 각각 산출하고 두 값을 합산하여 하위 및 상위 분위의 예측성능을 함께 평가하였다. 이는 두 분위 가운데 한쪽만을 잘 추정하는 모형이 아니라, 흉고직경 조건부 분포의 하위와 상위 위치를 종합적으로 잘 추정하는 모형을 선정하기 위한 것이다.

최종모형은 합산 pinball loss에 1 표준오차 규칙을 적용하여 선정하였다. 먼저 평균 교차검증 손실이 가장 작은 모형을 확인하고, 해당 모형의 평균 손실에 1 표준오차를 더한 값을 허용 가능한 손실의 상한으로 설정하였다. 이후 이 범위에 포함되는 후보 중 선형모형, 이차모형, B-spline 모형의 순서로 가장 단순한 모형을 선택하였다. 이 절차는 복잡한 모형의 예측성능 개선이 교차검증 과정의 표본 변동을 고려한 이후에도 충분히 큰 경우에만 비선형모형을 채택하고, 성능 차이가 크지 않은 경우에는 해석이 용이한 단순모형을 선택하기 위해 적용하였다.

최종모형의 성능은 교차검증 pinball loss, 관측비율, q10–q90 범위의 포괄률 및 분위회귀선의 교차 여부를 이용하여 평가하였다. 관측비율은 검증자료의 실제 흉고직경이 q10 또는 q90 예측값 이하에 위치하는 비율이며, 각각 목표 분위인 0.10과 0.90에 가까운지를 확인하였다. 포괄률은 실제 흉고직경이 q10과 q90 예측값 사이에 포함되는 비율로, 두 분위식이 의도한 조건부 중앙 80% 범위를 어느 정도 재현하는지를 확인하기 위해 사용하였다. 관측비율과 포괄률은 모형의 채택 또는 제외를 위한 절대적인 임계값이 아니라, 목표 분위와 분포 범위의 재현 정도를 확인하는 진단지표로 활용하였다.

회귀계수의 불확실성은 표본점 단위 군집 부트스트랩으로 평가하였다. 동일한 표본점의 수목은 생육환경을 공유하여 서로 유사할 수 있으므로, 개별 수목을 독립적으로 재표집하면 회귀계수의 불확실성이 실제보다 작게 평가될 가능성이 있다. 이에 개별 수목이 아닌 표본점을 복원추출하고, 선택된 표본점에 속한 모든 관측값을 함께 포함하여 최종모형을 다시 추정하였다. 이러한 재표집을 1,000회 반복하여 부트스트랩 표준오차와 백분위수 방식의 95% 신뢰구간을 산출하였다. 군집 부트스트랩은 표본점의 환경효과를 직접 제거하는 방법이 아니라, 동일 표본점의 관측값을 하나의 군집으로 유지하면서 표본점 구성에 따른 회귀계수의 불확실성을 평가하기 위해 적용하였다.

마지막으로 q10과 q90 분위회귀선이 조건부 분포의 순서를 일관되게 나타내는지를 확인하였다. q10은 하위 분위값이고 q90은 상위 분위값이므로, 동일한 수고에서 q10 예측값이 q90 예측값보다 크게 나타나면 두 분위식이 논리적인 순서를 충족하지 못한 것으로 볼 수 있다. 이에 수종별 관측 실측수고의 최솟값과 최댓값 사이에 500개의 평가 지점을 생성하고, 각 지점에서 두 분위회귀선의 교차 여부를 점검하였다. Pseudo R²는 최종모형의 표본 내 상대적 적합도를 설명하기 위한 보조적 기술지표로만 산출하였으며, 특정 값을 모형의 채택 또는 제외 기준으로 사용하지 않았다.

본 연구에서는 분석자료 구축과 기초통계량 산정, 반복 관측구조 확인, 후보 분위회귀모형의 비교, 최종모형 추정과 평가 및 민감도 분석의 순서로 분석을 수행하였다. 먼저 2006–2020년 국가산림자원조사 원자료에서 수종명, 실측수고, 흉고직경, 집락번호, 표본점번호, 조사차기, 조사연도, 개체목 번호 및 개체목 구분코드를 추출하였다. 수종명은 표 1에 제시한 52개 후보 수종과 대조하여 일치하는 관측값만을 추출하였으며, 후보 수종 목록에 포함되지 않은 수종은 이후의 표본 수 기준 적용에 앞서 제외하였다. 실측수고 또는 흉고직경이 결측이거나 0 이하인 관측값도 분석에서 제외하였다. 이후 수종별 유효 관측값 수와 고유 표본점 수를 산정하고, 유효 관측값이 100개 이상이면서 고유 표본점이 30개 이상인 수종을 주 분석대상으로 선정하였다.

구축된 자료의 표본 규모와 기본적인 분포 특성을 파악하기 위해 전체 자료와 수종별 자료에서 실측수고와 흉고직경의 관측값 수, 평균, 표준편차, 중앙값, 최솟값 및 최댓값을 산출하였다. 또한 수종별 고유 표본점 수와 표본점당 관측값 수를 산정하여 관측값이 일부 표본점에 집중되어 있는지를 검토하였다. 최솟값과 최댓값은 입력 오류 또는 생물학적으로 해석하기 어려운 관측값의 존재 가능성을 확인하는 데 활용하였다. 이러한 기초통계량 검토를 거쳐 구축된 자료를 이후의 수종별 분위회귀분석에 사용하였다.

다음으로 제5차부터 제7차 국가산림자원조사 자료를 통합함에 따라 발생할 수 있는 관측값 간 종속성을 확인하였다. 표본점번호를 기준으로 둘 이상의 조사차기에 포함된 표본점을 반복 표본점으로 분류하였다. 동일 표본점에서 조사된 수목은 토양, 지형, 광환경 및 수분 조건 등의 생육환경을 공유할 수 있으므로, 교차검증에서는 같은 표본점의 관측값을 하나의 집단으로 배정하고 군집 부트스트랩에서는 표본점을 재표집 단위로 설정하였다. 또한 표본점번호와 개체목 번호를 결합하여 개체 후보 식별자를 구성하고, 동일한 식별자가 둘 이상의 조사차기에 나타난 경우 반복 개체 후보로 분류하였다. 반복 개체 후보는 동일한 개체가 통합자료에 여러 차례 포함되었을 가능성이 분석 결과에 미치는 영향을 확인하기 위한 민감도 분석에 활용하였다.

이후 흉고직경을 종속변수로, m 단위로 환산한 실측수고를 설명변수로 설정하여 수종별 q10 및 q90 분위회귀모형을 구축하였다. 각 분위에 대해 선형모형, 평균중심화 이차모형 및 cubic B-spline 모형을 후보로 설정하였다. 후보모형은 표본점번호를 그룹으로 하는 5-fold 교차검증으로 비교하였으며, 같은 표본점의 자료가 훈련자료와 검증자료에 동시에 포함되지 않도록 하였다. 각 fold에서 q10과 q90의 pinball loss를 산출하여 합산하고, 합산 손실의 평균과 표준오차를 이용한 1 표준오차 규칙을 적용하였다. 평균 합산 손실이 가장 작은 모형의 1 표준오차 범위에 포함되는 후보 중 가장 단순한 모형을 수종별 최종모형으로 선정하였다.

선정된 최종모형은 수종별 전체 분석자료에 다시 적합하여 q10 및 q90 분위회귀계수를 산출하였다. 회귀계수의 불확실성은 동일 표본점 내 관측값의 유사성을 고려할 수 있도록 표본점을 복원추출하는 군집 부트스트랩을 1,000회 수행하여 평가하였다. 부트스트랩 계수 분포를 이용하여 표준오차와 백분위수 방식의 95% 신뢰구간을 산출하였으며, 분석의 재현성을 위해 고정된 난수 시드를 적용하였다. 최종모형의 검증성능은 표본점 단위 교차검증에서 얻은 예측값을 이용하여 q10 및 q90의 관측비율과 두 분위 사이의 포괄률로 평가하였다. 관측비율은 각각 목표 분위인 0.10과 0.90에 가까운지를, 포괄률은 조건부 중앙 80% 범위를 적절하게 재현하는지를 확인하기 위해 사용하였다. 또한 수종별 관측 실측수고 범위에 500개의 평가 지점을 생성하여 q10 예측값이 q90 예측값보다 커지는 분위회귀선 교차 여부를 확인하였다.

마지막으로 최소 관측값 기준의 변화가 분석대상 수종에 미치는 영향과 반복 관측 가능성이 결과에 미치는 영향을 확인하기 위해 민감도 분석을 수행하였다. 최소 관측값 기준은 50개, 75개, 100개 및 150개로 변경하되, 고유 표본점 30개 이상의 기준은 동일하게 유지하였다. 고유 표본점 기준을 유지한 이유는 많은 개체목이 소수의 표본점에 집중되는 것을 방지하고, 표본점 단위 5-fold 교차검증과 군집 부트스트랩에 필요한 집단 수를 확보하는 동시에 최소 관측값 기준의 변화에 따른 영향만을 구분하여 확인하기 위함이다. 이에 따라 각 최소 관측값 기준에서 분석대상으로 포함되는 수종의 구성을 비교하였다. 이 분석은 최소 관측값 기준의 변화에 따른 분석대상 수종 구성의 민감도를 평가한 것으로, 개별 수종의 표본 수를 단계적으로 변화시키면서 q10 및 q90 회귀계수나 예측선의 안정성을 평가한 분석과는 구분하였다.

이와 별도로 반복 관측 가능성에 대해서는 반복 개체 후보별 최초 관측값 한 건만 남긴 자료로 모형을 다시 추정하고, 전체 통합자료에서 산출된 회귀계수와 방향 및 크기를 비교하였다. 아울러 제5차, 제6차 및 제7차 자료를 각각 분리하여 조사차기별 회귀계수를 추정하고, 이를 전체 통합자료의 회귀계수와 비교하였다.

표 1에 제시된 52개 1차 후보수종과 국가산림자원조사 원자료의 수종명을 대조하고, 실측수고와 흉고직경이 모두 유효한 관측값을 추출한 결과, 40개 수종의 64,382개 관측값이 확보되었다. 이 자료는 8,889개의 고유 표본점에서 조사되었으며, 제5차 22,878개, 제6차 20,752개 및 제7차 20,752개의 관측값으로 구성되었다.

수종별 유효 관측값 100개 이상과 고유 표본점 30개 이상의 기준을 적용한 결과, 25개 수종의 64,072개 관측값이 최종 분석자료에 포함되었다. 최종 분석자료는 8,856개 고유 표본점으로 구성되었다. 분석대상 수종 가운데 상수리나무의 관측값이 14,868개로 가장 많았고, 왕벚나무가 174개로 가장 적었다. 고유 표본점 수는 밤나무가 2,683개로 가장 많았으며, 백합나무가 31개로 가장 적었다.

최종 분석자료의 실측수고는 평균 11.7 ± 4.2m, 중앙값 11.4m였으며, 관측범위는 1.5–34.1m였다. 흉고직경은 평균 18.8 ± 9.0cm, 중앙값 17.0cm였으며, 관측범위는 6.0–104.0cm였다. 수종별 평균 실측수고는 함박꽃나무가 6.6m로 가장 낮고 은사시나무가 16.1m로 가장 높았다. 평균 흉고직경은 함박꽃나무가 9.7cm로 가장 작고 가래나무가 26.7cm로 가장 크게 나타났다. 수종별 최솟값과 최댓값을 검토한 결과, 양수가 아닌 값의 제외 이후 추가로 제거한 관측값은 없었다(표 2 참조).

최종 분석자료에 포함된 8,856개 고유 표본점의 조사차기별 반복 여부를 확인한 결과, 2개 이상의 조사차기에 포함된 반복 표본점은 6,529개로 전체 표본점의 73.7%를 차지하였다. 이 가운데 1,664개 표본점은 2개 조사차기에, 4,865개 표본점은 제5차부터 제7차까지 3개 조사차기에 모두 포함되었다. 한 개 조사차기에만 포함된 표본점은 2,327개로 나타났다.

분석대상 수종 전체에서 표본점번호와 개체목 번호를 결합하여 구성한 개체 후보 식별자는 총 56,651개였다. 이 가운데 2개 이상의 조사차기에 동일한 식별자가 나타난 반복 개체 후보는 6,464개로 전체 개체 후보의 11.4%를 차지하였다. 반복 개체 후보 중 5,646개는 2개 조사차기에, 818개는 3개 조사차기에 나타났다. 다만 이 결과는 공개 원자료에서 동일한 표본점번호와 개체목 번호가 반복된 경우를 집계한 것으로, 조사차기 사이에 동일한 개체목임을 확정한 결과가 아니라 반복 가능성이 있는 개체 후보의 분포를 나타낸다(표 3 참조).

최종 분석대상 25개 수종에 대해 선형모형, 평균중심화 이차모형 및 cubic B-spline 모형의 표본점 단위 5-fold 교차검증을 수행한 결과, 모든 수종과 후보모형에서 실패한 fold 없이 q10 및 q90 합산 pinball loss가 산출되었다.

평균 합산 pinball loss가 가장 작은 모형만을 기준으로 비교하면 선형모형이 10개 수종, 이차모형이 9개 수종, B-spline 모형이 6개 수종에서 가장 낮은 손실을 나타냈다. 그러나 최소 손실 모형의 1 표준오차 범위에 포함되는 후보 중 가장 단순한 모형을 선택한 결과, 최종적으로 23개 수종에서 선형모형이 선정되었고 가죽나무와 갈참나무에서는 이차모형이 선정되었다. B-spline 모형이 최종모형으로 선정된 수종은 없었다.

가죽나무의 평균 합산 pinball loss는 이차모형에서 1.3456으로 가장 낮았으며, 1 표준오차 허용 상한은 1.4311이었다. 선형모형의 평균 손실은 1.4401로 허용 상한보다 높게 나타나 이차모형이 최종모형으로 선정되었다. 갈참나무에서도 이차모형의 평균 손실이 1.7328로 가장 낮았으며, 허용 상한은 1.7482였다. 선형모형의 평균 손실은 1.7494로 허용 상한을 초과하여 이차모형이 선정되었다.

느티나무, 마가목, 밤나무, 상수리나무, 은사시나무 및 층층나무에서는 B-spline 모형의 평균 손실이 가장 작게 나타났으나, 선형모형의 평균 손실도 최소 손실 모형의 1 표준오차 범위에 포함되었다. 이에 따라 해당 수종에서는 예측성능의 차이가 모형의 복잡성 증가를 뒷받침할 정도로 크지 않은 것으로 판단하여 선형모형을 선정하였다. 가래나무, 고로쇠나무, 당단풍나무, 백합나무, 쪽동백나무, 팽나무 및 함박꽃나무에서도 이차모형의 평균 손실이 가장 작았으나 선형모형이 1 표준오차 범위에 포함되어 최종적으로 선형모형이 선정되었다(표 4 참조).

1 표준오차 규칙에 따라 선정된 수종별 최종모형을 각 수종의 전체 분석자료에 다시 적합하여 q10 및 q90 분위회귀계수를 추정하였다(표 5 참조). H는 현장에서 측정한 수고(m), DBH는 흉고직경(cm)을 의미한다.

선형모형은 DBH = β0 + β1H로, 이차모형은 $DBH=\beta 0+\beta 1\left(H-\overline{H}\right)+\beta 2{\left(H-\overline{H}\right)}^2$로 표현하였다. 표 5에 제시한 수고 범위는 각 수종의 최종 분석자료에서 관측된 실측수고의 최솟값과 최댓값이며, 이를 벗어난 수고값에 대한 예측은 외삽에 해당하므로 회귀식은 해당 범위 내에서만 적용하였다. 수종별 최종 분위회귀선은 그림 1에 제시하였다. 회귀계수의 불확실성은 표본점을 재표집 단위로 하는 군집 부트스트랩을 1,000회 수행하여 평가하였으며, 부트스트랩 표준오차와 백분위수 방식의 95% 신뢰구간을 산출하였다(표 6 참조).

그림 1. 최종 수종별 분위회귀선 범례:

Download Original Figure Download Original Figure Download Original Figure Download Original Figure

최종모형의 1차항 추정치는 모든 수종의 q10과 q90에서 양수로 나타났다. 이는 관측된 실측수고 범위에서 수고가 큰 개체일수록 흉고직경의 조건부 하위 및 상위 분위값도 전반적으로 크게 나타났음을 의미한다. 또한 모든 수종에서 q90의 1차항 추정치가 q10보다 컸다. 따라서 개체 간 실측수고 차이에 대응하는 흉고직경 분위값의 차이는 조건부 분포의 하위보다 상위에서 상대적으로 크게 나타났다. 다만 함박꽃나무는 q10과 q90 모두 1차항의 군집 부트스트랩 신뢰구간 하한이 0에 근접하여 다른 수종보다 수고–흉고직경 관계의 불확실성이 컸다.

이차모형이 선정된 가죽나무와 갈참나무에서는 q10과 q90의 이차항이 모두 양수로 추정되었으며, 군집 부트스트랩 95% 신뢰구간도 0보다 큰 범위로 나타났다. 가죽나무의 이차항은 q10에서 0.065, q90에서 0.143이었고, 갈참나무는 각각 0.033과 0.049였다. 이는 두 수종에서 수고가 큰 개체군으로 갈수록 흉고직경 분위값의 개체 간 차이가 일정하게 유지되기보다 커지는 곡선 관계가 나타났음을 의미한다. 이차모형의 계수는 수종별 평균 실측수고를 기준으로 중심화하여 산출하였으며, 중심수고는 가죽나무 10.856m, 갈참나무 11.637m였다.

수종별 최종모형의 검증성능은 표본점 단위 5-fold 교차검증에서 산출된 예측값을 이용하여 평가하였다(표 7 참조). 모든 수종에서 5개의 교차검증 fold가 정상적으로 분석되어 모형 추정에 실패한 fold는 나타나지 않았다.

q10의 관측비율은 0.095–0.132의 범위로 나타났으며, 수종 평균은 0.105였다. q90의 관측비율은 0.885–0.906의 범위였고, 전체 평균은 0.896으로 나타났다. 따라서 검증자료에서 q10과 q90의 관측비율은 각각 목표 분위인 0.10과 0.90에 전반적으로 근접하였다.

q10과 q90 예측값 사이에 실제 흉고직경이 포함된 비율인 포괄률은 0.759–0.806의 범위로 나타났으며, 수종 평균은 0.791이었다. 이는 이론적인 중앙 80% 포괄률인 0.800에 근접한 결과이다. 자작나무의 포괄률이 0.806으로 가장 높았으며, 왕벚나무는 0.759로 가장 낮게 나타났다. 특히 왕벚나무는 q10 관측비율이 0.132로 목표 분위보다 높게 나타나 다른 수종보다 하위 분위의 재현 정도가 상대적으로 낮았다. 교차검증에서 산출된 q10과 q90의 합산 pinball loss는 수종별로 0.903–2.839의 범위로 나타났다.

다만 pinball loss의 크기는 수종별 흉고직경의 분포 범위와 변동성에 영향을 받으므로, 서로 다른 수종 사이의 절대적인 모형 우수성을 판정하는 기준으로 사용하지 않았다. 해당 값은 동일 수종 내에서 후보모형을 비교하고 최종모형을 선정하기 위한 상대적 예측오차로 해석하였다.

마지막으로 수종별로 관측된 실측수고의 최솟값과 최댓값 사이에 500개의 평가 지점을 생성하여 q10과 q90 분위회귀선의 교차 여부를 확인하였다. 분석대상 25개 수종 모두에서 q10 예측값이 q90 예측값보다 높아지는 지점은 확인되지 않았다. 따라서 최종모형은 각 수종에서 관측된 실측수고 범위 내에서 조건부 하위 분위와 상위 분위의 순서를 일관되게 유지하였다.

분석자료 선정 기준이 최종 분석대상에 미치는 영향을 확인하기 위해 고유 표본점 30개 이상 기준을 유지하면서 최소 관측값 기준을 50개, 75개, 100개 및 150개로 변경하였다. 모든 기준에서 분석대상은 동일한 25개 수종으로 나타났다. 이는 고유 표본점 기준을 충족한 수종 가운데 관측값이 가장 적은 왕벚나무도 174개의 유효 관측값을 보유하고 있었기 때문이다. 따라서 최소 관측값 기준을 50–150개 범위에서 변경하더라도 분석대상 수종의 구성은 달라지지 않았다. 해당 결과는 설정한 최소 관측값 기준의 변화에 대해 분석대상 수종 구성이 동일하게 유지되었음을 의미한다. 다만 개별 수종의 표본 수를 단계적으로 변화시키면서 q10 및 q90 회귀계수나 예측선의 안정성을 평가한 결과는 아니다.

반복 관측 가능성의 영향을 확인하기 위해 표본점번호와 개체목 번호가 동일한 반복 개체 후보에서 최초 관측값 한 건만 남기고 최종모형을 다시 추정하였다. 이에 따라 분석 관측값은 전체 통합자료의 64,072개에서 56,651개로 감소하였다. 그러나 25개 수종에서 산출된 q10과 q90의 1차항 추정치 50개는 모두 전체 통합자료와 동일한 양의 방향을 유지하였다. 전체 통합자료 대비 1차항 상대차이 절댓값의 중앙값은 1.61%로 나타났다. 다만 함박꽃나무 q90의 1차항은 전체 통합자료의 0.571에서 반복 개체 후보의 최초 관측값만 사용한 경우 0.784로 증가하여 가장 큰 상대차이인 37.25%를 보였다(표 8 참조).

조사차기에 따른 결과의 변화를 확인하기 위해 제5차, 제6차 및 제7차 자료를 각각 분리하여 전체 통합자료에서 선정된 것과 동일한 형태의 모형을 다시 추정하였다. 조사차기별 관측값이 50개 이상인 수종을 대상으로 분석한 결과, 각 조사차기에서 23개 수종의 q10 및 q90 계수 46개가 산출되었다. 모든 조사차기에서 1차항은 전체 통합자료와 동일한 양의 방향으로 나타났다(표 8 참조).

전체 통합자료 대비 1차항 상대차이 절댓값의 중앙값은 제5차에서 13.78%, 제6차에서 5.16%, 제7차에서 4.52%로 나타났다. 조사차기를 분리한 경우에도 실측수고가 큰 개체일수록 흉고직경의 조건부 분위값이 크게 나타나는 관계의 방향은 유지되었다. 다만 표본 수가 상대적으로 적은 일부 수종에서는 계수 크기의 차이가 컸다. 특히 함박꽃나무 q90의 1차항은 제6차에서 전체 통합자료보다 84.21% 크게 나타났으며, 제7차에서는 0에 가까운 값으로 추정되었다(표 8 참조).

종합하면, 최소 관측값 기준의 변화에 대해서는 분석대상 수종의 구성이 동일하게 유지되었다. 이 결과는 수종 선정기준의 민감도에 관한 것이며, 개별 수종에서 표본 수 변화에 따른 회귀계수나 예측선의 안정성을 의미하지는 않는다. 이와 별도로 반복 개체 후보의 최초 관측값만 사용한 분석과 조사차기별 분석에서는 추정 가능한 1차항의 방향이 전체 통합자료와 동일하게 나타났다.