1 15 R Part : website:

Transcription

1 1 15 R Part : website: kkarato@eco.u-toyama.ac.jp

2 1 2 1 R R R 3. r 4. data 11.csv data<- read.csv("data_11.csv") read.csv CSV(Comma Separated Value ) data R ctrl + R attach(data) obs<- length(data[,1]) str(data) # summary(data) 4 4 attach length(data[,1]) 1 str data # R summary data str(data)

3 : curve, lines, abline x-y curve y = 3 2x curve(3-2*x) lines 2 lines(c(0.2,0.8),c(2,2.5)) 1 (0.2, 2) 2 (0.8, 2.5) 2 abline 1 2 (y = 1 + 2x) abline(1,2) abline(v=0.5) x = 0.5 (vertical line) abline(h=1.5) y = 1.5 (horizontal line) 2 2 abline [18]age.sf [25]age.sp plot(age.sf, age.sp) reg<- lm(age.sf~age.sp) abline(reg) lm (fitting linear models) 2 age.sf~age.sp formula age.sf = α + β age.sp lm reg abline abline(lm(age.sf age.sp)) lm reg

4 2 4 reg Call: lm(formula = age.sf ~ age.sp) Coefficients: (Intercept) age.sp summary lm summary summary(reg) summary(lm(age.sf~age.sp)) summary(reg) Call: lm(formula = age.sf ~ age.sp) Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) e-08 *** age.sp < 2e-16 *** --- Signif. codes: 0 *** ** 0.01 * Residual standard error: on 2160 degrees of freedom Multiple R-squared: , Adjusted R-squared: F-statistic: 1.877e+04 on 1 and 2160 DF, p-value: < 2.2e-16 α, β u i age.sf = α + β age.sp + u i Residuals 2 Coefficients α, β (Estimate) (Std. Error) t (t value) p (Pr(> t )) (Intercept) ˆα age.sp

5 2 5 age.sp ˆβ * Signif. codes *** 0.1% ** 1% * 5%. 10% Residual standard error ˆσ = û 2 i / degrees of freedom Multiple R-squared R 2 Adjusted R-squared adj.r 2 F-statistic 0 H 0 : β = 0 F DF p summary summary names(lm(formula)) names(reg) [1] "coefficients" "residuals" "effects" "rank" [5] "fitted.values" "assign" "qr" "df.residual" [9] "xlevels" "call" "terms" "model" summary Value summary(reg)[6] summary(reg)$sig summary(reg)$sigma 1 1 summary(lm(formula)) Value [3] $res [4] $coef t p [6] $sig [8] $r.sq [9] $adj.r.sq [10] $fstat F [11] $cov

6 [21] hour.sf [18] age.sf [5] sex D<- ifelse(sex==2,1,0) sex==2 D hour.sf!=888 y<- hour.sf[hour.sf!=888] x<- age.sf[hour.sf!=888] D<- D[hour.sf!=888] y = β 1 + β 2 x + β 3 D + u (1) reg<- lm(y~x + D) summary(reg) y x2, x3, x4 y~x2+x3+x4 formula subset subset pref==16 [21] hour.sf subset(hour.sf,pref==16) pref==16 sex==1 [21] hour.sf subset(hour.sf,pref==16 & sex==1)

7 2 7 data subset(data,pref==16 & sex==1) hour.sf==888 data subset(data,hour.sf!=888) data.frame factor csv data<- read.csv("data_11.csv") data data.frame subset data.frame data2 data2<- data.frame(subset(data, hour.sf!=888)) (1) reg<- lm(hour.sf~age.sf + factor(sex), data=data2) summary(reg) lm 2 data=data2 factor sex 1 2 factor 1 2 lm formula factor Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) <2e-16 *** age.sf <2e-16 *** factor(sex) <2e-16 *** lm(y~ x + D) factor(sex)2 sex==2 hour.sf

8 hour.sf = β 1 age.sf + β 2 age.sf 2 + β 3 factor(sex)2 + u (2) 2 formula I (Inhibit Interpretation/Conversion of Objects) I(age.sf^2) reg<- lm(hour.sf~age.sf + I(age.sf^2) + factor(sex), data=data2) summary(reg) (2) Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) e-06 *** age.sf e-07 *** I(age.sf^2) e-11 *** factor(sex) < 2e-16 *** 2 summary(reg)$coef β 2 b2 β 3 b3 b2<- summary(reg)$coef[2,1] b3<- summary(reg)$coef[3,1] b2 = b3 = Ŷ Ŷ = 1.170x 0.015x 2 + extras (3) x extras curve(b2*x+b3*x^2, 20,80, axes=f) axis(1) axes=f axis(1) 30 (3) x 0 x = R -b2/(2*b3) [1]

9 2 9 1 data2 hour.sf 1. 2 (job.sf) (age.sf10) (job.sf) factor 1 1 reg<- lm(hour.sf~age.sf + I(age.sf^2) + factor(sex) + factor(job.sf), data=data2) summary(reg) Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) e-11 *** age.sf e-06 *** I(age.sf^2) e-09 *** factor(sex) < 2e-16 *** factor(job.sf) e-05 *** factor(job.sf) < 2e-16 *** factor(job.sf) ** factor(job.sf) e-10 *** factor(job.sf) ** factor(job.sf)2 factor(job.sf)6 job.sf== reg<- lm(hour.sf~factor(age.sf10) + factor(sex) + factor(job.sf), data=data2) summary(reg) age.sf10 10 factor(age.sf10)30 factor(age.sf10)80 20 age.sf10== , 40,

10 3 10 Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) < 2e-16 *** factor(age.sf10) factor(age.sf10) factor(age.sf10) factor(age.sf10) e-05 *** factor(age.sf10) e-08 *** factor(age.sf10) e-05 *** factor(sex) < 2e-16 *** factor(job.sf) e-05 *** factor(job.sf) < 2e-16 *** factor(job.sf) *** factor(job.sf) e-10 *** factor(job.sf) ** work.sf work.sf==1 work.sf==2 work.sf==3 { 1 Y i = 0 i = 1, 2,, n (4) 2 Y i = α + βx i + u i (5) (5) Y i 0 1 Y i 1 p i Y i = 0 Y i = 1 1 E(Y i ) = p i V (Y i ) = p i (1 p i )

11 3 11 (5) 2 Y i = α + βx i + u i (6) Y i < 0 Y i = 0 Y i 0 Y i = 1 Yi (latent variable) Yi 2 Y i α + βx i + u i Yi < 0 Y i = 0 Y i 0 Y i = 1 0 u i (5) 2 2 (generalized linear model; GLM) 3 Y i X i α + βx i Y i E(Y i ) = p i g g(e(y i )) = α + βx i (7) g u i N(0, σ 2 ) E(Y i ) = α + βx i (4) Pr(Y i = 1) = p i p i = eα+βx i 1+e α+βx i ln *1 E(Y i ) = p i g(p i ) p i 1 p i ln p i 1 p i = α + βx i (8) *1 10

12 (maximum likelihood method; ML) (likelihood) * 2 L = Pr(Y 1 = y 1 ) Pr(Y 2 = y 2 ) Pr(Y n = y n ) (9) y i 0 1 n p 2 B(n, p) Y y Pr(Y = y) = n C y p y (1 p) n y n = 1 Y = 0, 1 Pr(Y = 0) = p 0 (1 p) 1 0 = 1 p, Pr(Y = 1) = p 1 (1 p) 1 1 = p i (Y i = 1) (Y i = 0) Pr(Y i = y i ) = p y i i (1 p i) 1 y i (10) (7) X i Y i = 1 Pr(Y i = 1) = Pr(Y i 0) (11) = Pr(α + βx i + u i 0) = Pr(u i α βx i ) = p i Y i = 0 Pr(Y i = 0) = Pr(Y i < 0) (12) = Pr(α + βx i + u i < 0) = Pr(u i < α βx i ) = 1 p i (12) (13) (10) L = Pr(Y 1 = y 1 ) Pr(Y 2 = y 2 ) Pr(Y n = y n ) = Pr(u 1 α βx 1 ) y 1 Pr(u 1 < α βx 1 ) 1 y 1 Pr(u 2 α βx 2 ) y 2 Pr(u 2 < α βx 2 ) 1 y 2 Pr(u n α βx n ) y n Pr(u n < α βx n ) 1 y n X 1,, X n L α, β L α, β LL = ln L α, β ˆα, ˆβ *2

13 β = 0 0 (likelihood ratio test; LR) LR = 2(LL0 LL) χ 2 (k) (13) LL0 0 k k 2 F (Akaike s Information Criterion; AIC) AIC = 2LL + 2(k + 1) (14) AIC 2 (residual deviance) Dev Dev = 2LL (15) 2 LR AIC LR = Dev0 Dev, AIC = Dev + 2k (McFadden s pseudo R-squared) pr 2 = 1 LL LL0 (16) LR glm R glm y<- ifelse(work.sf==1 1,0) y.glm<- glm(y~ age.sf + factor(sex), family=binomial(link="logit")) summary(y.glm) 2 (4) work.sf y y age.sf factor(sex) glm family 2 binomial link="logit"

14 3 14 Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) <2e-16 *** age.sf <2e-16 *** factor(sex) <2e-16 *** --- Signif. codes: 0 *** ** 0.01 * (Dispersion parameter for binomial family taken to be 1) Null deviance: on 2161 degrees of freedom Residual deviance: on 2159 degrees of freedom AIC: Number of Fisher Scoring iterations: 4 z value 2 k = 2 Null deviance 2 0 (Dev0 = 2LL0) n = 2162 Null deviance n 1 Residual deviance LL Dev = 2LL n (k + 1) AIC = Dev + 2(k + 1) = sumamry deviance LL<- y.glm$dev/(-2) LL0<- y.glm$null.dev/(-2) LR<- -2*(LL0-LL) bhat<- y.glm$coef k<- length(bhat)-1 LRp<- 1-pchisq(LR,df=k) AIC<- y.glm$aic pr2<- 1 - LL/LL0 $dev (15) y.glm$dev/(-2) LL $coef length 1 k = 2 (13) p k 2 $aic (16)

15 ME *** *** *** LL LR LR p AIC pr X 1 Pr(Y = 1) (marginal effects) 10 ME = Pr(Y i = 1) X i = ˆβp (1 p ) (17) p (1 p ) X i ˆp i = Pr(u i ˆα ˆβX i ) p (1 p ) = i ˆp i(1 ˆp i ) n y.glm$fit phat<-y.glm$fit ME<- bhat*mean(phat*(1-phat)) ME (Intercept) age.sf D bhat phat mean(phat*(1-phat)) age.sf (1.43%) D (25.44%) ln p i 1 p i = β 1 + β 2 age.sf i + β 3 D i + u i 2

16 ownhouse==1 y = 1 y = 0 2 age.sf factor(income.hh) factor(pref) 3.5 (12) p i (12) p i N(0, σ 2 ) (probit model) p i = Pr(u i α + βx i ) = α+βxi f(u i )du i (18) f(u i ) (18) f(u i ) F (α + βx i ) α + βx i F 1 (p i ) F 1 (p i ) = α + βx i (19) glm y 2 y.glm<- glm(y~ age.sf + factor(sex), family=binomial(link="probit")) summary(y.glm) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) <2e-16 *** age.sf <2e-16 *** factor(sex) <2e-16 *** AIC R ME = Pr(Y i = 1) X i = ˆβ i f(α + βx i) n (20)

17 *** *** *** LL LR LR p AIC pr phat<-y.glm$fit ME<- bhat*mean(dnorm(qnorm(phat))) ME (Intercept) age.sf factor(sex) op56: Y i = {1, 2, 3, 4, 5} Y i 2 (7) Y i = βx i + u i (21) ζ j 1 < Y i < ζ j Y i = j (j = 1, 2, 3, 4, 5) (breakpoints) ζ 0, ζ 1, ζ 2, ζ 3, ζ 4, ζ 5 j

18 3 18 Y i = 1 ζ 0 < Yi < ζ 1 ζ 0 βx i < u i < ζ 1 βx i Y i = 2 ζ 1 < Yi < ζ 2 ζ 1 βx i < u i < ζ 2 βx i Y i = 3 ζ 2 < Yi < ζ 3 ζ 2 βx i < u i < ζ 3 βx i Y i = 4 ζ 3 < Yi < ζ 4 ζ 3 βx i < u i < ζ 4 βx i Y i = 5 ζ 4 < Yi < ζ 5 ζ 4 βx i < u i < ζ 5 βx i ζ 0 =, ζ 5 = ζ 1, ζ 2, ζ 3, ζ 4 Λ Y i = j X i p ij = Pr(Y i = j) = Λ(ζ j βx i ) Λ(ζ j 1 βx i ) (22) F p ij = Pr(Y i = j) = F (ζ j βx i ) F (ζ j 1 βx i ) (23) 2 (p y i1 i1 ) (py i2 i2 ) (py i3 i3 ) (py i4 i4 ) (py i5 i5 ) (22) (23) AIC = 2LL + 2(k + k ζ ) (24) k k ζ 5 ζ 1, ζ 2, ζ 3, ζ 4 : : Pr(Y i = j) = [λ(ˆζ j 1 X ˆβX i ) λ(ˆζ j ˆβX i )] ˆβ i (25) Pr(Y i = j) = [f(ˆζ j 1 X ˆβX i ) f(ˆζ j ˆβX i )] ˆβ i (26) λ f op56 age.sf facto(sex) 6 c c c1218 factor(edu.sf) pkg MASS R MASS library

19 3 19 library(mass) polr ologit<- polr(as.factor(op56)~ age.sf + factor(sex) + c06 + c c factor(edu.sf), method = "logistic") summary(ologit) op56 as.ordered method = "logistic" method = "probit" AIC LL<- ologit$deviance/(-2) ologit0<- polr(as.ordered(op56)~ 1, method = "logistic") LL0<- ologit0$deviance/(-2) LR<- -2*(LL0-LL); # LRp<- 1-pchisq(LR,df=k) # H0: 0 p bhat<- matrix(ologit$coef) # k<- length(bhat) # kz<- length(ologit$zeta) # - 1 AIC<- -2*LL + 2*(k+kz) # pr2<- 1-(LL/LL0) # polr 0 (null.deviance) polr(as.ordered(op56)~ 1, method = "logistic") LL0 u i ologit$fit (22) phat<- ologit$fit cphat<- array(0,c(ologit$n,4)) cphat[,1]<- phat[,1] cphat[,2]<- phat[,1] + phat[,2] cphat[,3]<- phat[,1] + phat[,2] + phat[,3] cphat[,4]<- phat[,1] + phat[,2] + phat[,3] + phat[,4] lambda<- dlogis(qlogis(cphat)) lambda 4 (25) λ(ζ 1 βx i ) λ(ζ 2 βx i ) λ(ζ 3 βx i ) λ(ζ 4 βx i )

20 3 20 dl1<- mean(-lambda[,1]) dl2<- mean(lambda[,1]-lambda[,2]) dl3<- mean(lambda[,2]-lambda[,3]) dl4<- mean(lambda[,3]-lambda[,4]) dl5<- mean(lambda[,4]) (25) ME ME<- cbind(bhat*dl1, bhat*dl2, bhat*dl3, bhat*dl4, bhat*dl5) ME [,1] [,2] [,3] [,4] [,5] [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] % 0.10% 0.11% 0.10% % 2.20% 2.33% 2.25%

21 Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 age.sf ** factor(sex) c c *** c *** factor(edu.sf) factor(edu.sf) factor(edu.sf) factor(edu.sf) : 1 2 ζ ζ ζ ζ LL LR LR p AIC pr % 2.47% 2.62% 2.52% 0.38% 0.28% 0.30% 0.28% 4 job.sati data3<- data.frame(subset(data, work.sf == 1)) 5 1 as.ordered(6-job.sati) age.sf factor(sex) factor(income.sf) hour.sf library(mass)

22 3 22 > data3<- data.frame(subset(data, work.sf == 1)) > ologit<- polr(as.ordered(6-job.sati)~ age.sf + factor(sex) + + factor(income.sf) + hour.sf +, method = "logistic", data=data3) > summary(ologit) Re-fitting to get Hessian Call: polr(formula = as.ordered(6 - job.sati) ~ age.sf + factor(sex) + factor(income.sf) + hour.sf, data = data3, method = "logistic") Coefficients: Value Std. Error t value age.sf factor(sex) factor(income.sf) factor(income.sf) factor(income.sf) hour.sf Intercepts: Value Std. Error t value Residual Deviance: AIC: