diff --git a/DESCRIPTION b/DESCRIPTION index c39c8e0c..f7370f39 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -32,6 +32,7 @@ Imports: reshape2, lattice, grid, + gridExtra, tweedie, utils, systemfit, diff --git a/NAMESPACE b/NAMESPACE index 62385ad4..05e173a3 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,7 @@ import(stats) -importFrom(ggplot2, ggplot, geom_density, facet_wrap, aes, "scale_x_continuous") +import(ggplot2) +#importFrom(ggplot2, ggplot, geom_density, facet_wrap, aes, "scale_x_continuous") +importFrom(gridExtra, marrangeGrob) importFrom(tweedie, rtweedie, "tweedie.profile") importFrom(statmod, tweedie) importFrom(systemfit, systemfit, systemfit.control) @@ -20,7 +22,8 @@ importFrom(cplm, cpglm) importFrom(MASS, glm.nb) importClassesFrom(cplm, NullNum, NullList) -importMethodsFrom(cplm, show, predict, vcov, residuals, resid, coef, fitted, plot, terms) +importMethodsFrom(cplm, show, predict, vcov, residuals, resid, coef, fitted, + plot, terms) importMethodsFrom(Matrix, summary) export(MackChainLadder, MunichChainLadder, BootChainLadder, @@ -59,7 +62,6 @@ S3method(print, ClarkLDF) S3method(print, ClarkCapeCod) S3method(print, ata) - S3method(as.triangle, matrix) S3method(as.triangle, data.frame) S3method(as.data.frame, triangle) @@ -89,8 +91,10 @@ S3method(mean, BootChainLadder) S3method(vcov, clark) exportMethods(predict, Mse, summary, show, coerce, - "[", "$", "[[", "[<-", names, coef, vcov, residCov, residCor, - residuals, resid, rstandard, fitted, plot, cbind2, rbind2, dim) + "[", "$", "[[", "[<-", names, coef, + vcov, residCov, residCor, + residuals, resid, rstandard, fitted, + plot, cbind2, rbind2, dim) S3method(CDR, MackChainLadder) S3method(CDR, BootChainLadder) @@ -101,4 +105,8 @@ export(tweedieReserve) S3method(print, tweedieReserve) S3method(summary, tweedieReserve) +export(plotParms) +S3method(plotParms, ChainLadder) +S3method(plotParms, MackChainLadder) + export(PaidIncurredChain) diff --git a/R/ChainLadder.R b/R/ChainLadder.R index 228d4358..d5ba4041 100644 --- a/R/ChainLadder.R +++ b/R/ChainLadder.R @@ -5,6 +5,8 @@ chainladder <- function(Triangle, weights=1, delta=1){ + + TriangleName <- deparse(substitute(Triangle)) Triangle <- checkTriangle(Triangle) n <- dim(Triangle)[2] @@ -27,7 +29,9 @@ chainladder <- function(Triangle, weights=1, } myModel <- lapply(c(1:(n-1)), lmCL, Triangle) - output <- list(Models=myModel, Triangle=Triangle, delta=delta, weights=weights) + output <- list(Models=myModel, Triangle=Triangle, delta=delta, + weights=weights, TriangleName = TriangleName) + output[["call"]] <- match.call(expand.dots = FALSE) class(output) <- c("ChainLadder", "TriangleModel", class(output)) return(output) } diff --git a/R/MackChainLadderFunctions.R b/R/MackChainLadderFunctions.R index ba310ecd..be1131a5 100755 --- a/R/MackChainLadderFunctions.R +++ b/R/MackChainLadderFunctions.R @@ -16,12 +16,22 @@ MackChainLadder <- function( tail.sigma=NULL, mse.method = "Mack") { + TriangleName <- deparse(substitute(Triangle)) ## idea: have a list for tail factor ## tail=list(f=FALSE, f.se=NULL, sigma=NULL, F.se=NULL) ## - # 2013-02-25 Parameter risk recursive formula may have a third term per + tail.input <- tail + tail.se.input <- tail.se + tail.sigma.input <- tail.sigma + mse.method.input <- mse.method + # 2013-02-25 Parameter risk recursive formula may have a third term per # Murphy and BBMW if (! mse.method %in% c("Mack", "Independence")) stop("mse.method must be 'Mack' or 'Independence'") + if (! est.sigma[1] %in% c("Mack", "log-linear")) { + if (!is.numeric(est.sigma)) + stop("est.sigma must be 'Mack' or 'log-linear' or numeric") + } + if (!is.logical(tail) & !is.numeric(tail)) stop("tail must be logical or numeric") Triangle <- checkTriangle(Triangle) m <- dim(Triangle)[1] @@ -118,6 +128,13 @@ MackChainLadder <- function( output[["Total.ProcessRisk"]] <- attr(Total.SE, "processrisk") output[["Total.ParameterRisk"]] <- attr(Total.SE, "paramrisk") output[["tail"]] <- tail + output[["TriangleName"]] <- TriangleName + output[["est.sigma"]] <- est.sigma + output[["tail.input"]] <- tail.input + output[["tail.se.input"]] <- tail.se.input + output[["tail.sigma.input"]] <- tail.sigma.input + output[["mse.method"]] <- mse.method + output[["mse.method.input"]] <- mse.method.input class(output) <- c("MackChainLadder", "TriangleModel", "list") return(output) } @@ -347,15 +364,41 @@ tail.SE <- function(FullTriangle, StdErr, Total.SE, tail.factor, tail.se = NULL, if(is.null(tail.se)){ .fse <- StdErr$f.se[start:(n-2)] - mse <- lm(log(.fse) ~ .dev) - tail.se <- exp(predict(mse, newdata=data.frame(.dev=tail.pos))) + ndx <- .fse > 0 + numpos <- sum(ndx) + if (numpos) { + .devse <- .dev[ndx] + .fse <- .fse[ndx] + mse <- lm(log(.fse) ~ .devse) + if (numpos > 1) tail.se <- exp(predict(mse, newdata=data.frame(.devse=tail.pos))) + else { + warning("Only one column for estimating tail.se") + tail.se <- suppressWarnings( + exp(predict(mse, newdata=data.frame(.dev=tail.pos)))) + } + } + else { + warning("No development variability in triangle. Setting tail.se = 0") + tail.se <- 0 + } } StdErr$f.se <- c(StdErr$f.se, tail.se = tail.se) if(is.null(tail.sigma)){ .sigma <- StdErr$sigma[start:(n-2)] - msig <- lm(log(.sigma) ~ .dev) - tail.sigma <- exp(predict(msig, newdata = data.frame(.dev = tail.pos))) + ndx <- .sigma > 0 + numpos <- sum(ndx) + if (numpos) { + .devsigma <- .dev[ndx] + .sigma <- .sigma[ndx] + msig <- lm(log(.sigma) ~ .devsigma) + if (numpos > 1) tail.sigma <- exp( + predict(msig, newdata = data.frame(.devsigma = tail.pos))) + else { + warning("No development variability in triangle. Setting tail.sigma = 0") + tail.sigma <- 0 + } + } } StdErr$sigma <- c(StdErr$sigma, tail.sigma = as.numeric(tail.sigma)) diff --git a/R/plotParms.R b/R/plotParms.R new file mode 100644 index 00000000..1a954fcf --- /dev/null +++ b/R/plotParms.R @@ -0,0 +1,414 @@ +#' plotParms Plot the estimated parameters of a model +#' +#' Methods to visualize the estimated parameters. +#' S3 methods currently exist for objects of class +#' ChainLadder and MackChainLadder. +#' +#' @param x object whose parameters are to be visualized +#' +#' @param title optional; character holding title of the plot; +#' defaults to something the class author deems appropriate. +#library(ggplot2) +#library(grid) +#library(gridExtra) +# In case need for aes_: see exchange and final solution at link below +# http://stackoverflow.com/questions/9439256/ +# how-can-i-handle-r-cmd-check-no-visible-binding-for- +# global-variable-notes-when + +plotParms <- function(x, which, ncol, nrow, title, ...) UseMethod("plotParms") +plotParms.default <- function(x, ...){ + stop("No 'plotParms' method exists for objects of class ", + class(x)) +} +plotParms.ChainLadder <- function(x, which = c(1L, 3L, 4L, 6L), + ncol = min(2L, length(which)), + nrow = ceiling(length(which) / ncol), + title = deparse(x$call), ...) { + grobs <- lapply(which, function(i) + switch(i, + plot.cl.f(x) + theme(axis.title.y=element_blank()), + plot.cl.f1(x) + theme(axis.title.y=element_blank()), + plot.cl.sigma(x) + theme(axis.title.y=element_blank()), + plot.cl.f.se(x) + theme(axis.title.y=element_blank()), + plot.cl.f.cv(x) + theme(axis.title.y=element_blank()), + plot.cl.f1.cv(x) + theme(axis.title.y=element_blank()) + ) + ) +# if (missing(title)) title <- paste0("chainladder(", +# x$TriangleName, +# ") parameter estimates") + + marrangeGrob(grobs = grobs, ncol = ncol, nrow = nrow, top = title, ...) +} +plot.cl.f <- function(x) { + n <- length(x$Models) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f <- sapply(smmry, function(x) x$coef["x","Estimate"]) + f.se <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + legend_captions <- c("lm", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[is.na(f)] <- src[2L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f.se, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f, colour = source)) + + geom_errorbar(aes(ymin=f-f.se, ymax=f+f.se), colour="black", width=.1, na.rm = TRUE) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(group = 1), na.rm = TRUE) + + geom_point() + + ggtitle("f estimates") + P +} +plot.cl.f1 <- function(x) { + n <- length(x$Models) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f <- sapply(smmry, function(x) x$coef["x","Estimate"]) - 1 + f.se <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + legend_captions <- c("lm", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[is.na(f)] <- src[2L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f.se, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f, colour = source)) + + geom_errorbar(aes(ymin=f-f.se, ymax=f+f.se), colour="black", width=.1, na.rm = TRUE) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(group = 1), na.rm = TRUE) + + geom_point() + + ggtitle("f-1 estimates") + P +} +plot.cl.sigma <- function(x) { + # require(ggplot2) + n <- length(x$Models) + smmry <- suppressWarnings(lapply(x$Models, summary)) + sigma <- sapply(smmry, function(x) x$sigma) + legend_captions <- c("lm", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[is.na(sigma)] <- src[2L] + sigmaNoNAs <- sigma + sigmaNoNAs[is.na(sigma)] <- 0 + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, sigma, sigmaNoNAs, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = sigma, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(colour = source, group = 1), na.rm = TRUE) + + geom_point(aes(y=sigmaNoNAs)) + + ggtitle("sigma estimates") + P +} +plot.cl.f.se <- function(x) { + # require(ggplot2) + n <- length(x$Models) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f.se <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + legend_captions <- c("lm", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[is.na(f.se)] <- src[2L] + f.seNoNAs <- f.se + f.seNoNAs[is.na(f.se)] <- 0 + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f.se, f.seNoNAs, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f.se, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(colour = source, group = 1), na.rm = TRUE) + + geom_point(aes(y = f.seNoNAs)) + + ggtitle("f.se estimates") + P +} +plot.cl.f.cv <- function(x) { + n <- length(x$Models) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f <- sapply(smmry, function(x) x$coef["x","Estimate"]) + f.se <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + f.cv <- f.se / f + legend_captions <- c("calc", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[is.na(f.cv)] <- src[2L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f.cv, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f.cv, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + ylab("cv(f)") + + geom_line(aes(group = 1), na.rm = TRUE) + + geom_point(na.rm = TRUE) + + ggtitle("cv(f) estimates") + P +} +plot.cl.f1.cv <- function(x) { + n <- length(x$Models) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f <- sapply(smmry, function(x) x$coef["x","Estimate"]) - 1 + f.se <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + f.cv <- f.se / f + legend_captions <- c("calc", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[is.na(f.cv)] <- src[2L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f.cv, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f.cv, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + ylab("cv(f)") + + geom_line(aes(group = 1), na.rm = TRUE) + + geom_point(na.rm = TRUE) + + ggtitle("cv(f-1) estimates") + P +} + +# Now MackCL +plotParms.MackChainLadder <- function(x, which = c(1L, 3L, 4L, 6L), + ncol = min(2L, length(which)), + nrow = ceiling(length(which) / ncol), + title = deparse(x$call), ...) { + grobs <- lapply(which, function(i) + switch(i, + plot.mackcl.f(x) + theme(axis.title.y=element_blank()), + plot.mackcl.f1(x) + theme(axis.title.y=element_blank()), + plot.mackcl.sigma(x) + theme(axis.title.y=element_blank()), + plot.mackcl.f.se(x) + theme(axis.title.y=element_blank()), + plot.mackcl.f.cv(x) + theme(axis.title.y=element_blank()), + plot.mackcl.f1.cv(x) + theme(axis.title.y=element_blank()) + ) + ) +# if (missing(title)) title <- paste0("MackChainLadder(", +# x$TriangleName, +# ") parameter estimates") + + marrangeGrob(grobs = grobs, ncol = ncol, nrow = nrow, top = title, ...) +} +plot.mackcl.f <- function(x) { + smmry <- suppressWarnings(lapply(x$Models, summary)) + f <- x$f + n <- length(f) + f.se <- x$f.se + if (length(f.se) < n) f.se <- c(f.se, NA) + # set display order + legend_captions <- c("lm", "log-lin", "default", "input", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[n] <- src[ifelse(is.logical(x$tail.input), + ifelse(x$tail.input, 2L, 3L), + 4L)] + source.se <- rep(src[1L], n) + source.se[is.na(f.se)] <- src[5L] + f.se[is.na(f.se)] <- 0 + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f, f.se, # sigma, na_sigma, f.cv, + source, source.se, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f, colour = source)) + + geom_errorbar(aes(ymin=f-f.se, ymax=f+f.se), colour="black", width=.1) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(group = 1), na.rm = TRUE) + + geom_point() + + ggtitle("f estimates") + P +} +plot.mackcl.f1 <- function(x) { + smmry <- suppressWarnings(lapply(x$Models, summary)) + f <- x$f - 1 + n <- length(f) + f.se <- x$f.se + if (length(f.se) < n) f.se <- c(f.se, NA) + # set display order + legend_captions <- c("lm", "log-lin", "default", "input", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + source[n] <- src[ifelse(is.logical(x$tail.input), + ifelse(x$tail.input, 2L, 3L), + 4L)] + source.se <- rep(src[1L], n) + source.se[is.na(f.se)] <- src[5L] + f.se[is.na(f.se)] <- 0 + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + df <- data.frame(xx, f, f.se, # sigma, na_sigma, f.cv, + source, source.se, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f, colour = source)) + + geom_errorbar(aes(ymin=f-f.se, ymax=f+f.se), colour="black", width=.1) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(group = 1), na.rm = TRUE) + + geom_point() + + ggtitle("f-1 estimates") + P +} +plot.mackcl.sigma <- function(x) { + sigma <- x$sigma + n <- length(x$f) + if (length(sigma) < n) sigma <- c(sigma, NA) + smmry <- suppressWarnings(lapply(x$Models, summary)) + sigmaregr <- sapply(smmry, function(x) x$sigma) + if (length(sigmaregr) < n) sigmaregr <- + c(sigmaregr, rep(NA, n - length(sigmaregr))) + notregr <- sigma != sigmaregr + notregr[is.na(notregr)] <- TRUE + # set display order + legend_captions <- + c("lm", "log-lin", "Mack", "input", "input-x", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + # For periods with insufficient data to estimate sigma, + # MackChainLadder has ways to inpute a value + if (x$est.sigma[1] %in% "log-linear") source[notregr] <- src[2L] + else + if (x$est.sigma[1] %in% "Mack") source[notregr] <- src[3L] + else source[notregr] <- src[4L] + # tail + if (is.null(x$tail.sigma.input)){ + if (x$f[n] > 1) source[n] <- src[2L] + else source[n] <- src[6L] + } else if (!x$tail.input) source[n] <- src[5L] + else source[n] <- src[4L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + sigmaNoNAs <- sigma + + sigmaNoNAs[is.na(sigma)] <- 0 + df <- data.frame(xx, sigma, sigmaNoNAs, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = sigma, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(colour = source, group = 1), na.rm = TRUE) + + geom_point(aes(y=sigmaNoNAs)) + + ggtitle("sigma estimates") + P +} +plot.mackcl.f.se <- function(x) { + f.se <- x$f.se + n <- length(x$f) + if (length(f.se) < n) f.se <- c(f.se, NA) + sigma <- x$sigma + if (length(sigma) < n) sigma <- c(sigma, NA) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f.seregr <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + if (length(f.seregr) < n) f.seregr <- + c(f.seregr, rep(NA, n - length(f.seregr))) + notregr <- f.se != f.seregr + notregr[is.na(notregr)] <- TRUE + # set display order + legend_captions <- + c("lm", "log-lin", "calc", "Mack", "input", "input-x", "NA") + src <- factor(legend_captions, levels = legend_captions) + source <- rep(src[1L], n) + # Currently, f.se (o/t tail) can get a value in two ways: + # As a result of the regression (notregr = FALSE) or as a result of + # a calculation based on an imputed sigma and f (notregr = TRUE). + source[notregr] <- src[3L] + # tail.se can get a value in two ways: input or estimated (log-linear) + # input: If argument tail=FALSE, then the user does not want a tail + # in the model, so although tail.se was provided, it will + # be ignored -- "input-x" + # estimated: + # It will be estimated when the tail > 1.000. + # Here, we cannot take the tail from x$tail b/c + # in MackChainLadder when argument tail = TRUE the output value + # named "tail" is the 'lm' object used to estimate it (the tail). + # Therefore, take actual tail value from last f to see if > unity. + if (is.null(x$tail.se.input)) { # tail.se not provided + if (x$f[n] > 1) source[n] <- src[2L] # tail.se est'd via log-linear + else source[n] <- src[7L] + } else if (!x$tail.input) source[n] <- src[6L] # tail not est'd so + # input tail.se ignored + else source[n] <- src[5L] # tail estimated therefore input tail.se used + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + f.seNoNAs <- f.se + + f.seNoNAs[is.na(f.se)] <- 0 + df <- data.frame(xx, f.se, f.seNoNAs, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f.se, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + geom_line(aes(colour = source, group = 1), na.rm = TRUE) + + geom_point(aes(y=f.seNoNAs)) + + ggtitle("f.se estimates") + P +} +plot.mackcl.f1.cv <- function(x) { + f <- x$f + n <- length(f) + f.se <- x$f.se + if (length(f.se) < n) f.se <- c(f.se, NA) + sigma <- x$sigma + if (length(sigma) < n) sigma <- c(sigma, NA) + f1.cv <- f.se/ (f - 1) + smmry <- suppressWarnings(lapply(x$Models, summary)) + f.seregr <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + if (length(f.seregr) < n) f.seregr <- + c(f.seregr, rep(NA, n - length(f.seregr))) + notregr <- f.se != f.seregr + notregr[is.na(notregr)] <- TRUE + # set display order + legend_captions <- c("calc", "NA") + src <- factor(legend_captions, levels = legend_captions) + # Currently, the cv is always the result of a calculation. + # Only exception (NA) is in the tail when argument tail=FALSE + source <- rep(src[1L], n) + if (is.logical(x$tail.input)) if (!x$tail.input) source[n] <- src[2L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + f1.cvNoNAs <- f1.cv + f1.cvNoNAs[is.na(f.se)] <- 0 + df <- data.frame(xx, f1.cv, f1.cvNoNAs, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f1.cv, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + ylab("cv(f-1)") + + geom_line(aes(colour = source, group = 1), na.rm = TRUE) + + geom_point(aes(y=f1.cvNoNAs)) + + ggtitle("cv(f-1) estimates") + P +} +plot.mackcl.f.cv <- function(x) { + f <- x$f + n <- length(f) + f.se <- x$f.se + if (length(f.se) < n) f.se <- c(f.se, NA) + sigma <- x$sigma + if (length(sigma) < n) sigma <- c(sigma, NA) + f.cv <- f.se/ f + smmry <- suppressWarnings(lapply(x$Models, summary)) + f.seregr <- sapply(smmry, function(x) x$coef["x","Std. Error"]) + if (length(f.seregr) < n) f.seregr <- + c(f.seregr, rep(NA, n - length(f.seregr))) + notregr <- f.se != f.seregr + notregr[is.na(notregr)] <- TRUE + + # set display order + legend_captions <- c("calc", "NA") + src <- factor(legend_captions, levels = legend_captions) + # Currently, the cv is always the result of a calculation. + # Only exception (NA) is in the tail when argument tail=FALSE + source <- rep(src[1L], n) + if (is.logical(x$tail.input)) if (!x$tail.input) source[n] <- src[2L] + xx <- factor(colnames(x$Triangle)[1:n], levels = colnames(x$Triangle)[1:n]) + f.cvNoNAs <- f.cv + f.cvNoNAs[is.na(f.se)] <- 0 + df <- data.frame(xx, f.cv, f.cvNoNAs, + source, + stringsAsFactors = FALSE) + P <- ggplot(df, aes(x = xx, y = f.cv, colour = source)) + + xlab(names(dimnames(x$Triangle))[2L]) + + ylab("cv(f)") + + geom_line(aes(colour = source, group = 1), na.rm = TRUE) + + geom_point(aes(y=f.cvNoNAs)) + + ggtitle("cv(f) estimates") + P +} diff --git a/man/plotParms.Rd b/man/plotParms.Rd new file mode 100644 index 00000000..8b76c7db --- /dev/null +++ b/man/plotParms.Rd @@ -0,0 +1,64 @@ +\name{plotParms} +\alias{plotParms} +\alias{plotParms.MackChainLadder} +\alias{plotParms.ChainLadder} +\alias{plotParms.default} +\title{ +Plot the estimated parameters of a model +} + +\encoding{UTF-8} + +\description{ +Generic function for visualizing the estimated parameters +of a model in the ChainLadder package. +S3 methods currently exist for objects of class +ChainLadder and MackChainLadder. +} + +\usage{ +plotParms(x, which, ncol, nrow, title, \dots) +\method{plotParms}{ChainLadder}(x, which = c(1, 3, 4, 6), ncol, nrow, title, \dots) +\method{plotParms}{MackChainLadder}(x, which = c(1, 3, 4, 6), ncol, nrow, title, \dots) +\method{plotParms}{default}(x, \dots) + +} +%- maybe also 'usage' for other objects documented here. +\arguments{ + + \item{x}{object whose parameters are to be visualized} + \item{which}{if a subset of the plots is required, + specify a subset of the numbers 1:6} + \item{ncol}{number of columns in the resulting plot} + \item{nrow}{number of rows in the resulting plot} + \item{title}{optional; character holding title of the plot; +defaults to method call. + } + \item{\dots}{arguments passed to other methods} +} +\details{ +Currently six plots are available for +\code{chainladder} and \code{MackChainLadder}: +\enumerate{ + \item {estimated link ratios f} + \item {estimated link ratios f less unity} + \item {estimated sigma values} + \item {estimated standard error of the link ratio f.se} + \item {coefficient of variation of the link ratio} + \item {coefficient of variation of the link ratio less unity} +} +} +\value{ +In the case of class ChainLadder or MackChainLadder, +list of class \code{arrangelist} returned by +\code{gridExtra::marrangeGrob}. +In all other cases, a \code{stop} error is executed at this time. +} +\author{ +dmm +} +\examples{ +plotParms(chainladder(GenIns)) +plotParms(MackChainLadder(GenIns)) +} +\keyword{ models } diff --git a/vignettes/plotParms.Rmd b/vignettes/plotParms.Rmd new file mode 100644 index 00000000..9faefb94 --- /dev/null +++ b/vignettes/plotParms.Rmd @@ -0,0 +1,229 @@ +--- +title: "plotParms" +author: "Dan Murphy" +date: "March 12, 2017" +output: + pdf_document: default + html_document: default +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +## Description + +*plotParms* is a generic^[S3 generic] +function in the *ChainLadder* package +that plots key parameters estimated by models in the package. +The goal is to -- at a glance -- get information about the +maginitude and potential variability of the parameter, +as well as its source: +estimated, specified, NA, etc. + +Currently two models are supported: + +* chainladder +* MackChainLadder + +For those models, plots of six parameters +are available to be rendered by *ggplot2*: + +1. estimated link ratios **f**, +including one-standard-error "error bars" +1. estimated link ratios **f** less unity, +including one-standard-error "error bars" +1. estimated **sigma** values +1. estimated standard error of the link ratios **f.se** +1. coefficient of variation of the link ratios +1. coefficient of variation of the link ratios less unity + +By default, +plots 1, 3, 4, and 6 are displayed in a 2x2 grid +(using *gridExtra::marrangeGrob*). +Some display flexibility is enabled; +see **plotParms Arguments** section below. + +## Model-Specific Details + +### chainladder + +*chainladder* uses *lm* to estimate +average age-to-age factors for the +given actuarial *Triangle* argument. +Therefore when a development period has only one observation -- +as can often occur in the tail -- +there is insufficient data for *lm* to estimate the **sigma** +parameter of that period's regression model. +When the **sigma** estimate +is technically **NA**, +*plotParms* indicates such with: + +* a broken line at that point +* a point value at zero (for want of a better choice) +* a color-coded legend entry indicating +that the source of the parameter is **NA** + +For example, Figure 1 below shows the default *plotParms* +of the *chainladder* call for the **GenIns** triangle. + +```{r echo=FALSE} +suppressPackageStartupMessages(library(ChainLadder)) +``` +```{r, fig.cap="chainladder(GenIns)"} +plotParms(chainladder(GenIns)) +``` + +Both the *chainladder* function and +*plotParms* borrow the *f* notation of Mack^[Thomas Mack. +Distribution-free calculation of the standard error of +chain ladder reserve estimates. +Astin Bulletin. Vol. 23. No 2. 1993. pp.213:225.]. +The legend indicates the source of the parameter's value: + +* lm: produced by the underlying regression +* NA: not available +* calc: the coefficient of variation is a calculated value +(cv(f) = f/f.se) or **NA** + +### MackChainLadder + +*MackChainLadder* uses +*chainladder* to estimate most parameters, +and extends the latter's capabilities in two situations: + +* sigma = NA +* a tail + +#### sigma = NA + +When *Triangle* has insufficient data to estimate **sigma** +for a given development period, +*MackChainLadder* has two methods for estimating ("imputing") a value: + +* log-linear regression (the default) +* the method of *Mack* + +*plotParms*' legend indicates which method is used. + +#### tail + +Suppose your triangle is 10x10. +In that case there are (typically) nine link ratio factors +to select, +and *chainladder* will generate nine such averages +depending on the value of arguments **weights** and **delta**. +In contrast, +*MackChainLadder* generates 10 factors because +*that method always includes a tail factor* as follows: + +* By default the tail value is 1.000, or +* The user can provide a tail (e.g., with **tail = 1.10**), or +* The user can ask *MackChainLadder* to estimate a tail +by specifying **tail = TRUE**. + +The **sigma** and **se** parameters for the tail are determined by +*MackChainLadder* arguments +**tail.sigma** and **tail.se**. +If the user does not explicitly provide those arguments +*MackChainLadder* will estimate their values. +(Refer to the *MackChainLadder* documentation for more information.) +*plotParms* will graph their value and their source. + +For example, +Figure 2 below is the default plot of the parameters of the +*MackChainLadder* model for the **GenIns** triangle. +Ten factors are plotted versus only nine in Figure 1. +The **NA** value of **sigma** for the ninth link ratio +in Figure 1 is replaced with a value +imputed via log-linear regression (the default). +The associated **f.se** +and **cv** values may now be calculated and are also no longer **NA**. +The default unity tail value +(the tenth factor) +is shown but, +since a unity value is assumed to be deterministic, +**tail.sigma** and **tail.se** values are not available. + +```{r, fig.cap="MackChainLadder(GenIns)"} +plotParms(MackChainLadder(GenIns)) +``` + +## plotParms Arguments + +*plotParms* for the two methods in this document is invoked as follows: + +``` +plotParms(x, which, ncol, nrow, title) +``` + +* x + + a *chainladder* or *MackChainLadder* object + + no default +* which + + a subset of integers 1 through 6 + + default: 1, 3, 4, 6 -- see list in **Description** section above + + follows in the tradition of *ChainLadder::plot.glmReserve* + by Wayne Zhang as well as *stats::plot.lm* +* ncol + + number of columns in the multi-plot output + + default: either 2, or 1 if only one plot is requested +* nrow + + number of rows in the multi-plot output + + default: sufficient number to display all plots on one page + given the number of plots selected and the **ncol** value +* title + + overall title to display at the top of the output page + + default: method's call + +## Final Example + +In this last example +we implement the *Mack Method* on the **GenIns** triangle +by specifying two key arguments: + +* mse.method = "Mack" (the default, but here specified +explicitly) +tells *MackChainLadder*'s recursion steps to duplicate Mack's +closed form formula +* est.sigma = "Mack" (not the default) +tells *MackChainLadder* to use Mack's +heuristic formula for imputing **sigma** values +when not available from a *chainladder* regression + +We also ask for a tail and associated uncertainty parameters +to be estimated via + +* tail = TRUE, +* tail.sigma = NULL, and +* tail.se = NULL. + +Finally, +we tell *plotParms* to display all six plots +and specify our own more succinct title +because the call +is too wide to fit in the display. + +```{r, fig.cap="Mack Method on GenIns"} +x <- MackChainLadder(GenIns, est.sigma = "Mack", mse.method = "Mack", + tail = TRUE, tail.sigma = NULL, tail.se = NULL) +plotParms(x, which = 1:6, title = "Mack Method on GenIns: Estimated Parameters") +``` + +Here all **cv** parameter values +are calculated. + +Also of note: +by comparing **cv(f-1)~9~** points in Figures 2 and 3, +the difference between +**mse.method = "log-linear"** (Figure 2) +and **mse.method = "Mack"** (Figure 3) +does not appear to be significant for this triangle. + +## Conclusion + +Let the author(s) know if you have any problems with plotParms +or wish to see an implementation +for another ChainLadder model. + +Cheers. \ No newline at end of file