Tobit Regression은 Censored Data에 적합한 Regression이다.

Tobit Regression

• Tobit 회귀(Tobit Regression)는 종속 변수가 특정 값에서 절단(Censored)된 상황에서 데이터를 분석하기 위해 사용되는 통계 기법입니다.
• James Tobin이 처음 제안한 이 모델은 경제학과 사회과학 분야에서 자주 사용되며, 일반 선형 회귀로는 설명할 수 없는 상황에서 효과적으로 적용할 수 있습니다.

Tobit Regression 수식

1. 관측된 종속 변수

$$y_i^* = X_i\beta + \varepsilon_i \quad \text{where} \quad \varepsilon_i \sim N(0, \sigma^2)$$

2. 실제 관측된 값

$$y_i =
\begin{cases}
y_i^* & \text{if } y_i^* > c \
c & \text{if } y_i^* \leq c
\end{cases}$$

3. 우도 함수

$$L(\beta, \sigma) = \prod_{y_i > c} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(y_i - X_i\beta)^2}{2\sigma2}\right)
\prod_{y_i \leq c} \Phi\left(\frac{c - X_i\beta}{\sigma}\right)$$

Motivation

• Tobit 회귀는 종속 변수가 특정 범위 내에서만 관찰될 수 있는 상황을 처리하기 위해 개발되었습니다
• . 예를 들어, 설문 조사에서 가구 소득을 조사할 때, 일부 응답자는 최대값(예: 10억 원 이상)으로 응답할 수 있습니다.
  • 이 경우 데이터는 "절단"되었다고 말하며, 일반적인 회귀 분석 기법으로는 이 절단 효과를 제대로 반영하기 어렵습니다.
• 이 기술은 이러한 절단(censoring) 데이터를 다루기 위해 설계되었으며, 절단 데이터의 특성을 모델에 반영하여 더 정확한 결과를 제공합니다.

장단점

장점

1. 절단된 데이터에 최적화: 일반 선형 회귀와 달리, Tobit 회귀는 절단 데이터를 효과적으로 처리할 수 있습니다.
2. 정확한 추정: 절단 데이터의 분포를 고려하기 때문에 종속 변수의 평균이나 분산을 더 잘 추정할 수 있습니다.
3. 다양한 응용 가능: 경제학, 마케팅, 의료 연구 등에서 다양한 상황에 적용할 수 있습니다.

단점

1. 모델 가정 의존성: Tobit 모델은 데이터가 정규분포를 따른다는 가정을 기반으로 하므로, 이 가정이 위반될 경우 결과가 왜곡될 수 있습니다.
2. 복잡한 해석: 모델 출력이 일반 선형 회귀보다 해석하기 어렵습니다.
3. 민감성: 관측치가 적거나 이상치(outliers)가 포함된 경우 민감하게 반응할 수 있습니다.

대안

Tobit 회귀 대신 사용할 수 있는 대안으로는 다음과 같은 방법이 있습니다:

1. Truncated Regression: 절단(censoring) 대신 절단된 데이터만 모델링합니다.
2. Heckman Selection Model: 표본 선택 편향을 교정하는 데 적합합니다.
3. Quantile Regression: 종속 변수의 조건부 분포의 특정 분위수를 모델링합니다.

Python으로 구현한 샘플 코드

[import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel

# 데이터 생성
np.random.seed(0)
n = 1000
x = np.random.normal(size=n)
y = 2.0 + 1.5 * x + np.random.normal(size=n)

# 절단 처리
y_censored = np.clip(y, 0, None)

# 데이터를 데이터프레임으로 정리
data = pd.DataFrame({"x": x, "y": y_censored})

# 사용자 정의 Tobit 모델 (GenericLikelihoodModel 활용)
class TobitModel(GenericLikelihoodModel):
    def loglike(self, params):
        beta = params[:-1]
        sigma = params[-1]
        xb = np.dot(self.exog, beta)
        ll = np.where(self.endog > 0,
                      -0.5 * np.log(2 * np.pi * sigma**2) - 0.5 * ((self.endog - xb) / sigma)**2,
                      np.log(1 - sm.distributions.norm.cdf(xb / sigma)))
        return ll.sum()

# 모형 피팅
exog = sm.add_constant(data['x'])
endog = data['y']
model = TobitModel(endog, exog)
result = model.fit()

print(result.summary())](<import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize
import seaborn as sns

# Tobit Regression 클래스 정의
class TobitRegression:
    def __init__(self, lower_bound=0):
        self.lower_bound = lower_bound
        self.coefficients = None
        self.sigma = None
    
    def _log_likelihood(self, params, X, y):
        beta = params[:-1]
        sigma = params[-1]
        mu = np.dot(X, beta)
        
        censored = (y %3C= self.lower_bound)
        uncensored = ~censored
        
        ll = 0
        if np.any(uncensored):
            ll += np.sum(stats.norm.logpdf(y[uncensored], mu[uncensored], sigma))
        if np.any(censored):
            ll += np.sum(stats.norm.logcdf((self.lower_bound - mu[censored]) / sigma))
        
        return -ll
    
    def fit(self, X, y):
        n_features = X.shape[1]
        initial_params = np.zeros(n_features + 1)
        initial_params[-1] = 1.0
        
        result = minimize(
            self._log_likelihood,
            initial_params,
            args=(X, y),
            method='BFGS'
        )
        
        self.coefficients = result.x[:-1]
        self.sigma = result.x[-1]
        return self
    
    def predict(self, X):
        if self.coefficients is None:
            raise ValueError("Model has not been fitted.")
        
        mu = np.dot(X, self.coefficients)
        lambda_term = stats.norm.pdf((self.lower_bound - mu) / self.sigma) / \
                     stats.norm.cdf((mu - self.lower_bound) / self.sigma)
        return mu + self.sigma * lambda_term

# Generate example data
np.random.seed(42)
n_samples = 1000

# Generate explanatory variables
X1 = np.random.normal(0, 1, n_samples)
X2 = np.random.normal(0, 1, n_samples)
X = np.column_stack([np.ones(n_samples), X1, X2])

# Set true coefficients
true_beta = [2.0, 1.5, -0.5]

# Generate latent variable
noise = np.random.normal(0, 1, n_samples)
y_latent = np.dot(X, true_beta) + noise

# Generate observed dependent variable (censored at 0)
y_observed = np.maximum(y_latent, 0)

# Create DataFrame
df = pd.DataFrame({
    'X1': X1,
    'X2': X2,
    'y_latent': y_latent,
    'y_observed': y_observed
})

# Fit Tobit model
model = TobitRegression(lower_bound=0)
model.fit(X, y_observed)

# Make predictions
y_pred = model.predict(X)

# Print results
print("True coefficients:", true_beta)
print("Estimated coefficients:", model.coefficients)
print("Estimated sigma:", model.sigma)

# Visualization
plt.figure(figsize=(15, 5))

# 1. Latent vs Observed Variables
plt.subplot(131)
plt.scatter(y_latent, y_observed, alpha=0.5)
plt.plot([min(y_latent), max(y_latent)], [0, max(y_latent)], 'r--')
plt.xlabel('Latent Variable')
plt.ylabel('Observed Variable')
plt.title('Latent vs Observed Variables')

# 2. Predicted vs Actual Values
plt.subplot(132)
plt.scatter(y_observed, y_pred, alpha=0.5)
plt.plot([0, max(y_observed)], [0, max(y_observed)], 'r--')
plt.xlabel('Actual Values')
plt.ylabel('Predicted Values')
plt.title('Predicted vs Actual Values')

# 3. Residual Distribution
plt.subplot(133)
residuals = y_observed - y_pred
sns.histplot(residuals[y_observed %3E 0], bins=30)
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.title('Residual Distribution (Uncensored)')

plt.tight_layout()
plt.show()

# Additional analysis: Proportion of censored observations
censored_ratio = (y_observed == 0).mean()
print(f"\nProportion of censored observations: {censored_ratio:.2%}")

# Calculate MSE (for uncensored observations only)
uncensored_mask = y_observed > 0
mse = np.mean((y_observed[uncensored_mask] - y_pred[uncensored_mask]) ** 2)
print(f"MSE (uncensored observations): {mse:.4f}")>)