R/DDPstar.f.R

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

f <-
function(..., by = NULL, nseg = 5, bdeg = 3, pord = 2, atau = 1, btau = 0.005, prior.2D = c("anisotropic", "isotropic")) {
	prior.2D <- match.arg(prior.2D)
	vars <- as.list(substitute(list(...)))[-1] # List
	args <- match.call()

	d <- length(vars)
	if(d == 0 | d > 2) {
		stop("Incorrect number of covariates")
	}

	if(d > 1 & !is.null(args$by)) {
		stop("by variables are not allowed with two-dimensional P-splines")
	}
	if(d > 1 & any(pord != 2)) {
		stop("For two-dimensional P-splines only second order differences are allowed")
	}

	if (length(nseg) < d) nseg = rep(nseg, d)
	if (length(pord) < d) pord = rep(pord, d)
	if (length(bdeg) < d) bdeg = rep(bdeg, d)

	term <- deparse(vars[[1]], backtick = TRUE, width.cutoff = 500)
	if (term[1] == ".") {
		stop("f(.) not yet supported.")
	}	
	if (d > 1) { 
		for (i in 2:d) {
			term[i] <- deparse(vars[[i]], backtick = TRUE, width.cutoff = 500)
			if (term[i] == ".") { 
				stop("f(.) not yet supported.")
			}
		}
	}
	for (i in 1:d){
		term[i] <- attr(terms(reformulate(term[i])), "term.labels")
	} 
	nseg.new <- round(nseg)
	if (all.equal(nseg.new,nseg) != TRUE) {
		warning("argument nseg of f() should be integer and has been rounded")
	}
	nseg <- nseg.new
	pord.new <-round(pord)
	if (all.equal(pord.new,pord) != TRUE) {
		warning("argument pord of f() should be integer and has been rounded")
	}
	pord <- pord.new
	
	if (length(unique(term)) != d) { 
		stop("Repeated variables as arguments of a smooth are not permitted")
	}

	if(d == 1 & is.null(args$by)) { # 1D Smooth effect			   
		cov <- c("-1", term[1])
		vars <- term[1]
	} else if (d == 1 & !is.null(args$by)) { # Factor by curve 	  
		cov <- c(deparse(args$by, backtick = TRUE, width.cutoff = 500), term[1])
		vars <- cov
	} else if (d == 2) { # 2D PS-ANOVA	  
		cov <- c(term[1], term[2])
		vars <- cov
	} else {
		stop("Invalid expression")
	}
	res <- list(cov = cov, vars = vars, term = term, by.var = !is.null(args$by), nseg = nseg, bdeg = bdeg, pord = pord, 
		atau = atau, btau = btau, prior.2D = prior.2D)
	res

}