Tukey's median polish

1. 개요

 Median polish는 John Tukey에 의해 제안된 간단하지만 강건한 (robust) 데이터 분석 기법으로 행과 열로 이루어진 테이블에서 행 영향 (row effect)와 열 영향 (column effect) 그리고 전체적인 중앙값 (overall median, 문헌에 따라서는 mean으로 표기)을 찾는 기법이다.

이것은 행렬에서 각 행과 열의 중앙값을 얻은 후 빼는 방식을 여러 번 반복해서 행 영향과 열 영향을 알아낸다. 중앙값을 이용하기 때문에 outlier에 평균값에 비해 덜 민감하다.

R에서는 medpolish 기능이 기본 제공되는 stat 패키지에 있지만 python에는 찾기가 어렵다. 그래서, stat 패키지의 코드에서 직접 numpy를 이용해 만들었다.

읽으면 좋은 자료 : https://mgimond.github.io/ES218/two_way.html

 

2. 코드

파라미터는 아래와 같다.

• x : 숫자를 포함하는 행렬
• eps : error tolerance로 어느정도의 오차를 허용할 것인지를 의미함.
• maxiter : error tolerance와 상관없이 몇 번 polish과정을 반복할 것인지를 의미한다. 단, error tolerance가 eps 밑으로 떨어지면 maxiter에 도달하지 않더라도 멈춘다.
• trace_iter : median polish과정의 오차값을 추적할 것인지를 의미함.
• na_rm : na값을 제거할지를 의미함. 사실 이것 그냥 아무 의미는 없다.

def medpolish(x, eps=0.01, maxiter=10, trace_iter=True, na_rm=False):
    z = np.array(x)
    nr, nc = z.shape
    t = 0
    r = np.zeros(nr)
    c = np.zeros(nc)
    oldsum = 0
    for iter in range(1, maxiter + 1):
        rdelta = np.median(z, axis=1, where=~np.isnan(z)) if na_rm else np.nanmedian(z, axis=1)
        z = z - np.tile(rdelta, (nc, 1)).T
        r = r + rdelta
        delta = np.median(c, where=~np.isnan(c)) if na_rm else np.nanmedian(c)
        c = c - delta
        t = t + delta
        cdelta = np.median(z, axis=0, where=~np.isnan(z)) if na_rm else np.nanmedian(z, axis=0)
        z = z - np.tile(cdelta, (nr, 1))
        c = c + cdelta
        delta = np.median(r, where=~np.isnan(r)) if na_rm else np.nanmedian(r)
        r = r - delta
        t = t + delta
        newsum = np.nansum(np.abs(z))
        converged = newsum == 0 or abs(newsum - oldsum) < eps * newsum
        if converged:
            break
        oldsum = newsum
        if trace_iter:
            print(f"{iter}: {newsum}")
    if converged:
        if trace_iter:
            print(f"Final: {newsum}")
    else:
        print(f"medpolish() did not converge in {maxiter} iterations")
    rownames = getattr(x, "index", None)
    colnames = getattr(x, "columns", None)
    ans = {
        "overall": t,
        "row": r,
        "col": c,
        "residuals": z,
        "name": repr(x),
    }
    if rownames is not None:
        ans["names_row"] = rownames
    if colnames is not None:
        ans["names_col"] = colnames
    return ans