Kwangmin Kim - FDA 5.3~5.4 — 함수-on-함수 회귀와 refund 패키지 통합 구현

1 두 절의 역할

이 포스트의 범위

절	주제	핵심 도구
5.3	함수-on-함수 회귀의 모형·식별성·벌점 추정	이변량 기저 + Riemann sum + 라플라시안 벌점
5.4	refund 패키지: fosr 와 pffr 사용법	OLS/GLS, pffrSim, ff() helper

5.3 은 “반응도 함수, 회귀자도 함수” 인 가장 일반적 함수 회귀 — 회귀계수가 두 시간 축을 잇는 이변량 핵 \(\psi(t, s)\). 5.1~5.2 의 텐서 곱 정칙화가 양방향으로 확장된다.

5.4 는 5.2 와 5.3 의 추정을 R 패키지 refund 가 어떻게 자동화하는지를 보여준다. 두 핵심 함수 — fosr (함수-on-스칼라 전용, OLS/GLS/순열검정) 와 pffr (일반 함수 회귀, mixed model framework) — 의 사용 패턴을 캐나다 기상 데이터와 시뮬레이션으로 익힌다.

두 절을 합쳐서 보면, “이론(5.3) 과 실무 도구(5.4) 의 다리” 가 만들어진다.

2 함수-on-함수 모형의 정의

2.1 가장 단순한 형태

정의: 함수-on-함수 회귀 (단순 형태, 식 5.10)

\(N\) 개의 곡선-곡선 쌍 \(\{(X_i, Y_i): i = 1, \ldots, N\}\) 에 대해:

\[ Y_i(t) = \alpha(t) + \int_{\mathcal{S}} \psi(t, s) X_i(s) \, ds + \varepsilon_i(t). \]

여기서:

\(Y_i: \mathcal{T} \to \mathbb{R}\) — 반응 곡선 (response curve)
\(X_i: \mathcal{S} \to \mathbb{R}\) — 회귀자 곡선 (regressor curve), 도메인 \(\mathcal{S}\) 가 반응의 \(\mathcal{T}\) 와 다를 수 있음
\(\alpha: \mathcal{T} \to \mathbb{R}\) — 절편 함수
\(\psi: \mathcal{T} \times \mathcal{S} \to \mathbb{R}\) — 회귀 핵 (regression kernel), 이변량 함수
\(\varepsilon_i\) — iid 평균 0 의 잡음 곡선

2.2 핵 ψ 의 의미

\(\psi(t, s)\) 는 “회귀자의 시점 \(s\) 의 값이 반응의 시점 \(t\) 의 값에 주는 영향력” 이다.

\(\psi(t_0, s_0) > 0\) → \(X(s_0)\) 의 증가가 \(Y(t_0)\) 를 증가시킴.
\(\psi(t_0, s_0) < 0\) → 반대 방향.
\(\psi(t_0, s_0) \approx 0\) → 영향 없음.

핵 자체는 2D 표면이며, persp 그림이나 히트맵으로 시각화한다.

2.3 직관 1: 무한차원 행렬-벡터 곱

다변량 행렬 곱 \(y_t = \sum_s \psi_{ts} x_s\) 에서 시간 인덱스 \(s, t\) 를 연속화하면 적분이 된다.

이산 (다변량)	연속 (함수)
\(y_t = \sum_{s=1}^p \psi_{ts} x_s\)	\(Y(t) = \int_{\mathcal{S}} \psi(t, s) X(s) \, ds\)
행렬 \(\boldsymbol{\Psi} \in \mathbb{R}^{m \times p}\)	핵 \(\psi(t, s)\)
행 인덱스 \(t\) × 열 인덱스 \(s\)	출력 시간 \(t\) × 입력 시간 \(s\)

핵 \(\psi\) 는 본질적으로 무한차원 행렬 — 입력 곡선을 적분 변환으로 출력 곡선에 매핑한다.

2.4 직관 2: 비유 — 음향 시스템의 임펄스 응답

스피커에 임펄스 신호(특정 순간의 짧은 펄스)를 주면 일정 시간 동안 잔향이 발생한다. 이 응답 패턴이 스피커의 “임펄스 응답 함수” \(h(t, s)\) — 시점 \(s\) 의 입력이 시점 \(t\) 의 출력에 주는 기여.

복잡한 입력 \(X(s)\) 는 무한히 많은 작은 임펄스의 중첩으로 볼 수 있고, 각 임펄스의 응답을 적분하면 출력 \(Y(t) = \int h(t, s) X(s) \, ds\) — 정확히 함수-on-함수 회귀의 형태.

스피커의 정체(임펄스 응답)를 입력-출력 쌍에서 추정하는 것이 “시스템 식별(system identification)” — 함수-on-함수 회귀의 공학적 동의어.

2.5 직관 3: 비유 — 빛의 굴절

프리즘은 백색광(여러 파장의 빛, \(X(s)\))을 굴절각별 강도 분포(\(Y(t)\))로 변환한다. 프리즘 내부의 굴절 특성은 \(\psi(\text{출력각}, \text{파장})\) — 두 변수의 함수.

여러 광원-스펙트럼 쌍을 관측하면 프리즘의 굴절 표면을 역추정할 수 있다. 함수-on-함수 회귀는 이런 “입력-출력 매핑의 표면 추정” 문제이다.

3 일반 모형 (스칼라 + 다중 함수 회귀자)

3.1 모형 식 (식 5.11)

실제 데이터에는 스칼라 회귀자와 함수 회귀자가 함께 등장한다. 함수 회귀자가 두 개인 경우:

\[ Y_i(t) = \alpha(t) + \mathbf{w}_i^T \boldsymbol{\gamma} + \int_{\mathcal{S}} \psi_1(t, s) X_{i1}(s) \, ds + \int_{\mathcal{R}} \psi_2(t, r) X_{i2}(r) \, dr + \varepsilon_i(t). \]

여기서:

\(\mathbf{w}_i = (w_{i1}, \ldots, w_{iq})^T\) — 스칼라 회귀자 벡터
\(\boldsymbol{\gamma} \in \mathbb{R}^q\) — 표준 회귀계수
\(X_{i1}, X_{i2}\) — 함수 회귀자 (서로 다른 도메인 \(\mathcal{S}, \mathcal{R}\) 가능)
\(\psi_1, \psi_2\) — 각 함수 회귀자에 대한 이변량 핵
\(\varepsilon_i\) — iid 가우스 random function

3.2 관측 데이터

각 단위 \(i\) 에 대해 다음을 관측:

\[ Y_i(t_{ij}), \ 1 \leq j \leq J_i; \quad X_{i1}(s_{ik}), \ 1 \leq k \leq K_i; \quad X_{i2}(r_{i\ell}), \ 1 \leq \ell \leq L_i; \quad \mathbf{w}_i. \]

본 절에서는 관측 시점이 밀집(dense) — 각 도메인을 충분히 잘 덮는 — 경우만 다룬다. 희소(sparse) 시점은 Ch.7 의 별도 방법론.

3.3 직관: 왜 다중 회귀자가 자연스러운가

실제 응용에서는 한 가지 정보만으로 반응을 설명하는 경우가 드물다.

응용	함수 반응	함수 회귀자들	스칼라 회귀자
의료 영상	한 뇌 트랙트 신호 곡선	다른 두 트랙트 신호 곡선	환자 나이, 성별, 질병 유무
환경	강수량 곡선	기온·습도 곡선	위도, 고도
분자생물학	RT-PCR 형광 곡선	표적 농도 곡선, 온도 곡선	시약 batch, 효소 활성도

식 (5.11) 의 일반 형태는 이런 다양한 정보를 통합한 모형이며, \(\psi_1, \psi_2, \boldsymbol{\gamma}, \alpha\) 를 동시에 추정해야 한다.

4 식별성 조건: 회귀자 중심화

4.1 문제

핵 \(\psi\) 와 절편 \(\alpha\) 가 데이터로부터 분리되어 식별 되려면 추가 제약이 필요하다. 식별성을 무시하면 다음과 같은 비고유성이 발생.

원 모형 \(Y(t) = \alpha(t) + \int \psi(t,s) X(s) ds\) 의 양변에서 \(\int \psi(t,s) \mu_X(s) ds\) 를 더하고 빼면

\[ Y(t) = \underbrace{[\alpha(t) + \int \psi(t,s) \mu_X(s) ds]}_{\widetilde{\alpha}(t)} + \int \psi(t,s) [X(s) - \mu_X(s)] \, ds. \]

이는 새 절편 \(\widetilde{\alpha}\) 와 같은 핵 \(\psi\) 의 모형. 즉 \((\alpha, \psi)\) 와 \((\widetilde{\alpha}, \psi)\) 가 같은 데이터 분포를 만든다 — 식별 불가능.

4.2 해결: 중심화

회귀자에서 표본 평균을 차감:

\[ X_{i1}^c(s) = X_{i1}(s) - \frac{1}{N} \sum_{i=1}^{N} X_{i1}(s), \quad X_{i2}^c(r) = X_{i2}(r) - \frac{1}{N} \sum_{i=1}^{N} X_{i2}(r). \]

중심화된 회귀자로 모형을 적합하면 절편과 핵이 분리되어 식별 가능. 이후 표기에서는 중심화가 완료된 것으로 가정하고 위첨자 \(c\) 를 생략한다.

4.3 직관: 평균 차감의 역할

다변량 회귀에서 절편을 추정할 때 회귀자의 평균을 차감하는 표준 기법과 같은 동기. 평균 차감 후의 회귀자는 원점을 지나며, 절편은 반응의 평균(계수와 회귀자 평균의 곱이 0이 되어)으로 분리 된다.

함수 회귀에서도 동일한 원리 — 다만 “평균” 이 함수 평균이고 “원점을 지난다” 는 적분 변환이 0 이 된다는 의미.

5 양방향 기저 전개

5.1 절편 함수의 기저 표현

절편 \(\alpha\) 는 일변량 기저로 전개:

\[ \alpha(t) = \sum_{m=1}^{M} \alpha_m B_m(t), \]

여기서 \(B_1, \ldots, B_M\) 은 일변량 기저, \(M\) 은 큰 수.

5.2 회귀 핵의 이변량 기저 표현

\(\psi_1(t, s)\) 는 두 인덱스를 가지므로 이변량 기저 가 필요:

\[ \psi_1(t, s) = \sum_{g=1}^{G} \psi_{1, g} B_{1, g}(t, s), \quad \psi_2(t, r) = \sum_{h=1}^{H} \psi_{2, h} B_{2, h}(t, r). \]

이변량 기저는 보통 두 일변량 기저의 텐서 곱 으로 구성:

\[ B_{1, g}(t, s) = B_g^t(t) \cdot B_g^s(s). \]

5.3 직관: 텐서 곱 기저의 의미

\(M_t\) 개의 \(t\)-기저와 \(M_s\) 개의 \(s\)-기저를 곱하면 \(G = M_t \times M_s\) 개의 이변량 기저가 생성된다. 예를 들어 \(t\) 에 30 개 B-spline, \(s\) 에 30 개 B-spline 이면 이변량 기저 \(30 \times 30 = 900\) 개.

이 큰 기저 수는 함수-on-함수가 함수-on-스칼라보다 본질적으로 어려운 이유 를 설명한다 — 모수의 개수가 차원에 곱셈적으로 폭발한다.

5.4 비유: 격자 직물

\(t\)-방향으로 \(M_t\) 줄, \(s\)-방향으로 \(M_s\) 줄이 짜여 있는 직물을 상상한다. 각 격자 점이 한 이변량 기저 \(B_{g_t}^t(t) \cdot B_{g_s}^s(s)\) 에 대응. 핵 \(\psi_1\) 의 추정은 “직물 위에 낚싯바늘을 박아 표면을 올려 데이터에 맞추는 작업” 에 비유할 수 있다 — \(G\) 개 격자 점의 높이(\(\psi_{1, g}\)) 를 결정한다.

5.5 Riemann sum 으로 적분 이산화

\(\psi_1\) 와 \(X_{i1}\) 의 적분을 격자 위의 합으로 근사:

\[ \int_{\mathcal{S}} \psi_1(t, s) X_{i1}(s) \, ds \approx \sum_{k=1}^{K} (s_k - s_{k-1}) \psi_1(t, s_k) X_{i1}(s_k). \]

기저 전개를 대입:

\[ = \sum_{g=1}^{G} \underbrace{\left[\sum_{k=1}^{K} (s_k - s_{k-1}) B_{1, g}(t, s_k) X_{i1}(s_k)\right]}_{B_{1, g, i}^*(t)} \psi_{1, g}. \]

새 일변량 함수 \(B_{1, g, i}^*(t)\) 는 이변량 기저와 단위 \(i\) 의 회귀자 곡선의 합성 — 한번 계산하면 알려진 함수.

5.6 직관: \(B^*\) 가 등장하는 이유

원 모형은 핵 \(\psi_1\) 과 회귀자 \(X_{i1}\) 의 적분 — 두 객체의 결합. 이를 추정 가능한 형태로 변환하려면 한 객체를 “안다고 가정한 함수” 로 처리해야 한다.

\(B^*_{1, g, i}(t)\) 는 본질적으로 “기저 \(B_{1, g}\) 를 회귀자 \(X_{i1}\) 으로 가중평균한 것” — 미지수가 아닌 데이터 의존적 함수. 이 변환을 통해 추정 문제는 “\(\psi_{1, g}\) 라는 미지수 벡터에 대한 선형 회귀” 로 환원된다.

6 이산 근사 모형

6.1 통합된 모형 식

이산화된 적분과 절편 기저를 모두 대입한 모형:

\[ Y_i(t_j) = \mathbf{w}_i^T \boldsymbol{\gamma} + \sum_{m=1}^{M} B_m(t_j) \alpha_m + \sum_{g=1}^{G} B_{1, g, i}^*(t_j) \psi_{1, g} + \sum_{h=1}^{H} B_{2, h, i}^*(t_j) \psi_{2, h} + \varepsilon_i(t_j). \]

6.2 미지 모수 벡터

추정해야 할 모수를 모두 모으면:

\[ \boldsymbol{\alpha} = (\alpha_1, \ldots, \alpha_M)^T, \quad \boldsymbol{\psi}_1 = (\psi_{1, 1}, \ldots, \psi_{1, G})^T, \quad \boldsymbol{\psi}_2 = (\psi_{2, 1}, \ldots, \psi_{2, H})^T, \]

그리고 스칼라 계수 \(\boldsymbol{\gamma} = (\gamma_1, \ldots, \gamma_q)^T\).

6.3 잔차 제곱합

식 (5.12):

\[ S(\boldsymbol{\gamma}, \boldsymbol{\alpha}, \boldsymbol{\psi}_1, \boldsymbol{\psi}_2) = \sum_{i=1}^{N} \sum_{j=1}^{J} \left( Y_i(t_j) - \mu_i(t_j; \boldsymbol{\gamma}, \boldsymbol{\alpha}, \boldsymbol{\psi}_1, \boldsymbol{\psi}_2) \right)^2, \]

여기서

\[ \mu_i(t_j) = \mathbf{w}_i^T \boldsymbol{\gamma} + \sum_{m=1}^{M} B_m(t_j) \alpha_m + \sum_{g=1}^{G} B_{1, g, i}^*(t_j) \psi_{1, g} + \sum_{h=1}^{H} B_{2, h, i}^*(t_j) \psi_{2, h}. \]

이는 표준 다변량 선형 회귀의 형태 — 미지수 벡터가 \((\boldsymbol{\gamma}, \boldsymbol{\alpha}, \boldsymbol{\psi}_1, \boldsymbol{\psi}_2)\) 의 결합이고 설계 행렬은 \(\mathbf{w}_i, B_m(t_j), B_{1,g,i}^*(t_j), B_{2,h,i}^*(t_j)\) 로 구성. 그러나 단순 LS 는 \(G, H\) 가 매우 크기 때문에 과적합이 발생한다 — 거칠기 벌점이 필요한 이유.

7 거칠기 벌점

7.1 절편 함수의 거칠기

\(\alpha\) 의 2차 도함수 제곱적분:

\[ P_0(\alpha) = \int_{\mathcal{T}} \{\alpha''(t)\}^2 \, dt. \]

기저 전개를 대입:

\[ P_0(\alpha) = \sum_{m, m' = 1}^{M} \alpha_m \alpha_{m'} \int_{\mathcal{T}} B_m''(t) B_{m'}''(t) \, dt = \boldsymbol{\alpha}^T \mathbf{D}_0 \boldsymbol{\alpha}, \]

여기서 \(\mathbf{D}_0\) 는 절편 기저의 거칠기 행렬 — 5.2 의 \(\mathbf{R}\) 과 같은 구조.

7.2 핵 ψ_1 의 거칠기: 라플라시안 벌점

이변량 함수 \(\psi_1(t, s)\) 의 거칠기는 다양한 방식으로 정의 가능. 흔한 선택은 라플라시안(Laplacian):

\[ (L \psi_1)(t, s) = \frac{\partial^2}{\partial t^2} \psi_1(t, s) + \frac{\partial^2}{\partial s^2} \psi_1(t, s). \]

벌점 함수:

\[ P_1(\psi_1) = \int_{\mathcal{T}} \int_{\mathcal{S}} \{(L \psi_1)(t, s)\}^2 \, dt \, ds = \boldsymbol{\psi}_1^T \mathbf{D}_1 \boldsymbol{\psi}_1, \]

여기서

\[ D_1(g, g') = \int_{\mathcal{T}} \int_{\mathcal{S}} (L B_{1, g})(t, s) (L B_{1, g'})(t, s) \, dt \, ds. \]

7.3 직관: 라플라시안이 왜 자연스러운가

라플라시안은 두 방향의 곡률을 같이 측정한다.

\(\psi_{tt}\) 만 벌점 → \(t\)-방향 곡률만 매끄럽게 강제. \(s\)-방향은 자유.
\(\psi_{ss}\) 만 벌점 → 반대.
라플라시안 \(\psi_{tt} + \psi_{ss}\) 벌점 → 두 방향 모두 매끄럽게.

도메인 지식이 두 방향의 매끄러움이 비슷하다고 시사하면 라플라시안이 적절. 한 방향이 본질적으로 거칠 것이 예상되면 다른 형태(혼합 도함수, 비대칭 가중치 등) 를 쓴다.

7.4 비유: 종이 위의 곡면

종이를 양손으로 잡고 곡면을 만들 때, 종이가 두 방향 모두 부드럽게 휘는 것이 라플라시안 벌점의 효과이다. 한 방향만 부드럽게 강제하면 다른 방향에 날카로운 주름이 생길 수 있다.

핵 \(\psi(t, s)\) 도 마찬가지로, 두 시간 축 모두에서 부드러운 표면 이 자연스러우므로 라플라시안 벌점이 표준 선택이다.

7.5 두 함수 회귀자의 결합 벌점

회귀자가 둘이면 각 핵에 별도 벌점을 부과:

\[ P_2(\psi_2) = \int \int \{(L \psi_2)(t, r)\}^2 \, dt \, dr = \boldsymbol{\psi}_2^T \mathbf{D}_2 \boldsymbol{\psi}_2. \]

8 벌점 최소제곱 추정량

8.1 통합 벌점 손실 (식 5.13)

\[ S(\boldsymbol{\gamma}, \boldsymbol{\alpha}, \boldsymbol{\psi}_1, \boldsymbol{\psi}_2) + \lambda_0 \boldsymbol{\alpha}^T \mathbf{D}_0 \boldsymbol{\alpha} + \lambda_1 \boldsymbol{\psi}_1^T \mathbf{D}_1 \boldsymbol{\psi}_1 + \lambda_2 \boldsymbol{\psi}_2^T \mathbf{D}_2 \boldsymbol{\psi}_2. \]

여기서:

\(S\) — 식 (5.12) 의 잔차 제곱합.
\(\lambda_0, \lambda_1, \lambda_2\) — 절편·핵 1·핵 2 의 매끄러움 모수 (각각 독립 선택).
\(\mathbf{D}_0, \mathbf{D}_1, \mathbf{D}_2\) — 비음정 거칠기 행렬.

8.2 직관: 모수마다 별도 정칙화

함수 회귀의 일관된 패턴 — 각 함수 모수에 별도의 매끄러움 모수.

모수	매끄러움 모수
절편 \(\alpha\)	\(\lambda_0\)
첫 핵 \(\psi_1\)	\(\lambda_1\)
둘째 핵 \(\psi_2\)	\(\lambda_2\)
스칼라 계수 \(\gamma\)	(없음 — 표준 ridge 가능하지만 보통 부과 안 함)

벌점 모수의 개수가 \(\geq 3\) 이 되면 그리드 검색이 비현실적이고, REML 같은 자동 방법이 필수.

8.3 추정의 실무적 도구: pffr

식 (5.13) 의 최소화는 refund::pffr 함수로 자동 수행된다. 사용자는 모형의 형태(어떤 변수가 함수 회귀자인지)만 명시하면 되고, 기저 선택·차원·벌점·모수 선택은 모두 기본값으로 처리된다 (필요 시 사용자가 조정 가능).

9 refund 패키지의 두 핵심 함수

9.1 fosr 와 pffr 의 비교

refund 의 두 함수

함수	모형	강점	약점
`fosr`	함수-on-스칼라 전용	OLS·GLS 지원, 순열 검정 (`fosr.perm`), `fd` 객체 직접 입력	함수-on-함수 불가
`pffr`	일반 함수 회귀 (스칼라/함수 혼합)	더 유연·강력, 함수-on-함수 가능, mixed model framework	GLS 자동화 부족

실무에서는 단순 함수-on-스칼라에는 fosr 를, 그 외 모든 경우에는 pffr 를 쓰는 것이 일반적이다. pffr 은 mgcv::gam 위에 구축되어 REML 매끄러움 선택이 자동.

9.2 두 함수의 인터페이스 차이

# fosr 사용 패턴
fit_fosr <- fosr(fdobj = Y_fd,    # 함수 반응을 fd 객체로 입력
                 X = design_mat,    # 설계 행렬
                 con = constraints, # 식별성 제약
                 method = "OLS")    # 또는 "GLS"

# pffr 사용 패턴
fit_pffr <- pffr(Y ~ x1 + ff(X_func1) + ff(X_func2),  # 수식 인터페이스
                 data = df_with_columns)

pffr 의 수식 인터페이스가 직관적 — 표준 R 의 lm/glm 과 비슷하며 함수 회귀자는 ff() (function-on-function), lf() (linear function-on-scalar), af() (additive function) 같은 헬퍼로 표시.

10 함수-on-스칼라 예시: 캐나다 기온 데이터

10.1 데이터 배경

fda 패키지의 CanadianWeather 데이터셋 — 캐나다 35 개 관측소의 일평균 기온/강수량 곡선 (1960~1994 평균).

\(Y_i(t)\) — 관측소 \(i\) 의 \(t\) 일째 기온 (\(t = 1, \ldots, 365\), 단위 \(^\circ\)C).
각 관측소는 4 개 기후권(Arctic, Atlantic, Continental, Pacific) 중 하나에 속함.

10.2 모형

\[ Y_i(t) = \mu(t) + \alpha_{r(i)}(t) + \varepsilon_i(t), \]

여기서 \(\mu(t)\) 는 캐나다 전체 평균 기온 곡선, \(\alpha_r(t)\) 는 기후권 \(r\) 의 평균과의 편차. 식별성을 위해 \(\sum_r \alpha_r(t) = 0\).

이 모형은 5.1 의 함수-on-스칼라 (식 5.1) 의 특수 경우로, \(\beta_1(t) = \mu(t)\), \(\beta_r(t) = \alpha_r(t)\) (\(r = 2, 3, 4, 5\)) 와 5x35 설계 행렬 \(\mathbf{X}\).

10.3 직관: 더미 코딩의 함수 일반화

다변량 ANOVA 의 더미 코딩과 같은 구조 — 각 관측소가 한 기후권에 속하므로 \(x_{i1} = 1\) (절편), \(x_{ir} = \mathbb{1}[\text{지역} = r]\). 단, 회귀계수가 시간에 따라 변하는 함수 — 4 월 vs 10 월에 기후권 차이의 패턴이 다르다.

10.4 OLS vs GLS

이웃한 관측소의 기온 곡선은 공간 의존성 을 가진다 — Vancouver 와 Victoria 의 잔차 곡선은 서로 상관됨.

OLS: 잔차가 비상관임을 가정. 점추정은 비편향이지만 표준오차 과소평가.
GLS: 잔차의 공분산 구조를 추정하여 가중 최소제곱. 의존 구조를 활용하므로 더 효율적.

fosr 는 두 옵션을 지원 — method = "OLS" 또는 "GLS". GLS 의 자세한 메커니즘은 Reiss et al. (2010) 참조.

10.5 코드: fosr 적합 (OLS + GLS 비교)

library(fda)
library(refund)

# 1. 함수 데이터 객체 생성 — 25 개 Fourier 기저로 평활
daybasis25 <- create.fourier.basis(rangeval = c(0, 365), nbasis = 25,
                                   axes = list("axesIntervals"))
Temp.fd <- with(CanadianWeather,
                smooth.basisPar(day.5,
                                dailyAv[, , "Temperature.C"],
                                daybasis25)$fd)

# 2. 설계 행렬 — 절편 + 4 개 기후권 더미
modmat <- cbind(1, model.matrix(~ factor(CanadianWeather$region) - 1))

# 3. 식별성 제약 — 4 개 기후권 효과의 합 = 0
constraints <- matrix(c(0, 1, 1, 1, 1), 1)

# 4. OLS 적합 (그리드 검색으로 lambda 선택)
olsmod <- fosr(fdobj = Temp.fd, X = modmat, con = constraints,
               method = "OLS", lambda = 100 * 10:30)

# 5. GLS 적합 (자동 lambda)
glsmod <- fosr(fdobj = Temp.fd, X = modmat, con = constraints,
               method = "GLS")

# 6. 결과 시각화 — 두 방법의 추정 효과 함수 비교
par(mfrow = c(2, 5), mar = c(5, 2, 4, 1))
plot(olsmod, split = 1, set.mfrow = FALSE,
     titles = c("OLS: Intercept",
                levels(factor(CanadianWeather$region))),
     ylab = "", xlab = "Day")
plot(glsmod, split = 1, set.mfrow = FALSE,
     titles = c("GLS: Intercept",
                levels(factor(CanadianWeather$region))),
     ylab = "", xlab = "Day")

10.6 결과 해석

전형적 출력:

\(\widehat{\mu}(t)\) — 캐나다 전체 평균 기온의 연주기. 1월 최저, 7월 최고 (대략 -10\(^\circ\)C ~ +15\(^\circ\)C).
\(\widehat{\alpha}_{Arctic}(t)\) — 매우 음의 값, 특히 겨울에 더 강하게 음. 북극권의 혹독한 겨울 반영.
\(\widehat{\alpha}_{Atlantic}(t), \widehat{\alpha}_{Continental}(t)\) — 비교적 0 에 가까운 변동.
\(\widehat{\alpha}_{Pacific}(t)\) — 약간 양의 값, 겨울에 더 따뜻함 (해양성 기후 효과).

OLS 와 GLS 의 점추정은 매우 비슷하지만, GLS 가 약간 더 매끄러운 경향이 있다 — 공간 의존성을 모형에 포함하여 잡음 중 일부가 “구조” 임을 인식하기 때문.

10.7 직관: GLS 의 추가 매끄러움

OLS 가 잔차의 시점 간 의존성과 관측소 간 의존성을 모두 무작위 잡음으로 처리하는 반면, GLS 는 관측소 간 공분산을 모형에 흡수한다. 이는 효과적으로 유효 표본 크기를 약간 줄이는 효과 — 의존된 정보는 독립 정보보다 가치가 작으므로, GLS 가 기온 변동의 “신호” 를 약간 더 보수적으로 추정한다.

GLS 가 항상 더 좋다는 보장은 없지만, 공간·시간 의존성이 명확한 경우 표준 도구이다.

11 함수-on-함수 예시: pffrSim 시뮬레이션

11.1 시뮬레이션 설정

refund::pffrSim 함수는 함수-on-함수 모형 (5.10) 을 따르는 데이터를 생성한다.

library(refund)

set.seed(9312)
data_ff <- pffrSim(scenario = c("int", "ff"), n = 100)

scenario = c("int", "ff") 옵션:

"int" — 절편 함수 \(\alpha(t)\) 를 포함.
"ff" — 함수-on-함수 항 \(\int \psi(t, s) X(s) ds\) 를 포함.

11.2 진짜 핵의 형태

pffrSim 의 기본 진짜 핵 (정의):

\[ \psi(s, t) = s \cos(\pi |s - t|) - 0.19. \]

코드 형태:

psi_st <- function(s, t) { s * cos(pi * abs(s - t)) - 0.19 }
s <- seq(0, 1, length = 40)
t <- seq(0, 1, length = 60)
psi_true <- outer(s, t, psi_st)

회귀자 \(X_i\) 는 40 개 등간격 점에서, 반응 \(Y_i\) 는 60 개 등간격 점에서 정의 — 도메인이 다른 일반 형태.

11.3 직관: 진짜 핵의 의미

\(\psi(s, t) = s \cos(\pi |s - t|) - 0.19\) 의 형태:

\(s\) 의 가중치 — \(s\) 가 커질수록 회귀자의 영향력이 커짐.
\(\cos(\pi |s - t|)\) — \(s\) 와 \(t\) 가 가까우면 양, 멀어지면 음으로 진동.
\(-0.19\) — 약간의 음 편향.

이는 “가까운 시점은 강한 양의 영향, 먼 시점은 음의 영향” 패턴 — 시간적 lag 효과의 반전을 모형화한다.

11.4 pffr 적합

# 수식 인터페이스로 함수-on-함수 모형 적합
m_ff <- pffr(Y ~ ff(X1), data = data_ff)

# ff() 헬퍼: X1 을 함수 회귀자로 표시
# 기본값: 텐서 곱 B-spline 기저, REML 로 lambda 자동 선택

11.5 추정 핵의 시각화

# pffr 객체에서 추정 핵의 표면 정보 추출
psi_plot <- plot(m_ff, select = 2, pers = TRUE)[[2]]

# 진짜 핵 시각화
opar <- par(no.readonly = TRUE)
layout(t(1:2))
par(mar = c(0, 1, opar$mar[3], 1))
persp(s, t, psi_true, xlab = "s", ylab = "t",
      main = expression(psi(t, s)),
      phi = 40, theta = 30, ticktype = "detailed",
      zlab = "", border = NA, col = "grey", shade = 0.7,
      zlim = range(psi_true))

# 추정 핵 시각화
persp(psi_plot$x, psi_plot$y,
      matrix(psi_plot$fit, 40, 40),
      xlab = "s", ylab = "t",
      phi = 40, theta = 30, ticktype = "detailed",
      main = expression(hat(psi)(t, s)),
      zlim = range(psi_true), zlab = "", border = NA,
      col = "grey", shade = 0.7)

전형적 결과: 추정 표면 \(\widehat{\psi}(t, s)\) 가 진짜 표면 \(\psi(t, s)\) 의 주요 패턴(양/음 영역, 진동) 을 잘 재현. 약간 더 매끄럽고 극단값이 약화된 형태 — 벌점의 효과.

11.6 적합값과 진짜 예측값의 비교

# 관측·진짜 예측·적합 곡선 비교
par(mar = opar$mar)
layout(t(1:3))

# 1. 관측 곡선
fplot <- function(x, y, col = rgb(0, 0, 0, max(0.1, sqrt(1/nrow(y)))),
                  lwd = 0.5, ...) {
  matplot(x, t(y), type = "l", lty = 1, col = col, lwd = lwd, ...)
}

t_grid <- seq(0, 1, length = 60)

fplot(t_grid, data_ff$Y, xlab = "t", ylab = "",
      main = "Observations", ylim = range(data_ff$Y))
matlines(t_grid, t(data_ff$Y[1:4, ]), col = 1, lty = 1:4, lwd = 2)

# 2. 진짜 예측 곡선 — alpha(t) + integral psi(t,s) X(s) ds
fplot(t_grid, attr(data_ff, "truth")$eta, xlab = "t", ylab = "",
      main = "True predictor", ylim = range(data_ff$Y))
matlines(t_grid, t(attr(data_ff, "truth")$eta[1:4, ]), col = 1, lty = 1:4, lwd = 2)

# 3. 적합 곡선
fplot(t_grid, fitted(m_ff), xlab = "t", ylab = "",
      main = "Estimates", ylim = range(data_ff$Y))
matlines(t_grid, t(fitted(m_ff)[1:4, ]), col = 1, lty = 1:4, lwd = 2)

세 그림의 비교가 모형의 성능을 직관적으로 보여준다.

패널	의미
Observations	잡음 있는 원시 데이터 — 진짜 신호 + \(\varepsilon_i(t)\)
True predictor	잡음 없는 진짜 신호 — 추정해야 할 목표
Estimates	모형 적합 결과

좋은 적합이라면 Estimates 가 True predictor 와 비슷한 형태이며, Observations 의 잡음을 적절히 평활화한 모습.

11.7 직관: 시뮬레이션이 왜 중요한가

실제 데이터에서는 진짜 핵 \(\psi\) 를 모르므로 추정 품질을 직접 평가할 수 없다. 시뮬레이션은 알려진 진리와 추정값을 비교 함으로써 다음을 검증한다.

핵의 패턴(양/음, 진동, 강도) 이 잘 재현되는가.
매끄러움 모수 선택(REML) 이 적절한가.
표본 크기 \(N\) 이 충분한가.

pffrSim 은 간편한 시뮬레이션 도구이며, 사용자는 psi_st 함수를 자신의 도메인에 맞게 수정하여 다양한 시나리오를 검증할 수 있다.

12 핵 추정 방법의 비교

12.1 5.3 의 양방향 기저 + 라플라시안 벌점 vs 5.5 의 FPCA 기반 추정

5.3 은 결정적 기저(B-spline, Fourier 등) 의 텐서 곱을 사용. 5.5 (다음 포스트) 는 데이터 기반 기저(EFPC) 를 사용. 두 접근의 비교:

측면	5.3 양방향 기저 + 벌점	5.5 FPCA 기반
기저	결정적 (B-spline 등)	데이터 기반 (EFPC)
모수 수	\(G = M_t \times M_s\) (큰 수)	\(p \times q\) (작은 수)
정칙화	라플라시안 벌점 (연속)	EFPC 절단 (이산)
계산	mgcv 로 효율적	KL 전개 한 번 + LS
강점	일반적, 자동화 (`pffr`)	데이터 적응적, 빠름
약점	큰 모수, 벌점 선택 민감	EFPC 가 \(Y\) 와 무관할 수 있음

실무에서는 두 방법을 모두 시도하고 결과를 비교하는 것이 바람직하다. pffr 가 기본 도구이지만, 데이터에 강한 주성분 구조가 있으면 5.5 의 FPCA 접근이 더 효율적일 수 있다.

13 두 절의 통합 시각

13.1 한 줄 요약

함수-on-함수 회귀의 핵 \(\psi(t, s)\) 는 두 시간 축 모두에서 무한차원인 모수이며, 양방향 텐서 곱 기저 전개 + 라플라시안 벌점으로 추정한다. refund 패키지의 fosr (함수-on-스칼라 전용, OLS/GLS) 와 pffr (일반 함수 회귀, mixed model framework) 가 이 추정을 자동화하며, 캐나다 기온/pffrSim 데이터로 사용법을 익힐 수 있다.

13.2 5.1~5.2 와의 비교

측면	5.1~5.2 (함수-on-스칼라)	5.3~5.4 (함수-on-함수)
회귀자	스칼라 (\(\mathbf{X} \in \mathbb{R}^{N \times q}\))	함수 (\(X_i: \mathcal{S} \to \mathbb{R}\))
모수	\(q\) 개 효과 함수 \(\beta_k(t)\)	이변량 핵 \(\psi(t, s)\)
자유도 (벌점 후)	\(qK\)	\(G = M_t \times M_s\) — 곱셈적 폭발
벌점	효과별 1차원 거칠기	라플라시안 (양방향 결합)
식별성	자동	회귀자 중심화 필요
R 함수	`fosr` (또는 `pffr` + `lf()`)	`pffr` + `ff()`
적용	시간에 따라 변하는 효과	두 함수 사이의 전체 변환

함수-on-함수는 함수-on-스칼라의 양방향 일반화 — 같은 도구·같은 직관·같은 정칙화 원리를 두 시간 축으로 확장한 것이다.

13.3 Chapter 5 후속 절과의 연결

후속 절	5.3~5.4 의 도구를 어떻게 확장하는가
5.5 FPCA 기반 핵 추정	결정적 기저 대신 \(X, Y\) 의 EFPC 를 기저로 사용 — 핵 분해 식 (5.19)
5.6 효과 없음 검정	추정된 핵 \(\psi\) 의 영가설 \(\psi = 0\) 에 대한 \(\chi^2_{pq}\) 검정
5.7 선형성 진단	EFPC 점수 산점도로 함수 선형 모형의 타당성 점검
Ch.6 함수 GLM	\(\psi\) 추정 + 링크 함수 — 비정규 함수 반응
Ch.8 FAR(1)	시간 시리즈 자기회귀 핵 \(\phi(t, s)\) — 5.3 의 도구 적용

5.3~5.4 의 framework 가 함수-on-함수 회귀의 표준 — 후속 챕터들이 모두 이 toolset 위에서 추가 확장을 한다.

14 관련 주제

선행 지식

후속 주제

관련 개념

텐서 곱 기저 — 이변량 함수 모형화의 표준 도구
라플라시안 연산자 — 다변수 곡률의 표준 측정
Mixed Model 과 REML — pffr 의 추정 framework
GLS (일반화 최소제곱) — 잔차 의존 구조를 가진 회귀
시스템 식별 (system identification) — 함수-on-함수 회귀의 공학 응용