原文链接：http://tecdat.cn/?p=6761

在逻辑回归中，我们将二元响应\（Y_i \）回归到协变量\（X_i \）上。下面的代码使用Metropolis采样来探索\（\ beta_1 \）和\（\ beta_2 \）的后验YiYi到协变量XiXi。

定义expit和分对数链接函数

1. logit<-function(x){log(x/(1-x))} 此函数计算\（（\ beta_1，\ beta_2）\）的联合后验。它返回后验的对数以获得数值稳定性。(β1,β2)(β1,β2)。它返回后验的对数以获得数值稳定性。

2. log_post<-function(Y,X,beta){

3. prob1  <- expit(beta[1] + beta[2]*X)

4. like+prior}

这是MCMC的主要功能.can.sd是候选标准偏差。

1. Bayes.logistic<-function(y,X,

2. n.samples=10000,

3. can.sd=0.1){

4. 


5. keep.beta     <- matrix(0,n.samples,2)

6. keep.beta[1,] <- beta

7. 


8. acc   <- att <- rep(0,2)

9. 


10. for(i in 2:n.samples){

11. 


12. for(j in 1:2){

13. 


14. att[j] <- att[j] + 1

15. 


16. # Draw candidate:

17. 


18. canbeta    <- beta

19. canbeta[j] <- rnorm(1,beta[j],can.sd)

20. canlp      <- log_post(Y,X,canbeta)

21. 


22. # Compute acceptance ratio:

23. 


24. R <- exp(canlp-curlp)

25. U <- runif(1)

26. if(U<R){

27. acc[j] <- acc[j]+1

28. }

29. }

30. keep.beta[i,]<-beta

31. 


32. }

33. # Return the posterior samples of beta and

34. # the Metropolis acceptance rates

35. list(beta=keep.beta,acc.rate=acc/att)}

生成一些模拟数据

1. set.seed(2008)

2. n         <- 100

3. X         <- rnorm(n)

4. true.p    <- expit(true.beta[1]+true.beta[2]*X)

5. Y         <- rbinom(n,1,true.p)

拟合模型

1. burn      <- 10000

2. n.samples <- 50000

3. fit  <- Bayes.logistic(Y,X,n.samples=n.samples,can.sd=0.5)

4. tock <- proc.time()[3]

5. 


6. tock-tick

结果

 abline(true.beta[1],0,lwd=2,col=2)

 abline(true.beta[2],0,lwd=2,col=2)

1. hist(fit$beta[,1],main="Intercept",xlab=expression(beta[1]),breaks=50)

2. 


1. hist(fit$beta[,2],main="Slope",xlab=expression(beta[2]),breaks=50)

2. abline(v=true.beta[2],lwd=2,col=2)

1. print("Posterior mean/sd")

2. ## [1] "Posterior mean/sd"

3. print(round(apply(fit$beta[burn:n.samples,],2,mean),3))

4. ## [1] -0.076  0.798

5. print(round(apply(fit$beta[burn:n.samples,],2,sd),3))

6. ## [1] 0.214 0.268

7. print(summary(glm(Y~X,family="binomial")))

8. ##

9. ## Call:

10. ## glm(formula = Y ~ X, family = "binomial")

11. ##

12. ## Deviance Residuals:

13. ##     Min       1Q   Median       3Q      Max

14. ## -1.6990  -1.1039  -0.6138   1.0955   1.8275

15. ##

16. ## Coefficients:

17. ##             Estimate Std. Error z value Pr(>|z|)

18. ## (Intercept) -0.07393    0.21034  -0.352  0.72521

19. ## X            0.76807    0.26370   2.913  0.00358 **

20. ## ---

21. ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

22. ##

23. ## (Dispersion parameter for binomial family taken to be 1)

24. ##

25. ##     Null deviance: 138.47  on 99  degrees of freedom

26. ## Residual deviance: 128.57  on 98  degrees of freedom

27. ## AIC: 132.57

28. ##

29. ## Number of Fisher Scoring iterations: 4

非常感谢您阅读本文，有任何问题请在下面留言！