Table Of ContentChapman & Hall/CRC
Computer Science and Data Analysis Series
BAYESIAN REGRESSION
MODELING WITH INLA
XIAOFENG WANG
Cleveland Clinic, Cleveland, Ohio
YU RYAN YUE
Baruch College, The City University of New York
JULIAN J. FARAWAY
University of Bath, UK
CRC Press
Taylor & Francis Group
6000 Broken Sound Parkway NW, Suite 300
Boca Raton, FL 33487-2742
© 2018 by Taylor & Francis Group, LLC
CRC Press is an imprint of Taylor & Francis Group, an Informa business
Version Date: 20180111
International Standard Book Number-13: 978-1-4987-2725-9 (Hardback)
Visit the Taylor & Francis Web site at
http://www.taylorandfrancis.com
and the CRC Press Web site at
http://www.crcpress.com
Contents
1 Introduction 1
1.1 QuickStart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Hubble’sLaw . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 StandardAnalysis . . . . . . . . . . . . . . . . . . . . . . 2
1.1.3 BayesianAnalysis . . . . . . . . . . . . . . . . . . . . . . 3
1.1.4 INLA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 BayesTheory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 PriorandPosteriorDistributions . . . . . . . . . . . . . . . . . . . 9
1.4 ModelChecking . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.5 ModelSelection . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6 HypothesisTesting . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.7 BayesianComputation . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7.1 Exact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7.2 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.7.3 Approximation . . . . . . . . . . . . . . . . . . . . . . . . 17
2 TheoryofINLA 19
2.1 LatentGaussianModels(LGMs) . . . . . . . . . . . . . . . . . . 19
2.2 GaussianMarkovRandomFields(GMRFs) . . . . . . . . . . . . . 21
2.3 LaplaceApproximationandINLA . . . . . . . . . . . . . . . . . . 23
2.4 INLAProblems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3 BayesianLinearRegression 39
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 BayesianInferenceforLinearRegression . . . . . . . . . . . . . . 40
3.3 Prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4 ModelSelectionandChecking . . . . . . . . . . . . . . . . . . . . 49
3.4.1 ModelSelectionbyDIC . . . . . . . . . . . . . . . . . . . 49
3.4.2 PosteriorPredictiveModelChecking . . . . . . . . . . . . 50
3.4.3 Cross-ValidationModelChecking . . . . . . . . . . . . . . 52
3.4.4 BayesianResidualAnalysis . . . . . . . . . . . . . . . . . 54
3.5 RobustRegression . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.6 AnalysisofVariance . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.7 RidgeRegressionforMulticollinearity . . . . . . . . . . . . . . . 59
3.8 RegressionwithAutoregressiveErrors . . . . . . . . . . . . . . . 63
4 GeneralizedLinearModels 71
4.1 GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2 BinaryResponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.3 CountResponses . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.1 PoissonRegression . . . . . . . . . . . . . . . . . . . . . . 77
4.3.2 NegativeBinomialRegression . . . . . . . . . . . . . . . . 79
4.4 ModelingRates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5 GammaRegressionforSkewedData . . . . . . . . . . . . . . . . 87
4.6 ProportionalResponses . . . . . . . . . . . . . . . . . . . . . . . 91
4.7 ModelingZero-InflatedData . . . . . . . . . . . . . . . . . . . . . 96
5 LinearMixedandGeneralizedLinearMixedModels 101
5.1 LinearMixedModels . . . . . . . . . . . . . . . . . . . . . . . . 101
5.2 SingleRandomEffect . . . . . . . . . . . . . . . . . . . . . . . . 102
5.2.1 ChoiceofPriors . . . . . . . . . . . . . . . . . . . . . . . 106
5.2.2 RandomEffects . . . . . . . . . . . . . . . . . . . . . . . . 109
5.3 LongitudinalData . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.3.1 RandomIntercept . . . . . . . . . . . . . . . . . . . . . . . 112
5.3.2 RandomSlopeandIntercept . . . . . . . . . . . . . . . . . 113
5.3.3 Prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.4 ClassicalZ-MatrixModel . . . . . . . . . . . . . . . . . . . . . . 119
5.4.1 RidgeRegressionRevisited . . . . . . . . . . . . . . . . . 121
5.5 GeneralizedLinearMixedModels . . . . . . . . . . . . . . . . . . 124
5.6 PoissonGLMM . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.7 BinaryGLMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
5.7.1 ImprovingtheApproximation . . . . . . . . . . . . . . . . 139
6 SurvivalAnalysis 141
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
6.2 SemiparametricModels . . . . . . . . . . . . . . . . . . . . . . . 143
6.2.1 PiecewiseConstantBaselineHazardModels . . . . . . . . 143
6.2.2 StratifiedProportionalHazardsModels . . . . . . . . . . . 146
6.3 AcceleratedFailureTimeModels . . . . . . . . . . . . . . . . . . 148
6.4 ModelDiagnosis . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
6.5 IntervalCensoredData . . . . . . . . . . . . . . . . . . . . . . . . 157
6.6 FrailtyModels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
6.7 JointModelingofLongitudinalandTime-to-EventData . . . . . . 164
7 RandomWalkModelsforSmoothingMethods 169
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
7.2 SmoothingSplines . . . . . . . . . . . . . . . . . . . . . . . . . . 170
7.2.1 RandomWalk(RW)PriorsforEqually-SpacedLocations . 170
7.2.2 ChoiceofPriorsonσ2andσ2 . . . . . . . . . . . . . . . . 176
ε f
7.2.3 RandomWalkModelsforNon-EquallySpacedLocations . 179
7.3 Thin-PlateSplines . . . . . . . . . . . . . . . . . . . . . . . . . . 185
7.3.1 Thin-PlateSplinesonRegularLattices. . . . . . . . . . . . 185
7.3.2 Thin-PlateSplinesatIrregularly-SpacedLocations . . . . . 188
7.4 BesagSpatialModel . . . . . . . . . . . . . . . . . . . . . . . . . 192
7.5 PenalizedRegressionSplines(P-Splines) . . . . . . . . . . . . . . 195
7.6 AdaptiveSplineSmoothing . . . . . . . . . . . . . . . . . . . . . 198
7.7 GeneralizedNonparametricRegressionModels . . . . . . . . . . . 201
7.8 ExcursionSetwithUncertainty . . . . . . . . . . . . . . . . . . . 206
8 GaussianProcessRegression 211
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
8.2 PenalizedComplexityPriors . . . . . . . . . . . . . . . . . . . . . 216
8.3 CredibleBandsforSmoothness . . . . . . . . . . . . . . . . . . . 217
8.4 Non-StationaryFields . . . . . . . . . . . . . . . . . . . . . . . . 220
8.5 InterpolationwithUncertainty . . . . . . . . . . . . . . . . . . . . 222
8.6 SurvivalResponse . . . . . . . . . . . . . . . . . . . . . . . . . . 226
9 AdditiveandGeneralizedAdditiveModels 229
9.1 AdditiveModels . . . . . . . . . . . . . . . . . . . . . . . . . . . 229
9.2 GeneralizedAdditiveModels . . . . . . . . . . . . . . . . . . . . 236
9.2.1 BinaryResponse . . . . . . . . . . . . . . . . . . . . . . . 237
9.2.2 CountResponse . . . . . . . . . . . . . . . . . . . . . . . 240
9.3 GeneralizedAdditiveMixedModels . . . . . . . . . . . . . . . . . 245
10 Errors-in-VariablesRegression 251
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251
10.2 ClassicalErrors-in-VariablesModels . . . . . . . . . . . . . . . . 254
10.2.1 A Simple Linear Model with Heteroscedastic Errors-in-
Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . 254
10.2.2 AGeneralExposureModelwithReplicatedMeasurements . 257
10.3 BerksonErrors-in-VariablesModels . . . . . . . . . . . . . . . . . 263
11 MiscellaneousTopicsinINLA 267
11.1 SplinesasaMixedModel . . . . . . . . . . . . . . . . . . . . . . 267
11.1.1 TruncatedPowerBasisSplines . . . . . . . . . . . . . . . . 267
11.1.2 O’SullivanSplines . . . . . . . . . . . . . . . . . . . . . . 268
11.1.3 Example:CanadianIncomeData. . . . . . . . . . . . . . . 269
11.2 AnalysisofVarianceforFunctionalData . . . . . . . . . . . . . . 272
11.3 ExtremeValues . . . . . . . . . . . . . . . . . . . . . . . . . . . . 278
11.4 DensityEstimationUsingINLA . . . . . . . . . . . . . . . . . . . 283
AppendixA Installation 289
AppendixB UninformativePriorsinLinearRegression 291
Bibliography 297
Index 309
Preface
INLAstandsforIntegratedNestedLaplaceApproximations.Itisamethodforfitting
abroadclassofBayesianmodels.Historically,itwasdifficulttofitanythingbutthe
mostsimpleBayesianmodels.Overthelasttwentyyears,aclassofBayesiancom-
putationalmethodsbasedonasimulationmethodcalledMarkovchainMonteCarlo
(MCMC) has been developed and has seen wide acceptance in statistics. Popular
packages using these methods include BUGS, JAGS and STAN. Despite impres-
sive improvements, these packages suffer from two problems. First, they are slow.
Forsomemorecomplexand/orlargerdataproblems,MCMCcanbeinfeasiblyslow,
evenifyouarepreparedtowaitdays.Butevenformoremodestproblems,theability
tofitmodelsquicklyiscrucialtoexploratorydataanalyses.INLAisanapproxima-
tion method. Typically, the approximations are more than adequate and remember
that simulation methods are inevitably approximations also if they are to finish in
finite time. INLA takes no more than a few seconds to fit the models found in this
book.
The other practical difficulty with MCMC methods is they take substantial ex-
pertise to use. You need to learn a specialized programming language to specify
themodels.Youalsoneedtounderstandthediagnosticsthatdeterminewhetherthe
modelhasbeenfitcorrectlyorwhetherthesimulationprocesshasfailed.Someskill
isnecessarytomakeasuccessofthis.Despiteongoingimprovementinthisarea,this
hasbeenanobstacletowideradoptionofthesemethodsinthescientificcommunity.
There are some drawbacks to INLA too. Although you do not need to learn a
specific programming language as the models can be specified and analyzed using
R,itisstillnotthatstraightforward.Indeed,thatiswhyyoushouldreadthisbook.
Furthermore, INLA only applies to a class called latent Gaussian models. If you
browsethroughthetableofcontentsforthisbook,youwillseethatthisclassisvery
broad.Nevertheless,therewillbesomemodelsthatcannotbefitwithINLAbutcan
bedonewithMCMC.
Wemakesomeassumptionsaboutyou,thereader.Weexpectthatyouhaveaba-
sicknowledgeofstatisticaltheoryandpractice.Weexpectyoualreadyknowsome-
thing about Bayesian methodology. This is not a theoretical book as we focus our
presentationaroundexamples.Wehopethatscientists,whoalreadyusesomestatis-
tics,willfindthisbookaccessible.WealsoexpectyouhavesomeknowledgeofR.
Wedoprovidefullyworkingcodeforalltheexamplessoyoumaybeabletoadapt
this to your own needs without proficiency in R. Even so, you will need some ex-
periencewithRtodrawthefullbenefit.Ifyouneedmoreintroductorymaterial,we
directyoutomorecomprehensiveandaccessibletexts.
This book is about regression models and we have not presented the extensive
spatialdataanalysiscapabilitiesofINLA.Forthosespecificallyinterestedinspatial
data,wedirectyoutoBlangiardoandCameletti(2015).
We have gathered the data and additional functions we use in this text as an R
packagewhichyoumayfindcurrentlyat:
https://github.com/julianfaraway/brinla
Our first thanks go to Håvard Rue and his coworkers for developing the theory
andproducingthesoftwareforINLA.ThanksalsotoFinnLindgren,DanielSimpson
andEgilFerkingstadforhelpfuladvice.
1
Introduction
INLA can be quite complex so we warm up with a simple regression example to
illustrate the most basic features. We do not pretend this as a comprehensive intro-
ductiontoBayesianmethods.InsteadwereviewtheBayesiantheorythatwillarise
laterinthebook.WealsodiscusswhyyoumightwanttouseINLAinpreferenceto
thealternatives.
1.1 QuickStart
1.1.1 Hubble’sLaw
Let’sgetstartedwithINLAusingasimpleexample.Supposewewanttoknowthe
ageoftheuniverse.Inthebeginning,therewasthebigbang.Galaxiesmovedaway
fromthecenteroftheuniverseataconstantvelocity.Ifwemeasurethevelocityof
other galaxies relative to us and their distance from us, we can estimate the age of
theuniverse.WeuseHubble’slawwhichstatesthat:
y=βx
whereyistherelativevelocitybetweenusandanothergalaxyandxisthedistance
fromustotheothergalaxy.βiscalledHubble’sconstant.Wecanestimatetheage
oftheuniverseasβ−1.
The Hubble Space Telescope has allowed the collection of the data vital to an-
sweringthisquestion.Weusedataon24galaxiesreportedinFreedmanetal.(2001).
ThesamedataisalsoanalyzedinWood(2006).NowrefertoAppendixAforhow
to install INLA and the R package for this book. The data can be found in our R
package,whichiscalledbrinla.
data(hubble, package = "brinla")
WehighlighttheRcodethroughoutthisbookinthemannerabove.Youshouldtype
this in at the R console or you may find it easier to copy and paste the code from
R scripts for each chapter. See the appendix for how to obtain these scripts. Al-
ternatively, you can gain access to the data and functions in our package with the
command:library(brinla).Wehaveusedthedatacommandinsteadsothatthe
sourceofthedataisclear.
WeplotthedataasseeninFigure1.1.
2 BayesianRegressionModelingwithINLA
plot(y ~ x, xlab = "Distance(Mpc)", ylab = "Velocity(km/s)",
data = hubble)
•
•
0
150 ••• • • ••
s) •
m/0
city(k100 • •
Velo • • •
• • • •
0
0
5
•
•
•
•
5 10 15 20
Distance(Mpc)
FIGURE1.1
Distanceandvelocityof24galaxiesrelativetoourlocation.
We see that the observations do not all lie on a line and so there must be some
measurement error. This error is hardly surprising given the difficulty of observing
suchdistantobjects.Evenso,wedoseeanapproximatelylinearrelationshipbetween
the variables, consistent with Hubble’s law. We allow for the error by specifying a
model:
y =βx +ε, i=1,...,24.
i i i
Fornow,weassumeonlythattheerrorshavemeanzeroandvarianceσ2.
1.1.2 StandardAnalysis
Themethodofleastsquaresisnowalmost200yearsoldandcanbeappliedtoesti-
mateβ.Wedonothaveaninterceptterminthismodelsoweputa-1inthemodel
formula.
lmod <- lm(y ~ x - 1, data = hubble)
coef(lmod)
x
76.581
Thisisourestimateofβ,whichisHubble’sconstant.Weneedtoconverttheunitsof
the data to something more convenient. We have 60 seconds, 60 minutes, 24 hours
Introduction 3
and 365.25 days in a year and one megaparsec (Mpc) is 3.09×1019 km. Here is a
function to transform Hubble’s constant into the age of the universe in billions of
years.Weapplyittoourestimate:
hubtoage <- function(x) 3.09e+19/(x * 60^2 * 24 * 365.25 * 1e+09)
hubtoage(coef(lmod))
x
12.786
Our point estimate of the age of the universe is 12.8 billion years old. We know
thereissomeuncertaintyinthisestimateandsowemaywantaconfidenceinterval
to express this. Now we need to make a distributional assumption on the errors.
WeassumetheyfollowaGaussiandistribution.Therearesomeotherassumptions.
We assume the structural form of the model is correct. We take this from Hubble’s
law and the data appear to support this. We have also assumed that the errors have
constantvariance.Thismaynotbequitetruebutweshallletthispassintheinterest
ofsimplicity.
Weconstructa95%confidenceintervalforHubble’sconstant:
(bci <- confint(lmod))
2.5 % 97.5 %
x 68.379 84.783
Theintervalisquitewide.Wecanconvertthistoanintervalfortheageoftheuni-
verse:
hubtoage(bci)
2.5 % 97.5 %
x 14.32 11.549
Theinversionneededtocalculatetheagechangestheorderofthelimitstotheinter-
val.Weseethatthe95%confidenceintervalrunsfrom11.5to14.3billionyears.
1.1.3 BayesianAnalysis
Wemaybesatisfiedwiththepreviousanalysis.Wehaveapointestimateandwehave
anexpressionoftheuncertaintyinthisestimate.Yetifweconsidertheresultsmore
closely, we may feel some dissatisfaction and be willing to consider an alternative
methodofanalysis.
Itistemptingtoviewtheconfidenceintervalasclaiminga95%probabilitythat
theageoftheuniverseliesintheinterval11.5to14.3billionyears.Butaccordingto
thetheoryusedtomakeconfidenceintervalsthisisnotatallwhatisclaimed.This
theory,sometimescalledFrequentistandsometimescalledclassicalinference,views
parametersasfixed,butusuallyunknown,quantities.Soitmakesnosensetoassign
probabilitiestoparametersaccordingtothistheory.The95%confidenceappliesto
theintervalandnottheparameter.Itclaimsa95%chancethattheintervalcoversthe
truevalueoftheparameter.
Wemightobjectthatthedatahavenowbeenobservedandareknownsothein-
terval is no longer random. To answer this objection, we have to construct a rather
