R/bbase.psanova.bs.R

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

bbase.psanova.bs <-
function(x1, x2, ndx, bdeg = c(3,3), pord = c(2,2), eps = 1e-5, prior.2D) {
	# Smooth part
	terms <- list()
	psanova.penalty.part <- list()

	# x1
	MM1 <- bbase.bs(x1, ndx[1], bdeg = bdeg[1], pord = 2, eps = eps, intercept = FALSE)
	terms[[1]] <- MM1
	# x2
	MM2 <- bbase.bs(x2, ndx[2], bdeg = bdeg[2], pord = 2, eps = eps, intercept = FALSE)
	terms[[2]] <- MM2

	X1 <- MM1$X; Z1 <- MM1$Z; d1 <- MM1$d; c1 <- ncol(X1) + ncol(Z1)
	X2 <- MM2$X; Z2 <- MM2$Z; d2 <- MM2$d; c2 <- ncol(X2) + ncol(Z2)

	one1. <- X1[,1, drop = FALSE]
	one2. <- X2[,1, drop = FALSE]
		
	x1. <- X1[,-1, drop = FALSE]
	x2. <- X2[,-1, drop = FALSE]
		
	# Fixed and random matrices
	X <- Rten2(X2, X1)
	# Delete the intercept
	X <- X[,-1,drop = FALSE]

	Z <- cbind(Rten2(one2., Z1), Rten2(Z2, one1.), Rten2(x2., Z1), Rten2(Z2, x1.), Rten2(Z2, Z1))
		
	# Variance/Covariance components
	if(prior.2D == "anisotropic") {
		psanova.penalty.part[[1]] <- diag(x = rep(0L, ncol(X)), ncol = ncol(X), nrow = ncol(X))		
		psanova.penalty.part[[2]] <- diag(d1)
		psanova.penalty.part[[3]] <- diag(d2)
		psanova.penalty.part[[4]] <- diag(rep(1, pord[2] - 1)%x%d1)
		psanova.penalty.part[[5]] <- diag(d2%x%rep(1, pord[1] - 1))
		psanova.penalty.part[[6]] <- diag(rep(1, c2 - pord[2])%x%d1 + d2%x%rep(1, c1 - pord[1]))
	} else {
		psanova.penalty.part[[1]] <- diag(x = rep(0L, ncol(X)), ncol = ncol(X), nrow = ncol(X))
		
		# All smooth component together!
		grand.K <- list()
		grand.K[[1]] <- diag(d1)
		grand.K[[2]] <- diag(d2)
		grand.K[[3]] <- diag(rep(1, pord[2] - 1)%x%d1)
		grand.K[[4]] <- diag(d2%x%rep(1, pord[1] - 1))
		grand.K[[5]] <- diag(rep(1, c2 - pord[2])%x%d1 + d2%x%rep(1, c1 - pord[1]))
		psanova.penalty.part[[2]] <- as.matrix(Matrix::bdiag(grand.K))
	}

	B <- cbind(X,Z)
	# Center the matrix
	# cm <- colMeans(B)
	cm <- rep(0, ncol(B))
	B <- sweep(B, 2, cm)

	res <- list()
	res$B <- B
	res$K <- psanova.penalty.part
	res$terms <- terms
	attr(res$terms,"cm") <- cm
	res
}