alt text

Introduction

In This Unit, I will show how to use Beyesian inference in a simple way. For more details you can see http://en.wikipedia.org/wiki/Bayesianinference.

A pharmaceutical company wants to test the efficiency of a drug in a group of patients affected by a certain ailment. To do this, it carries out a experiment, carried out on 15 patients, and note the results obtained, which can be of two types: positive (P) or negative (N). It is assumed that the probability that the test gives positive is p (the same for each patient). It is known that the number of positives observed is a random variable Y ⇠ Bin (n = 15, p). The laboratory experts, taking into takes into account the experience acquired in similar studies, they consider an a priori distribution p ⇠ Be (2, 4).

Suppose that, of the 15 patients in whom the study was done, in 6 of them positive results have been observed.

To carry out this exercise, I will use the R Bolstad library. The function binogcp within the Bolstad library allows us to evaluate and represent the density at posteriori of the probability of success in a Bernouilli experiment, p, with likelihood coming from a binomial and a continuous a priori distribution, which can be chosen specifying the density parameter, with possible values “beta”, “exp”, “normal”, or “uniform” (by default, if nothing is specified, it is considered an a priori distribution uniform - equivalent to a beta of parameters (1,1)). Depending on the distribution a priori chosen, the corresponding parameters must also be specified in the argument params: beta: a, b; exp: rate; normal: mean, sd; uniform: min, max. The values returned by the binogcp function are: + likelihood the scaled likelihood function of p, given x and n + posterior the posterior probability distribution of p, given x and n + pi the vector of possible values of p used in the prior distribution + pi.prior the associated mass probability function for the p values

1. I will Represent the a priori and a posteriori density of p.

library(Bolstad)
binogcp(6,15,density="beta", params = c(2,4),n.pi = 1e4)
results<-binogcp(6,15,density="beta", params = c(2,4),n.pi = 1e4)

mean(results$likelihood)
## [1] 1.248751e-05
head(results$posterior,5)
## [1] 7.868402e-25 1.718756e-21 6.132452e-20 6.456677e-19 3.745407e-18
mean(results$posterior)
## [1] 1
head(results$pi,10)
##  [1] 0.00005 0.00015 0.00025 0.00035 0.00045 0.00055 0.00065 0.00075 0.00085
## [10] 0.00095
mean(results$pi)
## [1] 0.5
mean(results$pi.prior)
## [1] 1

A priori Be (2, 4) and a posteriori Be (8, 13) distribution.

2. I will Calculate the posterior mean and posterior standard deviation.

To calculate the posterior mean and posterior standard deviation, we use the sintegral function, also from the Bolstad library, which calculates numerical integrals using Simpson’s rule.

p1<-results$pi
p2<-results$posterior
dens<-p1*p2
post.mean<-sintegral(p1,dens)$value
post.mean
## [1] 0.3809524
round(post.mean,4)
## [1] 0.381

posterior standard deviation.

dens2<-(p1-post.mean)^2*p2
post.var<-sintegral(p1,dens2)$value
post.var
## [1] 0.01071944
round(post.var,4)
## [1] 0.0107
post.sd<-sqrt(post.var)
post.sd
## [1] 0.1035347
round(post.sd,4)
## [1] 0.1035
post.sd<-sqrt(post.var)
round(post.sd,4)
## [1] 0.1035

3. I will Calculate a 95% probability interval for p.

cdf<-sintegral(p1,p2,n.pts = length(p1))$cdf


#Low limit 
lcb<-qbeta(0.025,8,13)
lcb
## [1] 0.1911901
round(lcb,4)
## [1] 0.1912
#Upper limit
ucb<-qbeta(0.975,8,13)
round(ucb,4)
## [1] 0.5922

4. I will Contrast the hypothesis H0: p <= 0.4 against H1: p> 0.4.

p0=0.4
x0=p1[p1<=p0]
x1=p1[p1>p0]
P0<-sintegral(x0,p2[p1<=p0])$value
P1<-sintegral(x1,p2[p1>p0])$value
round(P0,4)
## [1] 0.5839
round(P1,4)
## [1] 0.4157