JuliaCon 2021: DataFrames.jl 1.2 tutorial¶
- Julia 1.6.1 기반
이 튜토리얼은 다음과 같은 통계개념을 사용하는데 미리 알고 가면 좋습니다.
- confidence intervals - 신뢰구간
- density estimators - 밀도추정
- probit model - 프로빗모형
- 선택확률 $\pi(x)$가 [0,1]의 구간에 놓이도록 하기 위해서는 설명변수와 선택확률 사이에 선형이 아닌 비선형의 관계가 요구 됩니다. 이러한 비선형 관계를 표준정규분포의 누적확률 분포함수를 이용하여 표현하는 경우 이를 프로빗 모형이라고 하며, 이 때 사용되는 링크함수를 프로빗 함수하고 합니다. 예를 들어 구매금액에 따라 재구매 고객의 여부를 분석할 때, 구매금액이 많으면 재구매를 한다는 가정을 한다면, 여기서 단순한 구매금액의 분포는 정규분포의 형태를 띠게 될 것입니다. 하지만 가정에 따르면 어느 정도 구매금액이면 재구매를 할 것이므로 재구매의 분포는 누적정규분포의 형태일 것입니다.
- bootstrapping)
- 여기서 사용하는 데이터는 경제활동참가(Labor force participation) 데이터 입니다.
- 데이제공정보
- 국가 : 스위스
경제활동참여 (Labor force participation)
- 데이터 설명
| 컬럼명 | 영문 | 한글 |
|---|---|---|
| lfp | labour force participation ? | 경제활동참여 여부 |
| lnnlinc | the log of nonlabour incom | 비 노동소득의 로그값 |
| age | age in years devided by 10 | 나이를 10으로 나눈값 |
| educ | years of formal education | 정규교육연수 |
| nyc | the number of young children (younger than 7) | 7세 미만 자녀수 |
| noc | number of older children | 7세 이상 자녀수 |
| foreign | foreigner ? | 외국인 여부 |
using Bootstrap
using CSV
using CategoricalArrays
using Chain
using DataFrames
import Downloads
using GLM
using Plots
using Random
using StatsPlots
using Statistics
DataFrame 출력시 열(column)은 1000개, 행(row)은 20개 까지 표시 되도록 환경변수를 설정 한다.
ENV["LINES"] = 20
ENV["COLUMNS"] = 1000
Data preprocessing¶
- 데이터 다운로드
data_file = "data/participation.csv"
Downloads.download("https://vincentarelbundock.github.io/Rdatasets/csv/Ecdat/Participation.csv", data_file)
readlines(data_file)
데이터 분석의 목적은 경제활동참여여부(lfp) 예측모델을 만드는 것입니다
df_raw = CSV.read(data_file, DataFrame)
describe(df_raw)
- Data Transformation
selectfunction을 사용하여 아래와 같은 변환을 합니다.- recode :
lfp를 text에서 binary로 변환 - add square of :
age - change :
foreign컬럼을 categorical로 변환 - 그외 다른 모든 컬럼들은 그대로 둔다.
- 컬럼변환 옵션의 일반적인 syntax
source columns => transformation => target columns nameByRawwrapper는 row방향으로select를 수행하도록 합니다.(기본은 colum 방향으로 수행)
df = select(df_raw,
:lfp => (x -> recode(x,"yes" => 1, "no" => 0)) => :lfp,
:lnnlinc,
:age,
:age => ByRow(x -> x^2) => :age²,
Between(:educ, :noc),
:foreign => categorical => :foreign
)
컬럼명을 변경하지 않고 그대로 둔다면 renamecols=false를 인자로 준다.
df = select(df_raw,
:lfp => x->recode(x,"yes"=>1,"no"=>0),
:lnnlinc,
:age,
:age=>ByRow(x->x^2) => :age²,
Between(:educ,:noc),
:foreign=>categorical,
renamecols=false
)
describe(df)
데이터 분석¶
Exploratory data analysis(탐색적 데이터 분석)¶
“ '탐색적 데이터 분석(EDA)’은 우리가 존재한다고 믿는 것들은 물론이고 존재하지 않는다고 믿는 것들을 발견하려는 태도, 유연성, 그리고 자발성이다. “ - 존 튜키 (도서 Doing Data Science 중)
탐색적 데이터 분석이란 벨 연구소의 수학자 존 튜키가 제안한 데이터 분석 방법으로 통계적 가설 검정 등에 의존한 기존 통계학으로는 새롭게 나오는 많은 양의 데이터의 핵심 의미를 파악하는 데 어려움이 있다고 생각하여 이를 보완한 탐색적 데이터 분석을 도입했다고 합니다. 데이터를 분석하고 결과를 내는 과정에서 원 데이터에 대한 탐색과 이해를 기본으로 가지는 것이 가장 중요합니다. 이에 따라 탐색적 데이터 분석은 데이터의 분포와 값을 다양한 각도에서 관찰하며 데이터가 표현하는 현상을 더 잘 이해할 수 있도록 도와주고 데이터를 다양한 기준에서 살펴보는 과정을 통해 문제 정의 단계에서 미처 발견하지 못한 다양한 패턴을 발견하고 이를 바탕으로 기존의 가설을 수정하거나 새로운 가설을 추가할 수 있도록 합니다. 데이터에 대한 관찰과 지식이 이후에 통계적 추론이나 예측 모델 구축 시에도 사용되므로 데이터 분석 단계 중 중요한 단계라고 할 수 있습니다. EDA의 목표는 관측된 현상의 원인에 대한 가설을 제시하고, 적절한 통계 도구 및 기법의 선택을 위한 가이드가 되며, 통계 분석의 기초가 될 가정을 평가하고 추가 자료수집을 위한 기반을 제공합니다.
- 목표에 영향을 미치는 방향을 초기에 확인하기 위해 :lfp로 숫자 열의 평균을 계산하려고 합니다.
@chain df begin
groupby(:lfp)
combine([:lnnlinc, :age, :educ, :nyc, :noc] .=> mean)
end
[:lnnlinc, :age, :educ, :nyc, :noc].=>mean
- 위와 같은 동작을 하는 또 다른 방법
- 실수타입의 모든 컬럼을 추출 : names(df,Real)
@chain df begin
groupby(:lfp)
combine(names(df,Real) .=> mean)
end
- Categorical .변수 :
foreign취급 하기 nrow를combine에 넘겨주어 각 group의 row의 갯수를 return하게 한다.
@chain df begin
groupby([:lfp, :foreign])
combine(nrow)
end
- Cross-Tabulation (교차표)
- 교차분석 : 보통은 크로스탭(crosstab) 분석이라 하고, 2개의 명목 혹은 서열형 척도를 변수로 분석에 활용합니다. 즉 한 변수의 범주를 다른 변수의 범주와 교차 시키고 각각 교차된 경우에 해당하는 셀의 빈도를 분석하는 방법입니다.
- 아래 예는 외국인 여부(foreign)와 경제활동참여(lfp)여부의 빈도(nrow)를 교차표로 나타 냅니다.
- 경제할동에 참여하는(1) 사람들중 외국인(yes)은 147명, 내국인(no)는 254명
- 경제활동에 참여하지 않는(0) 사람들중 외국인(yes)은 69명, 내국인(no)는 402명
@chain df begin
groupby([:lfp, :foreign])
combine(nrow)
unstack(:lfp,:foreign,:nrow)
end
- 경제활동참여 여부에 따른 외국인의 비율 조사
@chain df begin
groupby([:lfp, :foreign])
combine(nrow)
unstack(:lfp,:foreign,:nrow)
select(:lfp,[:no,:yes] => ByRow((x,y)->y/(x+y))=>:foreign_yes)
end
위의 방법을 아래와 같이 간단하세 할 수 있다.
@chain df begin
groupby(:lfp)
combine(:foreign => (x -> mean(x .== "yes")) => :foreign_yes)
end
groupby 함수에 의해 생성된 GroupedDataFrame은 그 자체로 작업하기에 유용한 객체가 될 수 있습니다.
gd = groupby(df,:lfp)
gd를 인덱싱하는 여러 방법을 아래에 보여 준다
gd[1]
gd[(lfp=0,)]
gd[Dict(:lfp => 0)]
gd[(0,)]
- 경제할동참여 여부에 따른 나이별 확률 밀도를 그래프로 나타 냅니다.
- 아래 그래프에서 age는 10으로 나눈값을 사용하고 있으므로 실제는 10대 부터 70대 까지 경제활동에 참여 하는 나이별 확률 밀도를 나타 낸다.
- 아래 그래프에서 경제활동에 참여하는(빨간색 그래프) 가장 많은 나이대는 30대 후반
- 경제활동에 참여하지 않는 가장 많은 나이대는 30대 초반
@df df density(:age, group=:lfp)
예측모델 생성¶
GLM.jl 패키지를 사용한 probit model - 프로빗모형 생성¶
참조: R에서 프로빗 회귀분석(Probit Regression) 실시하기
프로빗 회귀분석(Probit Regression)은 종속변수가 이항형문제(즉, 유효한 범주의 개수가 두개인 경우)를 분류하는 모델로 일반화 선형 회귀모형(Generalized Linear Regression, GLM)중 하나입니다.
프로빗 회귀분석이 일반적인 회귀분석과 가장 크게 차이나는 부분은 종속변수를 단순히 Y로 두는 대신에 정규 누적함수($\Phi$)를 이용합니다. 정규누적함수를 이용하는 이유는 p(0 ~ 1)를 단순하게 종속변수로 둔다면 선형회귀모형을 만들면 다음과 같은 식이 만들어 집니다. $$ p = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_q x_q $$ 하지만 위의 식은 p가 0 ~ 1사이의 값을 갖는다는 것을 보장할 수 없게 되고, 따라서 다음과 같은 비선형함수를 사용합니다. $$ p = \Phi\left(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_q x_q\right) $$ 여기서 $\Phi$는 표준정규분포의 누적분포함수를 의미하고 식은 아래와 같습니다. $$ \Phi(z) = \frac{1}{\sqrt{2\pi}}\int^z_{-\infty}e^{-\frac{t^2}{2}}dt $$ 두번째 식에서 $\Phi(x)$함수의 특성상 $x_1,x_2,...,x_q$가 어떤값을 가지더라도 우변은 항상 0과 1의 값을 가지게 됩니다.. 프로빗 회귀모형은 아래와 같이 데이터 갯수가 n개에 대한 log likelihood 함수를 최대화 시키는 방향으로 추정모수 $\hat{\beta}$를 구하게 됩니다. Probit model 참조
- 다음은 비노동수익(lnnlinc),나이,교육정도(educ)등의 데이터가 주어 졌을 때 경제활동에 참여(lfp=1) 할 확률$Pr(lfp=1)$을 예측하는 것입니다
probit = glm(@formula(lfp ~ lnnlinc + age + age²+ educ + nyc + noc + foreign),
df, Binomial(), ProbitLink())
위 예에서 @formula를 손으로 넣었는데 아래와 같이 프로그램적으로 생성 할 수 있다.
probit = glm(Term(:lfp) ~ sum(Term.(propertynames(df)[2:end])),
df, Binomial(), ProbitLink() )
Note the following:
Term(:lfp) ~ sum(Term.(propertynames(df)[2:end]))
vs.
@formula(lfp ~ lnnlinc + age + age²+educ + nyc + noc + foreign)
마지막으로 @formula가 :age의 제곱을 자동으로 계산할 만큼 충분히 강력하다는 것을 봅시다.
probit = glm(@formula(lfp ~ lnnlinc + age + age^2 + educ + nyc + noc + foreign),
df, Binomial(), ProbitLink())
다시 한번 formula를 체크 하면 아래와 같습니다.
@formula(lfp ~ lnnlinc + age + age^2 + educ + nyc + noc + foreign)
다음으로 우리는 다른 모든 변수를 고정하고 :age를 수정함에 따라 모델의 예측이 어떻게 변경되는지 확인하기위해 새 데이터 프레임을 준비합니다.
예) 비노동소득이 22,026 ($log(22026) \approx 10.0$)이고 연령대가 20 ~ 62, 정규교육을 9년 정도 받고 7세 이상의 자녀를 1명둔 외국인 노동자들의 데이터
df_pred = DataFrame(lnnlinc=10.0, age=2.0:0.01:6.2, educ =9, nyc = 0, noc=1,
foreign="yes")
신뢰 구간과 함께 예측을 수행합니다.
probit_pred = predict(probit,df_pred,interval=:confidence)
데이터가 주어졌을 때 나이에 따른 경제활동에 참여할 확율을 plot으로 나타냅니다.
데이터 프레임에서 행렬(matrix)을 쉽게 만들기 위해 Matrix를 사용합니다.
plot(df_pred.age, Matrix(probit_pred),labels=["lfp" "lower" "upper"],
xlabel="age", ylabel="Pr(lfp=1)")
Advanced DataFrames.jl functionalities: bootstrapping example¶
프로빗 객체를 다시 조사해 보겠습니다.
probit
매개변수에 대한 매개변수 신뢰구간을 얻었음을 알 수 있습니다. 그러나 샘플이 그리 크지 않았기 때문에 부트스트랩을 사용하여 검증하고자 합니다.
먼저 수동으로 부트스트랩을 수행하고 다음으로 결과를 Bootstrap.jl 패키지가 생성하는 것과 비교할 것입니다.
통상 DataFrames.jl의 몇 가지 새로운 기능을 배우려고 할 것입니다. 데이터 프레임을 취하는 함수로 시작하고
1. 콘텐츠의 부트스트랩 샘플 하나를 생성합니다.
2. 프로빗 모델을 부트스트랩된 데이터에 적용 합니다.
3. 계산된 계수와 함께 NamedTuple을 반환합니다. function boot_sample(df)
# df에서 샘플을 df의 크기만큼 무작위로 뽑는다. (복원추출)
df_boot = df[rand(1:nrow(df),nrow(df)),:]
probit_boot = glm(@formula(lfp ~ lnnlinc + age + age^2 + educ + nyc + noc + foreign),
df_boot, Binomial(), ProbitLink())
return (;(Symbol.(coefnames(probit_boot)) .=> coef(probit_boot))...)
end
boot_sample 함수를 여러 번 실행해야 합니다. 결과를 coef_boot 데이터 프레임에 저장합니다.
function run_boot(df,reps)
coef_boot = DataFrame()
for _ in 1:reps
push!(coef_boot, boot_sample(df))
end
return coef_boot
end
Bootstrap.jl이 생성하는 것과 유사한 결과를 원하기 때문에 난수 생성기에 시드를 사용합니다(수동 코드에서 동일한 방식으로 부트스트랩을 위한 샘플 행을 확인했습니다).
Random.seed!(1234)
@time coef_boot = run_boot(df,1000)
이 데이터를 사용하여 백분위수 부트스트랩을 사용하여 95% 신뢰 구간을 계산합니다.
conf_boot = mapcols(x->quantile(x,[0.025,0.975]),coef_boot)
다음은 GLM.jl에 의해 계산된 매개변수 신뢰구간입니다. 부트스트랩 결과와 비교하고 싶습니다.
confint(probit)
conf_param = DataFrame(permutedims(confint(probit)),names(conf_boot))
그리고 conf_boot 데이터 프레임에 'append!' 합니다.
append!(conf_boot,conf_param)
데이터의 각 행이 보유하고 있는 내용을 추적하는 것은 좋습니다. 따라서 데이터 프레임에 새 열을을 삽입합니다. 앞에 놓고 싶을 때 insertcols! 함수를 사용합니다.
insertcols!(conf_boot,1,:statistic=>["boot lo","boot hi","parametric lo", "parametric hi"])
데이터 프레임도 전치될 수 있습니다. 그러나 대상 데이터 프레임에서 열 이름으로 사용할 열을 제공해야 합니다(데이터 프레임 개체에는 열 이름이 있어야 함).
conf_boot_t = permutedims(conf_boot, :statistic)
계수의 추정치를 표에 추가해 보겠습니다.
coef(probit)
insertcols!(conf_boot_t,2,:estimate => coef(probit))
이제 좀 더 발전된 것들을 위한 시간입니다. 신뢰 구간의 끝을 유지하는 열(3에서 6까지의 열)을 추정치의 절대 편차로 변환하려고 합니다. 이러한 변환은 플로팅에 유용합니다.
select!(conf_boot_t,:statistic, :estimate,3:6 .=> x->abs.(x .- conf_boot_t.estimate),
renamecols=false)
scatter(0.05 .+ (1:8),conf_boot_t.estimate,
yerror=(conf_boot_t."boot lo",conf_boot_t."boot hi"),
label="bootstrap",
xticks=(1:8,conf_boot_t.statistic), xrotation=45)
scatter!(-0.05 .+ (1:8), conf_boot_t.estimate,
yerror=(conf_boot_t."parametric lo",conf_boot_t."parametric hi"),
label="parametric")
- 보시다시피 이 경우 두 유형의 간격이 매우 가깝습니다. 마치기 전에 Bootstrap.jl 패키지를 사용하여 동일한 작업을 수행할 수 있었던 방법에 대한 샘플을 살펴보겠습니다. 이번에는 통계를 컴파일하는 함수가 Bootstrap.jl 패키지에 의해 처리되므로 샘플링을 수행할 필요가 없습니다.
function boot_probit(df_boot)
probit_boot = glm(@formula(lfp ~ lnnlinc + age + age^2 + educ + nyc + noc + foreign),
df_boot, Binomial(), ProbitLink())
return (;(Symbol.(coefnames(probit_boot)) .=> coef(probit_boot))...)
end
Random.seed!(1234)
@time bs = bootstrap(boot_probit,df,BasicSampling(1000))
다음으로 95% 백분위수 신뢰 구간을 계산합니다.
bs_ci = confint(bs, PercentileConfInt(0.95))
결과가 수동 간격과 일치하는지 확인하겠습니다. 먼저 튜플이 포함된 데이터 프레임에 새 열을 만듭니다. 이것은 DataFrames.jl에서 문제 없이 처리되는 중첩 데이터 구조의 일반적인 예입니다.
이전 계산과 일치하도록 추정치에서 신뢰 구간의 하한 및 상한 편차를 계산합니다.
conf_boot_t.bootstrap = [(ci[1],ci[1]-ci[2],ci[3]-ci[1]) for ci in bs_ci]
conf_boot_t
하지만 조금 불편합니다.
먼저 :bootstrap 열을 3개의 열로 중첩 해제합니다.
select!(conf_boot_t, Not(:bootstrap), :bootstrap=>["estimate 2","boot lo 2","boot hi 2"])
Next reorder the columns using regular expressions.
select(conf_boot_t,:statistic,r"estimate",r"lo",r"hi")
지원되는 다양한 colummn selector에 대한 자세한 설명은 여기에서 찾을 수 있습니다.
이제 수동 계산이 Bootstrap.jl 패키지와 정확히 동일한 결과를 생성한다는 것을 더 쉽게 알 수 있습니다.
완료하기 전에 :estimate로 데이터 프레임을 정렬합니다.
sort(conf_boot_t,:estimate)
다음은 신뢰 구간의 너비를 기준으로 행을 정렬하는 고급 예입니다.
conf_boot_t[sortperm(conf_boot_t."boot hi" + conf_boot_t."boot lo"),:]
데이터 프레임 정렬의 더 많은 예는 여기에서 찾을 수 있습니다.