library(knitr) opts_chunk$set(cache=FALSE,comment=NA, fig.path="figs/", warning=FALSE, message=FALSE, pngquant='--speed=1 --quality=0-50') options(digits=5,show.signif.stars=FALSE,width=120) knit_hooks$set(pngquant = hook_pngquant) knitr::knit_hooks$set(setPch = function(before, options, envir) { if(before) par(pch = 20) }) opts_chunk$set(setPch = TRUE)
See the introduction for an overview.
Load the libraries:
library(ggplot2) library(INLA) library(dplyr)
Load in and plot the data:
data(psid, package="faraway") head(psid) summary(psid) psid$cyear <- psid$year-78 psid20 <- filter(psid, person <= 20) ggplot(psid20, aes(x=year, y=income))+geom_line()+facet_wrap(~ person) ggplot(psid20, aes(x=year, y=income+100, group=person)) +geom_line()+facet_wrap(~ sex)+scale_y_log10()
psid$slperson <- psid$person formula <- log(income) ~ cyear*sex+age+educ + f(person, model="iid") + f(slperson, cyear , model="iid") result <- inla(formula, family="gaussian", data=psid) result <- inla.hyperpar(result) summary(result)
Seems like the default priors might be adequate. Compute the transforms to an SD scale for the random effect terms. Make a table of summary statistics for the posteriors:
sigmaint <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[2]]) sigmaslope <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[3]]) sigmaepsilon <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[1]]) restab=sapply(result$marginals.fixed, function(x) inla.zmarginal(x,silent=TRUE)) restab=cbind(restab, inla.zmarginal(sigmaint,silent=TRUE)) restab=cbind(restab, inla.zmarginal(sigmaslope,silent=TRUE)) restab=cbind(restab, inla.zmarginal(sigmaepsilon,silent=TRUE)) mm <- model.matrix(~ cyear*sex+age+educ,psid) colnames(restab) = c(colnames(mm),"intSD","slopeSD","epsilon") data.frame(restab)
Also construct a plot the SD posteriors:
ddf <- data.frame(rbind(sigmaint,sigmaslope,sigmaepsilon),errterm=gl(3,nrow(sigmaint),labels = c("slope","intercept","epsilon"))) ggplot(ddf, aes(x,y, linetype=errterm))+geom_line()+xlab("log(income)")+ylab("density")
Posteriors look OK. Slope SD is on a different scale so not so comparable.
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.