10月9日（木）
傾向スコア分析
傾向スコアについてＲでやっているサイトがあったので、それをほぼもろパクリして自分でやってみた。
install.packages("Matching") library(Matching) data(lalonde) logi <- glm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,family=binomial,data=lalonde) Y78 <- lalonde$re78 Tre <- lalonde$treat lm(Y78~Tre) mout1 <- Match(Y=Y78,Tr=Tre,X=logi$fitted) summary(mout1) > lm(Y78~Tre) Call: lm(formula = Y78 ~ Tre) Coefficients: (Intercept) Tre 4555 1794 > summary(mout1) Estimate... 2138.6 AI SE...... 797.76 T-stat..... 2.6807 p.val...... 0.0073468 Original number of observations.............. 445 Original number of treated obs............... 185 Matched number of observations............... 185 Matched number of observations (unweighted). 322
ちなみに
mout1 <- Match(Y=Y78,Tr=Tre,X=logi$fitted)
には感激した。
これで層化した分析をやってくれているようである。（どうやっているのか厳密なことはよく分からないが…）
これによって、「だいたい同じ条件を持っている人」同士で比較ができるので、バイアスの影響を取り除くことができる、ということらしい。
ただし、「アウトカムに影響しないが曝露率には影響する因子」「アウトカムに影響するが曝露率には影響しない因子」が含まれていたり、あるいは重要な交絡因子が含まれていなかったりすると不都合が起こるはずである。
ちょっと直観的にわかりにくい部分があるので、Ｒでシミュレーションして確かめてみる。
まずは理想的な場合においてのデータ作成＆効果の測定
kouraku1,kouraku2によって曝露率e、outcomeが影響を受けているモデル
effect <- 2 kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) summary(mout1) > summary(mout1) Estimate... 2.0176 AI SE...... 0.018846 T-stat..... 107.06 p.val...... < 2.22e-16 Original number of observations.............. 100 Original number of treated obs............... 77 Matched number of observations............... 77 Matched number of observations (unweighted). 77
これは申し分ない結果。
感動的である。ちなみにこのモデルで交絡を無視すると
> lm(outcome~bakuro) Call: lm(formula = outcome ~ bakuro) Coefficients: (Intercept) bakuro 0.6888 2.3023
こんなかんじ。
それぞれ何度も繰り返して箱ひげ図を書いたものが下の通り。
ここで、「アウトカムに影響しないが曝露率には影響する因子」あり（Ａ）、「アウトカムに影響するが曝露率には影響しない因子」あり（Ｂ）、交絡因子が抜けている（Ｃ）のそれぞれの場合についても、同じようにつくって箱ひげ図を描いて並べてみた。
effect <- 2 #理想的 hako1 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2+insi1,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } #Ａ hako3 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) insi1 <- runif(100,0,1) x <- kouraku1+kouraku2+insi1 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,insi1,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2+insi1,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } #Ｂ hako4 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) insi2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+insi2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,insi2,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2+insi2,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } #Ｃ hako5 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+insi2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,bakuro,outcome)) logi <- glm(bakuro~kouraku1,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } x <- c() for (i in 1:1000) x <- c(x,hako1()) A <- c() for (i in 1:1000) A <- c(A,hako3()) B <- c() for (i in 1:1000) B <- c(B,hako4()) C <- c() for (i in 1:1000) C <- c(C,hako5()) boxplot(x,A,B,C,names=c("理想的","Ａ","Ｂ","Ｃ"))
やはり、Ａ，Ｂはいずれも理想的な場合よりも分散の値が大きくなっているし、Ｃには交絡が残っている。
なお、確認として、Ａのデータについてinsi1を除いて傾向スコアを計算してやると、（結果をＡ.とする）
hako3. <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) insi1 <- runif(100,0,1) x <- kouraku1+kouraku2+insi1 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,insi1,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } A. <- c() for (i in 1:1000) A. <- c(A.,hako3.()) boxplot(x,A,A.,names=c("理想的","Ａ","Ａ."))
このように、理想的なときと同じ結果が得られる。
まとめると、傾向スコアの算出に用いる因子には、
「曝露率とアウトカムの両方と関係している交絡因子だけを、できるだけ多く加える」ことが重要であることがシミュレーションからも明らかとなった。
逆に言えば、その部分さえ十分正確に出来れば、そこそこ安心して結果を受け入れることができそうに思える。
（回帰モデルは、入れるべき因子の項はそれでいいのか？とか、そのモデル式で本当に正しいのか？とか、気になる点がたくさんあるが…。）
もちろん厳密な値ではないかもしれないけど、かといって大きく外れてさえなければそれでいいかな、というとき（実際にそういうことはよくあるように思われる）に、このようなセミパラメトリックな方法は相当な威力を発揮すると感じた。
以上
10月8日（水）
検出力について
数理統計の教科書を読んでいたら、検定における「検出力」に関する内容があって、すぐにはピンと来なかったので、Ｒで検出力を求める関数をつくってみた。
（ここへ来てより道しすぎな感は否めないが…。）
正規母集団 N(mu,1)からn個の無作為標本をとってくる場合を考える。
いま、mu>mu0のとき、mu=mu0を帰無仮説として有意水準95％で片側検定を行って帰無仮説が棄却される確率を求める。
とりあえずモンテカルロ法が使えると思ったのでそのバージョン。
mu0 <- 0 ken.p <- function(x){ a <- qnorm(0.975,mu0,1/x[2]) cnt <- 0 for (i in 1:10000){ if (rnorm(1,x[1],1/x[2])>a) cnt <- cnt+1 else cnt <- cnt } return(cnt/10000) } mu <- seq(0,1,length=50) n <- 1:10 mu.n <- as.matrix(expand.grid(mu,n)) ken <- apply(mu.n,1,ken.p) library(rgl) plot3d(mu.n[,1],mu.n[,2],ken,xlab="mu",ylab="n")
mu,nが大きくなるほど検出力が大きくなることが確認できた。
そういえばモンテカルロ・シミュレーションを使ったのは初めてなので、感慨深い。
ただ、モンテカルロを使う必要はなくて、
ken.p. <- function(x){ a <- qnorm(0.975,mu0,1/x[2]) b <- pnorm(a,x[1],1/x[2],lower=F) return(b) }
としても可能。
とりあえず検出力の概念については理解できた。
要するに、標準正規分布の表しか与えられていなくて紙と鉛筆だけでやる場合（教科書のやり方）、式を同値変形して解かなくてはいかないのでイメージとしてつかみにくかったということらしい。
確率分布とかについて考える場合、直観的に処理できるという意味でもコンピュータは優れていると感じた。
以上
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になって、それなりに妥当な回帰が出来た。
来週の予定としては、まず「回帰の妥当性」についてもう一度考えた後、「傾向スコア分析」なるものについて勉強しようと思う。
それが終わったら論文をいくつか読んで、あとは最終レポート作成…という流れになりそうである。
以上
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
Ｒって便利だな、と改めて思った。
以上
10月1日（水）
ロジスティック回帰モデル
３Ｄプロットは僕の理解の範疇を超えているようである。
意味が分からなさすぎて発狂しそうになったので、諦めてperspでグラフをつくってみた。
age <- c(18,21,22,25,26,28,33,34,35,37,42,47,55,56,58,61,65,68,75,77) out <- c(0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,1,1,1,1,1) d <- cbind(age,out) age.n <- d[which(d[,2]==0),1] age.o <- d[which(d[,2]==1),1] p.n <- function(a,x0){ 1-1/(1+exp(-a*(age.n-x0))) } p.o <- function(a,x0){ 1/(1+exp(-a*(age.o-x0))) } p <- function(a,x0){ prod(p.n(a,x0))*prod(p.o(a,x0)) } a. <- seq(0,0.5,0.01) x0. <- 20:70 z <- c() for(i in a.) for(j in x0.) z <- c(z,p(i,j)) z. <- matrix(z,ncol=length(x0.),byrow=T) persp(a.,x0.,z.,theta = 30)
また、max(z.)を調べたところ、a=0.3,x0=34のときに最も尤度が高くなることが分かった。
実際に係数をこれにしてグラフを描くと、
a <- 0.3 x0 <- 34 f <- function(x){ 1/(1+exp(-a*(x-x0))) } plot(age,out,xlim=c(15,77),ylim=c(0,1),ylab="リスク",main="年齢とリスクの関係") par(new=T) curve(f,15,77,xlim=c(15,77),ylim=c(0,1),ann=F)
なるほど確かに、それらしい曲線が引けている。
なお、対数をとって微分する方法もやってみたのだが、うまく解が求められなかったので挫折した。
一応そういう方法もある、ということだけ頭に入れておこうと思う。
回帰モデルによる交絡の制御
回帰モデルの強みは、例えば曝露群と非曝露群とで年齢層が全く異なっていて、層化が使えないような場合でも統計量の推定が行えることである。
（ただし、モデルが合っている保証はないが）
曝露群、非曝露群で全く年齢層が異なっているデータを作って、曝露による効果の推定を行ってみた。
ここで、アウトカム＝ k0 + k1*（年齢）+ k2*（曝露の有無1or0）
（k0,k1,k2は定数）
というモデル式を設定した。
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) d <- data.frame(D) d.e <- d[sample(which(d[,1]<=50&d[,2]==1),n),] d.n <- d[sample(which(d[,1]>50&d[,2]==0),n),] plot(d.e[,1],d.e[,3],xlim=c(20,80),ylim=c(0,30)) par(new=T) plot(d.n[,1],d.n[,3],xlim=c(20,80),ylim=c(0,30)) result.e <- lm(d.e[,3]~d.e[,1],data=d) result.n <- lm(d.n[,3]~d.n[,1],data=d) E <- result.e$coefficients[1]-result.n$coefficients[1] E > E (Intercept) 1.736189
このように、Rに装備されているlm()の利用によって推定値を出すことは出来た。
ただ、この場合年齢による効果（k2）は曝露によらず一定とするモデルで考えているので、その条件を加えて推定を行った方がよいと思われる。
そこで、自分で最小二乗法を利用してやってみた。
（もしかしたらRに装備されている関数で出来るのかもしれないが）
k0 <- 0.2
k1 <- 0.3
k2 <- 3
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]-f1(X[1],X[3]))^2)
}
P(c(0.3,0.2,3.2))
op <- optim(c(0.3,0.2,3.2),P)
E <- op$par[3]-op$par[2]
E
> E
[1] 6.1551理論的には合ってるような気がするのだが、値がうまく合わない。
ちょっと理由が分からないので、明日また検討することにする。
以上
9月30日（火）
ロジスティックモデル
「ロスマンの疫学」p.303のデータにおいて、モデル式を
リスク＝
とし、最尤法による係数推定を行った。
age <- c(18,21,22,25,26,28,33,34,35,37,42,47,55,56,58,61,65,68,75,77) out <- c(0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,1,1,1,1,1) d <- cbind(age,out) plot(age,out) age.n <- deta[which(d[,2]==0),1] age.o <- deta[which(d[,2]==1),1] p.n <- function(a,x0){ 1-1/(1+exp(-a*(age.n-x0))) } p.o <- function(a,x0){ 1/(1+exp(-a*(age.o-x0))) } p <- function(a,x0){ prod(p.n(a,x0))*prod(p.o(a,x0)) } a. <- seq(0,0.5,0.01) x0. <- 20:70 z <- c() for(i in a.) for(j in x0.) z <- c(z,p(i,j)) z. <- matrix(z,ncol=length(x0.),byrow=T) plot3d(a.,x0.,z.)
しかし、なんだか変な図が出来た上に、ブログへのデータの貼り方がよくわからない。
明日ちゃんとまとめようと思う。
以上
9月29日（月）
今日で、「ロスマンの疫学」を最後まで読み切った。
とりあえず、Ｒでやってみたい題材には困らなさそうである。
ただ、統計の教科書で回帰分析について読んでから始めたいと思ったので、あとは統計の教科書を読んでいた。
徐々に残りの期間も少なくなってきているので、明日～明後日ぐらいには教科書を読み終えたいところである。
以上
