R/interpret.DDPstarformula.R

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

interpret.DDPstarformula <-
function(formula, data) {
	env <- environment(formula) 
	if(inherits(formula, "character"))		  
		formula <- as.formula(formula)

	tf <- terms.formula(formula, specials = c("f", "ns", "bs", "rae"))
    if(!is.null(attr(tf,"specials")$ns) | !is.null(attr(tf,"specials")$bs)) {
        stop("'ns' (natural splines) or 'bs' (B-splines) are not allowed in the formula. Please use 'f' instead.")
    }

	if (attr(tf, "response") > 0) {
		response <- as.character(attr(tf, "variables")[2])
	} else {
		stop("The formula should include the response variable (left hand side)")
	}

	# Number of tems
	terms <- attr(tf, "term.labels")
	nt <- length(terms)

	# Smooth terms
	ns <- sort(attr(tf,"specials")$f) - 1 # Response is in the formula
	
	# Random terms
	nre <- sort(attr(tf,"specials")$rae) - 1 # Response is in the formula
   
	II <- list()
	h  <- list()
	nseg <- list()
	bdeg <- list()
	pord <- list()
	atau <- list()
	btau <- list()
	by.var <- list()
	prior.2D <- list()
	partial <- vector()
	k <- 0
	vars.formula <- NULL
	if(nt) {
		for (i in 1:nt) {
			if (i %in% ns) { # Smooth components (1D and 2D)        
                k <- k + 1                 
                #st <- eval(parse(text = paste("DDPstar.",terms[i], sep = "")), envir = env)
                st <- eval(parse(text = terms[i]), envir = env)
                II[[k]] <- st$cov
                h[[k]] <- -1
                nseg[[k]] <- st$nseg
                bdeg[[k]] <- st$bdeg
                pord[[k]] <- st$pord
                atau[[k]] <- st$atau
                btau[[k]] <- st$btau
                by.var[[k]] <- st$by.var
                prior.2D[[k]] <- st$prior.2D
                partial[k] <- terms[i]
                vars.formula <- c(vars.formula, st$vars)
            } else if (i %in% nre) { # Random effects
            	k <- k + 1
            	#st <- eval(parse(text = paste("DDPstar.",terms[i], sep = "")), envir = env)                 
                st <- eval(parse(text = terms[i]), envir = env)
                II[[k]] <- st$cov
                h[[k]] <- -2
                nseg[[k]] <- 0
                bdeg[[k]] <- 0
                pord[[k]] <- 0
                atau[[k]] <- st$atau
                btau[[k]] <- st$btau
                by.var[[k]] <- FALSE
                prior.2D[[k]] <- NULL
                partial[k] <- terms[i]
                vars.formula <- c(vars.formula, st$vars)
        	} else {	# Parametric components
				k <- k + 1
				II[[k]]<- c(-1, terms[i])
				h[[k]] <- 0
				nseg[[k]] <- 0
				bdeg[[k]] <- 0
				pord[[k]] <- 0
				atau[[k]] <- 0
                btau[[k]] <- 0
				by.var[[k]] <- FALSE
				prior.2D[[k]] <- NULL
				partial[k] <- terms[i]
				vars.formula <- c(vars.formula, all.vars(formula(paste("~", terms[i]))))
			}
        }
        #names.cov <- all.vars(formula)[-1]
        names.cov <- vars.formula 
        data.cov <- data[, names(data) %in% names.cov, drop = FALSE]
        numeric.var <- names.cov[!sapply(names.cov, function(x, data) is.factor(data[,x]), data = data.cov)]
        if(length(numeric.var) != 0) {
            cov.std <- matrix(ncol = length(numeric.var), nrow = 2, dimnames = list(c("Mean", "sd"), numeric.var))
            data.cov.std <- data.cov
            for(i in 1:length(numeric.var)) {
                cov.std[1,i] <- mean(data.cov[,numeric.var[i]], na.rm = TRUE)
                cov.std[2,i] <- sd(data.cov[,numeric.var[i]], na.rm = TRUE)
                data.cov.std[,numeric.var[i]] <- (data.cov[,numeric.var[i]] - cov.std[1,i])/cov.std[2,i]
            }
        } else {
            cov.std <- NULL
            data.cov.std <- data.cov
        }
    } else { # Only the intercept
        names.cov <- NULL
        data.cov <- NULL
        cov.std <- NULL
        data.cov.std <- NULL
    }	   
	II <- if(length(II)) {
		matrix(unlist(II), nrow = 2)
	} else {
		matrix(0, nrow = 2)
	}	   
	res <- list(response = response, II = II, h = unlist(h), by.var = unlist(by.var), 
		nseg = nseg, bdeg = bdeg, pord = pord, atau = atau, btau = btau,
		npartial = k, partial = partial, prior.2D = unlist(prior.2D),
		data.cov = data.cov, cov.std = cov.std, data.cov.std = data.cov.std)
	res
}