LM

#install.packages(c("psych", "ppcor", "car"))
library(psych)
library(ppcor)

## Loading required package: MASS

library(car)

## Loading required package: carData

##
## Attaching package: 'car'

## The following object is masked from 'package:psych':
##
##     logit

データ生成

#a <- rnorm(100, 0, 1)
#b <- 1.5 * a + rnorm(100, 7, 1)
#c <- rnorm(100, 5, 0.3)
#y <- 2.0 * a + 3.5 * b + 0.8 * c + rnorm(100, 23, 1)
#d <- data.frame(a = a, b = b, c = c, y = y)
load("env.RData")
head(d)

##             a        b        c        y
## 1  1.75552469 9.698773 4.717665 63.19326
## 2  1.22438170 9.333587 5.376546 62.64162
## 3 -0.26704143 6.173867 4.685969 47.13082
## 4 -1.60687568 7.024980 4.909123 48.55534
## 5  0.07469729 7.123100 5.002922 52.19395
## 6  0.75229215 7.910237 5.373781 55.88746

単相関を確認

cor(d, use="complete.obs", method = "pearson")

##            a           b           c           y
## a  1.0000000  0.80664757 -0.02110610  0.87392995
## b  0.8066476  1.00000000 -0.04288845  0.98449714
## c -0.0211061 -0.04288845  1.00000000 -0.01027933
## y  0.8739300  0.98449714 -0.01027933  1.00000000

こっちのほうが見やすい。

psych::pairs.panels(d, rug = F, ellipses = F, lm = T, method = "pearson")

aとbの生成メカニズムのせいで、yに対する相関係数が高くなっている。また、本来cは正の相関があるべきだが、yからcをみると負の相関と判断される。

偏相関係数

偏相関係数を計算してみると、別の見方ができる。

ppcor::pcor(d, method = "pearson")

## $estimate
##            a          b          c         y
## a  1.0000000 -0.6465238 -0.1870975 0.7787333
## b -0.6465238  1.0000000 -0.2609780 0.9749029
## c -0.1870975 -0.2609780  1.0000000 0.2580397
## y  0.7787333  0.9749029  0.2580397 1.0000000
##
## $p.value
##              a            b           c            y
## a 0.000000e+00 6.451696e-13 0.065076464 3.773494e-21
## b 6.451696e-13 0.000000e+00 0.009444017 1.943756e-64
## c 6.507646e-02 9.444017e-03 0.000000000 1.030903e-02
## y 3.773494e-21 1.943756e-64 0.010309028 0.000000e+00
##
## $statistic
##           a         b         c         y
## a  0.000000 -8.303401 -1.866127 12.162166
## b -8.303401  0.000000 -2.648849 42.905400
## c -1.866127 -2.648849  0.000000  2.616885
## y 12.162166 42.905400  2.616885  0.000000
##
## $n
## [1] 100
##
## $gp
## [1] 2
##
## $method
## [1] "pearson"

aの係数が大幅に低くなり、bの係数は微減、cの係数は正となった。

重回帰分析

さらに重回帰分析を掛けてみる。aとbの間に多重共線性があるはずだが、試してみる。

d.lm <- lm(y ~ ., data = d)
summary(d.lm)

##
## Call:
## lm(formula = y ~ ., data = d)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -2.78876 -0.56731  0.03203  0.42872  2.55007
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.7195     1.6041  14.163   <2e-16 ***
## a             1.8765     0.1543  12.162   <2e-16 ***
## b             3.5485     0.0827  42.905   <2e-16 ***
## c             0.7697     0.2941   2.617   0.0103 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8813 on 96 degrees of freedom
## Multiple R-squared:  0.9883, Adjusted R-squared:  0.9879
## F-statistic:  2702 on 3 and 96 DF,  p-value: < 2.2e-16

うまく言っているようにみえる。

多重共線性を確認

a,bが6。cは1。a,bが高めだが、vifの数値だけみると取り除くか迷う。

car::vif(d.lm)

##        a        b        c
## 2.864202 2.868202 1.002366