R/predictive.checks.DDPstar.R

In DDPstar: Density Regression via Dirichlet Process Mixtures of Normal Structured Additive Regression Models

predictive.checks.DDPstar <-
function(object, statistics = c("min","max","kurtosis","skewness"), devnew = TRUE, ndensity = 512, trim = 0.025) {
	if(class(object)[1] != "DDPstar") {
		stop(paste0("This function can not be used for this object class: ", class(object)[object]))
	}

	y <- object$data_model$y
	n <- length(y)
	nsim <- nrow(object$fit$probs)
	yrep <- matrix(0, nrow = n, ncol = nsim)
	
	aux <- t(apply(object$fit$probs[1:nsim,,drop = FALSE], 1, function(x, n)  sample(1:length(x), n, replace = TRUE, prob = x), n = n))
	for(l in 1:nsim) {
		yrep[,l] <- rnorm(n = n, mean = colSums(t(object$data_model$X)*t(object$fit$beta[l,aux[l,],])), sd = object$fit$sd[l,aux[l,]])
	}

	i <- 1
	for(stat in statistics) {
		if(i != 1 & devnew) dev.new()
		yrepstat <- apply(yrep, 2, function(y, stat) {do.call(stat, list(y))}, stat = stat)
		yrepstat <- trim.out(yrepstat, trim = trim)
		ystat <- do.call(stat, list(y))
		xlim <- range(c(yrepstat,ystat))
		hist(yrepstat, col = "gray60", main = stat, xlim = xlim, xlab = "Statistic")
		abline(v = ystat, col = "red", lwd = 3)
		i <- i + 1
	}
	# Density
	if(devnew) dev.new()
	ylim <- c(0, max(density(y)$y) + 0.2*max(density(y)$y))
	xlim <- c(min(density(y)$x) - 0.2, max(density(y)$x) - 0.2)
	plot(density(yrep[,1], n = ndensity),col = "lightskyblue1", ylim = ylim, xlim = xlim, main = "Density", xlab = "Response variable")
	# Only a sample
	s <- sample(1:nsim, ifelse(nsim < 500, nsim, 500))
	#s <- 1:nsim 
	for(i in s){
		lines(density(yrep[,i]), col="lightskyblue1")
	}
	lines(density(y), col = "black", lwd = 4)

	res <- list()
	res$yrep <- yrep
	res$y <- y
	invisible(res)
}