Applying Nonlinear Mixed-effects Models to Taper Equations: A Case Study of Pinus densiflora in Gangwon Province, Republic of Korea

비선형 혼합효과 모형의 수간곡선 적용: 강원지방 소나무를 대상으로

Abstract

In this study, the performance of a nonlinear mixed-effects (NLME) model used to estimate the stem taper of Pinus densiflora in Gangwon Province was compared with that of a nonlinear fixed-effects (NLFE) model using several performance measures. For the diameters of whole tree stems, the NLME model improved on the performance of the NLFE model by 26.4%, 42.9%, 43.1%, and 0.9% in terms of BIAS, MAB, RMSE, and FI, respectively. For the cross-section areas of whole tree stems, the NLME model improved on the performance of the NLFE model by 67.7%, 44.7%, 45.8%, and 1.0% in terms of BIAS, MAB, RMSE, and FI, respectively. Based on the analysis of 12 relative height classes of tree stems, stem taper estimation performance was also reasonably improved by the NLME model, which showed better MAB, RMSE, and FI at every relative height class compared with those of the NLFE model. In some classes, the NLFE model had better BIAS than the NLME model (stem diameter: 0.05, 0.2, 0.3, and 0.8; stem cross-section area: 0.05, 0.3, 0.5, 0.6, and 1.0). However, the NLME model enhanced the performance of stem diameter and cross-section area estimations at the lowest stem part (0.2 m from the ground). Improvements for stem diameter in terms of BIAS, MAB, RMSE, and FI were 84.2%, 69.8%, 68.7%, and 3.1%, respectively. For stem cross-section areas, the improvements in BIAS, MAB, RMSE, and FI were 98.5%, 70.1%, 68.7%, and 3.1%, respectively. The cross-section area at 0.2 m from the ground occupied 22.7% of total cross-section area. Improvements in estimation of cross-section area at the lowest stem part indicate that stem volume estimation performance could also be enhanced. Although NLME models are more difficult to fit than NLFE models, the use of NLME models as a standard method for the estimating the parameters of stem taper equations should be considered.

강원지방소나무의 수간곡선 추정에 비선형 혼합효과 모형(nonlinear mixed-effects models: NLME)을 적용하고 몇 가지 성능 척도를 이용하여 비선형 고정효과 모형과 비교하였다. 혼합효과 모형이 고정효과 모형을 개선한 정도를 전체 성능의 측면에서 설명하면, 수간직경에 대해서는 BIAS 26.4%, MAB 42.9%, RMSE 43.1%, FI 0.9%만큼이었고, 수간단면적에 대해서는 BIAS 67.7%, MAB 44.7%, RMSE 45.8%, FI 1.0%만큼이었다. 수간을 12개의 상대수고 구간으로 세분화하여 분석한 결과에서도 수간곡선의 성능은 혼합효과 모형에 의해 크게 개선된 것으로 나타났다. 혼합효과 모형은 모든 상대수고 구간에서 수간직경 및 수간단면적에 대한 성능이 고정효과 모형보다 더 나은 MAB, RMSE, FI를 나타내었고, BIAS의 경우 일부 구간(수간직경: 0.05, 0.2, 0.3, 0.8, 수간단면적: 0.05, 0.3, 0.5, 0.6, 1.0)에서만 고정효과 모형보다 뒤떨어지는 것으로 확인되었다. 특히 지상에 근접한 수간 최하부(수고 0.2 m 지점)에서 수간직경 및 수간단면적 추정 성능이 혼합효과 모형에 의해 크게 향상되었다. 수간직경의 경우 BIAS 84.2%, MAB 69.8%, RMSE 68.7%, FI 3.1%, 수간단면적의 경우 BIAS 98.5%, MAB 70.1%, RMSE 68.7%, FI 3.1%만큼 향상된 것으로 분석되었다. 지상으로부터 0.2 m 높이 지점의 수간단면적이 수간단면적 전체에서 차지하는 비중은 22.7%에 달하였다. 이렇게 수간재적에서 가장 큰 비중을 차지하는 수간 최하부에서의 추정 성능이 크게 향상되었다는 것은 전체 수간재적의 추정 성능 역시 큰 폭으로 향상될 수 있음을 시사한다. 비록 모형 적합 과정이 고정효과 모형보다 까다롭다는 단점이 있지만, 추정 성능의 개선 효과를 고려하면 NLME를 수간곡선 추정의 표준적인 방법으로 사용하는 것을 검토할 필요가 있다.

서론

최근 기후변화에 따라 발생하고 있는 여러 가지 문제들은 인류사회가 탄소 중립을 달성해야만 하는 이유를 말해준다. 탄소 중립이라는 시대적 과제는 당장 변화와 행동을 요구하고 있으며, 이는 산림 분야도 예외는 아니다. 산림은 탄소흡수원으로서 그 가치를 인정받고 있으며, 탄소 중립 시나리오에서도 대기 중 탄소 저감의 수단으로 유용하게 사용될 수 있다. 따라서 산림의 탄소흡수량을 정확하게 추정하는 일은 산림자원의 측면에서 뿐만이 아니라 탄소 중립 계획에 매우 기초적이며 중요한 일이다. 산림의 탄소 흡수량은 산림에 존재하는 입목의 재적을 추정하고, 이 추정된 재적을 바이오매스 및 탄소 전환 계수를 이용하여 탄소 저장량으로 전환하게 된다. 이를 바탕으로 두 시점간의 탄소저장량 차이로서 탄소흡수량을 산출하게 된다. 바이오매스 및 탄소 전환 계수는 고정된 상수임을 고려하면, 산림 탄소흡수량의 정확성은 입목의 재적량을 얼마나 정확하게 추정하느냐에 큰 영향을 받는다고 할 수 있다.

국내에서 입목의 재적을 추정하기 위해 사용되는 방법은 수간곡선식(taper equations)이다. 수간곡선식이란 지상으로부터 입목의 초두부에 이르기까지 수간의 직경 변화를 나타내는 함수식으로서, 이를 이용하면 입목의 형상을 입체화할 수 있을 뿐만 아니라 적분을 통해 입목의 재적을 원하는 구간에서 계산할 수 있고 이용수고(merchantable height) 역시 구할 수 있는 등 다양한 장점을 갖는다(Burkhart and Tomé, 2012). 국내에서는 Lee(1994)가 최초로 수간곡선식을 소개한 이래로, 많은 연구가 이루어졌다(Son et al., 2002; Lee et al., 2003; Kang et al., 2014). 산림청과 국립산림과학원은 입목의 수간 재적(stem volume)에 관한 정보를 제공하기 위하여 재적·중량표 및 임분수확표(Korea Forest Service, 2009)를 작성한 이래로 이를 꾸준히 보완하여 입목재적·바이오매스 및 임분수확표(Korea Forest Service and National Institute of Forest Science, 2020)를 발표하기에 이르렀다.

위에서 인용된 입목의 수간 재적 정보를 작성하기 위해 사용된 수간곡선식은 Kozak(1988)이 개발한 변량지수(variable-exponent) 수간곡선식이다. 이는 입목의 부위에 따라 수간곡선의 형태가 변화한다는 가정을 바탕으로 그 형태의 변화에 맞추어 지수를 변화시켜 수간곡선을 정확하게 적합할 목표로 개발되었다. 그런데 모든 수간곡선식은 식을 추정하는데 사용되는 자료들이 위계적(hierarchical)이거나 혹은 그룹화(grouped)되어 있다는 공통점을 가지고 있다. 즉, 한 개체목마다 수간 상의 직경 측정값이 다수 포함되는 수간곡선 자료의 이러한 구조적인 특징은 필연적으로 자기상관(autocorrelation)의 문제를 수반하게 된다. 게다가 Kozak(1988)의 수간곡선식에 사용된 예측변수들에는 총수고, 흉고직경, 상대수고 외에도 이들에 약간의 변형을 가한 변수들이 추가로 포함되어 있어서 다중공선성(multicollinearity)의 가능성 역시 크다. Kozak(1997) 역시 이 문제를 인지하여 관련된 연구를 수행하였다. 해당 연구에 따르면 극심한 다중공선성과 자기상관이 수간곡선식에 존재하더라도 그 예측값은 불편(unbiased)하다고 보고하였으며, 특히 개체목마다 한 개의 측정값만을 사용하여 자기상관을 제거했음에도 예측성능은 향상되지 않았다고 했다.

혼합효과 모형(mixed-effects models)은 한 가지 이상의 분류 요인에 따라 그룹화된 예측변수와 반응변수 사이의 관계를 설명하기 위해 주로 사용된다(Pinheiro and Bates, 2000). 혼합효과 모형은 같은 수준의 분류인자를 공유하는 관측값들을 공통의 임의효과(random effect)로 연관시킴으로써 자료의 그룹화(혹은 위계적 군집화)로 초래된 공분산 구조를 유연하게 표현할 수 있다. 이러한 특징을 가진 혼합효과모형은 수간곡선 자료에 매우 적절한 적용 사례가 될 수 있으며, 실제로 비선형 혼합효과 모형(nonlinear mixed-effects models: NLME)을 적용하여 수간 곡선식의 예측성능을 향상시킨 사례가 다수 있다(Garber and Maguire, 2003; Trincado and Burkhart, 2006; Li et al., 2012; Özçelik, et al., 2019). 최근 수간곡선식의 개발 역사와 활용을 다뤘던 McTague and Weiskittel(2021)은 지난 10년에 걸쳐 NLME가 수간곡선식을 적합하는 대세적인 방법이 되었다고 하였다.

그런데 국내에서 수행된 수간곡선 관련 연구 중에서 NLME를 적용한 사례는 아직 없다. 이것이 의미하는 바는 국내에서 제작되어 사용되는 수간의 재적 정보는 수간곡선 자료의 구조에서 유래하는 공분산을 설명하지 않은 채 추정된 수간직경을 통해 작성되었다는 것이다. 그러므로 수간곡선 자료의 특성을 더욱 잘 표현할 수 있는 NLME를 적용하여 수간곡선을 추정하면 더 정확한 수간직경값을 추정할 가능성이 있다. 이를 통하여 더 정확한 수간재적을 추정하고 나아가 산림 탄소흡수량의 추정값이 더욱 정확해질 수 있다는 것이다. 본 연구의 목적은 국내의 대표수종 중 하나인 강원지방소나무를 대상으로 하여 추정된 Kozak(1988)의 수간곡선식에 NLME를 적용하고 그 성능의 변화를 분석하는 것이다.

재료 및 방법

1. 자료

본 연구는 입목재적·바이오매스 및 임분수확표(Korea Forest Service and National Institute of Forest Science, 2020)에서 강원지방소나무의 입목 수간재적표의 작성에 이용한 표본목 600본을 대상으로 하였다. 표본목마다 총수고, 흉고(1.2 m)직경 및 수간 여러 지점의 수피포함 직경(diameter outside bark: DOB)을 측정하였다. 수간직경을 측정한 지점은 흉고를 포함하여 흉고 아래에서는 지상으로부터 0.2, 0.7 m 높이 지점에서, 흉고 위에서는 대부분 2 m 간격(3.2, 5.2, 7.2, ⋯, 25.2) 지점에서 수피포함 직경을 측정하였다. 표본목과 수간직경 측정치의 기술 통계량은 Table 1과 같다.

Table 1. Descriptive statistics of sample trees and diameter measurements.

^a SD: standard deviation.

2. 수간곡선식

분석에 사용된 수간곡선식은 Kozak(1988)의 변량지수 수간곡선식으로서 그 형태는 아래의 식 1과 같다.

\(\begin{aligned}\hat{d}_{i}=\beta_{0} D^{\beta_{1}} \beta_{2}^{D} X_{i}^{\beta_{3} Z^{2}+\beta_{4} \ln (Z+0.001)+\beta_{5} \sqrt{Z}+\beta_{6} e^{Z}+\beta_{7}\left(\frac{D}{H}\right)}\end{aligned}\)       (1)

여기서 \(\begin{aligned}\hat{d}_{i}\end{aligned}\)는 h_i에서 수피제외 직경 추정치, h_i는 지상으로부터 높이(0≤ h_i ≤H), H는 총수고, HI는 지상으로부터 변곡점의 높이, \(\begin{aligned}Z=\frac{h_{i}}{H}\end{aligned}\), \(\begin{aligned}X=\frac{(1-\sqrt{Z})}{(1-\sqrt{p})}\end{aligned}\), \(\begin{aligned}p=\left(\frac{H I}{H}\right) \times 100\end{aligned}\), D는 흉고직경(수피포함), p는 0.1∼0.3 사이의 값을 갖는 변곡점의 상대높이이며 (본 연구에서 강원지방소나무의 p는 0.3), β_i는 매개변수이다(직경의 단위는 cm, 수고의 단위는 m).

하지만, 본 연구에서 사용된 자료는 수피를 제외한 수간직경(diameter inside bark: DIB)이 아니라 DOB이다. 이것은 Kozak(1988) 수간곡선식의 개발 조건에는 맞지 않으나 일반적으로 DOB와 DIB 사이에는 선형관계가 있으므로 DOB를 예측하기 위하여 해당 수간곡선식을 사용하였다(Li et al., 2012).

3. 비선형 혼합효과 모형

고정효과(fixed effects)는 전체 모집단 혹은 실험요인의 반복 가능한 수준과 연관된 매개변수(parameters)이며, 임의효과(random effects)는 모집단에서 무작위로 선발된 개별 실험 단위와 관련되어 있다(Pinheiro and Bates, 2000). 본 연구의 수간곡선 추정에 관하여 다시 설명하면, 고정효과는 강원지방소나무 모집단의 일반적 경향을, 임의효과는 각 표본목의 변이를 설명한다고 할 수 있다. 이 두 효과를 함께 포함하는 모형을 혼합효과 모형이라 하며, 이는 모집단의 일반적 경향에 대한 설명력을 잃지 않으면서도 표본목 간의 작은 차이를 설명하는 유연함을 제공한다(Li and Weiskitell, 2010).

고정효과와 임의효과의 일부 혹은 전부가 모형함수에서 비선형적으로 나타나는 혼합효과 모형을 NLME라 한다(Pinheiro and Bates, 2000). 어떤 높이 값을 갖는 수간 지점에서 그 높이 값과 해당 지점의 수간직경 사이의 관계를 추정하는 수간곡선은 일반적으로 비선형의 형태를 취한다. 그리고 수간곡선 자료는 한 개체목에서 여러 수간직경 값을 측정하는 구조적인 형태를 보이기 때문에 수간곡선의 추정에 NLME는 매우 적합한 방법이라 할 수 있다. 그룹이 단일 수준인 경우(single-level of grouping), 비선형 혼합효과 모형은 다음과 같은 형태를 취한다(Pinheiro and Bates, 2000).

y_ij = f(𝜙_ij, v_ij) + 𝜖_ij, i = 1, ⋯, M, j = 1, ⋯, n_i       (2)

𝜙_ij = A_ijβ + B_ijb_i, b_i ∼ N(0, ψ)       (3)

여기서 y_ij는 i번째 그룹의 j번째 관측값, M은 그룹의 개수, n_i는 i번째 그룹의 관측값 개수, f는 그룹 특정적(group-specific) 매개변수 벡터 𝜙_ij와 공변량(예측변수) 벡터 v_ij의 미분 가능한 실가 일반 함수(a general, real-valued, differential function), 𝜖_ij는 정규분포하는 그룹 내(withingroup) 오차이며, 함수 f는 𝜙_ij 중 최소 하나의 성분에서 비선형이다. β는 고정효과의 p차원 벡터이고, b_i는 i번째 그룹과 연관된 임의효과의 q차원 벡터로서 그 공분산 행렬로 ψ를 갖는다. A_ij와 B_ij는 적절한 차원을 가진 행렬로서 그룹과 j번째 관측값의 공변량 값에 의해 결정된다. 그리고 서로 다른 그룹에 해당하는 관측값들은 서로 독립이며, 𝜖_ij는 N(0, σ²)을 따라 서로 독립적으로 분포하며 b_i와도 독립인 것으로 가정한다.

혼합효과 수간곡선 모형의 𝜖_ij에서 확인되는 이분산성, 즉 직경이나 수고가 증가함에 따라 𝜖_ij의 분산 역시 증가하는 경향이 여러 관련 문헌에서 보고되었다(Trincado and Burkhart, 2006; Valentine and Gregoire, 2001). 이러한 𝜖_ij의 이분산적 분산구조를 설명하기 위해서 아래의 분산함수가 고려되었다(Pinheiro and Bates, 2000; Li and Weiskitell, 2010).

Var(𝜖_i)=σ²|v_i|^2δ       (4)

여기서 σ²는 잔차제곱합(residual sum of squares), υ_i는 가중치 변수(본 연구에서는 상대수고 값을 사용), δ는 분산함수의 매개변수이다.

한편 Garber and Magiure(2003)는 수간곡선 모형에 임의효과를 추가하는 것만으로는 오차의 자기상관을 완전히 제거할 수 없었으며, 1차 연속 자기상관(a first-order continuous autoregressive: CAR(1)) 오차 구조를 모형에 추가하여 자기상관을 제거할 수 있었다고 보고하였다. 본 연구에서도 𝜖_ij 사이의 상관관계, 즉 자기상관을 설명하기 위해 CAR(1) 오차 구조를 모형 추정에 도입하였다. CAR(1)의 구조는 아래와 같다(Pinheiro and Bates, 2000; Li and Weiskitell, 2010):

Corr(𝜖_t, 𝜖_s) = θ^{|t - s|}       (5)

여기서 𝜖_t와 𝜖_s는 두 관측값에 대한 모형의 오차, θ는 한 단위만큼 떨어진 두 개의 관측값 사이의 상관관계, |t - s|은 두 관측값 사이의 높이 차이이다.

요약하면 임의효과를 추가한 혼합효과 모형을 적용하여 수간직경의 전체적인 경향과 더불어 각 개체목 내의 수간직경 변이를 설명하고자 하였으며, 모형의 성능을 더 개선하기 위해 분산함수로 개체목 내 이분산성(within-tree variance heterogeneity)을, CAR(1) 오차 구조로 개체목 내 상관관계(within-tree correlation)를 모형에 포함하고자 하였다.

4. 모형 평가 기준

추정된 모형의 수간직경 추정능력을 평가 및 비교하고자 몇 가지 성능 척도를 사용하였다: 평균편향(mean bias: BIAS), 평균편향 백분율(BIAS%), 평균절대편향(mean absolute bias: MAB), 평균절대편향 백분율(MAB%), 평균제곱근오차(root mean square error: RMSE), 평균제곱근오차 백분율(RMSE %), 적합도지수(fitness index: FI).

\(\begin{aligned}B I A S=\frac{1}{n} \sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right), \quad B I A S \%=\frac{B I A S}{\bar{y}} \times 100\end{aligned}\)         (6)

\(\begin{aligned}M A B=\frac{1}{n} \sum_{i=1}^{n}\left|y_{i}-\hat{y}_{i}\right|, M A B \%=\frac{M A B}{\bar{y}} \times 100\end{aligned}\)       (7)

\(\begin{aligned}\mathrm{SE}=\sqrt{\frac{1}{\mathrm{n}} \sum_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{y}_{\mathrm{i}}-\hat{\mathrm{y}}_{\mathrm{i}}\right)^{2}}, \mathrm{SE} \%=\frac{\mathrm{SE}}{\overline{\mathrm{y}}} \times 100\end{aligned}\)       (8)

\(\begin{aligned}F I=1-\left[\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2}\right] /\left[\sum_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}\right]\end{aligned}\)       (9)

여기서 y_i는 i번째 수간직경의 관측값이며, \(\begin{aligned}\hat{y}_{i}\end{aligned}\)은 i번째 수간직경의 예측값, \(\begin{aligned}\bar{y}=\frac{1}{n} \sum_{i=1}^{n} y_{i}\end{aligned}\), n은 수간직경 관측값의 개수이다.

수간의 재적을 산출하는 것은 수간직경 추정의 주요한 목적 중 하나이기 때문에, 수간의 재적과 비례하는 수간 단면적의 추정능력도 함께 평가하였다(수간단면적의 단위는 cm²). 수간단면적 추정의 성능 척도는 식 6～9와 동일하지만, y_i 대신 \(\begin{aligned}y_{i}^{*}=\frac{\pi}{4} y_{i}^{2}\end{aligned}\), \(\begin{aligned}\hat{y}_{i}\end{aligned}\) 대신 , \(\begin{aligned}\hat{y}_{i}^{*}=\frac{\pi}{4} \hat{y}_{i}^{2}\end{aligned}\), \(\begin{aligned}\bar y\end{aligned}\) 대신 \(\begin{aligned}\bar{y}^{*}=\frac{\pi}{4 n} \sum_{i=1}^{n} y_{i}^{2}\end{aligned}\)를 사용하였다. 하지만, 이처럼 추정한 수간단 면적의 추정능력이 수간재적의 추정능력을 편향적으로 반영한다는 점에 유의해야 한다. 본 연구에서 사용한 수간곡선식(Kozak, 1988)은 수간직경을 추정하는 식이지 수간직경의 제곱이나 수간단면적을 추정하는 식이 아니다. Mehtätalo and Lappi(2020)는 수간직경의 추정값을 제곱한 것은 수간직경 제곱의 불편추정량(unbiased estimator)이 아니며, 이는 아래의 사실에서 확인된다고 하였다.

E(X²) = [E(X)]² + Var(X)       (10)

여기서 X는 수간직경이다.

Burkhart(2016)가 지적하였듯이 수간직경을 불편추정하기 위해서는 수간곡선의 종속변수가 수간직경, 수간단면적이나 수간재적을 불편추정하기 위해서는 수간곡선의 종속변수가 수간직경의 제곱 혹은 수간단면적이어야만 한다. 그렇지 않다면 종속변수의 변수변환에 따른 오차를 보정해 주어야 한다. 비록 본 연구에서 수간단면적 분석은 이러한 종속변수 변환에 따른 편향을 보정하지 않았으나, Czaplewski and Bruce(1990)가 제시한 변수 변환에 따른 편향의 산점도에서 보여준 것과 같은 어떤 특정한 패턴이 본 연구의 수간단면적 잔차도에서는 발견되지 않았다(Figure 2). 이는 아마도 변수변환에 따른 편향이 상대적으로 크지 않기 때문으로 판단되며, 본 연구의 목적과 이러한 상황을 고려하여 편향된 수간단면적 추정치를 분석에 사용하였다.

수간곡선 추정의 국부적인 성능을 살펴보기 위하여 잔차의 산점도(잔차도)와 상대수고 구간별 성능 척도를 사용하였다. Brack(1999)은 수간곡선을 추정하는 여러 방법은 전체적인 편향이 작더라도 국부적으로는 큰 편향이 존재할 수 있다고 하였다. Son(2000)도 수간의 1/3 이하 지점이 수간에서 차지하는 비율이 상대적으로 크고 그 가치 또한 높다는 것을 지적하며, 수간곡선 모형을 다룰 때는 수간 하부의 편향에 주목해야 한다고 주장하였다.

국내 수종을 대상으로 수간곡선의 성능을 국부적으로 살펴본 연구는 다수 존재한다(Lee et al., 2003; Doyog et al., 2017). 특히 이들은 상대수고를 0.1 단위로 총 10개 구간으로 나누어 성능을 비교하였다. 하지만 본 연구에서 가장 큰 수간직경을 갖는 흉고 아래 지점들(0.2, 0.7 m 높이)의 상대수고 값은 대부분 0.05 이하의 값을 갖고 있어서(Table 2), 위와 같은 상대수고 구획으로는 수간 하부에서 추정의 성능을 제대로 구분하여 살펴보기 어려웠다. 따라서 상대수고 0.1 이하의 구간을 더 세분화하고자 하였다. Table 2에 따르면 0.2 m 높이에 있는 수간 지점의 상대수고는 최대 0.025를 넘지 않았으며, 다른 어떤 측정지점의 상대수고도 0.025 이하에 위치하지 않았다. 0.7 m 높이지점은 상대수고 최솟값이 0.0271에서 최대 0.0854였고 평균은 0.0445이다. 이러한 상대수고의 분포를 고려하여 상대수고 0.1 이하의 구간을 0 초과 0.025 이하, 0.025 초과 0.05 이하, 0.05 초과 0.1 이하 총 3개 구간으로 세분화하고, 나머지 구간은 0.1 간격으로 하여 총 12개 상대수고 구간을 설정하여 분석하였다. 이런 상대수고 구간 설정은 본 연구에서 사용된 수간곡선식 개발 연구(Kozak, 1988)에서 흉고 이하 수간 지점에 대하여 0.5 m 이하 구간과 0.5∼1.3 m 사이 구간을 설정한 것과 유사하다고 할 수 있다. Kozak(1988)은 구간 설정 원리에 대해서 따로 설명하지는 않았지만, 수간 하부에 이런 구간을 추가로 설정해 놓은 것으로 보아 수간 하부에서의 성능이 매우 중요함을 파악한 것으로 판단된다.

Table 2. Distribution of relative height for measurement locations sorted by absolute height around breast height on tree stem.

5. 모형의 추정

비선형 혼합효과 모형을 추정하는 데 있어 중요한 부분은 어떤 매개변수에 임의효과를 추가할지를 결정하는 것이다. 본 연구에서 사용한 수간곡선식에는 β₀에서 β₇까지 총 8개의 매개변수가 포함되어 있는데, 이 중 최대 3개의 매개변수를 포함하는 것을 조건으로 고정효과 모형을 가장 크게 개선하는 최상의 조합을 찾기 위해 가능한 모든 조합에 대하여 분석하였으며, 그 기준으로 AIC(Akaike Information Criterion)를 사용하였다.

선형모형과 달리 비선형 모형을 추정하는 과정에서 흔히 겪게 되는 어려움은 바로 수렴 문제(convergence problems)이다. 특히 비선형 혼합효과 모형에서 임의효과가 포함되는 경우 이 수렴 문제는 흔히 발생한다. 이러한 수렴 문제에 대처하는 방법에는 매개변수 추정의 반복 횟수 증가, 새로운 매개변수 초깃값 제공, 수렴기준 완화, 모형의 단순화(임의효과의 일부 공분산을 고정하거나 일부 임의효과를 제거) 등이 있다(Mehtätalo and Lappi, 2020). 또 Pinheiro and Bates(2000)는 임의효과의 개수가 관측대상의 개수보다 상대적으로 큰 경우에는 과대 매개변수화(overparameterized)에 따른 수렴 문제를 피하고자 임의효과 벡터의 공분산행렬 ψ를 대각행렬로 하는 것(임의효과들이 서로 독립적이라고 조정)을 일반적으로 추천한다고 하였다.

최적의 임의효과 조합을 찾는 본 연구의 분석과정에서도 위에서 언급한 수렴 문제는 자주 발생하였으며, 수렴문제가 발생한 모형에 대해서는 위에서 언급한 방법들을 모두 적용해 보았다. 하지만, 해당 방법들을 적용하였음에도 불구하고 수렴 문제를 해결하지 못하는 모형이 다수 나타났다(수렴문제가 해결되지 않은 모형은 임의효과 매개변수 선발에서 제외하였다). 특히 임의효과 4개 이상부터는 수렴문제가 발생하지 않는 조합을 찾기 어려웠으며, 개별모형의 추정 시간 역시 매우 길어졌기 때문에 임의효과의 개수를 최대 3개로 제한하였다. 또 본 연구에서 사용된 수간곡선식은 아니지만 동일 저자가 개량·수정한 수간곡선식(Kozak, 2004)에 NLME를 적용한 연구들에서는 2개의 임의효과를 고려한 사실 역시 고려하였다(Garber and Maguire, 2003; Li and Weiskitell, 2010; Li et al., 2012; Poudel et al., 2018). 본 연구와 관련된 모든 분석은 통계분석 프로그램인 R version 4.0.5(R Core Team, 2021) 및 R 통합개발환경 RStudio version 1.3.1093(RStudio Team, 2020)에서 수행하였다. 혼합효과 모형의 분석을 위해 R 패키지 nlme(Pinheiro et al., 2021)를, 그림의 출력을 위해 R 패키지 ggplot2(Wickham, 2016), cowplot(Wilke, 2020), dplyr(Wickham et al., 2021) 및 함수 ggplot.corr(Liu, 2018)를 사용하였다.

결과 및 고찰

1. 혼합효과 모형 선발

총 8개의 매개변수 중 최대 3개의 매개변수에 대해 임의 효과를 추가하는 모든 경우에 대하여 모형을 적합하였다. 매개변수 추정이 불가하거나 추정 시 수렴 경고가 발생한 경우를 제외한 모형 중에서 최상의 개선 효과를 나타낸 5개 모형의 AIC와 수간직경에 대한 성능 척도(백분율 척도 제외)를 Table 3에 표시하였다. 그중에서 매개변수 β₁, β₅, β₆에 임의효과를 추가한 모형이 가장 낮은 AIC 값과 더불어 가장 좋은 성능 척도를 나타내었기에, 이 모형을 혼합효과 모형으로 선발하였다.

Table 3. AIC and some performance measures of the taper equation for the five best combinations of random-effects parameters.

^a Fixed-effects model

2. 고정효과 모형과 선발된 혼합효과 모형의 성능 비교

1) 전체 성능 비교

모형 선발 이후, CAR(1) 오차 구조와 분산함수를 추가하여 모형의 매개변수를 추정하였는데, 분산함수를 적용한 모형과 CAR(1) 및 분산구조를 동시에 적용한 모형은 수렴문제가 발생하여 분석에서 제외하였다. 최종적으로 고정효과 모형, 기본 혼합효과 모형 및 CAR(1)을 적용한 혼합효과 모형 총 3개 모형에 대해 비교를 수행하였으며, 각 모형의 매개변수 추정치는 Table 4와 같다.

Table 4. Parameter estimates and their standard error for Kozak(1988) taper model.

^a θ is the estimate of first-order continuous autoregressive parameter from the mixed-effects model with CAR(1).

^b Statistically not significant at α = 0.05.

선정된 모형들의 추정 성능을 분석한 결과 수간직경 추정에서는 기본 혼합효과 모형이 모든 성능 척도에서 가장 우수한 것으로 확인되었다(Table 5). 기본 혼합효과 모형은 고정효과 모형을 BIAS 26.4%, MAB 42.9%, RMSE 43.1%, FI 0.9%만큼 개선하였으며, CAR(1) 혼합효과 모형은 같은 척도에 대해 각각 –268.2%, 11.9%, 13.2%, 0.3%를 개선하였다. CAR(1) 모형이 BIAS에서 성능이 더 떨어지는 결과를 보였으나 실제 그 차이는 –0.0295에 불과했다. 수간단면적 추정에서도 기본 혼합효과 모형이 가장 우수하였으며, 고정효과 모형 대비 개선의 정도는 BIAS 67.7%, MAB 44.7%, RMSE 45.8%, FI 1.0%였다(Table 5).

Table 5. Performance measures for diameter and cross-section area by model type.

^a The number in the parenthesis represents the percentage by which the corresponding model improves the fixed-effects one.

^b Biased determined by t-test at α = 0.05.

유의수준 0.05에서 BIAS에 대해 실시한 t-검정 결과, 수간직경에서는 CAR(1) 혼합효과 모형이 수간단면적에서는 고정효과 모형이 편향된 것을 확인하였다. 고정효과 모형으로 추정한 수간단면적이 편향된 것은 시사하는 바가큰데, 왜냐하면 입목재적·바이오매스 및 임분수확표(Korea Forest Service and National Institute of Forest Science, 2020)에서 입목수간재적표를 작성하는데 사용된 추정 방법이 위에서 설명한 편향이 발생하는 방법과 같기 때문이다. 이에 따라 추정된 수간재적의 정확도도 낮아지게 되므로 현재 사용 중인 입목수간재적표의 편향성에 대한 검토가 필요할 것으로 판단된다.

CAR(1) 오차 구조를 적용한 모형의 성능이 기본 혼합효과 모형의 성능보다 떨어졌는데, 이것에 대해서는 아래에서 다룰 자기상관 부분에서 더 자세히 언급하지만, CAR(1) 구조가 개체목 내 오차의 상관관계를 적절하게 설명하지 못하는 것이 아닌가 추측된다. 하지만 이에 대한 정확한 원인은 추가 연구를 통해서 다룰 수 있을 것이다. 또 Mehtätalo and Lappi(2020)는 임의효과를 추가한 혼합효과 모형을 이용하여 잔차의 변동을 설명한 이후 남아있는 오차는 아주 작아지므로, 오차를 설명하기 위한 상관모형은 크게 중요하지 않다고도 하였다.

2) 수간고별 성능 비교

기존의 고정효과 모형을 가장 크게 개선하도록 임의효과를 부여할 매개변수를 결정한 후 오차의 분산과 자기상관을 추가로 설명하도록 분산함수와 오차 구조를 적용하였으나 이를 통해 기본 혼합효과 모형의 성능이 개선되지 않았기 때문에, 수간고별 성능 비교에서는 고정효과 모형과 기본 혼합효과 모형만을 비교하였다.

앞서 설명한 대로 설정된 12개의 상대수고 구간에서 각 모형의 성능 척도를 계산한 결과를 Table 6에 표시하였다. 기본 혼합효과 모형은 수간직경과 수간단면적 추정 모두 전체 상대수고 구간에 걸쳐 고정효과 모형보다 더 나은 MAB, RMSE, FI를 나타내었으며, BIAS의 일부 상대수고구간(수간직경: 0.05, 0.2, 0.3, 0.8; 수간단면적: 0.05, 0.3, 0.5, 0.6, 1.0)에서만 고정효과 모형보다 뒤떨어지는 것으로 확인되었다. 특히 고정효과 모형을 개선한 측면에서 가장 큰 효과를 나타낸 구간은 0.025 구간, 즉 상대수고가 0.025 이하인 수간 최하부 지점(수고 0.2 m)이었다. 수간직경의 경우 BIAS 84.2%, MAB 69.8%, RMSE 68.7%, FI 3.1%의 개선을 가져왔다. 수간단면적도 기본 혼합효과 모형은 수간 최하부에서 추정 성능을 BIAS 98.5%, MAB 70.1%, RMSE 68.7%, FI 3.1%만큼 개선하였다.

Table 6. Performance measures for diameter and cross-section area by model type and relative height class.

(Note) The blue number in bold indicates that the relative height class including that blue number has the poorest value in the corresponding performance measure. Likewise, the red number in bold indicates that the relative height class including that red number has the best value in the corresponding performance measure.

^a Relative height classes: 0.025 = (0, 0.025], 0.05 = (0.025, 0.05], 0.1 = (0.05, 0.1], 0.2 = (0.1, 0.2], …, 0.9 = (0.8, 0.9], 1.0 = (0.9, 1.0]

^b The number in the parenthesis represents the percentage by which the mixed-effects model improves the fixed-effects one.

FI는 고정효과 모형에서는 지상에서 흉고까지는 약간 증가하다가 흉고 이상 지점에서부터는 초두부에 가까워질수록 지속해서 감소하였다. 기본 혼합효과 모형에서도 대체로 비슷한 경향을 보였으나 가장 아래쪽에 있는 상대수고 구간 0.025에서 가장 높은 FI 값을 나타낸다는 점에서 차이를 보였다. 고정효과 모형에서는 FI가 상대수고 구간 0.6(수간직경)이나 0.7(수간단면적)에서부터 0.9 이하의 값으로 내려갔지만 기본 혼합효과 모형에서는 상대수고 구간 1.0을 제외한 모든 구간에서 0.9 이상의 값을 유지하였다. 또 상대수고 구간 1.0에서 고정효과 모형은 약 0.6의 FI 값을 보였지만 기본 혼합효과 모형은 0.7 이상의 값을 나타내었다.

기본 혼합효과 모형의 수간직경과 수간단면적 추정 성능에서 주목할 것은 MAB와 RMSE에서 최고 및 최저 성능을 나타낸 구간이 상반된다는 점이다. 수간직경의 경우에는 가장 좋은 성능은 구간 0.025, 가장 뒤떨어지는 성능은 구간 1.0에서 관찰되었다. 수간단면적의 경우에는 그 반대로서 가장 좋은 성능은 구간 1.0, 가장 뒤떨어지는 성능은 구간 0.05에서 나타났다. 고정효과 모형에서는 두 변수에서 가장 뒤떨어지는 성능을 보인 구간은 0.025로 같았고, 가장 뛰어난 성능을 보인 구간은 수간직경에서는 0.1, 수간단면적에서는 1.0으로 혼합효과 모형과 같은 방식은 아니지만 두 변수의 평가 사이에 차이점이 확인되었다. 비록 본 연구에서는 수간직경과 수간단면적 양쪽 모두에서 혼합효과 모형이 가장 우수한 것으로 나타났지만, 수간곡선의 추정 성능이 수간직경과 수간단면적 사이에서 이런 변동성을 나타내었다는 것은 시사하는 바가 적지 않다. 추정 변수에 따른 추정 성능의 이러한 변동성은 사용하는 수간곡선 모형이나 대상 수종에 따라서도 크게 영향을 받을 가능성이 있다고 판단되며, 이는 곧 수간직경과 수간재적의 추정 평가는 독립적으로 이루어져야 한다는 Burkhart (2016)의 주장을 더욱 뒷받침한다.

백분율 성능 척도 측면에서 기본 혼합효과 모형은 수간직경과 수간단면적에서 모두 가장 뛰어난 성능을 상대수고 구간 0.025, 즉 수간 최하단부에서 나타내었고, 가장 뒤떨어지는 성능을 수간 최상단부인 상대수고 구간 1.0에서 나타내었다. 이는 고정효과 모형에서도 유사하게 나타났는데, 가장 뒤떨어지는 성능이 상대수고 구간 1.0에서 나타난 것은 같았고, 성능이 가장 뛰어난 곳은 역시 수간하부라고 간주할 수 있는 상대수고 구간 0.1에 집중되었다.

위에서 설명한 혼합효과 모형의 개선 효과는 잔차도와 잔차의 상자그림을 통해서 시각적으로 확인할 수 있다. Figure 1(수간직경 잔차도)과 2(수간단면적 잔차도)를 살펴보면, 전체 상대수고 구간에서 혼합효과 모형의 잔차가 고정효과 모형의 잔차보다 0 근처에 더 집중적으로 분포하고 있음을 확인할 수 있다. Figure 1과 2는 수간 하부에서 실제 수고에 따른 혼합효과 모형의 개선 효과를 살펴보기 위해 측정지점의 실제 높이(0.2 m, 0.7 m, 1.2 m)에 따라 잔차를 구분하여 출력하고 있다. 1.2 m에서는 눈에 띄는 분산의 감소를 확인하기 어려웠으나, 0.2 m(상대수고구간 0.025)와 0.7 m에서는 잔차의 분산 감소를 쉽게 확인할 수 있었다.

Figure 1. Residual plots for diameter estimation along relative height by model type.

The color of residual point represents the absolute height from ground of the measurement location along tree stem. The horizontal blue lines of above and below indicate residual values 2 and -2, respectively.

Figure 2. Residual plots for cross-section area estimation along relative height by model type.

The color of residual point represents the absolute height from ground of the measurement location along tree stem. The horizontal blue lines of above and below represent residual values 100 and -100, respectively.

상대수고 구간별 성능 척도와 잔차도에서 혼합효과 모형에 의한 잔차의 분산 감소를 확인하였으나 잔차의 분포를 좀 더 상세히 살펴보기 위해 잔차의 구간별 상자그림을 출력하였다(Figure 3, 4). 대부분 구간에서 잔차의 전체 범위뿐만 아니라 사분위수 범위가 줄어든 것을 확인할 수 있었으며, 0.025 구간에서는 수간직경과 수간단면적 모두에서 그 감소의 폭이 역시 매우 큰 것으로 나타났다. 그리고 수간직경의 상대수고 구간 1.0에서 혼합효과 모형의 이상점이 증가하여 전체적인 범위가 고정효과 모형보다 넓어진 것이 관찰되었다.

Figure 3. Boxplots of residuals for diameter estimation along relative height classa by model type.

The white box represents the interquartile range, and the horizontal line in the box the median.

^a Relative height classes: 0.025 = (0, 0.025], 0.05 = (0.025, 0.05], 0.1 = (0.05, 0.1], 0.2 = (0.1, 0.2], …, 0.9 = (0.8, 0.9], 1.0 = (0.9, 1.0]

Figure 4. Boxplots of residuals for cross-section area estimation along relative height classa by model type.

The white box represents the interquartile range, and the horizontal line in the box the median.

^a Relative height classes: 0.025 = (0, 0.025], 0.05 = (0.025, 0.05], 0.1 = (0.05, 0.1], 0.2 = (0.1, 0.2], …, 0.9 = (0.8, 0.9], 1.0 = (0.9, 1.0)

이상 상대수고 구간에 따른 수간곡선의 추정 성능을 비교한 결과, 전체 구간에 대한 평가에서는 알 수 없던 사실, 수간재적에서 가장 큰 비중을 차지하는 수간 최하부에서의 추정 성능이 크게 향상된 것을 확인할 수 있었다. 수간하부가 재적에서 차지하는 비율이 어느 정도인지 가늠하기 위해 지상에서 가까운 수간직경 측점 지점 몇 곳을 대상으로 전체 수간단면적에서 차지하는 비율을 계산해 보았다(Table 7). 0.2 m 지점은 전체 수간단면적의 22.7%를, 0.7 m 지점은 18.3%를 차지하고 있었으며, 이 두 지점은 수간단면적에서 합계 41.0%의 비율을 차지하고 있었으며, 이는 수간직경에서 차지하는 비율 31.6%보다 약 10%가량 큰 수치였다. 즉, NLME를 사용하여 수간 하부에서의 추정 성능이 대폭 개선된 효과는 수간직경보다 수간단면적에서 더 크게 나타난다는 것이다. 결국, NLME를 적용하게 되면 수간재적의 추정능력 또한 큰 폭으로 향상될 것으로 판단된다.

Table 7. Percentage of sum of stem diameter and cross-section area at some measurement locations closer to ground.

3) 개체목 수준에서 수간곡선의 비교

앞서 기술한 전체 및 구간별 성능의 향상이 개체목 수준에서 어떻게 반영되는지를 시각적으로 살펴보기 위해 일부 표본목을 대상으로 수간직경 측정값과 고정효과 모형과 혼합효과 모형에 의한 수간직경 추정값을 이용하여 수간곡선을 그려 비교해보았다. 선택된 표본목들은 2그루로서 그 흉고직경이 표본목 중에서 최소(8.7 cm)와 최대(46.8 cm)인 것으로 하였다(Figure 5). 수간곡선을 살펴본 결과, 수간고별 성능 비교에서 기술한 바와 같이 지상부에 근접한 근주부에서 혼합효과 모형이 고정효과 모형에 비하여 실제 수간직경을 더 근사하게 추정하는 것을 양쪽 표본목 모두에서 확인할 수 있었다.

Figure 5. Stem profiles at two selected sample trees, which were drawn through diameters measured, ones estimated by fixed-effects model, and ones estimated by mixed-effects model, respectively.

4) 자기상관

비교 대상 3개 모형의 자기상관함수(autocorrelation function: ACF) 그림을 분석해보면 기본 혼합효과 모형이 고정효과 모형의 자기상관을 어느 정도는 해소하였으나 여전히 자기상관이 남아있는 것으로 확인되었다(Figure 6). CAR(1)을 적용한 혼합효과 모형은 오차의 자기상관을 추가로 설명하기 위한 목적에도 불구하고 고정효과 모형의 자기상관을 거의 해소하지 못한 모습을 보인다. Mehtätalo and Lappi(2020)에 따르면 CAR(1) 오차 구조는 수간직경 사이의 상관관계가 일반적으로 음의 값을 띈다는 사실을 고려하지 않고 있다고 언급한 것을 미루어볼 때, CAR(1) 구조가 강원지방소나무의 개체목 내 오차 간의 상관관계를 설명하는데 부적절한 것이 아닌가 추측하였으나, 이에 대한 정확한 원인은 추가적인 연구가 필요한 것으로 판단된다. 이런 상관관계의 조건에 영향을 받지 않고 오차의 공분산 구조를 설명할 수 있는 준모수적 방법(Kublin et al., 2008)이나 비모수적 방법(Lappi, 2006)이 제시되었으나, 본 연구에서는 고려하지 않았다.

Figure 6. ACF plots for fixed-effects, mixed-effects and mixed-effects with CAR(1) models. The two horizontal dotted lines represent 95% confidence interval for autocorrelation.

결론

본 연구는 수간곡선식 추정에 NLME를 적용할 때 어떤 효과를 얻을 수 있는지를 강원도지방소나무를 대상으로 분석해보았다. 국내에서 가장 널리 사용되는 Kozak(1988)의 수간곡선식에 NLME를 적용하여 비선형 회귀분석보다 더 개선된 수간곡선 추정 결과를 얻을 수 있었다. 특히 수간재적에서 높은 비율을 차지하는 지상부 근처 수간 하부 지점에서 수간직경 및 수간단면적 추정 성능이 크게 개선되었다. 이것은 수간곡선 추정의 주요한 목적 중 하나인 수간재적 추정의 정확도를 크게 향상시킬 것으로 기대된다. 한편 수간직경 추정값을 이용하여 수간단면적을 추정하는 방식은 고정효과 모형에서 수간단면적을 편향 추정(biased estimation)하였으며, 이것은 수간재적의 편향 추정으로 이어지기 때문에, 향후 현재 사용 중인 입목수간재적표의 편향에 대한 면밀한 검토가 필요하다.

임의효과를 부여할 적절한 매개변수를 결정해야 하고, 결정된 매개변수에 임의효과를 부여하여 전체 NLME 모형의 매개변수를 추정하는 과정에서 수렴문제가 발생하는 등 NLME 모형을 적합(model fitting)하는 것이 비선형 회귀분석보다 훨씬 까다롭다는 단점도 존재한다. 하지만 위계적 특성이 있는 수간곡선 자료에 내재한 자기상관성을 줄이는 성질, 수간 하부의 수간곡선 추정능력을 크게 향상시킨 연구 결과, 그리고 지난 10년에 걸쳐 NLME가 수간곡선식을 적합하는 대세적인 방법이 되었다는 사실(McTague and Weiskittel, 2021)을 생각할 때 앞으로 수간곡선 추정에 NLME를 표준적으로 적용하는 것을 검토해야 할 것으로 판단된다.