R/rpar.R

Defines functions drpar.mlogit prpar.mlogit qrpar.mlogit drpar.rpar prpar.rpar qrpar.rpar drpar prpar qrpar rg.mlogit stdev.mlogit med.mlogit mean.mlogit rg.rpar stdev.rpar med.rpar mean.rpar med rg stdev s2norm m2norm cov.mlogit cor.mlogit plot.rpar plot.mlogit summary.rpar print.rpar rpar

Documented in cor.mlogit cov.mlogit drpar drpar.mlogit drpar.rpar mean.mlogit mean.rpar med med.mlogit med.rpar plot.mlogit plot.rpar print.rpar prpar prpar.mlogit prpar.rpar qrpar qrpar.mlogit qrpar.rpar rg rg.mlogit rg.rpar rpar stdev stdev.mlogit stdev.rpar summary.rpar

#' random parameter objects
#' 
#' `rpar` objects contain the relevant information about estimated
#' random parameters. The homonymous function extract on `rpar` object
#' from a `mlogit` object.
#' 
#' `mlogit` objects contain an element called `rpar` which contain a
#' list of `rpar` objects, one for each estimated random
#' parameter. The `print` method prints the name of the distribution
#' and the parameter, the `summary` behave like the one for numeric
#' vectors.
#' 
#' @name rpar
#' @aliases rpar print.rpar summary.rpar
#' @param x,object a `mlogit` object,
#' @param par the name or the index of the parameters to be extracted
#'     ; if `NULL`, all the parameters are selected,
#' @param norm the coefficient used for normalization if any,
#' @param ... further arguments.
#' @param digits the number of digits
#' @param width the width of the printed output
#' @return a `rpar` object, which contains:
#'
#' - dist: the name of the distribution,
#' - mean: the first parameter of the distribution,
#' - sigma: the second parameter of the distribution,
#' - name: the name of the parameter.
#' 
#' @export
#' @author Yves Croissant
#' @seealso [mlogit()] for the estimation of a random parameters logit
#'     model.
#' @keywords regression
rpar <- function(x, par = NULL, norm = NULL, ...){
    if (is.null(par)) par <- names(x$rpar)
    if (length(par) == 1){
        result <- x$rpar[[par]]
        if (!is.null(norm)) result$norm <- abs(coef(x)[norm])
    }    
    else{
        result <- x$rpar[par]
        if (!is.null(norm))
            lapply(result, function(x){x$norm <- abs(coef(x)[norm]); x})
    }
    result
}

#' @rdname rpar
#' @export
print.rpar <- function(x, digits = max(3, getOption("digits") - 2),
                       width = getOption("width"), ...){
    dist <- switch(x$dist,
                   "n"  = "normal",
                   "ln" = "log-normal",
                   "cn" = "censored normal",
                   "t"  = "triangular",
                   "u"  = "uniform",
                   "zbu" = "uniform",
                   "zbt" = "triangular"
                   )
    npar1 <- switch(x$dist,
                    "n"  = "mean",
                    "ln" = "meanlog",
                    "cn" = "mean",
                    "t"  = "center",
                    "u"  = "center",
                    "zbu" = "center",
                    "zbt" = "center"
                    )
    
    npar2 <- switch(x$dist,
                    "n"  = "sd",
                    "ln" = "sdlog",
                    "cn" = "sd",
                    "t"  = "span",
                    "u"  = "span",
                    "zbu" = NA,
                    "zbt" = NA
                    )
    par1 <- x$mean
    par2 <- x$sigma
    cat(paste(dist, " distribution with parameters ",round(par1, 3),
              " (", npar1, ")", " and ", round(par2, 3),
              " (",npar2,")", "\n",sep = ""))
}

#' @rdname rpar
#' @export
summary.rpar <- function(object, ...){
    norm <- object$norm
    rg <- rg.rpar(object, norm)
    Q1 <- qrpar(object, norm)(0.25)
    M <-  qrpar(object, norm)(0.5)
    Q3 <- qrpar(object, norm)(0.75)
    m <- mean(object, norm)
    c('Min.' = as.numeric(rg[1]), '1st Qu.' = as.numeric(Q1), 'Median' = as.numeric(M),
      'Mean' = as.numeric(m), '3rd Qu.' = as.numeric(Q3), 'Max.' = as.numeric(rg[2]))
}

#' Plot of the distribution of estimated random parameters
#' 
#' Methods for `rpar` and `mlogit` objects which provide a plot of the
#' distribution of one or all of the estimated random parameters
#' 
#' For the `rpar` method, one plot is drawn. For the `mlogit` method,
#' one plot for each selected random parameter is drawn.
#'
#' @name plot.mlogit
#' @aliases plot.mlogit plot.rpar
#' @importFrom graphics lines plot polygon segments title
#' @param x a `mlogit` or a `rpar` object,
#' @param type the function to be plotted, whether the density or the
#'     probability density function,
#' @param par a subset of the random parameters ; if `NULL`, all the
#'     parameters are selected,
#' @param norm the coefficient's name for the `mlogit` method or the
#'     coefficient's value for the `rpar` method used for
#'     normalization,
#' @param ... further arguments, passed to `plot.rpar` for the
#'     `mlogit` method and to `plot` for the `rpar` method.
#' @export
#' @author Yves Croissant
#' @seealso [mlogit()] the estimation of random parameters logit
#'     models and [rpar()] for the description of `rpar` objects and
#'     [distribution] for functions which return informations about
#'     the distribution of random parameters.
#' @keywords regression
plot.mlogit <- function(x, par = NULL, norm = NULL, type = c("density", "probability"), ...){
    if (is.null(x$rpar)) stop("the plot method is only relevant for random parameters")
    if (is.null(par)) par <- names(x$rpar)
    if (!is.null(norm)) norm = abs(coef(x)[norm])
    rpar <- x$rpar[par]
    K <- length(rpar)
    if (K > 1){
        nrow <- 1 + (K > 2) + (K > 6)
        ncol <- 1 + (K > 1) + (K > 4)
        opar <- par(mfrow = c(nrow, ncol))
        for (i in names(rpar)){
            plot(rpar(x, i), norm = norm, type = type, ...)
        }
        par(opar)
    }
    else{
        plot(rpar(x, 1), norm = norm, type = type, ...)
    }
}

#' @rdname plot.mlogit
#' @export
plot.rpar <- function(x, norm = NULL, type = c("density", "probability"), ...){
    type <- match.arg(type)
    if (type == "density") f <- drpar
    if (type == "probability") f <- prpar
    marg <- .05
    law <- x$dist
    rg <- rg(x)
    low <- rg[1]
    np <- x$name
    neg.values <- ifelse(low < 0,TRUE,FALSE)
    up <- rg[2]
    if (!is.finite(low)) low <- qrpar(x, norm = norm)(0.005)
    if (!is.finite(up)) up <- qrpar(x, norm = norm)(0.995)
    ptstot <- gnrpoints(low, up, 1000)
    ytot <- do.call(f, list(x = x, norm = norm))(ptstot)
    ymax <- max(ytot)*(1 + marg)
    plot(ptstot, ytot, type = "n", ann = F, xaxs = "i", yaxs = "i",
         las = 1, ylim = c(0, ymax),
         xlim = c(low - marg *(up - low), up + marg * (up - low)))
    ma <- paste("Distribution of", np)
    if (neg.values){
        pourc0 <- prpar(x)(0)
        ma <- paste(ma, ":", round(pourc0 * 100, 0),"% of 0")
        if (type == "density"){
            if (low < 0){
                ptsneg <- gnrpoints(low, 0, 10)
                yneg <- do.call(f, list(x = x, norm = norm))(ptsneg)
                polygon(c(low, ptsneg, 0),c(0, yneg, 0), col = "lightblue", border = NA)
            }
        }
        else{
            segments(low - marg * (up - low), pourc0, 0, pourc0, lty = "dotted")
            segments(0, 0, 0, pourc0, lty = "dotted")
        }
    }
    lines(ptstot, ytot)
    if (law == "u" && type == "density"){
        segments(up, 0, up, drpar(x)(up))
        segments(low, 0, low, drpar(x)(low))
    }
    title(main = ma)
}

#' Correlation structure of the random parameters
#' 
#' Functions that extract the correlation structure of a mlogit object
#' 
#' @name cor.mlogit
#' @aliases cor.mlogit cov.mlogit
#' @param x an `mlogit` object with random parameters and `correlation
#'     = TRUE`.
#' @details These functions are deprecated, use [vcov][vcov.mlogit].
#'     instead.
#' @return A numerical matrix which returns either the correlation or
#'     the covariance matrix of the random parameters.
#' @export
#' @author Yves Croissant
#' @keywords regression
cor.mlogit <- function(x){
    if (is.null(x$rpar) || is.null(attr(x$rpar, 'covariance')))
        stop('cor.mlogit only relevant for random models with correlation')
    cor.mlogit <- cov.mlogit(x)
    K <- nrow(cor.mlogit)
    sd.mlogit <- sqrt(diag(cor.mlogit))
    for (i in 1:K){
        for (j in 1:K){
            cor.mlogit[i,j] <- cor.mlogit[i,j] / sd.mlogit[i] / sd.mlogit[j]
        }
    }
    cor.mlogit
}

#' @rdname cor.mlogit
#' @export
cov.mlogit <- function(x){
    if (is.null(x$rpar) || is.null(attr(x$rpar, 'covariance')))
        stop('cov.mlogit only relevant for random models with correlation')
    attr(x$rpar, 'covariance')
}

# if a normalization coefficient is used, m2norm and s2norm transform
# the two parameters of the distribution
m2norm <- function(m, dist, norm){
    switch(dist,
           "n"  = m / norm,
           "ln" = m - log(norm),
           "t"  = m / norm,
           "cn" = m / norm,
           "u"  = m / norm,
           "zbu" = m / norm,
           "zbt" = m / norm
           )
}

s2norm <- function(s, dist, norm){
    switch(dist,
           "n"  = s / norm,
           "ln" = s,
           "t"  = s / norm,
           "cn" = s / norm,
           "u"  = s / norm,
           "zbu" = s / norm,
           "zbt" = s / norm
           )
}

#' Functions used to describe the characteristics of estimated random
#' parameters
#'
#' @name distribution
#' 
#' @aliases distribution med rg stdev mean.rpar med.rpar stdev.rpar
#'     rg.mlogit mean.mlogit med.mlogit stdev.mlogit rg.rpar qrpar
#'     prpar drpar qrpar.rpar prpar.rpar drpar.rpar qrpar.mlogit
#'     prpar.mlogit drpar.mlogit
#' @param x a `mlogit` or a `rpar` object,
#' @param norm the variable used for normalization if any : for the
#'     `mlogit` method, this should be the name of the parameter, for
#'     the `rpar` method the absolute value of the parameter,
#' @param par the required parameter(s) for the `mlogit` methods
#'     (either the name or the position of the parameter(s). If
#'     `NULL`, all the random parameters are used.
#' @param y values for which the function has to be evaluated,
#' @param ... further arguments.
#'
#' @details `rpar` objects contain all the relevant information about
#'     the distribution of random parameters. These functions enables
#'     to obtain easily descriptive statistics, density, probability
#'     and quantiles of the distribution.
#' 
#' `mean`, `med`, `stdev` and `rg` compute respectively the mean, the
#' median, the standard deviation and the range of the random
#' parameter. `qrpar`, `prpar`, `drpar` return functions that compute
#' the quantiles, the probability and the density of the random
#' parameters (note that `sd` and `range` are not generic function in
#' `R` and that `median` is, but without `...`).
#' 
#' @return a numeric vector for `qrpar`, `drpar` and `prpar`, a
#'     numeric vector for `mean`, `stdev` and `med` and a numeric
#'     matrix for `rg`.
#' @author Yves Croissant
#' @seealso [mlogit()] for the estimation of random parameters logit
#'     models and [rpar()] for the description of `rpar` objects.
#' @keywords regression


# mean, med, rg and stdev methods for rpar and mlogit objects (sd and
# range are not generic, so create a stdev and a rg generic ; median
# is generic, but without ..., so create a med generic

#' @rdname distribution
#' @export
stdev <- function(x, ...){
    UseMethod("stdev")
}

#' @rdname distribution
#' @export
rg <- function(x, ...){
    UseMethod("rg")
}

#' @rdname distribution
#' @export
med <- function(x, ...){
    UseMethod("med")
}

#' @rdname distribution
#' @export
mean.rpar <- function(x, norm = NULL, ...){
    if (is.null(norm) & ! is.null(x$norm)) norm <- as.numeric(x$norm)
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (! is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm)
    }
    switch(dist,
           "n"  = m,
           "ln" = exp(m + 0.5 * s^2),
           "u"  = m,
           "t"  = m,
           "cn" = s * dnorm(- m / s) + m * (1 - pnorm(- m / s)),
           "zbu" = m,
           "zbt" = m
           )
}

#' @rdname distribution
#' @export
med.rpar <- function(x, norm = NULL, ...){
    if (is.null(norm) & ! is.null(x$norm)) norm <- as.numeric(x$norm)
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (! is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm) 
    }
    switch(dist,
           "n"  = m,
           "ln" = exp(m),
           "u"  = m,
           "t"  = m,
           "cn" = max(0, m),
           "zbu" = m,
           "zbt" = m
           )
}

#' @rdname distribution
#' @export
stdev.rpar <- function(x, norm = NULL, ...){
    if (is.null(norm) & ! is.null(x$norm)) norm <- as.numeric(x$norm)
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (! is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm)
    }
    switch(dist,
           "n"  = s,
           "ln" = sqrt(exp(s ^ 2) - 1) * exp(m + 0.5 * s ^ 2),
           "u"  = s ^ 2 / 3,
           "t"  = s,
           "cn" = sqrt( s ^ 2 * (1 - pnorm(- m / s)) + m * (s * dnorm(- m / s) + m * (1 - pnorm(- m / s))) -
                        (s * dnorm(- m / s) + m * (1 - pnorm(- m / s))) ^ 2),
           "zbu" = m ^ 2 / 3,
           "zbt" = m
           )
}


#' @rdname distribution
#' @export
rg.rpar <- function(x, norm = NULL, ...){
    if (is.null(norm) & ! is.null(x$norm)) norm <- as.numeric(x$norm)
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (! is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm)
    }
    result <- switch(dist,
                     "n"  = c(-Inf, +Inf),
                     "ln" = c(0, +Inf),
                     "u"  = c(m - s, m + s),
                     "t"  = c(m - s, m + s),
                     "cn" = c(0, +Inf),
                     "zbu" = c(0, 2 * m),
                     "zbt" = c(0, 2 * m)
                     )
    names(result) <- c('Min.', 'Max.')
    result
}

#' @rdname distribution
#' @export
mean.mlogit <- function(x, par = NULL, norm = NULL, ...){
    if (!is.null(norm)) norm <- abs(coef(x)[norm])
    if (is.null(par)) par <- names(x$rpar)
    sapply(x$rpar[par], function(x) mean(x, norm = norm, ...))
}

#' @rdname distribution
#' @export
med.mlogit <- function(x, par = NULL, norm = NULL, ...){
    if (!is.null(norm)) norm <- abs(coef(x)[norm])
    if (is.null(par)) par <- names(x$rpar)
    sapply(x$rpar[par], function(x) med(x, norm = norm, ...))
}

#' @rdname distribution
#' @export
stdev.mlogit <- function(x, par = NULL, norm = NULL, ...){
    if (!is.null(norm)) norm <- abs(coef(x)[norm])
    if (is.null(par)) par <- names(x$rpar)
    sapply(x$rpar[par], function(x) stdev(x, norm = norm, ...))
}

#' @rdname distribution
#' @export
rg.mlogit <- function(x, par = NULL, norm = NULL, ...){
    if (!is.null(norm)) norm <- abs(coef(x)[norm])
    if (is.null(par)) par <- names(x$rpar)
    if (length(par) > 1)
        sapply(x$rpar[par], function(x) rg(x, norm = norm, ...))
    else rg(x$rpar[[par]])
}

#' @rdname distribution
#' @export
qrpar <- function(x, ...){
    UseMethod("qrpar")
}

#' @rdname distribution
#' @export
prpar <- function(x, ...){
    UseMethod("prpar")
}

#' @rdname distribution
#' @export
drpar <- function(x, ...){
    UseMethod("drpar")
}

#' @rdname distribution
#' @export
qrpar.rpar <- function(x, norm = NULL, ...){
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (!is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm)
    }
    switch(dist,
           "n"  = function(x = (1:9) / 10) qnorm(x, m, s),
           "ln" = function(x = (1:9) / 10) qlnorm(x, m, s),
           "u"  = function(x = (1:9) / 10) qunif(x, m - s, m + s),
           "t"  = function(x = (1:9) / 10) (m - s + sqrt(2 * s^2 * x)) *
                                       (x <= 0.5) + (m + s - sqrt(2 * s^2 *(1 - x))) *(x > 0.5),
           "cn" = function(x=(1:9)/10) max(0, qnorm(x, m, s)),
           "zbu" = function(x = (1:9) / 10) qunif(x, 0, 2 * m),
           "zbt" = function(x = (1:9) / 10) sqrt(2 * m ^ 2 * x) *
                                            (x <= 0.5) + (2 * m - sqrt(2 * m ^ 2 *(1 - x))) * (x > 0.5)
           )
}

#' @rdname distribution
#' @export
prpar.rpar <- function(x, norm = NULL, ...){
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (!is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm)
    }
    switch(dist,
           "n"  = function(x) pnorm(x, m, s),
           "ln" = function(x) plnorm(x, m, s),
           "u"  = function(x) punif(x, m - s, m + s),
           "t"  = function(x) (x >= (m - s) & x < m) * (x - m + s) ^ 2 / (2 * s ^ 2) +
                              (x >= m & x <= (m + s)) * (1 - (m + s - x) ^ 2 / (2 * s ^ 2)) + (x > (m + s)) * 1 + 0,
           "cn" = function(x) pnorm(x, m, s),
           "zbu" = function(x) punif(x, 0, 2 * m),
           "zbt"  = function(x) (x >= 0 & x < m) * x ^ 2 / (2 * m ^ 2) +
                                (x >= m & x <= 2 * m) * (1 - (2 * m) ^ 2 / (2 * m ^ 2)) + (x > 2 * m) * 1 + 0

           )
}

#' @rdname distribution
#' @export
drpar.rpar <- function(x, norm = NULL, ...){
    dist <- x$dist
    m <- x$mean
    s <- abs(x$sigma)
    if (!is.null(norm)){
        s <- s2norm(s, dist, norm)
        m <- m2norm(m, dist, norm)
    }
    switch(dist,
           "n"  = function(x) dnorm(x, m, s),
           "ln" = function(x) dlnorm(x, m, s),
           "u"  = function(x) (1 / s + x * 0) * (x >= m - s & x <= m + s) + 0,
           "t"  = function(x) (x >= (m - s) & x < m) * (x - m + s) / s ^ 2 +
                              (x >= m & x <= (m + s)) * (s + m - x) / s ^ 2 + 0,
           "cn" = function(x) dnorm(x, m, s),
           "zbu" = function(x) dunif(x, 0, 2 * m),
           "zbt" = function(x) (x >= 0 & x < m) * x / m ^ 2 +
                               (x >= m & x <= 2 * m) * (2 * m - x) / m ^ 2 + 0,
           )
}

#' @rdname distribution
#' @export
qrpar.mlogit <- function(x, par = 1, y = NULL, norm = NULL, ...){
    if (is.null(rpar))
        stop("qrpar function only relevant for random parameters models")
    if (length(par) > 1)
        stop("only one parameter should be selected")
    x <- x$rpar[[par]]
    if (is.null(norm)) norm <- abs(coef(x)[norm])
    if (is.null(y)){
        qrpar(x, norm = norm, ...)
    }
    else{
        qrpar(x, norm = norm, ...)(y)
    }
}

#' @rdname distribution
#' @export
prpar.mlogit <- function(x, par = 1, y = NULL, norm = NULL, ...){
    if (is.null(rpar))
        stop("prpar function only relevant for random parameters models")
    if (length(par) > 1)
        stop("only one parameter should be selected")
    x <- x$rpar[[par]]
    if (!is.null(norm)) norm <- coef(x)[norm]
    if (is.null(y)){
        prpar(x, norm = norm, ...)
    }
    else{
        prpar(x, norm = norm, ...)(y)
    }
}

#' @rdname distribution
#' @export
drpar.mlogit <- function(x, par = 1, y = NULL, norm = NULL, ...){
    if (is.null(rpar))
        stop("drpar function only relevant for random parameters models")
    if (length(par) > 1)
        stop("only one parameter should be selected")
    x <- x$rpar[[par]]
    if (!is.null(norm)) norm <- coef(x)[norm]
    if (is.null(y)){
        drpar(x, norm = norm, ...)
    }
    else{
        drpar(x, norm = norm, ...)(y)
    }
}

Try the mlogit package in your browser

Any scripts or data that you put into this service are public.

mlogit documentation built on Oct. 23, 2020, 5:29 p.m.