mcmc: Obtain posterior samples for all model parameters

Description

This is the function to actually run the MCMC machinery to produce posterior samples from all model parameters and required derived values. It is a generic function, so that customized versions may be conveniently defined for specific subclasses of GeneralData, GeneralModel, and McmcOptions input.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
mcmc(data, model, options, ...)

## S4 method for signature 'GeneralData,GeneralModel,McmcOptions'
mcmc(data, model, options,
  program = c("JAGS", "OpenBUGS", "WinBUGS"), verbose = FALSE,
  fromPrior = data@nObs == 0L, ...)

## S4 method for signature 'Data,LogisticNormal,McmcOptions'
mcmc(data, model, options,
  verbose = FALSE, ...)

## S4 method for signature 'Data,LogisticIndepBeta,McmcOptions'
mcmc(data, model, options, ...)

## S4 method for signature 'DataDual,Effloglog,McmcOptions'
mcmc(data, model, options, ...)

## S4 method for signature 'DataDual,EffFlexi,McmcOptions'
mcmc(data, model, options, ...)

Arguments

data

The data input, an object of class GeneralData

model

The model input, an object of class GeneralModel

options

MCMC options, an object of class McmcOptions

program

the program which shall be used: either “JAGS” (default), “OpenBUGS” or “WinBUGS”

verbose

shall progress bar and messages be printed? (not default)

fromPrior

sample from the prior only? Defaults to checking if nObs is 0. For some models it might be necessary to specify it manually here though.

...

unused

Details

Reproducible samples can be obtained by setting the seed via set.seed before in the user code as usual. However, note that because the RNG sampler used is external to R, running this MCMC function will not change the seed position – that is, the repeated call to this function will then result in exactly the same output.

Value

The posterior samples, an object of class Samples.

Methods (by class)

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
# create some data from the class 'Data'
myData <- Data(x=c(0.1,0.5,1.5,3,6,10,10,10),
               y=c(0,0,0,0,0,0,1,0),
               doseGrid=c(0.1,0.5,1.5,3,6,
                          seq(from=10,to=80,by=2)))

# Initialize the CRM model 
model <- LogisticLogNormal(mean=c(-0.85, 1),
                           cov=
                             matrix(c(1, -0.5, -0.5, 1),
                                    nrow=2),
                           refDose=56)


# Sample from the posterior distribution
options <- McmcOptions(burnin=100,
                       step=2,
                       samples=1000)

samples <- mcmc(data = myData, model = model, options=options)



# create some data from the class 'Data'
myData <- Data(x=c(0.1,0.5,1.5,3,6,10,10,10),
               y=c(0,0,0,0,0,0,1,0),
               doseGrid=c(0.1,0.5,1.5,3,6,
                          seq(from=10,to=80,by=2)))

# Initialize the CRM model 
model <- LogisticLogNormal(mean=c(-0.85, 1),
                           cov=
                             matrix(c(1, -0.5, -0.5, 1),
                                    nrow=2),
                           refDose=56)


# Sample from the posterior distribution
options <- McmcOptions(burnin=100,
                       step=2,
                       samples=1000)

samples <- mcmc(data = myData, model = model, options=options)


##obtain mcmc DLE samples given the data, LogisticIndepBeta (DLE model) and mcmc simulations options
## data must be of 'Data' class
data<-Data(x=c(25,50,50,75,100,100,225,300),y=c(0,0,0,0,1,1,1,1),
           doseGrid=seq(25,300,25))
## model must be of 'LogisticIndepBeta' class
model<-LogisticIndepBeta(binDLE=c(1.05,1.8),DLEweights=c(3,3),DLEdose=c(25,300),data=data)
## options must be ''McmcOptions' class
options<-McmcOptions(burnin=100,step=2,samples=200)
set.seed(94)
samples<-mcmc(data=data,model=model,options=options)
##obtain mcmc efficacy samples given the data, 'Effloglog' model (efficacy model) and
## mcmc simulations options data must be of 'DataDual' class
data<-DataDual(x=c(25,50,25,50,75,300,250,150),
              y=c(0,0,0,0,0,1,1,0),
              w=c(0.31,0.42,0.59,0.45,0.6,0.7,0.6,0.52),
              doseGrid=seq(25,300,25))
## model must be of 'Effloglog' class
Effmodel<-Effloglog(Eff=c(1.223,2.513),Effdose=c(25,300),nu=c(a=1,b=0.025),data=data)

## options must be ''McmcOptions' class
options<-McmcOptions(burnin=100,step=2,samples=200)
set.seed(94)
samples<-mcmc(data=data,model=Effmodel,options=options)
##obtain mcmc efficacy samples given the data, 'EffFlexi' model (efficacy model) and 
## mcmc simulations options
## data must be of 'DataDual' class
data<-DataDual(x=c(25,50,25,50,75,300,250,150),
               y=c(0,0,0,0,0,1,1,0),
               w=c(0.31,0.42,0.59,0.45,0.6,0.7,0.6,0.52),
               doseGrid=seq(25,300,25))
## model must be of 'EffFlexi' class

Effmodel<- EffFlexi(Eff=c(1.223, 2.513),Effdose=c(25,300),
                    sigma2=c(a=0.1,b=0.1),sigma2betaW=c(a=20,b=50),smooth="RW2",data=data)

## options must be ''McmcOptions' class
options<-McmcOptions(burnin=100,step=2,samples=200)
set.seed(94)
samples<-mcmc(data=data,model=Effmodel,options=options)

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.