통계 모형 - lm 과 glm 으로 판매 예측

이 장에서 배우는 것

앞 장에서 지점별 판매 기록을 분석하기 좋은 형태로 정리했다. 이제 정리한 표에서 판매량과 관련된 변수를 골라 관계를 모형으로 표현한다. 평균 판매량을 구하는 것에서 한 걸음 더 나아가, 지점과 행사 여부, 기온이 주어졌을 때 예상 판매량을 계산한다. 재고 부족처럼 두 가지 결과 중 하나가 발생하는 사건은 발생 확률로 예측한다.

통계 모형은 자료의 규칙성을 간결하게 표현하는 도구다. 계산된 계수만 읽어서는 모형이 잘 작동하는지 알기 어렵다. 어떤 자료로 적합했는지, 잔차에 무엇이 남았는지, 새 자료에서 오차가 얼마나 커지는지 함께 살펴야 한다. 이 장에서는 작은 편의점 체인의 가상 자료를 만들고 두 종류의 모형을 적합한 뒤 예측과 진단을 한 파일에서 실행한다.

• 모형 공식(formula)으로 반응변수와 설명변수의 관계를 지정한다.
• lm() 결과에서 계수, 잔차 표준편차, 결정계수를 구분해 읽는다.
• 잔차 그림으로 선형 관계와 오차에 관한 가정을 점검한다.
• glm()으로 재고 부족의 확률을 예측하고 출력 척도를 구분한다.
• predict()를 새 자료에 적용하고 과적합을 피하기 위한 평가 기준을 세운다.

문제 상황

편의점 체인의 발주 담당자는 다음 주 음료 판매량을 예상하려 한다. A지점과 B지점은 유동 인구가 다르고, 할인 행사가 있는 날에는 판매량이 늘어난다. 기온도 음료 수요에 영향을 줄 수 있다. 전체 평균 하나로 발주량을 정하면 지점마다 부족하거나 남는 날이 반복된다.

담당자에게는 서로 다른 두 질문이 있다. 첫째는 “이 조건에서 몇 개가 팔릴 것으로 예상하는가”다. 둘째는 “수요에 비해 재고 여유가 적을 때 재고 부족이 발생할 확률은 얼마인가”다. 첫 질문의 결과는 수량이고, 둘째 질문의 결과는 0과 1 사이의 확률이다. 결과의 형태가 다르므로 같은 모형을 그대로 적용하지 않는다.

실제 기록에서는 행사일이 더운 날에 몰리거나 지점별 기록 기간이 다를 수 있다. 그러면 여러 설명변수의 효과를 구분하기 어려워진다. 이 장의 자료는 관계를 읽는 데 집중하도록 균형 있게 만든 가상 자료다. 판매 자료에는 지점, 행사 여부, 두 기온을 같은 횟수로 넣는다. 재고 자료에는 수요 압력별로 성공과 실패가 모두 나타나게 한다.

예제의 계수와 오차는 자료를 만드는 규칙에서 의도적으로 정한 값이다. 실행 결과는 코드가 설명대로 동작하는지 확인하는 기준이며, 실제 편의점에서 같은 관계가 성립한다는 근거는 아니다. 모형의 계산을 이해한 뒤에는 실제 자료의 수집 방식과 예측 시점을 함께 검토해야 한다.

공식으로 관계를 쓰고 lm 결과를 읽기

왼쪽은 결과, 오른쪽은 설명변수다

sales ~ branch + promo + temp_gap에서 물결표 왼쪽은 예측할 판매량이고 오른쪽은 설명변수다. 오른쪽의 덧셈 기호는 값을 미리 더하라는 뜻이 아니다. 각 변수를 모형의 항으로 넣으라는 뜻이다. 기본적으로 절편도 포함된다. 이 모형은 지점 차이, 행사 차이, 기온 차이를 각각 더해 예상 판매량을 계산한다.

branch는 범주를 나타내는 요인이다. 수준을 c("A", "B")로 지정하고 기본 처리 대비를 사용하면 A가 기준이 된다. branchB 계수는 다른 설명변수가 같은 조건에서 B와 A의 예상 판매량 차이다. promo는 행사가 없으면 0, 있으면 1인 수치 변수다. 이 경우 계수는 0에서 1로 바뀔 때의 예상 판매량 차이를 나타낸다.

temp_gap은 기온에서 20을 뺀 값이다. 이렇게 기준을 옮기면 절편이 A지점, 행사 없음, 기온 20도에서의 예상 판매량이 된다. 기온을 그대로 넣어도 적합값은 같지만 절편은 기온 0도에서의 예상값이 된다. 관측 범위와 가까운 기준을 쓰면 절편을 설명하기 편하다.

공식의 기호는 모형에 들어갈 항을 정한다

표현	의미	편의점 예시
y ~ x + z	절편과 두 변수의 항	행사와 기온의 효과를 더한다
y ~ x * z	두 변수의 항과 상호작용 항	지점마다 행사 효과가 다르다
y ~ x:z	상호작용 항	다른 항은 별도로 지정해야 한다
y ~ x + I(x^2)	원래 값과 실제 제곱값의 항	기온과 판매량의 곡률을 표현한다

공식 안의 ^는 일반 계산식의 제곱과 다르게 해석될 수 있다. 실제 제곱값을 항으로 넣으려면 I()로 감싼다. 또 branch * promo는 지점과 행사 항에 상호작용을 추가한다. 이때 행사 계수는 기준 지점에서의 효과가 되고, 상호작용 계수는 다른 지점에서 더해지는 효과가 된다. 계수의 의미는 공식에 어떤 항이 들어갔는지에 따라 달라진다.

계수와 적합 정도는 다른 질문에 답한다

lm()은 잔차 제곱합을 최소화하는 계수를 구한다. 잔차는 관측 판매량에서 모형의 적합값을 뺀 값이다. 잔차가 양수면 실제 판매량이 적합값보다 많았고, 음수면 적었다. 절편이 있는 최소제곱 적합에서는 잔차 합이 수치 오차 범위에서 0이 된다. 그러나 합이 0이라는 사실만으로 개별 예측이 정확하다고 판단할 수는 없다.

summary()의 계수 표에는 추정값, 표준오차, 검정통계량, 유의확률이 들어 있다. 추정값은 관계의 크기를 나타내고 표준오차는 그 추정의 불확실성을 나타낸다. 유의확률은 해당 계수가 0이라는 가설과 자료의 양립 정도를 정해진 가정 아래에서 계산한 값이다. 효과의 크기나 예측 성능을 대신하는 숫자는 아니다.

잔차 표준편차는 남은 오차의 규모를 판매량과 같은 단위로 나타낸다. 계산할 때 잔차 제곱합을 잔차 자유도로 나눈다. 학습 자료의 평균제곱근오차는 같은 제곱합을 관측 수로 나눈 뒤 제곱근을 취하므로 두 값은 대개 다르다. 예제에서는 계수 네 개를 관측 16개로 추정하므로 잔차 자유도는 12다.

결정계수는 절편이 있는 이 예제에서 전체 판매량 변동 중 모형이 설명한 비율이다. 높은 값은 학습 자료를 잘 설명한다는 뜻이지만 새 자료에서도 잘 예측한다는 보장은 아니다. 설명변수를 추가하면 학습 자료의 결정계수는 낮아지지 않는다. 변수를 늘렸다는 사실만으로 분석이 개선됐다고 말하기 어려운 이유다.

또한 계수는 다른 항을 일정하게 두었을 때의 모형상 차이다. 실제 행사 효과를 인과적으로 설명하려면 행사 선정 과정과 누락 변수 등을 검토해야 한다. 매출이 높은 날에만 행사를 배치한 자료에서는 행사 계수에 원래 수요 차이가 섞일 수 있다. 회귀 계산 자체가 이런 자료 수집 문제를 해소하지는 않는다.

잔차에서 모형이 놓친 구조를 찾기

모형을 적합하면 설명된 부분과 설명되지 않은 부분으로 관측값을 나눌 수 있다. 설명되지 않은 부분인 잔차에 규칙이 남아 있는지 살펴보는 것이 진단의 출발점이다. 한 숫자로 오차를 요약하기 전에 그림을 그리면 곡선, 퍼지는 폭, 특정 관측의 큰 오차를 구분하기 쉽다.

적합값을 가로축, 잔차를 세로축에 두었을 때 잔차가 0 주변에 특별한 구조 없이 분포하는지 확인한다. 곡선 모양이 보이면 관계의 형태를 더 살펴야 한다. 오른쪽으로 갈수록 폭이 넓어지면 오차 분산이 일정하다는 가정이 맞지 않을 수 있다. 지점별로 양수와 음수가 갈리면 지점 효과나 지점별 관계를 충분히 표현하지 못했을 가능성이 있다.

정규 분위수 그림은 잔차 분포의 모양을 정규분포와 비교하는 도구다. 특히 작은 표본에서 계수 검정과 구간 추론을 해석할 때 참고한다. 정규성이 계수를 계산하기 위한 필수 조건인 것은 아니다. 또한 잔차 그림이 보기 좋다는 이유만으로 독립성이나 예측 성능까지 확인했다고 볼 수는 없다.

이 장의 판매 잔차는 설명을 위해 의도적으로 -2와 2만 사용한다. 따라서 잔차 그림에는 두 개의 수평 띠가 나타나고 정규 분위수 그림에서도 정규분포의 모습과 다르다. 이는 프로그램 오류가 아니라 자료 생성 규칙의 결과다. 계수와 오차 계산은 확인할 수 있지만 이 자료로 현실적인 검정 결과를 논의하는 데에는 제한이 있다.

일별 기록이라면 시간 순서로 잔차를 그리는 것도 필요하다. 며칠 연속 같은 방향의 잔차가 나오면 날씨 변화나 요일 구조가 남았을 수 있다. 같은 날 여러 지점의 오차가 함께 커지면 공통 사건의 영향도 생각할 수 있다. 관측끼리 독립인지 여부는 수집 구조와 시간 순서를 함께 봐야 한다.

큰 잔차를 발견했다고 해당 행을 곧바로 삭제하지 않는다. 입력 오류인지, 휴점이나 납품 지연 같은 특별한 상황인지 먼저 확인한다. 실제로 발생할 수 있는 특수 상황이라면 예측 도구가 그 상황을 어떻게 다룰지도 결정해야 한다. 분석하기 편한 행만 남기면 운영 중 마주칠 오차를 과소평가할 수 있다.

glm으로 확률을 예측하고 새 자료에서 평가하기

재고 부족 여부에는 확률 모형을 쓴다

일반화 선형 모형(generalized linear model)은 반응변수의 형태에 맞는 분포와 연결함수를 지정한다. 재고 부족 여부를 0과 1로 기록한 예제에서는 이항 분포와 기본 로짓 연결함수를 사용한다. 이를 로지스틱 회귀(logistic regression)라고 부른다. glm(..., family = binomial())로 적합한다.

모형의 선형 계산 결과를 η, 재고 부족 확률을 p라고 하면 로짓 연결에서는 η가 log(p / (1 - p))다. 확률은 1 / (1 + exp(-η))로 되돌린다. 이 변환 덕분에 설명변수의 선형 계산 결과가 어떤 값이든 예측 확률은 0과 1 사이에 놓인다.

예제의 pressure는 예측 시점에 알고 있는 수요 압력 지표다. 값이 클수록 예상 수요에 비해 재고 여유가 적다고 가정한다. -1, 0, 1에서 각각 12건을 만들고 재고 부족을 3건, 6건, 9건 넣는다. 세 집단의 관측 비율은 0.25, 0.50, 0.75이며 이 관계는 한 개의 기울기를 가진 로지스틱 모형으로 표현할 수 있다.

로지스틱 계수는 확률이 직접 얼마나 증가하는지를 나타내지 않는다. 기울기는 지표가 한 단위 증가할 때 로그 오즈가 얼마나 변하는지 나타낸다. 계수에 exp()를 적용하면 오즈의 배수가 된다. 이 예제에서는 약 3배이며, 같은 한 단위 변화라도 출발 확률에 따라 확률의 변화량은 달라질 수 있다.

predict()에서 type = "response"를 지정하면 확률을 받는다. 기본 출력은 연결함수 척도이므로 로짓 값이다. 확률에 0.5 같은 기준을 적용하면 이진 판단을 만들 수 있지만, 발주에서는 놓친 재고 부족과 과잉 발주의 비용이 다르다. 확률 계산과 운영 기준 결정은 나눠서 생각한다.

평가는 예측 시점의 정보로 한다

판매 모형은 적합에 사용한 자료와 따로 만든 평가 자료에서 오차를 계산한다. 평균제곱근오차(root mean squared error)는 sqrt(mean((actual - predicted)^2))다. 단위가 판매량과 같아서 오차 규모를 읽기 쉽고 큰 오차에 더 많은 무게를 준다. 예제에서는 학습 오차가 2, 평가 오차가 3이 되도록 잔차를 설계한다.

과적합(overfitting)은 학습 자료의 우연한 변동까지 따라가 새 자료에서 성능이 떨어지는 현상이다. 복잡한 항을 추가한 뒤 학습 오차가 줄었다는 사실만으로 선택하지 않는다. 예측에 쓰려는 자료와 비슷한 평가 자료에서 성능이 유지되는지 살핀다. 이 장의 작은 가상 자료는 평가 코드의 동작을 보여 줄 뿐 실제 일반화 성능을 입증하지 않는다.

실제 다음 주 판매량을 예측한다면 과거 기간으로 적합하고 더 나중 기간으로 평가하는 구성이 자연스럽다. 미래 날짜의 기록이 학습 쪽으로 섞이지 않도록 한다. 변수를 고르거나 설정을 바꾸는 데 평가 자료를 반복해서 사용했다면 그 자료도 선택 과정에 영향을 준 것이다. 최종 확인용 기간을 별도로 확보하는 편이 낫다.

설명변수도 예측 시점에 알 수 있어야 한다. 하루가 끝나고 확인한 실제 품절 시간을 오전 발주 예측에 넣을 수는 없다. 기온 역시 운영 시점에는 예보값을 쓰게 된다면 평가도 그 조건에 맞춰야 한다. 편리하게 구한 사후 정보가 예측에 섞이는 일을 정보 누출이라고 한다.

newdata에는 공식에서 사용한 변수 이름과 호환되는 자료형이 필요하다. 새 지점은 기존 모형이 학습한 요인 수준에 없을 수 있고, 새 기온은 학습 범위를 벗어날 수 있다. 예제의 평가 기온은 22도로 학습한 20도와 24도 사이에 둔다. 관측 범위 밖에서 관계를 연장하는 외삽은 별도의 근거가 필요한 예측이다.

완성 코드

아래 내용을 main.R로 저장한다. 외부 파일이나 패키지는 필요하지 않다. 난수를 사용하지 않으므로 실행할 때마다 같은 자료와 수치 출력이 나온다. 대비 방식과 출력 언어를 지정해 지점 계수의 이름과 그림의 글꼴 환경에 따른 차이를 줄인다. 두 그림은 현재 작업 디렉터리에 저장된다.

# main.R
# 1. 실행 설정과 검사 함수
options(contrasts = c("contr.treatment", "contr.poly"))

check_close <- function(actual, expected, tolerance = 1e-7) {
  stopifnot(
    length(actual) == length(expected),
    all(is.finite(actual)),
    all(abs(actual - expected) < tolerance)
  )
}

rmse <- function(actual, predicted) {
  stopifnot(
    length(actual) > 0L,
    length(actual) == length(predicted),
    all(is.finite(actual)),
    all(is.finite(predicted))
  )
  sqrt(mean((actual - predicted)^2))
}

# 2. 판매 학습 자료
train <- expand.grid(
  branch = c("A", "B"),
  promo = c(0, 1),
  temp = c(20, 24),
  replicate = c(1L, 2L),
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)
train$branch <- factor(train$branch, levels = c("A", "B"))
train$temp_gap <- train$temp - 20
train$error <- ifelse(train$replicate == 1L, -2, 2)
train$sales <- 100 +
  20 * as.integer(train$branch == "B") +
  30 * train$promo +
  5 * train$temp_gap +
  train$error

stopifnot(nrow(train) == 16L, !anyNA(train))

# 3. 販売モ형の適合
sales_fit <- lm(
  sales ~ branch + promo + temp_gap,
  data = train,
  na.action = na.fail
)
sales_summary <- summary(sales_fit)

check_close(
  unname(coef(sales_fit)),
  c(100, 20, 30, 5)
)
stopifnot(df.residual(sales_fit) == 12L)
check_close(sum(residuals(sales_fit)), 0)

# 4. 적합에 사용하지 않은 판매 평가 자료
test <- expand.grid(
  branch = c("A", "B"),
  promo = c(0, 1),
  replicate = c(1L, 2L),
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)
test$branch <- factor(test$branch, levels = levels(train$branch))
test$temp <- 22
test$temp_gap <- test$temp - 20
test$error <- ifelse(test$replicate == 1L, -3, 3)
test$sales <- 100 +
  20 * as.integer(test$branch == "B") +
  30 * test$promo +
  5 * test$temp_gap +
  test$error

stopifnot(nrow(test) == 8L, !anyNA(test))
train_pred <- predict(sales_fit, newdata = train)
test_pred <- predict(sales_fit, newdata = test)
train_rmse <- rmse(train$sales, train_pred)
test_rmse <- rmse(test$sales, test_pred)

check_close(train_rmse, 2)
check_close(test_rmse, 3)

# 5. 다음 발주 조건의 예상 판매량
next_day <- data.frame(
  branch = factor("B", levels = levels(train$branch)),
  promo = 1,
  temp_gap = 2
)
next_sales <- predict(sales_fit, newdata = next_day)
check_close(as.numeric(next_sales), 160)

# 6. 재고 부족 자료와 로지스틱 회귀
inventory <- data.frame(
  pressure = rep(c(-1, 0, 1), each = 12L),
  trial = rep(seq_len(12L), times = 3L)
)
inventory$stockout <- as.integer(
  inventory$trial <= rep(c(3L, 6L, 9L), each = 12L)
)
stopifnot(
  !anyNA(inventory),
  all(inventory$stockout %in% c(0L, 1L))
)

stock_fit <- glm(
  stockout ~ pressure,
  data = inventory,
  family = binomial(),
  na.action = na.fail
)
stopifnot(stock_fit$converged)
check_close(
  unname(coef(stock_fit)),
  c(0, log(3)),
  tolerance = 1e-6
)

new_inventory <- data.frame(pressure = c(-1, 0, 1))
stock_prob <- predict(
  stock_fit,
  newdata = new_inventory,
  type = "response"
)
check_close(as.numeric(stock_prob), c(0.25, 0.50, 0.75))
stopifnot(all(stock_prob > 0 & stock_prob < 1))

# 7. 잔차 진단 그림
png("sales-diagnostics.png", width = 1000, height = 480)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(
  fitted(sales_fit), residuals(sales_fit),
  xlab = "Fitted sales", ylab = "Residual",
  main = "Residuals vs fitted",
  pch = 19, col = "#315675", ylim = c(-3, 3)
)
abline(h = 0, col = "#c0392b", lty = 2)
qqnorm(
  residuals(sales_fit),
  main = "Normal Q-Q",
  pch = 19, col = "#315675"
)
qqline(residuals(sales_fit), col = "#c0392b")
invisible(dev.off())

# 8. 수요 압력과 재고 부족 확률 그림
pressure_grid <- data.frame(
  pressure = seq(-1, 1, length.out = 101L)
)
prob_grid <- predict(
  stock_fit,
  newdata = pressure_grid,
  type = "response"
)
png("stockout-probability.png", width = 720, height = 480)
plot(
  pressure_grid$pressure, prob_grid,
  type = "l", lwd = 2, col = "#315675",
  xlab = "Demand pressure", ylab = "Stockout probability",
  main = "Predicted stockout probability", ylim = c(0, 1)
)
points(new_inventory$pressure, stock_prob, pch = 19)
invisible(dev.off())

# 9. 결정적인 수치 출력
cat(sprintf(
  "판매 자료: 학습 %d행, 평가 %d행\n",
  nrow(train), nrow(test)
))
cat("판매 모형 계수:\n")
for (term in names(coef(sales_fit))) {
  cat(sprintf("  %s = %.3f\n", term, coef(sales_fit)[[term]]))
}
cat(sprintf("잔차 자유도: %d\n", df.residual(sales_fit)))
cat(sprintf("잔차 표준편차: %.3f\n", sales_summary$sigma))
cat(sprintf("결정계수: %.3f\n", sales_summary$r.squared))
cat(sprintf("학습 RMSE: %.3f\n", train_rmse))
cat(sprintf("평가 RMSE: %.3f\n", test_rmse))
cat(sprintf(
  "B지점, 행사 있음, 22도 예상 판매량: %.3f\n",
  as.numeric(next_sales)
))
cat(sprintf(
  "수요 압력 계수: %.3f\n",
  coef(stock_fit)[["pressure"]]
))
cat(sprintf(
  "수요 압력 1 증가의 오즈 배수: %.3f\n",
  exp(coef(stock_fit)[["pressure"]])
))
cat("재고 부족 확률:\n")
for (i in seq_len(nrow(new_inventory))) {
  cat(sprintf(
    "  pressure=%d: %.3f\n",
    new_inventory$pressure[i], stock_prob[i]
  ))
}
cat("그림 파일: sales-diagnostics.png\n")
cat("그림 파일: stockout-probability.png\n")
cat("검사: 모두 통과\n")

줄별 해설

실행 설정과 검사 함수. options()는 순서 없는 요인에 처리 대비를 적용한다. 따라서 A를 기준으로 B의 차이를 추정한다. check_close()는 소수 계산을 허용 오차 안에서 검사한다. 회귀 계산에는 부동소수점 오차가 있으므로 계수를 ==로 비교하지 않는다. rmse()는 길이가 같고 비어 있지 않으며 유한한 수치인지 확인한 뒤 오차를 계산한다.

판매 학습 자료. expand.grid()는 지정한 값의 모든 조합을 만든다. 지점 두 개, 행사 상태 두 개, 기온 두 개, 반복 두 개를 곱하면 16행이다. 반복 번호는 자료 생성용 열이며 모형에는 넣지 않는다. 같은 조건에서 오차 -2와 2를 한 번씩 부여하므로 조건별 평균에는 오차가 남지 않는다.

판매량 생성. as.integer(train$branch == "B")는 B에서 1, A에서 0이다. 기본 100개에 B지점이면 20개, 행사면 30개, 20도보다 한 도 높을 때마다 5개를 더한다. 이는 예제 자료의 생성 규칙이다. 현실 자료에서는 이 계수를 미리 알 수 없고 관측 자료에서 추정해야 한다.

모형 적합. 공식은 생성 규칙에 사용한 설명변수만 포함한다. na.action = na.fail은 결측 행을 조용히 제외하는 대신 적합을 중단하게 한다. 자료 크기가 작을 때 몇 행이 빠졌는지 모르고 모형을 비교하는 일을 줄일 수 있다. summary() 결과는 객체에 보관하고 필요한 수치만 꺼내 출력한다.

계수 검사. coef()의 순서는 절편, B지점, 행사, 기온 차이다. 이름을 제외하고 예상 계수와 비교한다. 잔차 자유도가 12인지 검사하면 자료 수와 추정한 계수 수의 관계도 확인할 수 있다. 이 검사는 자료 생성 규칙이 알려진 예제에서만 가능한 확인이며, 실제 모형의 계수 정답을 검사한다는 뜻은 아니다.

판매 평가 자료. 지점과 행사 조합은 유지하고 기온을 22도로 둔다. 같은 조건에 -3과 3의 오차를 넣어 학습 자료보다 오차 규모를 키운다. 요인 수준은 학습 자료에서 가져온다. 이 단계에서 계산한 실제 판매량은 평가를 위한 정답이며 predict()가 예측을 계산하는 데 사용하지 않는다.

판매 예측. predict(sales_fit, newdata = test)는 평가 자료의 설명변수로 예상값을 계산한다. newdata를 생략하면 학습 자료의 적합값을 돌려준다. 두 결과는 목적이 다르므로 평가 코드에서는 새 자료를 명시한다. next_day에는 반응변수인 판매량을 넣지 않아도 된다.

재고 자료. trial은 각 압력 집단 안에서 1부터 12까지의 번호다. 부족 건수를 정하는 데만 쓰며 설명변수가 아니다. 이 번호를 모형에 넣으면 자료를 만든 인위적인 순서까지 학습하게 된다. 예측 시점에 의미가 있는 변수와 예제 생성용 변수를 구분하는 사례다.

로지스틱 적합. 반응변수는 재고 부족이면 1, 아니면 0이다. 각 집단에 두 결과가 모두 있어서 설명변수만으로 모든 결과가 분리되지는 않는다. 적합 후에는 수렴 여부와 계수를 확인한다. 수렴은 반복 계산이 정지 기준을 만족했다는 뜻이며 모형의 유용성을 보장하는 판정은 아니다.

확률과 그림. 세 조건의 확률을 확인한 뒤 -1부터 1까지 촘촘한 격자에서 곡선을 계산한다. 곡선 사이의 모든 위치를 실제로 관측했다는 뜻은 아니다. 자료를 바탕으로 적합한 함수의 형태를 보여 준 것이다. 각 png() 다음의 dev.off()는 파일 기록을 마치며 invisible()은 장치 번호 출력을 숨긴다.

출력 형식. sprintf()로 소수점 자릿수를 고정한다. 기본 모형 출력은 표시 설정에 따라 달라질 수 있으므로 이 프로그램은 필요한 값만 직접 출력한다. 맨 마지막 통과 문장은 모든 stopifnot() 검사를 지나고 그림 파일 기록까지 끝났을 때 나타난다.

실행 결과

터미널에서 파일을 저장한 디렉터리로 이동한 뒤 다음 명령을 실행한다. R은 별도 컴파일 단계 없이 스크립트를 읽어 실행한다.

Rscript main.R

예상 표준 출력은 다음과 같다.

판매 자료: 학습 16행, 평가 8행
판매 모형 계수:
  (Intercept) = 100.000
  branchB = 20.000
  promo = 30.000
  temp_gap = 5.000
잔차 자유도: 12
잔차 표준편차: 2.309
결정계수: 0.991
학습 RMSE: 2.000
평가 RMSE: 3.000
B지점, 행사 있음, 22도 예상 판매량: 160.000
수요 압력 계수: 1.099
수요 압력 1 증가의 오즈 배수: 3.000
재고 부족 확률:
  pressure=-1: 0.250
  pressure=0: 0.500
  pressure=1: 0.750
그림 파일: sales-diagnostics.png
그림 파일: stockout-probability.png
검사: 모두 통과

B지점의 행사일, 기온 22도에서는 100 + 20 + 30 + 5 * 2로 160개를 예상한다. 이는 주어진 조건의 평균 판매량에 대한 예측이다. 특정 하루의 판매량이 반드시 160개라는 뜻은 아니다. 발주량을 정할 때에는 남은 재고, 납품 주기, 부족 비용 등도 고려해야 한다.

판매 학습 자료의 잔차 제곱합은 64다. 학습 오차는 sqrt(64 / 16)으로 2이고 잔차 표준편차는 sqrt(64 / 12)로 약 2.309다. 두 값의 차이는 오류가 아니라 분모의 차이다. 평가 자료에서는 각 잔차의 크기가 3이므로 평가 오차도 3이다.

sales-diagnostics.png에는 잔차 대 적합값 그림과 정규 분위수 그림이 저장된다. 같은 조건의 두 관측이 -2와 2로 갈리는 모습을 확인할 수 있다. stockout-probability.png에는 압력이 증가할수록 부족 확률이 증가하는 곡선과 세 조건의 예측 확률이 표시된다. PNG의 글꼴이나 렌더링은 환경에 따라 조금 달라질 수 있지만 계산한 수치와 파일 이름은 같다.

실무에서 자주 틀리는 것

학습 적합값을 새 자료 평가로 착각한다

평가하려는 자료를 지정하지 않으면 학습 자료의 적합값을 받는다. 학습 오차는 모형이 이미 본 자료에서 계산한 값이다. 이를 다음 기간의 예측 성능으로 보고하면 오차를 작게 보일 수 있다.

틀린 코드는 다음과 같다.

# 새 자료 성능이라고 보고하지만 학습 자료를 평가한다.
pred <- predict(sales_fit)
reported_error <- rmse(train$sales, pred)

평가 자료를 명시하고 실제값도 같은 자료에서 가져온다.

pred <- predict(sales_fit, newdata = test)
reported_error <- rmse(test$sales, pred)

행 수가 우연히 같으면 길이 검사만으로 자료의 잘못된 짝을 찾을 수 없다. 실제 업무에서는 날짜와 지점 식별자로 예측과 정답이 같은 관측에 대응하는지도 확인한다.

로지스틱 기본 출력을 확률로 읽는다

연결함수 척도의 값은 음수이거나 1보다 클 수 있다. 이 값을 확률로 출력하거나 확률 기준과 비교하면 판단이 달라진다.

틀린 코드는 다음과 같다.

prob <- predict(stock_fit, newdata = new_inventory)
needs_review <- prob >= 0.5

확률로 변환한 뒤 운영 기준을 적용한다.

prob <- predict(
  stock_fit,
  newdata = new_inventory,
  type = "response"
)
needs_review <- prob >= 0.5

연결함수 척도를 사용할 목적이라면 기본 출력도 유효하다. 문제는 계산값의 척도와 해석을 맞추지 않은 데 있다. 또한 0.5는 예시 기준이며 사업상의 비용을 반영해 검토해야 한다.

새 지점을 기존 기준 지점처럼 처리한다

학습한 적 없는 C지점의 차이를 기존 모형에서 알아낼 수는 없다. 알려진 수준으로 요인을 강제로 만들면 C가 결측값이 될 수 있다. 이를 확인하지 않으면 예측에 결측값이 섞인다.

틀린 코드는 다음과 같다.

incoming <- data.frame(branch = "C", promo = 0, temp_gap = 2)
incoming$branch <- factor(
  incoming$branch,
  levels = levels(train$branch)
)
pred <- predict(sales_fit, newdata = incoming)

수준을 변환하기 전에 알려진 지점인지 검사한다.

incoming <- data.frame(branch = "B", promo = 0, temp_gap = 2)
known_branches <- levels(train$branch)
stopifnot(all(incoming$branch %in% known_branches))
incoming$branch <- factor(
  incoming$branch,
  levels = known_branches
)
pred <- predict(sales_fit, newdata = incoming)

C지점을 지원하려면 자료를 확보해 모형을 다시 적합하거나 새 지점에 적용할 별도의 규칙을 설계해야 한다. 이름만 기존 지점으로 바꾸면 지점 차이를 잘못 전달하게 된다.

학습 결정계수가 높은 모형만 고른다

항을 늘리면 학습 자료를 더 세밀하게 따라갈 수 있다. 하지만 항이 많아질수록 자료가 뒷받침하지 못하는 관계까지 추정할 수 있다. 다음 코드는 학습 자료의 설명력만으로 모형을 고른다.

simple <- sales_fit
richer <- lm(
  sales ~ branch * promo + temp_gap,
  data = train,
  na.action = na.fail
)
fits <- list(simple = simple, richer = richer)
scores <- vapply(
  fits,
  function(fit) summary(fit)$r.squared,
  numeric(1)
)
chosen <- names(which.max(scores))

후보를 비교하는 평가 자료가 있다면 예측 오차를 계산한다. 이름을 바꿔 선택 기준을 분명하게 드러낸다.

validation_rmse <- vapply(
  fits,
  function(fit) {
    pred <- predict(fit, newdata = test)
    rmse(test$sales, pred)
  },
  numeric(1)
)
chosen <- names(which.min(validation_rmse))

이 가상 자료에서는 추가 상호작용이 필요한 구조가 없으므로 두 후보의 예측은 수치 오차 범위에서 같다. 비슷한 성능이라면 단순한 관계를 유지하는 판단도 가능하다. 실제 분석에서 이 자료로 후보를 골랐다면 최종 성능 확인에는 다시 분리한 자료가 필요하다.

한눈에 보기

판매량과 재고 부족 확률은 적합과 해석의 기준이 다르다

확인 항목	판매량 모형	재고 부족 모형
적합 함수	lm()	glm(family = binomial())
반응변수	판매 수량	부족이면 1, 아니면 0
계수 해석	조건별 예상 수량의 차이	조건별 로그 오즈의 차이
새 자료 예측	predict(fit, newdata)	predict(fit, newdata, type = "response")
핵심 점검	잔차 구조와 새 자료의 오차	수렴, 확률, 실제 발생률과의 관계

공식은 관계를 선언하고, 모형 객체는 추정한 계수와 적합에 관한 정보를 보관한다. summary()는 적합 결과를 요약하고 predict()는 조건을 예상값으로 바꾼다. 같은 이름의 함수가 모형 종류에 맞는 동작을 제공하므로 객체 종류와 출력 척도를 함께 확인한다.

이 프로그램은 base R만 사용한다. 대응 도구로는 tidyverse의 tidyr::expand_grid()가 조합 자료를 만들 때 쓰이고, dplyr::mutate()가 계산 열을 추가할 때 쓰인다. 모형 적합과 예측은 그 도구를 사용하더라도 lm(), glm(), predict()로 수행할 수 있다. 이 장의 실행에는 추가 설치가 필요하지 않다.

검정과 예측 구간은 모형의 가정에 의존한다. 판매 예측의 불확실성을 계산할 때에는 조건부 평균의 불확실성과 실제 하루 판매량의 변동을 구분해야 한다. 다음 장에서는 자료를 다시 뽑는 계산을 통해 추정 결과의 흔들림을 살펴본다.

연습 문제

1. 판매 모형으로 A지점, 행사 없음, 기온 24도인 날의 예상 판매량을 계산하라. newdata를 만들어 예측하고 check_close()로 검사하라.
2. 같은 판매 모형을 사용하되 평가 자료의 오차를 -6과 6으로 바꿔 판매량을 다시 만들라. 평가 오차를 구하고 계수를 다시 적합할 필요가 없는 이유를 설명하라.
3. 수요 압력이 0.5일 때 재고 부족 확률을 계산하라. predict()의 결과를 plogis()로 직접 계산한 값과 비교하라.
4. temp_gap 대신 원래 기온 temp를 넣어 판매 모형을 적합하라. 절편과 기온 계수는 어떻게 바뀌는지 확인하고 학습 자료의 예측값이 같은지 검사하라.

정답과 해설

1. 조건을 공식의 변수로 전달한다

case_a <- data.frame(
  branch = factor("A", levels = levels(train$branch)),
  promo = 0,
  temp_gap = 4
)
answer_a <- predict(sales_fit, newdata = case_a)
check_close(as.numeric(answer_a), 120)

예상 판매량은 120개다. A는 기준 지점이므로 지점 차이를 더하지 않고 행사 효과도 0이다. 기온 차이 4에 기온 계수 5를 곱해 20개를 더한다. temp만 넣고 temp_gap을 생략하면 기존 공식에 필요한 변수가 없으므로 예측할 수 없다.

2. 평가 자료의 오차 규모만 바뀐다

test_wider <- test
test_wider$error <- ifelse(
  test_wider$replicate == 1L, -6, 6
)
test_wider$sales <- 100 +
  20 * as.integer(test_wider$branch == "B") +
  30 * test_wider$promo +
  5 * test_wider$temp_gap +
  test_wider$error

wider_pred <- predict(sales_fit, newdata = test_wider)
wider_rmse <- rmse(test_wider$sales, wider_pred)
check_close(wider_rmse, 6)

평가 오차는 6이다. 설명변수가 같으므로 예측값은 이전과 같고 실제값의 변동만 커졌다. 학습 자료와 적합한 모형은 바꾸지 않았으므로 계수를 다시 추정하지 않는다. 평가 결과를 본 뒤 평가 자료까지 포함해 적합하면 원래의 분리 평가가 유지되지 않는다.

3. 선형 계산값을 확률로 되돌린다

case_pressure <- data.frame(pressure = 0.5)
answer_prob <- predict(
  stock_fit,
  newdata = case_pressure,
  type = "response"
)
manual_prob <- plogis(
  coef(stock_fit)[["(Intercept)"]] +
  coef(stock_fit)[["pressure"]] * 0.5
)
check_close(as.numeric(answer_prob), manual_prob)
check_close(as.numeric(answer_prob), plogis(log(3) * 0.5))

확률은 약 0.634다. plogis()는 로짓 값을 확률로 바꾸는 계산을 제공한다. 압력 0에서 0.5로 이동했다고 확률이 일정한 양만큼 증가하는 것은 아니다. 선형으로 변하는 것은 로그 오즈이며 확률 변환은 곡선이다.

4. 기준을 바꾸면 절편이 달라진다

raw_temp_fit <- lm(
  sales ~ branch + promo + temp,
  data = train,
  na.action = na.fail
)
check_close(
  unname(coef(raw_temp_fit)),
  c(0, 20, 30, 5)
)
check_close(
  as.numeric(predict(raw_temp_fit, newdata = train)),
  as.numeric(predict(sales_fit, newdata = train))
)

기온 계수는 5로 같고 절편은 100에서 0으로 바뀐다. 원래 모형의 100 + 5 * (temp - 20)을 풀면 5 * temp가 되기 때문이다. 두 모형은 같은 선을 다른 기준으로 표현한다. 기온 0도는 학습 범위 밖이므로 새 절편을 실제 관측 조건의 예상 판매량처럼 해석하는 데에는 주의가 필요하다.