#' @importFrom nlme gls
standardES <- function( sigma=NA, sigma_1=NA, sigma_2=NA, nest=1:3, standardize=TRUE) {
if (standardize) {
nest <- c(nest, c(5,5,5,5)) # add dummy slots so nest has length >= 4
nest[1:2] <- sort(nest[1:2])
nest[3:4] <- sort(nest[3:4])
nest <- nest[nest < 5] # remove dummy slots
}
if (length(nest) == 2L) {
res <- data.frame(a=sigma_1)
names(res) <- paste0("sigma_", nest[1], nest[2])
return(res)
}
if (length(nest) == 3L) {
res <- data.frame(a = sigma_1, b = sigma)
names(res) <- c( paste0("sigma_", nest[1], nest[2]),
paste0("sigma_", nest[1], nest[2], "_", nest[3])
)
return(res)
}
res <- data.frame(a = sigma_1, b=sigma_2, c=sigma)
names(res) <- c( paste0("sigma_", nest[1], nest[2]),
paste0("sigma_", nest[3], nest[4]),
paste0("sigma_", nest[1], nest[2], "_", nest[3], nest[4])
)
return(res)
}
# This function provides alpha_i values for the CES model.
# It converts from delta_1 and delta to alpha_i in a nest-aware way.
# The standardCoefs provide the importance of each parameter
# IN THE ORDER THEY APPEAR IN THE FORMULA,
# regardless of which nest you choose
# and regardless of how the values for delta and delta_1 turn out.
# The point here is that the values of delta and delta_1 are applicable to the permuted
# variable order during the fitting process.
# But, we want coefficients that apply to the parameters in the order they appeared
# in the formula.
# standardCoefs does just that.
# Example: A formula is specified as y ~ a + b + c + time,
# nest is given as c(3,1,2) such that the parameters (during the estimation process)
# are permuted to y ~ c + a + b + time.
# Results will be the following:
# delta delta_1 alpha_1 alpha_2 alpha_3
# 1 1 0 0 1
# 1 0 1 0 0
# 0 0 or 1 0 1 0
standardCoefs <- function(delta=NA, delta_1=NA, nest=NULL, digits=5) {
# convert to standard coefficents taking nest order into account
# basically we are just permuting things so that alpha_i can always
# refer to the same quantities, even when they show up in different parts of the model.
# when nest is NULL, we make nest be 1:3
if (is.null(nest)) nest <- 1:3
# if delta is NA, set it to 1 or local calculation.
# This collapses the 2-factor case to the 3-factor case and handles cases where delta can't be approximated.
ldelta <- if (is.na(delta)) 1 else delta
ldelta_1 <- if (!is.na(delta) && round(delta, digits) == 0.0) 0.5 else delta_1
# make sure we have exactly 3 slots
nest <- head(c(nest,3), 3)
if (length(nest) < 3) stop("bad nest")
# alternative ways to handle delta = NA
# second works well for 2-factor models
modprod <- function(a, b) { # want 0 * NA == NA * 0 == 0
if (!is.na(a) && a==0) return(0)
if (!is.na(b) && b==0) return(0)
a * b
}
# res <- switch(method,
# "1" = list(delta * ldelta_1, delta * (1.0 - ldelta_1), 1.0 - delta), # NA, NA, NA
# "2" = list(ldelta * delta_1, ldelta * (1.0 - delta_1), 1.0 - delta), # delta_1, 1-delta_1, NA
# "3" = list(ldelta * delta_1, ldelta * (1.0 - delta_1), 1.0 - ldelta), # delta_1, 1-delta_1, 0
# "4" = list(
# modprod(delta, delta_1),
# modprod(delta, (1.0 - delta_1)),
# 1.0 - delta)
# )
res <- list( modprod(delta, delta_1),
modprod(delta, (1.0 - delta_1)),
1.0 - delta)
# invert the nest permutation
res<- res[order(nest)]
names(res) <- c("alpha_1", "alpha_2", "alpha_3")
return(as.data.frame(res))
}
nestMatch <- function( n1, n2 ) {
(length(n1) == length(n2)) && all(n1==n2)
}
#' Natural coefficients of a model
#'
#' A convenience function that returns the "natural" coefficients of a model object.
#' @param object the model object from which you want to extract the natural coefficients.
#' @param ... additional arguments
#' @return the coefficients of the model on a "natural" scale.
#' @export
naturalCoef <- function(object, ...) {
UseMethod("naturalCoef")
}
#' @export
naturalCoef.default <- function(object, ...) {
attr(object, "naturalCoeffs")
}
#' @export
naturalCoef.SFmodel <- function(object, ...) {
as.data.frame(
dplyr::tibble(
logscale=coef(object)[1],
scale=exp(logscale),
lambda = if (object$winner == 2) coef(object)[2] else coef(object)[3],
m = if (object$winner == 2) 1.0 else coef(object)[2]
)
)
}
#' @export
naturalCoef.CDEmodel <- function(object, ...) {
leftCoef <- c(3,2,3,3,1,2,3)[object$winner]
cf <- coef(object)[c("alpha_1", "alpha_2", "alpha_3")]
names(cf) <- c("alpha_1", "alpha_2", "alpha_3")
cf[is.na(cf)] <- 0
cf[leftCoef] <- 1 - sum(cf)
as.data.frame(
dplyr::tibble(
lambda = coef(object)["lambda"],
logscale = coef(object)["logscale"],
scale = exp(logscale),
alpha_1 = cf["alpha_1"],
alpha_2 = cf["alpha_2"],
alpha_3 = cf["alpha_3"]
)
)
}
#' @export
naturalCoef.LINEXmodel <- function( object, ...) {
as.data.frame(
dplyr::tibble(
logscale = coef(object)[1],
scale = exp(logscale),
a_0 = coef(object)[2],
a_1 = coef(object)[3],
c_t = a_1 / a_0
)
)
}
#' @export
naturalCoef.plm <- function(object, ...) {
makeNatCoef(object, ...)
}
#' @export
naturalCoef.cesEst <- function(object, ...) {
makeNatCoef(object, ...)
}
makeNatCoef <- function(object, nest=object$nest, ...) {
coefList <- as.list(coef(object))
if (inherits(object, "cesEst") && length(coefList) == 4) {
# We disagree with the naming of the coefficients by cesEst when only two factors of production are involved.
# cesEst calls the coefficients gamma, lambda, delta, and rho.
# Depending on nest, we sometimes prefer gamma, lambda, delta_1 and rho_1,
# because this case provides, essentially, the inner nest.
# This code appends "_1" to the appropriate names and adds appropriate delta, rho, and sigma values.
# We chose to rename the coefficients based on position.
# We could have renamed by name to defend against position changes in cesEst.
# But, cesEst could also change names in the future, so there is no clear benefit to renaming by name.
# So, we'll stick with renaming by position.
if (all(nest <= 2)) { # nest = 1 2 or 2 1
names(coefList)[3:4] <- paste0(names(coefList)[3:4], "_1")
coefList[["delta"]] <- 1
coefList[["sigma"]] <- NA
coefList[["rho"]] <- NA
} else {
if (any(nest == 1)) { # nest = 1 3 or 3 1
coefList[["delta_1"]] <- 1
coefList[["sigma_1"]] <- NA
coefList[["rho_1"]] <- NA
} else { # nest = 2 3 or 3 2
coefList[["delta_1"]] <- 0
coefList[["sigma_1"]] <- NA
coefList[["rho_1"]] <- NA
}
}
# special treatment for boundary model problem child
# in this case alpha_1 and alpha_2 don't make sense because we are creating
# a new first factor as pmin(x1, x2)
if (!is.null(object$boundary) && grepl("^17:", object$boundary)) {
coefList[["delta_1"]] <- NA
coefList[["sigma_1"]] <- 0
}
}
# We also prefer the name "scale" for the pre-multiplier coefficient.
# cesEst gives this as "gamma"
scale <- ifelse(is.null(coefList$logscale), NA, exp(coefList$logscale))
if (is.na(scale)) scale <- coefList$gamma
if (is.null(scale)) scale <- NA
lambda <- ifelse(is.null(coefList$lambda), NA, coefList$lambda)
delta_1 <- ifelse(is.null(coefList$delta_1), NA, coefList$delta_1)
delta <- ifelse(is.null(coefList$delta), NA, coefList$delta)
rho_1 <- ifelse(is.null(coefList$rho_1), NA, coefList$rho_1)
sigma_1 <- ifelse(is.null(coefList$sigma_1), NA, coefList$sigma_1)
rho <- ifelse(is.null(coefList$rho), NA, coefList$rho)
sigma <- ifelse(is.null(coefList$sigma), NA, coefList$sigma)
# We no longer include sse in naturalCoef, because sse is not a coefficient
# sse <- sum(resid(object)^2)
if (is.null(nest)) { nest <- 1:3 }
sc <- standardCoefs(delta_1=delta_1, delta=delta, nest=nest)
res <-
as.data.frame(
dplyr::tibble(
scale = scale,
logscale = log(scale),
lambda = lambda,
delta = delta,
delta_1 = delta_1,
sigma_1 = if (is.na(sigma_1)) 1/(1 + rho_1) else sigma_1,
rho_1 = if (is.na(rho_1)) 1/sigma_1 - 1 else rho_1,
sigma = if (is.na(sigma)) 1/(1 + rho) else sigma,
rho = if (is.na(rho)) 1/sigma - 1 else rho,
alpha_1 = sc$alpha_1,
alpha_2 = sc$alpha_2,
alpha_3 = sc$alpha_3
)
)
bind_cols( res, standardES(sigma = res$sigma, sigma_1 = res$sigma_1, nest=nest))
}
#' Extract CDEmodel Residuals
#'
#' The residuals of the linear model are extracted.
#' This method is identical to `residuals.gls()` except that the
#' "std" and "label" attributes are not retained in the residuals object.
#'
#' @param object a CDEmodel object
#' @param type an optional character string specifying the type of residuals to
#' be used. If "response", the "raw" residuals (observed - fitted) are used;
#' else, if "pearson", the standardized residuals (raw residuals divided by
#' the corresponding standard errors) are used; else, if "normalized", the
#' normalized residuals (standardized residuals pre-multiplied by the inverse
#' square-root factor of the estimated error correlation matrix) are used.
#' Partial matching of arguments is used, so only the first character needs to
#' be provided. Defaults to "response".
#' @param ... some methods for this generic function require additional
#' arguments. None are used in this method.
#' @export
residuals.CDEmodel <-
function (object, type = c("response", "pearson", "normalized"),
...)
{
type <- match.arg(type)
val <- object$residuals
if (type != "response") {
val <- val/attr(val, "std")
lab <- "Standardized residuals"
if (type == "normalized") {
if (!is.null(cSt <- object$modelStruct$corStruct)) {
val <- recalc(cSt, list(Xy = as.matrix(val)))$Xy[,
1]
lab <- "Normalized residuals"
}
}
}
else {
lab <- "Residuals"
if (!is.null(aux <- attr(object, "units")$y))
lab <- paste(lab, aux)
}
attr(val, "label") <- lab
if (!is.null(object$na.action)) {
res <- naresid(object$na.action, val)
# remove std and label attributes
attr(res, "std") <- NULL # naresid(object$na.action, attr(val, "std"))
attr(res, "label") <- NULL # attr(val, "label")
res
}
else {
# remove std and label attributes
attr(val, "std") <- NULL
attr(val, "label") <- NULL
val
}
}
#' Extract fitting history of a model
#'
#' This convenience function returns the fitting history for a model.
#' @param model the model object from which you want to extract the `history`.
#' @return the fitting history for a model.
#' @export
getHistory <- function(model) {
model$history
}
#' Extract data associated with an object
#'
#' Extract data associated with an object
#' @param object an object from which you want to extract data.
#' @return data associated with `object`.
#' @export
getData <- function(object, ...) {
UseMethod("getData")
}
#' @export
getData.default <- function(object, ...) {
object$data
}
#' @export
getData.SFmodel <- function(object, ...) {
object$data
}
#' @export
getData.CDEmodel <- function(object, ...) {
object$data
}
#' @export
getData.LINEXmodel <- function(object, ...) {
object$data
}
#' @export
getData.cesModel <- function(object, ...) {
object$data
}
#' @export
getData.lm <- function(object, full = TRUE, ...) {
if (full) {
eval(object$call[["data"]], environment(object$call[["formula"]]))
} else {
model.frame(object, ...)
}
}
#' Compute SSE from a model object
#'
#' Compute SSE from a model object and also check that certain constraints are met
#' by the coefficients.
#'
#' @param object a model object that inherits from `cesModel`.
#' @return a data frame with with variables `sse`, `constrained`, `sse.contrained`.
#'
sse <- function(object) {
sse <- sum(resid(object)^2)
coefs <- naturalCoef(object)
# call = Reduce(paste, gdata::trim(deparse(object$call))
constrained <-
(is.na(coefs$delta) || (0 <= coefs$delta && coefs$delta <= 1)) &&
(is.na(coefs$delta_1) || (0 <= coefs$delta_1 && coefs$delta_1 <= 1)) &&
(is.na(coefs$rho) || coefs$rho >= -1) &&
(is.na(coefs$rho_1) || coefs$rho_1 >= -1)
as.data.frame(
dplyr::tibble(
sse = sse,
constrained = constrained,
sse.constrained =
if(constrained) sse else Inf
)
)
}
#' Extract the best model (least sse) from a list of models
#'
#' @param models the list of models
#' @param digits the number of digits of `sse` that is considered significant
#' @param orderOnly if `FALSE`, returns a reordered list of models.
#' If `TRUE`, returns an integer vector representing the order of the models.
#' @note The ordering process preserves the original order in the event of ties at the desired significance level.
#' @note This function relies upon the `sse` value being stored
#' in an attribute of the model called `naturalCoeffs`.
#' @return an integer vector representing the order of the models if `orderOnly` is `TRUE`.
#' A reordered list of models if `orderOnly` is `FALSE` (the default).
#' @export
bestModel <- function(models, digits=6, orderOnly=FALSE, constrained=FALSE) {
if (constrained) {
o <- order(sapply( models, function(model) {
if (is.null(model)) NA else round(sse(model)$sse.constrained, digits=digits) } ) )
} else {
o <- order(sapply( models, function(model) {
if (is.null(model)) NA else round(sse(model)$sse, digits=digits) } ) )
}
if (orderOnly) return(o)
models[[ o[1] ]]
}
#' Compute fitted values on natural scale
#'
#' This is similar to `fitted`, but will invert logarithmic
#' transformations of the response variable for certain models (e.g., LINEX, and Cobb-Douglas
#' models fit in this package).
#'
#' @param object An object returned by a fitting function.
#' @export
yhat <- function(object, ...) {
UseMethod('yhat')
}
# This works for cesEst objects because fitted() retunrs fits on the natural scale.
#' @export
yhat.default <- function(object,...) {
fitted(object,...)
}
# Note: fitted() returns on the natural scale for cesEst objects but on the log-scale for plm objects,
# but both return residuals on the log-scale, so this works for both of them.
#' @export
yhat.cesModel <- function(object, ...) {
object$response / exp(resid(object))
}
#' @export
yhat.CDEmodel <- function( object, ... ) {
# model has form log(y) - log(x_0) ~ iYear + I(log x_1 - log x_0) + ... + I(log(x_k) - log(x_0))
#lx0 <- eval( parse( text = gsub( ".* - ", "", names(object$model)[1]) ), getData(object))
# grabbing log(x_0) from from the call
# 2 -> formula/model argument
# 2 -> lhs of ~
# 3 -> rhs of (last) - ; use ( ) to group for handcrafted hacky kludge
coef_env <- as.environment(as.list(object$coefficients))
parent.env(coef_env) <- parent.frame()
# two ways to skin this cat, but first requires more careful use of parens.
# lx0 <- eval(object$call[[2]][[2]][[3]], getData(object), coef_env)
lx0 <- 0 - eval(object$call[[2]][[2]], getData(object), coef_env) + log(object$response)
exp( fitted(object,...) + lx0[!is.na(lx0)] )
}
#' @export
yhat.LINEXmodel <- function( object, ... ) {
# model has form log(y) - log(x_0) ~ iYear + I(log x_1 - log x_0) + ... + I(log(x_k) - log(x_0))
# lx0 <- eval( parse( text = gsub( ".* - ", "", names(object$model)[1]) ), getData(object))
lx0 <- eval(object$call[[2]][[2]][[3]], getData(object))
exp( fitted(object, ...) + lx0[!is.na(lx0)] )
}
#' @export
yhat.SFmodel <- function( object, ... ) {
# model has form log(y) - log(x_0) ~ iYear + I(log x_1 - log x_0) + ... + I(log(x_k) - log(x_0))
# lx0 <- eval( parse( text = gsub( ".* - ", "", names(object$model)[1]) ), getData(object))
lx0 <- eval(object$call[[2]][[2]], getData(object))
exp( fitted(object, ...) + lx0[!is.na(lx0)] )
}
#' Return response values from original data used to fit a model
#'
#' This function returns the values of the original response variable
#' The values are calculated from fits and residuals.
#' For models of class `"LINEXmodel"`
#' `"CDEmodel"` or `"cesEst"`, the logarthmic transformation, if it
#' was used, will be undone, returning the values to their natural scale.
#' @param object a model object from one of the model fitting functions.
#' @param ... additional arguments
#' @return a numeric vector
#' @export
response <- function(object, ...) {
UseMethod('response')
}
#' @export
response.LINEXmodel <- function(object, ...)
return( object$response )
#' @export
response.CDEmodel <- function(object, ...)
return( object$response )
#' @export
response.cesModel <- function(object, ...) {
return( object$response )
}
#' @export
response.default <- function(object, ...) {
return( fitted(object) + resid(object) )
}
#' @export
predict.CDEmodel <- function( object, newdata, ... ) {
# model has form log(y) - log(x_0) ~ iYear + I(log x_1 - log x_0) + ... + I(log(x_k) - log(x_0))
# 2 -> lhs of ~
coef_env <- as.environment(as.list(object$coefficients))
parent.env(coef_env) <- parent.frame()
if (missing(newdata)) {
lx0 <- 0 - eval(object$call$formula[[2]], getData(object), coef_env) + log(object$response)
lx0 <- lx0[!is.na(lx0)]
} else {
lx0 <- eval(object$call$formula[[2]][[3]], newdata, coef_env)
}
exp( NextMethod() + lx0 )
}
#' @export
predict.LINEXmodel <- function( object, ... ) {
# model has form log(y) - log(x_0) ~ iYear + I(log x_1 - log x_0) + ... + I(log(x_k) - log(x_0))
# lx0 <- eval( parse( text = gsub( ".* - ", "", names(object$model)[1]) ), getData(object))
lx0 <- eval(object$call[[2]][[2]][[3]], getData(object))
exp( NextMethod() + lx0[!is.na(lx0)] )
}
#' Fit single factor models
#'
#' @param formula a formula of the form ` y ~ factor + time `
#' @param data a data frame in which `formula` is evaluated
#' @param constrained a logical indicating whether the model parameters are constrained
#' @param correlation an optional [`nlme::corStruct`] object describing the
#' within-group correlation structure. See the documentation of [`nlme::corClasses`]
#' for a description of the available [`nlme::corStruct`] classes.
#' If a grouping variable is to be used,
#' it must be specified in the `form` argument to the [`nlme::corStruct`] constructor.
#' Defaults to `NULL`, corresponding to uncorrelated errors.
#' @return an lm object with some additional attributes
#' @examples
#' data(EconUK)
#' sfModel(iGDP ~ iK + iYear, data = EconUK)
#' sfModel(response = "iGDP", factor = "iK", time = "iYear", data = EconUK)
#'
#' @export
sfModel <- function(formula, data, response, factor, time, constrained=FALSE,
save.data=TRUE, correlation = NULL) {
if ( missing(formula) ) {
formula <- substitute( response ~ factor + time,
list( response = as_name_or_null(response),
factor = as_name_or_null(factor),
time = as_name_or_null(time)
)
)
}
formulas <- list( log(y) ~ x + time,
log(y) - log(x) ~ time)
formulas <- lapply( formulas,
function(x) do.call(substitute,
list(x,
list( y = formula[[2]],
x = formula[[3]][[2]],
time = formula[[3]][[3]]
)
)
)
)
if (constrained){
res <- eval( substitute(nlme::gls(f, data=data, correlation = C),
list(C = correlation, f=formulas[[2]])) )
res$winner <- 2
logscale <- coef(res)[1]
lambda <- coef(res)[2]
} else {
res <- eval( substitute(nlme::gls(f, data=data, correlation = C),
list(C = correlation, f=formulas[[1]])) )
res$winner <- 1
logscale <- coef(res)[1]
lambda <- as.vector(coef(res)[3])
}
sdata <- subset(data, select = all.vars(formula))
sdata <- data[complete.cases(sdata), unique(c(all.vars(formula), names(data)))]
if (save.data) { res$data <- sdata }
res$response <- eval( formula[[2]], sdata, parent.frame() )
res$formula <- formula
class(res) <- c("SFmodel", class(res))
return(res)
}
#' Fit Cobb-Douglas models
#'
#' @param formala a formula of the form
#' `response ~ x1 + x2 + time` or
#' `response ~ x1 + x2 + x3 + time`
#' @param data a data frame in which `formula` is evaluated
#' @param response instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param x1 instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param x2 instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param x3 instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param time instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param constrained a logical indicating whether the parameters should be constrained in the fitting process.
#' @return a CDEmodel object, which is an lm object with some additioanl attributes.
#' @examples
#' data(EconUK)
#' cdModel(response = "iGDP", x1 = "iK", x2 = "iL", time = "iYear", data = EconUK)
#' cdModel(iGDP ~ iK + iL + iYear, data = EconUK)
#'
#' @export
cdModel <- function(formula, data, response, x1, x2, x3, time,
constrained=TRUE, save.data=TRUE,
correlation = NULL, ...) {
if (missing(formula)) {
if (missing(x3)) {
return( cd2Model( data=data,
response=response,
x1=x1,
x2=x2,
time=time,
constrained = constrained, save.data=save.data,
correlation = correlation, ...) )
} else {
return( cd3Model( data=data,
response=response,
x1=x1,
x2=x2,
x3=x3,
time=time,
constrained = constrained, save.data=save.data,
correlation = correlation, ...) )
}
}
if (ncol( attr(terms(formula),"factors") ) == 3 ) {
# formula contains response, x1, and x2
cd2Model( formula=formula, data=data, constrained = constrained,
save.data=save.data, correlation = correlation, ... )
} else {
# formula contains response, x1, x2, and x3
cd3Model( formula=formula, data=data, constrained = constrained, save.data=save.data,
correlation = correlation, ... )
}
}
# These formulas provide C-D models at all boundaries.
# Also, these formulas all assume constant returns to scale.
CDformulas <- list(
log(y) - log(x3) ~
I(log(x1) - log(x3)) + I(log(x2) - log(x3)) + time, # [[1]] Full model
log(y) - log(x2) ~ I(log(x1) - log(x2)) + time, # [[2]] With x3 exponent = 0
log(y) - log(x3) ~ I(log(x1) - log(x3)) + time, # [[3]] With x2 exponent = 0
log(y) - log(x3) ~ I(log(x2) - log(x3)) + time, # [[4]] With x1 exponent = 0
log(y) - log(x1) ~ time, # [[5]] With x1 exponent = 1, x2 and x3 exponents = 0
log(y) - log(x2) ~ time, # [[6]] With x2 exponent = 1, x1 and x3 exponents = 0
log(y) - log(x3) ~ time # [[7]] With x3 exponent = 1, x1 and x2 exponents = 0
)
CDcoefNames <- list(
c("logscale", "alpha_1", "alpha_2", "lambda"), # [[1]] Full model
c("logscale", "alpha_1", "lambda"), # [[2]] With x3 exponent = 0
c("logscale", "alpha_1", "lambda"), # [[3]] With x2 exponent = 0
c("logscale", "alpha_2", "lambda"), # [[4]] With x1 exponent = 0
c("logscale", "lambda"), # [[5]] With x1 exponent = 1, x2 and x3 exponents = 0
c("logscale", "lambda"), # [[6]] With x2 exponent = 1, x1 and x3 exponents = 0
c("logscale", "lambda") # [[7]] With x3 exponent = 1, x1 and x2 exponents = 0
)
#' Fit Cobb-Douglas Models
#'
#' @param formala a formual of the form `response ~ x1 + x2 + time`
#' @param data a data fram in which `formula` is evaluated
#' @param response instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param x1 instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param x2 instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param time instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param constrained a logical indicating whether the parameters are constrained
#' @param correlation an optional [`nlme::corStruct`] object describing the
#' within-group correlation structure. See the documentation of [`nlme::corClasses`]
#' for a description of the available [`nlme::corStruct`] classes.
#' If a grouping variable is to be used,
#' it must be specified in the `form` argument to the [`nlme::corStruct`] constructor.
#' Defaults to `NULL`, corresponding to uncorrelated errors.
#' @return a CDEmodel object, which is an lm object with some additional attributes.
#' @keywords internal
#' @examples
#' data(EconUK)
#' cd2Model(response = iGDP, x1 = iK, x2 = iL, time = iYear, data = EconUK)
#' cd2Model(iGDP ~ iK + iL + iYear, data = EconUK)
#' cd2Model(response = "iGDP", x1="iK", x2="iL", time="iYear", data = EconUK)
#'
cd2Model <- function(formula, data, response, x1, x2, time, constrained=TRUE,
save.data=TRUE, correlation = NULL, ...) {
if ( missing(formula) ) {
formula <- substitute( response ~ x1 + x2 + time,
list( response = as_name_or_null(response),
x1 = as_name_or_null(x1),
x2 = as_name_or_null(x2),
time = as_name_or_null(time)
)
)
}
sdata <- subset(data, select=all.vars(formula))
sdata <- data[complete.cases(sdata), unique(c(all.vars(formula), names(data)))]
# we will only look at formulas[[c(2, 5, 6)]] below; but keep all for indexing convenience
formulas <- lapply( CDformulas,
function(x) do.call(substitute,
list(x,
list( y = formula[[2]],
x1 = formula[[3]][[2]][[2]],
x2 = formula[[3]][[2]][[3]],
time = formula[[3]][[3]]
)
)
)
)
res <- eval(substitute(nlme::gls( f, correlation = C, data=sdata ),
list(C = correlation, f=formulas[[2]])))
winner <- 2
names(res$coefficients) <- CDcoefNames[[2]]
alpha_1 <- coef(res)["alpha_1"]
if (constrained){
# Adjust alpha_1 if we are outside the range 0 <= alpha_1 <= 1
if (alpha_1 < 0.0){
alpha_1 <- 0.0
winner <- 6
} else if (alpha_1 > 1) {
alpha_1 <- 1.0
winner <- 5
}
if (winner > 2) {
res <- eval(substitute(nlme::gls( f, data=sdata, correlation = C ),
list(C = correlation, f=formulas[[winner]])))
names(res$coefficients) <- CDcoefNames[[winner]]
}
}
if (save.data) {
res$data <- sdata
}
res$response <- eval( formula[[2]], sdata, parent.frame() )
res$formula <- formula
res$winner <- winner
class(res) <- c("CDEmodel", class(res))
return(res)
}
respectsConstraints <- function( model ) {
# Tells whether we're fitting within the constraints for a Cobb-Douglas model.
# This assumes that the model has coefficients named alpha_i or that
# the paramters of interest are all but the first (intercept) and last (time).
cf <- coef(model)
if (any(c("alpha_1","alpha_2","alpha_3") %in% names(cf) ) ) {
cf <- cf[ names(cf) %in% c("alpha_1", "alpha_2", "alpha_3") ]
} else {
# get rid of first and last (intercept and time)
cf <- tail(head(cf,-1), -1)
}
all( cf >= 0 ) & ( sum(cf) <=1 )
}
#' Fit Cobb-Douglas models with 3 Factors
#'
#' @param formula a formula of the form `response ~ x1 + x2 + x3 + time`
#' @param data a data frame in which `formala` is evaluated
#' @param response instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param x1 instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param x2 instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param x3 instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param time instead of specifying a formula, expressions for
#' the components can be specified individually as character strings.
#' @param constrained a logical indicated whether the coefficents are constrained. See details
#' @param correlation an optional [`nlme::corStruct`] object describing the
#' within-group correlation structure. See the documentation of [`nlme::corClasses`]
#' for a description of the available [`nlme::corStruct`] classes.
#' If a grouping variable is to be used,
#' it must be specified in the `form` argument to the [`nlme::corStruct`] constructor.
#' Defaults to `NULL`, corresponding to uncorrelated errors.
#' @param \dots additional arguments; currently unused.
#' @details More about contranints TBA.
#' @keywords internal
#' @examples
#' data(EconUK)
#' cd3Model(response = iGDP, x1 = iK, x2 = iL, x3 = iQp, time = iYear, data = EconUK)
#' cd3Model(iGDP ~ iK + iL + iQp + iYear, data = EconUK)
#' cd3Model(response = "iGDP", x1="iK", x2="iL", x3="iQp", time="iYear", data = EconUK)
# y ~ x1 + x2 + x3 + time
cd3Model <- function( formula, data, response, x1, x2, x3, time,
constrained=TRUE, save.data=TRUE,
correlation = NULL,
...){
if ( missing(formula) ) {
formula <- substitute( response ~ x1 + x2 + x3 + time,
list( response = as_name_or_null(response),
x1 = as_name_or_null(x1),
x2 = as_name_or_null(x2),
x3 = as_name_or_null(x3),
time = as_name_or_null(time)
)
)
}
sdata <- subset(data, select = all.vars(formula))
sdata <- data[complete.cases(sdata), unique(c(all.vars(formula), names(data)))]
formulas <- lapply(CDformulas, function(x) do.call( substitute, list( x, list(
time = formula[[3]][[3]],
x3 = formula[[3]][[2]][[3]],
x2 = formula[[3]][[2]][[2]][[3]],
x1 = formula[[3]][[2]][[2]][[2]],
y = formula[[2]]
)
) ) )
models <- lapply( formulas, function(form)
eval(substitute(nlme::gls(f, correlation = C, data=sdata), list(f=form, C = correlation)))
)
sse <- sapply( models, function(m) sum( resid(m)^2 ) )
logLik <- sapply(models, function(m) m$logLik)
if ( constrained ) {
good <- sapply( models, respectsConstraints)
good[ !good ] <- NA
} else {
good <- TRUE
}
winner <- which.max( logLik * good )
res <- models[[winner]]
names( res$coefficients ) <- CDcoefNames[[winner]]
res$winner <- winner
res$formula <- formula
attr(res, "good") <- sapply( models, respectsConstraints )
res$response <- eval( formula[[2]], sdata, parent.frame() )
if (save.data) { res$data <- sdata }
class(res) <- c( "CDEmodel", class(res) )
return(res)
}
#' Fit LINEX models
#'
#' @param formula a formula of the form `response ~ x1 + x2 + x3 + time`
#' @param response instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param x1 instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param x2 instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param x3 instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param time instead of specifying a formula, expressions for
#' the components can be specified individually.
#' @param data a data frame in which `formula` is evaluated
#' @param correlation an optional [`nlme::corStruct`] object describing the
#' within-group correlation structure. See the documentation of [`nlme::corClasses`]
#' for a description of the available [`nlme::corStruct`] classes.
#' If a grouping variable is to be used,
#' it must be specified in the `form` argument to the [`nlme::corStruct`] constructor.
#' Defaults to `NULL`, corresponding to uncorrelated errors.
#' @examples
#' data(EconUK)
#' linexModel(iGDP ~ iK + iL + iXp + iYear, data = EconUK)
#' linexModel(response = "iGDP", x1="iK", x2="iL", x3="iXp", time="iYear", data = EconUK)
#' naturalCoef(linexModel(response = "iGDP", x1="iK", x2="iL", x3="iXp", time="iYear", data = EconUK))
#'
#' @export
linexModel <- function(formula, data, response, x1, x2, x3, time, save.data=TRUE,
correlation = NULL) {
if ( missing(formula) ) {
formula <- substitute( response ~ x1 + x2 + x3 + time,
list( response = as_name_or_null(response),
x1 = as_name_or_null(x1),
x2 = as_name_or_null(x2),
x3 = as_name_or_null(x3),
time = as_name_or_null(time)
)
)
}
formulas <- list( log(y) - log(x3) ~
I(2 * (1 - 1/(x1 / ( .5 * (x3 + x2) ) )) ) +
I( x2/x3 - 1 )
)
formulas <- lapply(formulas, function(x) do.call( substitute, list( x, list(
time = formula[[3]][[3]],
x3 = formula[[3]][[2]][[3]],
x2 = formula[[3]][[2]][[2]][[3]],
x1 = formula[[3]][[2]][[2]][[2]],
y = formula[[2]]
)
) ) )
res <- eval( substitute(nlme::gls(f, data=data, correlation = C),
list(C = correlation, f=formulas[[1]])) )
sdata <- subset(data, select = all.vars(formula))
sdata <- data[complete.cases(sdata), unique(c(all.vars(formula), names(data)))]
if (save.data) { res$data <- sdata }
res$response <- eval( formula[[2]], sdata, parent.frame() )
res$formula <- formula
class(res) <- c("LINEXmodel", class(res))
return(res)
}
nautralCoef.gls <- function(object, ...) {
EconModels:::makeNatCoef(object, ...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.