콘크리트 압축강도 회귀분석 — 계수를 뽑는 시간보다 그 계수를 믿어도 되는지 확인하는 데 더 썼다

출발점 — 균열은 배합에서 오는데, 기준은 한 줄이었다

시멘트·물·골재를 섞는 재료 특성에서 비롯되는 균열이 대부분인데, 표준이 말하는 배합 기준은 상한선 한 줄이었다.

주제를 고를 때 근거로 삼은 자료는 한국고속철도건설공단, 「콘크리트 구조물의 균열관리」(2001)였다. 구조물 균열의 대부분이 시멘트·물·골재를 혼합해 쓰는 콘크리트의 재료적 특성에서 비롯된다는 내용이다. 원인이 배합에 있다면, 배합 비율에 관한 기준이 얼마나 정밀한지가 곧 균열 관리의 정밀도가 된다.

그런데 건설공사 표준시방서(국가건설기준센터)에서 확인한 배합 관련 조항은 “콘크리트의 물-결합재비는 원칙적으로 60% 이하이어야 한다”(1.9.3 내구성) 정도였다. 강도에 관해서도 “콘크리트의 강도는 일반적으로 표준양생을 실시한 콘크리트 공시체의 재령 28일일 때 시험값을 기준으로 한다”(1.9.2 강도)는 정의가 있을 뿐, 어떤 재료를 얼마나 넣으면 강도가 얼마나 나오는지에 대한 관계식은 없었다.

발표 도입부에는 당시 기사 세 건을 붙였다. 시멘트를 덜 넣은 레미콘 제조사에 대한 경찰 수사 착수 보도(뉴스1, 2019.05.04), 외벽 낙하를 다룬 「마른 하늘에 ‘외벽 날벼락’… 머리 위가 불안하다」(조선일보), 「붕괴위험 아파트 속, 오늘도 불안한 삶」(2019.05.08)이다. “시멘트 혼합 비율에 대한 정확한 기준이 필요한 시점”이라는 문제 제기를 그 위에 얹었다.

목적을 예측이 아니라 설명으로 고정했다

"어떤 재료가 강도를 얼마나 좌우하는가"가 질문이었으므로, 계수를 읽을 수 있는 모형만 후보가 되었다.

이 프로젝트가 만들려 한 것은 강도를 잘 맞히는 모형이 아니라 배합 기준을 대신할 수 있는 관계식이었다. 결론 슬라이드에도 "'물과 혼합재의 비율을 60%로…'와 같은 기준 대신 이 프로젝트에서 도출해낸 회귀식에 따라 더 명확한 기준을 고려해야 할 것"이라고 적었다. 목적이 그렇다면 계수 하나하나가 곧 답이므로, 예측력이 더 높더라도 계수를 읽을 수 없는 모형은 후보에서 빠진다.

그래서 다중회귀를 택했고, 대신 계수를 믿어도 되는지를 확인하는 쪽에 작업의 무게를 두었다. 발표 목차 6단계 중 회귀분석은 세 번째 항목 하나였고, 나머지 절반(기본가정 검정 · 회귀모형·데이터 교정 · 교차검증)이 전부 "이 계수를 믿어도 되는가"에 배정되어 있다.

데이터 — 재료 7종과 경화 시간, 관측치 1,030개

독립변수 8개는 이미 정해져 있었다. 판단할 것은 그중 무엇을 남길지였다.

독립변수는 배합에 들어가는 재료 7종과 경화 시간을 합해 8개, 종속변수는 압축강도다. 관측치 수는 발표자료에 직접 적혀 있지 않지만, 첫 모형 요약의 잔차 자유도 1,021과 추정 계수 9개(절편 포함)로 역산하면 1,030행이다.

발표자료에 실린 데이터 일부는 다음과 같다. 재료 변수와 시간 변수가 한 모형에 함께 들어가므로, 배합 자체의 효과와 시간에 따른 강도 발현 효과를 계수 수준에서 나눠 읽을 수 있었다.

EDA — 상관행렬이 먼저 알려 준 두 가지

강도와 가장 강하게 붙어 있는 변수는 시멘트였고, 변수들 사이에서 가장 강한 관계는 물과 감수제였다.

변수별 BoxPlot으로 분포와 이상치 후보를 먼저 봤다. 시멘트는 약 100~540 범위에 중앙값 270 부근으로 비교적 고르게 퍼져 있었지만, 광석 찌꺼기는 1사분위가 0에 붙어 있고 중앙값이 22 부근이었다. 석탄재는 중앙값 자체가 0이다. 뒤에서 이상치 절단이 어긋나는 변수가 바로 이 분포들이다.

이어서 여덟 개 독립변수와 압축강도의 상관행렬을 구하고, 상관계수를 격자 위에 수치로 표시해 함께 봤다. 강도와의 상관은 다음과 같았다.

여기서 두 가지가 읽혔다. 첫째, 강도와 가장 강하게 붙어 있는 것은 시멘트(0.50)이고 골재 두 종류는 -0.16, -0.17로 약했다. 뒤에서 조골재·세골재를 제거 후보로 본 판단의 근거가 이미 여기에 있었다. 둘째, 독립변수끼리도 무시할 수 없는 관계가 있었다 — 물과 감수제 -0.66, 물과 세골재 -0.45, 시멘트와 석탄재 -0.40, 석탄재와 감수제 0.38이다.

둘째 항목은 다중공선성 검정을 그냥 넘길 수 없다는 뜻이었다. 특히 시멘트·광석 찌꺼기·석탄재는 서로를 대체하는 결합재이고 물과 감수제는 유동성을 두고 맞바꾸는 관계라, 이들 사이에 공선성이 남아 있으면 “어떤 재료가 얼마나”라는 애초의 질문이 성립하지 않는다.

첫 모형 — 8개 변수를 전부 넣었을 때

결정계수 0.6155. 골재 두 개만 유의수준 0.05를 아슬아슬하게 넘겼다.

8개 변수를 모두 넣어 다중회귀를 적합했다.

잔차 표준오차 10.4(자유도 1,021), 결정계수 0.6155, 조정 결정계수 0.6125, F 204.3(8, 1021), p < 2.2e-16이었다. 절편은 유의하지 않았고, 조골재(0.0544)와 세골재(0.0595)만 유의수준 0.05를 조금씩 넘겼다. 발표자료에는 이 두 변수 옆에 “제거 고려”라고 적어 두었다.

부호는 도메인 상식과 어긋나지 않았다. 시멘트·광석 찌꺼기·석탄재 같은 결합재와 경화 시간은 양(+), 물은 음(-)이다. 물이 많아지면 강도가 떨어진다는 물-결합재비의 방향이 계수로 그대로 나왔다.

변수 선택 — 자동 선택은 아무것도 빼지 않았다

단계적 회귀와 전진선택법은 8개 변수를 그대로 유지했고, 두 변수를 뺀 것은 결국 사람의 판단이었다.

단계적 회귀를 돌렸다. 시작 AIC는 4832.91이었고, 어떤 변수를 빼도 AIC가 올라갔다. 즉 단계적 회귀는 아무 변수도 제거하지 않는 것을 최적으로 판단했다.

전진선택법도 결과가 같았다. 그럼에도 p값을 근거로 조골재와 세골재를 제거해 6변수 모형을 만들었다. 자동 선택 기준(AIC)과 사람의 기준(p값 0.05)이 서로 다른 답을 냈고, 이 프로젝트는 후자를 택한 셈이다. 목적이 설명이라 유의하지 않은 계수를 남겨 두고 해석하고 싶지 않았다는 것이 이유였지만, 자동 선택 결과와 어긋난다는 점은 발표에서 짚지 않았다.

제거의 대가는 거의 없었다. 6변수 모형의 결정계수는 0.614로 전체 변수 모형의 0.6155에서 0.0015밖에 떨어지지 않았고, 조정 결정계수는 0.6125 → 0.6118, 잔차 표준오차는 10.4 → 10.41이었다. 대신 남은 6개 변수의 p값이 전부 0.005 아래로 내려가 계수 해석이 훨씬 단정해졌다. 이후의 가정 검정은 모두 이 6변수 모형을 대상으로 진행했다.

기본가정 검정 ① 등분산성 · 독립성 · 다중공선성

세 가정 중 둘은 통과했고, 등분산성만 완만하게 위배되는 모양이었다.

등분산성은 잔차가 0을 중심으로 하는 수평 밴드 형태로 흩어져 있는지를 본다. 6변수 모형의 잔차-적합값 그림은 대체로 그 형태였지만 적합값 50 부근부터 추세선이 0 아래로 내려갔고, 표준화 잔차의 산포를 보는 Scale-Location 그림의 추세선은 왼쪽에서 오른쪽으로 완만하게 올라갔다. 강도가 높은 구간에서 잔차의 산포가 더 크다는 뜻이다. 같은 판단을 수치로도 확인하려고 비일정분산 검정을 함께 돌렸고, spread-level plot의 추세선도 오른쪽으로 갈수록 올라가 진단 그림과 같은 방향을 가리켰다.

잔차 그림에는 382·384·225번 관측치가 반복해서 이름표를 달고 나타났다. 이 세 점은 뒤의 Q-Q plot에서도 직선을 벗어난 양 끝에 다시 등장한다. 잔차-지렛대 그림에서는 지렛대 값이 최대 0.033 정도에 그쳤고 쿡 거리 경계선이 그림 안에 그려지지 않았다. 지렛대가 큰 단일 관측치 때문에 회귀선이 끌려가는 상황은 아니었다.

오차의 독립성은 더빈-왓슨 검정으로 확인했다. 발표자료의 판단은 “p값이 의미가 없으므로 자기상관은 없다고 할 수 있다”였고, 시차 1은 각 자료를 바로 다음 자료와 비교했다는 뜻이라는 설명을 함께 적었다. 관측치 순서 자체에 시간적 의미가 있는 데이터가 아니므로, 통과가 예상된 항목이었다.

다중공선성은 VIF로 봤고 결과는 “모두 4가 안 넘음 → 다중공선성 없음”이었다. EDA에서 물과 감수제가 -0.66으로 얽혀 있는 것을 보고 가장 걱정했던 항목인데, 8개 변수를 그대로 넣은 선형 모형에서는 문제가 되지 않았다. 상관계수 -0.66은 두 변수 사이의 관계지만 VIF는 한 변수를 나머지 전부로 회귀했을 때의 설명력을 보기 때문에, 쌍 상관이 다소 높아도 VIF가 낮게 나올 수 있다. 결국 네 가정 중 문제가 된 것은 등분산성의 완만한 위배와 — 다음 절에서 볼 — 정규성이었다.

기본가정 검정 ② 정규성 — 여기서 전부 걸렸다

Q-Q plot은 그럴듯했지만, 샤피로-윌크 검정은 8개 변수 전부에서 정규성을 기각했다.

표준화 잔차의 normal Q-Q plot은 대부분 구간에서 45도 직선 위에 놓였다. 다만 양 끝이 직선에서 벗어났고, 여기서도 382·384·225번이 이름표를 달고 나왔다. 그림만 보면 “대체로 정규”라고 넘길 만한 모양이다.

그런데 각 변수에 샤피로-윌크 검정을 돌리자 결과가 달랐다. 8개 변수 전부에서 p값이 유의수준 아래로 나와 정규성 귀무가설이 기각되었다.

가장 심한 것은 굳기 시간(W=0.59071)과 석탄재(W=0.762)였다. 굳기 시간은 3·14·28·56일처럼 몇 개의 값에만 찍히는 이산적인 변수이고(정제 후 남은 데이터에서도 이 네 값만 보인다), 석탄재는 BoxPlot의 중앙값 자체가 0이다 — 절반가량의 배합에 석탄재가 아예 들어가지 않았다는 뜻이다. 분포 모양이 정규와 거리가 먼 변수들이었다.

여기서 남은 선택지는 두 가지였다 — 결과를 그대로 보고하거나, 데이터를 고쳐 다시 검정하거나. 다음 세 절이 후자를 택한 기록이다.

교정 ① 이상치 절단 — 여덟 변수 중 둘만 나아졌다

변수마다 상·하한을 정해 바깥값을 결측으로 바꿨는데, 정규성이 개선된 변수는 둘뿐이었고 조건 두 개는 잘못 적혀 있었다.

첫 번째 교정은 변수별로 상·하한을 정해 그 바깥의 값을 결측으로 바꾸는 방식이었다. BoxPlot에서 읽은 극단치 경계를 참고해 다음과 같이 잘랐다.

절단 규칙을 다시 읽으면서 잘못 적힌 조건 두 개를 확인했다. 첫째, 광석 찌꺼기의 상한 조건이 광석 찌꺼기가 아니라 시멘트 양을 기준으로 적혀 있었다. 광석 찌꺼기가 250을 넘는 배합이 아니라 시멘트가 250을 넘는 배합에서 광석 찌꺼기를 지운 셈이다. 둘째, 조골재의 하한 조건은 데이터에 존재하지 않는 변수를 가리키고 있어 애초에 성립하지 않았다. 둘 다 발표 당시에는 잡지 못했다.

여기에 더해 EDA의 BoxPlot을 같이 놓고 보면 광석 찌꺼기의 하한 조건 자체도 문제가 있다. 광석 찌꺼기는 1사분위가 0에 붙어 있는 변수라, ‘3 미만을 결측으로’는 이상치를 자르는 것이 아니라 하위 사분위 상당 부분을 통째로 지우는 조건이었다.

결과는 발표자료의 표현 그대로 “일부 향상”이었다. 개선된 것은 시멘트(p 2.2e-16 → 6.481e-15)와 광석 찌꺼기(2.2e-16 → 0.005078) 둘뿐이고, 물(1.463e-10 → 2.339e-14), 석탄재(W 0.762 → 0.66835), 세골재(1.842e-10 → 2.745e-12)는 오히려 나빠졌다. 감수제·조골재·굳기 시간의 W값은 절단 전과 동일한 값으로 표에 남았다(조골재는 위의 잘못된 하한 조건을 생각하면 실제로 아무것도 잘리지 않았을 가능성이 있다). 어떤 변수의 꼬리를 자르면 다른 변수의 분포가 함께 잘려 나가기 때문에, 변수별로 따로 자르는 방식은 여기서 한계가 분명했다.

교정 ② 마할라노비스 거리 — 결정계수는 올랐지만 1,030행이 86행이 되었다

다변량 거리로 한 번에 정제하니 설명력은 0.7611까지 올랐다. 남은 표본은 8%였다.

변수별로 자르는 대신 8개 변수를 한꺼번에 보는 거리 기준으로 바꿨다. 공분산행렬을 이용한 마할라노비스 거리를 계산하고, 카이제곱 분포에서 컷오프를 잡아 그보다 먼 관측치를 제외했다.

정제 결과로 남은 데이터를 보면 시멘트부터 압축강도까지 9개 열이 그대로 들어 있다. 즉 거리 계산에 종속변수인 압축강도가 함께 들어갔고, 카이제곱 컷오프의 자유도도 8이 아니라 9였다는 뜻이다. 정규성 p값은 발표자료 표현으로 “많이 향상 됨을 확인”이었고, 정제한 데이터로 다시 적합한 모형의 결정계수는 0.7611(조정 0.7363), 잔차 표준오차는 10.41에서 5.325로 절반 아래가 되었다.

문제는 남은 표본이었다. 이 모형의 잔차 자유도는 77이고 추정 계수는 9개이므로 남은 행은 86개다. 1,030행에서 86행, 8.3%만 남긴 것이다. 컷오프를 카이제곱 분포의 상위 70% 지점에서 잡았다는 것은 거리 기준으로 상위 70%를 전부 이상치로 본다는 뜻이라, 통상적인 이상치 기준(상위 5%나 1%)과 방향이 반대다. 결정계수가 0.6155에서 0.7611로 오른 것은 모형이 좋아져서라기보다 설명하기 쉬운 데이터만 남겨서일 가능성이 크다.

실제로 이 모형에서는 물(p 0.74927)·감수제(0.72250)·조골재(0.15696)가 전부 유의하지 않게 되었다. 표본이 86개로 줄면서 표준오차가 커진 결과다. 물의 계수는 -0.03687로, 6변수 모형의 -0.218088에서 거의 0으로 내려앉았다. 물-결합재비를 다시 쓰겠다는 이 프로젝트의 출발점을 생각하면 물 계수가 사라진 모형은 쓸 수 없다. 이 점을 발표에서는 지적하지 않았다.

교정 ③ 로그변환 — 정규성은 끝내 회복되지 않았다

여덟 중 일곱 변수의 p값이 올랐지만, 여덟 개 모두 여전히 기각이었다.

세 번째 교정은 로그변환이었다. 변환한 데이터로 정규성 검정을 다시 돌린 결과를 앞 두 단계와 나란히 놓으면 다음과 같다.

원본과 비교하면 여덟 변수 중 일곱 개의 p값이 올라갔고, 특히 물은 1.463e-10에서 0.0002237로 여섯 자릿수 이상 개선되었다. 그런데도 여덟 개 전부 여전히 유의수준 아래다. 정규성 가정은 끝내 만족시키지 못했고, 발표자료의 문구도 “정규성 p값이 일부 향상됨이 확인”에서 더 나아가지 않는다.

로그변환한 데이터로 적합한 모형은 결정계수 0.7616(조정 0.7528)이었다. 잔차 자유도 216, 관측치 225개다. 잔차 표준오차 0.2101은 로그 스케일 값이라 다른 모형의 10.4나 5.325와 직접 비교할 수 없다. 여기서도 물 계수(-0.13639, p 0.564908)는 유의하지 않았다.

세 번의 교정을 정리하면 이렇다.

결정계수만 보면 정제·변환 쪽이 낫지만, 둘 다 표본을 크게 줄인 대가로 얻은 값이고 물 계수를 잃었다. 최종 결론에 쓴 것은 6변수 모형 계열이다.

교차검증 — 새 데이터에서는 얼마나 떨어지는가

"이 회귀식이 실제 현실에서 얼마나 잘 들어맞을까"를 10겹으로 확인했다.

마지막 단계는 교차검증이었다. 자료를 10겹으로 나눠, 원본 적합값과 교차검증 적합값 각각에 대해 실측값과의 상관을 제곱해 결정계수를 구하고 그 차이를 보는 방식이었다.

발표자료의 해석은 “이 모형의 회귀식을 새로운 데이터에 적용했을 때 이 모형이 설명할 수 있는 결정계수 값은 0.645 정도일 것”이었고, 결론 슬라이드에도 “결정계수 0.6458487 정도의 설명력을 가진 회귀식”이라고 적었다. 즉 최종적으로 보고한 설명력은 적합 시점의 값이 아니라 교차검증 값이다.

다만 지금 다시 보면 이 절에 맞지 않는 숫자가 하나 있다. 교차검증에 쓴 회귀식은 6변수 모형과 같은 식인데, 6변수 모형의 결정계수는 0.614이지 0.7136이 아니다. 앞선 교정 과정에서 분석에 쓰던 데이터가 이미 바뀐 상태였을 가능성이 있으나, 남아 있는 자료로는 어느 데이터에 대해 돌린 값인지 확정할 수 없다.

결론과 남은 문제

회귀식은 얻었고, 얻는 과정에서 무엇을 못 했는지도 함께 남았다.

결론 슬라이드에 적은 문장은 이것이다 — 시멘트 양, 광석 찌꺼기 양, 석탄재, 물, 고성능 감수제, 굳히는 시간에 따라 콘크리트 강도가 결정계수 0.6458 정도의 설명력을 가진 회귀식으로 설명된다. 그리고 실제 건축 산업에서도 “물과 혼합재의 비율을 60%로” 같은 기준 대신 이런 회귀식에 따른 기준을 고려할 만하다는 제언을 덧붙였다.

지금 다시 읽으면 넘어간 것들이 보인다.

• 정규성은 끝내 만족시키지 못했다. 세 번의 교정 뒤에도 여덟 변수 전부 기각이었는데, 결론에서는 이 사실을 다루지 않았다. 표본이 1,030개면 중심극한정리로 계수 추정의 정규성은 어느 정도 확보되므로 회귀 자체가 무효는 아니지만, 그 논리를 발표에 적어 두지 않았다.
• 검정 대상이 잔차가 아니라 원변수였다. 회귀의 정규성 가정은 잔차에 대한 것인데 정규성 검정은 각 독립변수에 돌렸다. Q-Q plot만 잔차 기준이었다.
• 교정이 표본을 대량으로 버렸다. 마할라노비스 정제 후 남은 것은 86행(원본의 8%)이고, 로그변환 회귀에 쓰인 데이터는 225행(원본의 22%)이었다. 결정계수 상승분의 상당 부분이 여기서 왔을 것이다. 다만 225행이 로그변환 자체 때문에 줄어든 것인지 앞 단계의 결측 치환이 함께 반영된 것인지는 자료로 구분되지 않는다.
• 등분산성 위배를 수치로 확정하지 못했다. 비일정분산 검정을 실행한 기록은 남아 있지만 검정 통계량과 p값이 자료에 없다.

그래도 이 프로젝트에서 만들어진 순서 — 모형 적합 → 네 가지 기본가정 검정 → 교정 → 재검정 → 교차검증 — 는 이후 통계 작업의 기본이 되었다. 특히 “교정했더니 결정계수가 올랐다”를 성과로 쓰기 전에 무엇을 버려서 올랐는지를 먼저 세는 습관은 이 프로젝트를 다시 읽으면서 확실해졌다.

같은 폴더에 함께 남아 있던 참고 자료

2차항을 넣으면 VIF가 폭증하는데, 그것이 진짜 다중공선성인가를 다룬 다른 사람의 발표자료다.

이 프로젝트 폴더에는 「콘크리트 압축강도 회귀분석 PPT.pptx」라는 다른 발표자료가 함께 남아 있다. 표지의 이름은 최현호이고, 슬라이드 안의 실행 결과 캡처에는 2018년 1월 15일이라는 날짜가 찍혀 있으며, R이 아니라 Python으로 작업한 자료다. 데이터도 UCI의 Concrete Slump Test(관측치 103개)로 우리 팀이 쓴 것과 다르다. 우리 팀 발표(2019년 6월)와는 별개의 자료이고, 다중공선성 부분을 참고하려고 보관한 것으로 보인다.

내용 자체는 우리가 VIF를 “모두 4 미만이므로 문제없음”으로 넘긴 지점을 한 단계 더 파고든다. 물 변수에 2차항(Water²)을 추가하자 VIF가 다음처럼 튀었다.

자료의 결론은 “같은 변수로 만든 2차항이므로 다중공선성이 생기는 것이 당연하고, 다른 변수와 생긴 것이 아니므로 이 경우에는 다중공선성으로 보지 않는다”였다. 이어서 중심화(centering)로 VIF를 낮춘 모델 3을 만들고, 모델 2와 모델 3의 R²(0.938)와 MSE가 같다는 점을 근거로 두 모형의 설명력이 동일함을 보였다. VIF 수치가 크게 달라져도 설명력은 그대로라는 것 — VIF를 기준선 하나로만 읽으면 놓치는 부분이다.

목적이 달라지면 모델과 변수도 달라진다

같은 시기에 한 통계 프로젝트 세 건은 모두 표 형태의 데이터였지만, 목적이 달라 고른 모델과 설계한 변수가 전부 달랐다.