10月3日（金）

plot3dについて

先生のアドバイスを頂いて、plot3dについてようやく分かった。

あるサイトに

library(rgl)
 
z <- matrix(1:20, ncol = 4, nrow = 5)

x <- c(1, 2, 3, 4)       # z 行列の各行の座標
y <- c(1, 2, 3, 4, 5)    # z 行列の各列の座標
 
plot3d(x, y, z)

こんなコードがのっていて、これで書けたと思っていたが、よく見たらこのプロット間違っていた。
x、yのベクトルがループして書けているっぽく見えていただけらしい。
とりあえず正しい使い方が分かったのでよかった。

交互作用を含むモデルの重回帰分析

交互作用の最も単純なモデルとして、性別（1か0）によって年齢の係数がそれぞれalpha,betaに変化するようなモデルをつくって回帰分析をしてみた。
（従属変数が連続量のアウトカム、説明変数が二値的な性別と連続量の年齢の2つで、その2つが交互作用するようなモデル）

k <-　0.2
alpha <-　3
beta <-　2
N <- 100

out. <- function(x){
if (x[2]==1) return(k+alpha*x[1]+rnorm(1))
else return(k+beta*x[1]+rnorm(1))
}

age <- sample(20:80,N,replace=T)
sex <- sample(0:1,N,replace=T)
outcome <- apply(cbind(age,sex),1,out.)

d <- data.frame(cbind(age,sex,outcome))

result <- lm(outcome~age+sex+age*sex,d)

plot3d(age,sex,outcome)

r.c <- result$coefficient

x <- 20:80
y <- 0:1
xy <- as.matrix(expand.grid(x,y))

f <- function(x){
r.c[1]+r.c[2]*x[1]+r.c[3]*x[2]+r.c[4]*x[1]*x[2]
}

z <- apply(xy,1,f)

points3d(xy[,1],xy[,2],z,col=2)

黒がデータ点で、赤が回帰モデルによって求めた直線（面）をプロットしたものである。

> summary(result)

Call:
lm(formula = outcome ~ age + sex + age * sex, data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.68398 -0.73942 -0.00973  0.67727  2.69020 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  0.669753   0.437379    1.531    0.129    
age         -2.006617   0.008452 -237.420   <2e-16 ***
sex         -0.895740   0.668276   -1.340    0.183    
age:sex      5.015197   0.012611  397.697   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.056 on 96 degrees of freedom
Multiple R-squared:  0.9999,    Adjusted R-squared:  0.9999 
F-statistic: 4.962e+05 on 3 and 96 DF,  p-value: < 2.2e-16

おそらくここに回帰の妥当性にあたる部分が書かれていると思われるが、解読は不能。

試しにこのデータを、交互作用がないようなモデル式に当てはめてみた。

k <-　0.2
alpha <-　3
beta <-　-2
N <- 100

out. <- function(x){
if (x[2]==1) return(k+alpha*x[1]+rnorm(1))
else return(k+beta*x[1]+rnorm(1))
}

age <- sample(20:80,N,replace=T)
sex <- sample(0:1,N,replace=T)
outcome <- apply(cbind(age,sex),1,out.)

d <- data.frame(cbind(age,sex,outcome))

result <- lm(outcome~age+sex,d)

summary(result)

plot3d(age,sex,outcome)

r.c <- result$coefficient

x <- 20:80
y <- 0:1
xy <- as.matrix(expand.grid(x,y))

f <- function(x){
r.c[1]+r.c[2]*x[1]+r.c[3]*x[2]
}

z <- apply(xy,1,f)

points3d(xy[,1],xy[,2],z,col=2)

結果はこのように、無残なことに。

> summary(result)

Call:
lm(formula = outcome ~ age + sex, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-78.094 -37.588  -2.309  38.346  93.767 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -104.9162    13.8046  -7.600 1.88e-11 ***
age            0.1241     0.2511   0.494    0.622    
sex          239.7482     8.8818  26.993  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 44.21 on 97 degrees of freedom
Multiple R-squared:  0.8828,    Adjusted R-squared:  0.8804 
F-statistic: 365.4 on 2 and 97 DF,  p-value: < 2.2e-16

明らかに回帰の妥当性としては劣っているはずだが、どこにその情報があるのか…？
p値をガリガり計算してもよいが、かなり面倒である。

ちなみに、交互作用がないモデルを交互作用を含むようなモデル式に当てはめた場合は、

k0 <-　0.2
k1 <-　0.3
k2 <-　2
N <- 100

out. <- function(x){
k0+k1*x[1]+k2*x[2]+rnorm(1)
}

age <- sample(20:80,N,replace=T)
sex <- sample(0:1,N,replace=T)
outcome <- apply(cbind(age,sex),1,out.)

d <- data.frame(cbind(age,sex,outcome))

result <- lm(outcome~age+sex+age*sex,d)

plot3d(age,sex,outcome)

r.c <- result$coefficient

x <- 20:80
y <- 0:1
xy <- as.matrix(expand.grid(x,y))

f <- function(x){
r.c[1]+r.c[2]*x[1]+r.c[3]*x[2]+r.c[4]*x[1]*x[2]
}

z <- apply(xy,1,f)

points3d(xy[,1],xy[,2],z,col=2)

> result

Call:
lm(formula = outcome ~ age + sex + age * sex, data = d)

Coefficients:
(Intercept)          age          sex      age:sex  
  -0.195967     0.305939     1.750063     0.003695  

予想通りではあるが、age*sexの項の係数はほぼ0になって、それなりに妥当な回帰が出来た。

来週の予定としては、まず「回帰の妥当性」についてもう一度考えた後、「傾向スコア分析」なるものについて勉強しようと思う。
それが終わったら論文をいくつか読んで、あとは最終レポート作成…という流れになりそうである。

以上