10月2日（木）

回帰モデルによる交絡の制御　2

ただの入力ミスだった。

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

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

a1.d <- sample(20:80,N,replace=T)
a2.d <- sample(0:1,N,replace=T)
a3.d <- apply(cbind(a1.d,a2.d),1,out.)

d <- cbind(a1.d,a2.d,a3.d)

x <- sample(which(d[,1]<=50&d[,2]==1),n)
x. <- sample(which(d[,1]>50&d[,2]==0),n)

age <- d[x,1]
age. <- d[x.,1]

f1 <- function(k1,b1){
k1*age+b1
}

f2 <- function(k1,b2){
k1*age.+b2
}

P <- function(X){
sum((d[x,3]-f1(X[1],X[2]))^2)+sum((d[x.,3]-f2(X[1],X[3]))^2)
}
op <- optim(c(0.3,0.2,3.2),P)

E <- op$par[2]-op$par[3]
E

> E
[1] 2.22436

値をいろいろ変えてやってみたが、最初にやったやり方よりも安定して正確な値が出ているようである。
やはり傾きが一緒という条件は入れたほうがよいらしい。

理論が分かっていれば、このように状況に合わせて応用することができる。
やはり勉強は大事だと思った。

分散分析

分散分析について勉強したので、上でつくったモデルの年齢の項の係数を0として用いて、曝露によって差が出るかどうかの分散分析を実際にやってみることにした。
（サンプル数は100とした）

k0 <-　0.2
k1 <-　0
k2 <-　5
N <- 100

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

a1.d <- sample(20:80,N,replace=T)
a2.d <- sample(0:1,N,replace=T)
a3.d <- apply(cbind(a1.d,a2.d),1,out.)

d <- cbind(a1.d,a2.d,a3.d)

a <- mean(d[,3])
b <- mean(d[which(d[,2]==1),3])
c <- mean(d[which(d[,2]==0),3])

e.ef <- b-a
n.ef <- c-a

e.er <- d[which(d[,2]==1),3]-b
n.er <- d[which(d[,2]==0),3]-c

MS.ef <- sum((e.ef)^2)+sum((n.ef)^2)
MS.er <- (sum((e.er)^2)+sum((n.er)^2))/(N-2)

F <- MS.ef/MS.er

F>qf(0.95,1,N-2)

> F>qf(0.95,1,N-2)
[1] TRUE

これで一応、曝露による効果がありそうであることが分かった。
ただ、k2=2などの小さい値だと、有意差なしという結果になってしまうため、あくまで分散分析は効果があるかどうかの目安としてしか利用できないようである。

重回帰分析

昨日今日と、曝露・非曝露群で年齢に重なりがないときについていろいろとやったが、線形モデルに従うと考えているのだから、重回帰分析すればよいだけのことに気が付いた。

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

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

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

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

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

result$coefficients

> result$coefficients
(Intercept)         age    exposure 
  0.2570855   0.2991444   2.0300088 

Ｒって便利だな、と改めて思った。

以上