Nothing
#' @name zambia_graph
#' @title Prior precision matrix for spatial variable in Zambia data set
#' @description This is a 57x57 matrix containing row- and columwise the regions of Zambia, and the entries
#' define the neighbourhoodstructure. The corresponding map sambia.bnd can be downloaded from \url{http://www.stat.uni-muenchen.de/~kneib/regressionsbuch/daten_e.html}.
#' from the bnd file the prior precision matrix is obtained by
#' library(BayesX)
#' map <- read.bnd("zambia.bnd")
#' K <- bnd2gra(map)
#' @docType data
NULL
#' @name zambia_height92
#' @title Malnutrition in Zambia
#' @description The primary goal of a statistical analysis is to determine the effect of certain socioeconomic
#' variables of the child, the mother, and the household on the child`s nutritional
#' condition
#' \itemize{
#' \item zscore child`s Z-score
#' \item c _breastf duration of breastfeeding in months
#' \item c_age child`s age in months
#' \item m_agebirth mother`s age at birth in years
#' \item m_height mother`s height in centimeter
#' \item m_bmi mother`s body mass index
#' \item m_education mother`s level of education
#' \item m_work mother`s work status
#' \item region region of residence in Zambia
#' \item district district of residence in Zambia
#' }
#' @format A data frame with 4421 rows and 21 variables
#' @source \url{http://www.stat.uni-muenchen.de/~kneib/regressionsbuch/daten_e.html}
#' @docType data
NULL
#' Compute Density Function of Approximated (Differentiably) Uniform Distribution.
#'
#' @param x denotes the argument of the density function.
#' @param scale the scale parameter originally defining the upper bound of the uniform distribution.
#' @param tildec denotes the ratio between scale parameter \eqn{\theta} and \eqn{s}. The latter is responsible for the closeness of
#' the approximation to the uniform distribution. See also below for further details and the default value.
#' @return the density.
#' @author Nadja Klein
#' @details The density of the uniform distribution for \eqn{\tau} is approximated by
#' \deqn{p(\tau)=(1/(1+exp(\tau\tilde{c}/\theta-\tilde{c})))/(\theta(1+log(1+exp(-\tilde{c}))))}.
#' This results in \deqn{p(\tau^2)=0.5*(\tau^2)^(-1/2)(1/(1+exp((\tau^2)^(1/2)\tilde{c}/\theta-\tilde{c})))/(\theta(1+log(1+exp(-\tilde{c}))))} for \eqn{tau^2}.
#' \eqn{\tilde{c}} is chosen such that \eqn{P(\tau<=\theta)>=0.95}.
#'
#' @references Nadja Klein and Thomas Kneib (2015). Scale-Dependent Priors for Variance Parameters in Structured Additive Distributional Regression.
#' \emph{Working Paper}.
#'
#' @seealso \code{\link{rapprox_unif}},\code{\link{papprox_unif}}
#' @export
#'
dapprox_unif <- function(x,scale,tildec=13.86294)
{
arg <- x^(0.5)*tildec/scale - tildec
ret <- 0.5*x^(-0.5)*(1/(1+exp(arg)))/(scale*(1+log(1+exp(-tildec))/tildec))
return(ret)
}
NULL
#' Draw Random Numbers from Approximated (Differentiably) Uniform Distribution.
#'
#' @param n number of draws.
#' @param scale the scale parameter originally defining the upper bound of the uniform distribution.
#' @param tildec denotes the ratio between scale parameter \eqn{\theta} and \eqn{s}. The latter is responsible for the closeness of
#' the approximation to the uniform distribution. See also below for further details and the default value.
#' @param seed denotes the seed
#' @return n draws with density \code{\link{papprox_unif}}.
#' @author Nadja Klein
#' @details The method is based on the inversion method and the quantile function is computed numerically using \code{\link{uniroot}}.
#'
#' @references Nadja Klein and Thomas Kneib (2015). Scale-Dependent Priors for Variance Parameters in Structured Additive Distributional Regression.
#' \emph{Working Paper}.
#' @seealso \code{\link{rapprox_unif}},\code{\link{papprox_unif}}
#' @import stats
#' @export
#'
rapprox_unif<-function (n=100, scale,tildec=13.86294,seed=123)
{
set.seed(seed)
invF <- function(x,u)
{
arg <- x^(0.5)*tildec/scale
( tildec*x^(0.5) - scale*log(exp(arg)+exp(tildec))) / ( scale * (tildec + log(1+exp(-tildec)) ))+1-u
}
us <- runif(n)
vec <- rep(NA,length=n)
for(i in 1:length(us))
{
vec[i] <- uniroot(invF,lower=0,upper=10*scale^2,u=us[i])$root
}
return(vec)
}
NULL
#' Compute Cumulative Distribution Function of Approximated (Differentiably) Uniform Distribution.
#'
#' @param x denotes the argument of cumulative distribution function
#' @param scale the scale parameter originally defining the upper bound of the uniform distribution.
#' @param tildec denotes the ratio between scale parameter \eqn{\theta} and \eqn{s}. The latter is responsible for the closeness of
#' the approximation to the uniform distribution. See also below for further details and the default value.
#' @return the cumulative distribution function.
#' @author Nadja Klein
#' @details The cumulative distribution function of \code{\link{dapprox_unif}} is given by
#' \deqn{(1/(log(1+exp(-\tilde{c}))+\tilde{c}))*(\tilde{c}*(\tau^2)^(1/2)/\theta-log(exp((\tau^2)^(1/2)*\tilde{c}/\theta)+exp(\tilde{c})))}
#' \eqn{\tilde{c}} is chosen such that \eqn{P(\tau^2<=\theta)>=0.95}.
#'
#' @references Nadja Klein and Thomas Kneib (2015). Scale-Dependent Priors for Variance Parameters in Structured Additive Distributional Regression.
#' \emph{Working Paper}.
#'
#' @seealso \code{\link{rapprox_unif}},\code{\link{dapprox_unif}}
#' @export
#'
papprox_unif <- function(x,scale,tildec=13.86294)
{
arg <- x^(0.5)*tildec/scale
#ret <- (sqrt(x)*tildec - log(1+exp(arg-tildec)) + log(1+exp(-tildec))) / (scale*(1+log(exp(-tildec))+tildec))
ret <- ( tildec*x^(0.5)/scale - log(exp(arg)+exp(tildec))) / ((tildec + log(1+exp(-tildec)) ))+1
return(ret)
}
NULL
dapprox_unif_exp <- function(x,scale,tildec=13.86294,log=FALSE)
{
ex <- exp(x)
arg <- ex^(0.5)*tildec/scale - tildec
if(log)
ret <- 0.5*x-log(2)-log(scale)-log(1+log(1+exp(-tildec))/tildec)-log(1+exp(arg))
else
ret <- 0.5*sqrt(ex)*(1/(1+exp(arg)))/(scale*(1+log(1+exp(-tildec))/tildec))
return(ret)
}
NULL
#' Computing Designmatrix for Splines
#'
#' This function computes the design matrix for Bayesian P-splines as it would be done in
#' BayesX. The implementation currently on works properly for default values (knots=20, degree=3).
#'
#' @param x the covariate vector.
#' @param degree of the B-splines, default is 3.
#' @param m number of knots, default is 20.
#' @param min_x the left interval boundary, default is min(x).
#' @param max_x the right interval boundary, defalut is max(x).
#' @return a list with design matrix at distinct covariates, design matrix at all observations,
#' index of sorted observations, the difference matrix, precision matrix and the knots used.
#' @author Nadja Klein
#' @references Stefan Lang and Andy Brezger (2004). Bayesian P-Splines.
#' \emph{Journal of Computational and Graphical Statistics}, \bold{13}, 183--212.
#'
#' Belitz, C., Brezger, A., Klein, N., Kneib, T., Lang, S., Umlauf, N. (2015): BayesX - Software for Bayesian inference in structured additive regression models.
#' Version 3.0.1. Available from http://www.bayesx.org.
#' @import splines
#' @export
#'
DesignM <- function(x,degree=3,m=20,min_x=min(x),max_x=max(x)) {
##############################################################################################
#reduce establishment of design matrix to unique and sorted observations
x_work <- unique(sort.int(x,index.return=TRUE)$x)
min_x_work <- min_x
max_x_work <- max_x
#find indices to get back to whole covariate vector of length n(number of observations.)
ix <- c()
for(i in 1:length(x)) {
ix <- c(ix,which(x_work==x[i]))
}
##############################################################################################
#set degree of B-splines and number of inner knots
degree <- degree
l <- degree
hilfsgr <- (max_x_work-min_x_work)/100
step <- (max_x_work+hilfsgr-min_x_work+hilfsgr)/(m-1)
knoten_B=seq(min_x_work-hilfsgr-degree*step,max_x_work+hilfsgr+degree*step,by=step)
if(length(knoten_B)>1) {
Z_BayesX<- splineDesign(knoten_B,x_work,ord=degree+1,outer.ok=T) #designmatrix of x_work as in BayesX
} else {
knoten=seq(min_x_work-0.001-(degree)*((max_x_work+0.001-min_x_work+0.001)/(m-1)),max_x_work+0.001+(degree)*((max_x_work+0.002-min_x_work)/(m-1)),length=m+2*l)
Z_BayesX<- splineDesign(knoten,x_work,ord=degree+1,outer.ok=T) #designmatrix of x_work as in BayesX
}
Zcomp_BayesX <- Z_BayesX[ix,]
#Construction of matrix of second differences (beta_j-2*beta_{j-1}+beta_{j-2})
Dmatrix <- matrix(0,nrow=m+l-3,ncol=m+l-1)
for(i in 1:(m+l-3)){
Dmatrix[i,i] <- 1
Dmatrix[i,i+1] <- -2
Dmatrix[i,i+2] <- 1
}
#Kmatrix, product of the transposed Dmatrix and the Dmatrix (needed for penalizing term)
Kmatrix <- t(Dmatrix)%*%Dmatrix
return(list(Z_B=Z_BayesX,Zcomp_B=Zcomp_BayesX,ix=ix,Dmatrix=Dmatrix,Kmatrix=Kmatrix,knoten_B=knoten_B))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.