データの準備
AERパッケージのCPS1985データを利用する。
library(AER) data(CPS1985) d1 <- CPS1985
caretパッケージのdummyVars関数を用いて、factor型のものをすべてダミー変数化する。
library(caret) dummy <- dummyVars(~.,data=d1) d2 <- as.data.frame(predict(dummy, d1)) str(d2)
このようにダミー変数化される。
'data.frame': 534 obs. of 24 variables: $ wage : num 5.1 4.95 6.67 4 7.5 ... $ education : num 8 9 12 12 12 13 10 12 16 12 ... $ experience : num 21 42 1 4 17 9 27 9 11 9 ... $ age : num 35 57 19 22 35 28 43 27 33 27 ... $ ethnicity.cauc : num 0 1 1 1 1 1 1 1 1 1 ... $ ethnicity.hispanic : num 1 0 0 0 0 0 0 0 0 0 ... $ ethnicity.other : num 0 0 0 0 0 0 0 0 0 0 ... $ region.south : num 0 0 0 0 0 0 1 0 0 0 ... $ region.other : num 1 1 1 1 1 1 0 1 1 1 ... $ gender.male : num 0 0 1 1 1 1 1 1 1 1 ... $ gender.female : num 1 1 0 0 0 0 0 0 0 0 ... $ occupation.worker : num 1 1 1 1 1 1 1 1 1 1 ... $ occupation.technical : num 0 0 0 0 0 0 0 0 0 0 ... $ occupation.services : num 0 0 0 0 0 0 0 0 0 0 ... $ occupation.office : num 0 0 0 0 0 0 0 0 0 0 ... $ occupation.sales : num 0 0 0 0 0 0 0 0 0 0 ... $ occupation.management: num 0 0 0 0 0 0 0 0 0 0 ... $ sector.manufacturing : num 1 1 1 0 0 0 0 0 1 0 ... $ sector.construction : num 0 0 0 0 0 0 0 0 0 0 ... $ sector.other : num 0 0 0 1 1 1 1 1 0 1 ... $ union.no : num 1 1 1 1 1 0 1 1 1 1 ... $ union.yes : num 0 0 0 0 0 1 0 0 0 0 ... $ married.no : num 0 0 1 1 0 1 1 1 0 1 ... $ married.yes : num 1 1 0 0 1 0 0 0 1 0 ...
Mplusへのデータの出力
まずは変数名を書き出す。Mplusはコードに変数名を入れる必要があるので行う作業である。MplusAutomationパッケージを使って変数名を出力しておくと、コードを書くのが非常に楽になる。
library(MplusAutomation) variable.names(d2) # 変数名を書き出し
変数の書き出し結果。エディタなどを利用してMplusで使えるように半角スペース区切りにしていく。
[1] "wage" "education" "experience" "age" "ethnicity.cauc" "ethnicity.hispanic" [7] "ethnicity.other" "region.south" "region.other" "gender.male" "gender.female" "occupation.worker" [13] "occupation.technical" "occupation.services" "occupation.office" "occupation.sales" "occupation.management" "sector.manufacturing" [19] "sector.construction" "sector.other" "union.no" "union.yes" "married.no" "married.yes"
csvで書き出すというのでも良いかもしれない。エクセルを使った方が多少速くなる気もする。
df <- variable.names(d2) write.csv(df, file = "varnames.csv", row.names=TRUE, fileEncoding ="UTF-8")
最新版のExcelでは文字コードがUTF-8になっているので、ここではUTF-8で出力しているが、ExcelのエンコードがShift-JISの人はfileEncoding ="cp932"と書き換える必要がある。
Mplus用にデータ本体を出力する。
prepareMplusData(d2, filename="CPS1985.dat",
overwrite=T)
重回帰分析
Rでの実行例
分析例は、賃金が従属変数、性別(男)と年齢が独立変数になった重回帰分析である。
lm01 <- glm(wage ~ factor(gender.male) + age, data=d2) summary(lm01)
結果。
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.65425 0.76105 6.116 1.87e-09 ***
gender.male 2.27469 0.43029 5.286 1.82e-07 ***
age 0.08522 0.01830 4.656 4.08e-06 ***
Mplusでの実行例
同じことをMplusで行う。まずはコード。
TITLE: Linear Regression
DATA: file =CPS1985.dat;
VARIABLE:
NAMES = wage education experience age ethnicity_cauc ethnicity_hispanic
ethnicity_other region_south region_other gender_male gender_female
occupation_worker occupation_technical occupation_services occupation_office
occupation_sales occupation_management sector_manufacturing
sector_construction sector_other union_no union_yes married_no married_yes;
USEVARIABLES= wage gender_male age;
MODEL: wage on gender_male age;
ANALYSIS: estimator = ML;
OUTPUT: TECH1 SAMPSTAT STDYX RESIDUAL CINTERVAL;
PLOT: TYPE = PLOT1;
結果。
MODEL RESULTS
Two-Tailed
Estimate S.E. Est./S.E. P-Value
WAGE ON
GENDER_MAL 2.275 0.429 5.302 0.000
AGE 0.085 0.018 4.670 0.000
Intercepts
WAGE 4.654 0.759 6.133 0.000
Residual Variances
WAGE 24.257 1.484 16.341 0.000
Rでの推定結果と見比べると、推定値は全く同じである。標準誤差は微妙に違う。Rで通常の重回帰をすると最小二乗法を使うが、Mplusの方は最尤法だからだ。
ちなみに、ロバスト最尤法＝MLRを使うと下記のようになる。さきほどのコードをANALYSIS: estimator = MLR;と書き換えるだけである。
MODEL RESULTS
Two-Tailed
Estimate S.E. Est./S.E. P-Value
WAGE ON
GENDER_MAL 2.275 0.434 5.244 0.000
AGE 0.085 0.020 4.237 0.000
Intercepts
WAGE 4.654 0.847 5.494 0.000
Residual Variances
WAGE 24.257 3.195 7.591 0.000
推定値は同じだが、標準誤差が大きく違っている。MLRの方が標準誤差が大きい。標準誤差は小さい方がよいため、MLでの推定の方がよいととりあえず判断してよいだろう。
ロジスティック回帰分析
Rでの実行例
lm04 <- glm(factor(married.yes) ~ factor(gender.male) + age, family=binomial, data=d2) summary(lm04)
結果。
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.420497 0.353612 -4.017 5.89e-05 ***
factor(gender.male)1 0.049341 0.191166 0.258 0.796
age 0.057319 0.009245 6.200 5.65e-10 ***
Mplusでの実行例
コード。
title: Logistic Regression
data: file =CPS1985.dat;
variable:
names = wage education experience age ethnicity_cauc ethnicity_hispanic
ethnicity_other region_south region_other gender_male gender_female
occupation_worker occupation_technical occupation_services occupation_office
occupation_sales occupation_management sector_manufacturing
sector_construction sector_other union_no union_yes married_no married_yes;
USEVARIABLES= married_yes gender_male age;
categorical = married_yes;
model:
married_yes on gender_male age;
analysis: estimator = ml;
output: tech1 sampstat standardized;
結果。
MODEL RESULTS
Two-Tailed
Estimate S.E. Est./S.E. P-Value
MARRIED_YE ON
GENDER_MAL 0.049 0.191 0.258 0.796
AGE 0.057 0.009 6.199 0.000
Thresholds
MARRIED_$1 1.420 0.354 4.017 0.000
推定値が同じというのは先ほどの重回帰と同じである。 気をつけるべきところは、Mplusでは切片(Intercept)という表現ではなく閾値(Thresholds)になっていることである。値は、絶対値は同じだが、符号のプラスマイナスが逆である。もし、Mplusで分析した結果をレポートなり、論文なりに切片として表記する場合には、符号を逆転させる必要性がある。
