Gamma distribution in alm() (issue #13). Note how it is implemented i…

…n the function (see vignette)!
config-i1 · Jan 8, 2021 · 16eb934 · 16eb934
1 parent facc9ed
commit 16eb934
Show file tree

Hide file tree

Showing 10 changed files with 135 additions and 47 deletions.
diff --git a/CRAN-RELEASE b/CRAN-RELEASE
@@ -1,2 +1,2 @@
-This package was submitted to CRAN on 2021-01-04.
-Once it is accepted, delete this file and tag the release (commit e019feb).
+This package was submitted to CRAN on 2021-01-08.
+Once it is accepted, delete this file and tag the release (commit facc9ed).
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: greybox
 Type: Package
 Title: Toolbox for Model Building and Forecasting
-Version: 0.6.6
-Date: 2021-01-07
+Version: 0.6.7.41001
+Date: 2021-01-08
 Authors@R: c(person("Ivan", "Svetunkov", email = "[email protected]", role = c("aut", "cre"),
                   comment="Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK"),
              person("Yves R.", "Sagaert", role = c("ctb"),

diff --git a/NAMESPACE b/NAMESPACE
@@ -237,6 +237,7 @@ importFrom(stats,cov)
 importFrom(stats,dbeta)
 importFrom(stats,dchisq)
 importFrom(stats,deltat)
+importFrom(stats,dgamma)
 importFrom(stats,dlnorm)
 importFrom(stats,dlogis)
 importFrom(stats,dnbinom)

diff --git a/NEWS b/NEWS
@@ -1,3 +1,10 @@
+greybox v0.6.7 (Release data: 2021-01-08)
+==============
+
+Changes:
+* Gamma distribution in alm().
+
+
 greybox v0.6.6 (Release data: 2021-01-07)
 ==============
 

diff --git a/R/alm.R b/R/alm.R
@@ -20,6 +20,7 @@
 #' \item \link[greybox]{dbcnorm} - Box-Cox normal distribution,
 # \item \link[stats]{dchisq} - Chi-Squared Distribution,
 #' \item \link[statmod]{dinvgauss} - Inverse Gaussian distribution,
+#' \item \link[stats]{dgamma} - Gamma distribution,
 #' \item \link[greybox]{dlogitnorm} - Logit-normal distribution,
 #' \item \link[stats]{dbeta} - Beta distribution,
 #' \item \link[stats]{dpois} - Poisson Distribution,
@@ -244,14 +245,15 @@
 #' @importFrom pracma hessian
 #' @importFrom nloptr nloptr
 #' @importFrom stats model.frame sd terms model.matrix
-#' @importFrom stats dchisq dlnorm dnorm dlogis dpois dnbinom dt dbeta
+#' @importFrom stats dchisq dlnorm dnorm dlogis dpois dnbinom dt dbeta dgamma
 #' @importFrom stats plogis
 #' @importFrom statmod dinvgauss
 #' @importFrom forecast Arima
 #' @export alm
 alm <- function(formula, data, subset, na.action,
                 distribution=c("dnorm","dlaplace","ds","dgnorm","dlogis","dt","dalaplace",
-                               "dlnorm","dllaplace","dls","dlgnorm","dbcnorm","dfnorm","dinvgauss",
+                               "dlnorm","dllaplace","dls","dlgnorm","dbcnorm","dfnorm",
+                               "dinvgauss","dgamma",
                                "dpois","dnbinom",
                                "dbeta","dlogitnorm",
                                "plogis","pnorm"),
@@ -411,6 +413,7 @@ alm <- function(formula, data, subset, na.action,
 
         mu[] <- switch(distribution,
                        "dinvgauss"=,
+                       "dgamma"=,
                        "dpois" =,
                        "dnbinom" = exp(matrixXreg %*% B),
                        "dchisq" = ifelseFast(any(matrixXreg %*% B <0),1E+100,(matrixXreg %*% B)^2),
@@ -447,6 +450,7 @@ alm <- function(formula, data, subset, na.action,
                         "dlgnorm" = (other*sum(abs(log(y[otU])-mu[otU])^other)/obsInsample)^{1/other},
                         "dbcnorm" = sqrt(sum((bcTransform(y[otU],other)-mu[otU])^2)/obsInsample),
                         "dinvgauss" = sum((y[otU]/mu[otU]-1)^2 / (y[otU]/mu[otU]))/obsInsample,
+                        "dgamma" = sum((y[otU]/mu[otU]-1)^2)/obsInsample,
                         "dlogitnorm" = sqrt(sum((log(y[otU]/(1-y[otU]))-mu[otU])^2)/obsInsample),
                         "dfnorm" = abs(other),
                         "dt" = ,
@@ -512,6 +516,8 @@ alm <- function(formula, data, subset, na.action,
                                    "dfnorm" = dfnorm(y[otU], mu=fitterReturn$mu[otU], sigma=fitterReturn$scale, log=TRUE),
                                    "dinvgauss" = dinvgauss(y[otU], mean=fitterReturn$mu[otU],
                                                            dispersion=fitterReturn$scale/fitterReturn$mu[otU], log=TRUE),
+                                   "dgamma" = dgamma(y[otU], shape=1/fitterReturn$scale,
+                                                     scale=fitterReturn$scale*fitterReturn$mu[otU], log=TRUE),
                                    "dchisq" = dchisq(y[otU], df=fitterReturn$scale, ncp=fitterReturn$mu[otU], log=TRUE),
                                    "dpois" = dpois(y[otU], lambda=fitterReturn$mu[otU], log=TRUE),
                                    "dnbinom" = dnbinom(y[otU], mu=fitterReturn$mu[otU], size=fitterReturn$scale, log=TRUE),
@@ -538,6 +544,9 @@ alm <- function(formula, data, subset, na.action,
                                               # "dinvgauss" = 0.5*(obsZero*(log(pi/2)+1+suppressWarnings(log(fitterReturn$scale)))-
                                               #                                 sum(log(fitterReturn$mu[!otU]))),
                                               "dinvgauss" = obsZero*(0.5*(log(pi/2)+1+suppressWarnings(log(fitterReturn$scale)))),
+                                              "dgamma" = obsZero*(1/fitterReturn$scale + log(fitterReturn$scale) +
+                                                                  log(gamma(1/fitterReturn$scale)) +
+                                                                  (1-1/fitterReturn$scale)*digamma(1/fitterReturn$scale)),
                                               "dlaplace" =,
                                               "dllaplace" =,
                                               "ds" =,
@@ -572,6 +581,7 @@ alm <- function(formula, data, subset, na.action,
                                 "dnorm" =,
                                 "dgnorm" =,
                                 "dinvgauss" =,
+                                "dgamma" =,
                                 "dlaplace" =,
                                 "dalaplace" =,
                                 "dlogis" =,
@@ -1013,12 +1023,13 @@ alm <- function(formula, data, subset, na.action,
     errors <- vector("numeric", obsInsample);
     ot <- vector("logical", obsInsample);
 
-    if(any(y<0) & any(distribution==c("dfnorm","dlnorm","dllaplace","dls","dbcnorm","dinvgauss","dchisq","dpois","dnbinom"))){
+    if(any(y<0) & any(distribution==c("dfnorm","dlnorm","dllaplace","dls","dbcnorm","dchisq","dpois","dnbinom",
+                                      "dinvgauss","dgamma"))){
         stop(paste0("Negative values are not allowed in the response variable for the distribution '",distribution,"'"),
              call.=FALSE);
     }
 
-    if(any(y==0) & any(distribution==c("dinvgauss")) & !occurrenceModel){
+    if(any(y==0) & any(distribution==c("dinvgauss","dgamma")) & !occurrenceModel){
         stop(paste0("Zero values are not allowed in the response variable for the distribution '",distribution,"'"),
              call.=FALSE);
     }
@@ -1262,7 +1273,8 @@ alm <- function(formula, data, subset, na.action,
         if(is.null(B)){
             #### I(0) initialisation ####
             if(iOrder==0){
-                if(any(distribution==c("dlnorm","dllaplace","dls","dlgnorm","dpois","dnbinom","dinvgauss"))){
+                if(any(distribution==c("dlnorm","dllaplace","dls","dlgnorm","dpois","dnbinom",
+                                       "dinvgauss","dgamma"))){
                     if(any(y[otU]==0)){
                         # Use Box-Cox if there are zeroes
                         B <- .lm.fit(matrixXreg[otU,,drop=FALSE],bcTransform(y[otU],0.01))$coefficients;
@@ -1340,7 +1352,8 @@ alm <- function(formula, data, subset, na.action,
                     matrixXregForDiffs <- matrixXregForDiffs[-c(1:iOrder),,drop=FALSE];
                 }
 
-                if(any(distribution==c("dlnorm","dllaplace","dls","dlgnorm","dpois","dnbinom","dinvgauss"))){
+                if(any(distribution==c("dlnorm","dllaplace","dls","dlgnorm","dpois","dnbinom",
+                                       "dinvgauss","dgamma"))){
                     B <- .lm.fit(matrixXregForDiffs,diff(log(y[otU]),differences=iOrder))$coefficients;
                 }
                 else if(any(distribution==c("plogis","pnorm"))){
@@ -1665,6 +1678,7 @@ alm <- function(formula, data, subset, na.action,
                        "dnorm" =,
                        "dgnorm" =,
                        "dinvgauss" =,
+                       "dgamma" =,
                        "dlaplace" =,
                        "dalaplace" =,
                        "dlogis" =,
@@ -1697,7 +1711,8 @@ alm <- function(formula, data, subset, na.action,
                        "dgnorm" =,
                        "dnbinom" =,
                        "dpois" = y - mu,
-                       "dinvgauss" = y / mu,
+                       "dinvgauss" =,
+                       "dgamma" = y / mu,
                        "dchisq" = sqrt(y) - sqrt(mu),
                        "dlnorm" =,
                        "dllaplace" =,

diff --git a/R/methods.R b/R/methods.R
@@ -243,6 +243,7 @@ pointLik.alm <- function(object, ...){
                             "dbcnorm" = dbcnorm(y, mu=mu, sigma=scale, lambda=object$other$lambdaBC, log=TRUE),
                             "dlogitnorm" = dlogitnorm(y, mu=mu, sigma=scale, log=TRUE),
                             "dinvgauss" = dinvgauss(y, mean=mu, dispersion=scale/mu, log=TRUE),
+                            "dgamma" = dgamma(y, shape=1/scale, scale=scale*mu, log=TRUE),
                             "dlaplace" = dlaplace(y, mu=mu, scale=scale, log=TRUE),
                             "dllaplace" = dlaplace(log(y), mu=mu, scale=scale, log=TRUE),
                             "dalaplace" = dalaplace(y, mu=mu, scale=scale, alpha=object$other$alpha, log=TRUE),
@@ -275,6 +276,7 @@ pointLik.alm <- function(object, ...){
                                    "dlgnorm" = 1/object$other$shape -
                                        log(object$other$shape / (2*scale*gamma(1/object$other$shape))),
                                    "dinvgauss" = 0.5*(log(pi/2)+1+log(scale)),
+                                   "dgamma" = 1/scale + log(scale) + log(gamma(1/scale)) + (1-1/scale)*digamma(1/scale),
                                    "dlaplace" =,
                                    "dllaplace" =,
                                    "dalaplace" = (1 + log(2*scale)),
@@ -1133,8 +1135,8 @@ plot.greybox <- function(x, which=c(1,2,4,6), level=0.95, legend=FALSE,
         outliers <- outlierdummy(x, level=level, type=type);
         outliersID <- outliers$id;
         statistic <- outliers$statistic;
-        # Analyse stuff in logarithms if the distribution is dinvgauss
-        if(x$distribution=="dinvgauss"){
+        # Analyse stuff in logarithms if the distribution is dinvgauss / dgamma
+        if(any(x$distribution==c("dinvgauss","dgamma"))){
             ellipsis$y[] <- log(ellipsis$y);
             statistic[] <- log(statistic);
         }
@@ -1197,7 +1199,7 @@ plot.greybox <- function(x, which=c(1,2,4,6), level=0.95, legend=FALSE,
 
         ellipsis$x <- as.vector(fitted(x));
         ellipsis$y <- as.vector(residuals(x));
-        if(x$distribution=="dinvgauss"){
+        if(any(x$distribution==c("dinvgauss","dgamma"))){
             ellipsis$y[] <- log(ellipsis$y);
         }
         if(type=="abs"){
@@ -1332,6 +1334,17 @@ plot.greybox <- function(x, which=c(1,2,4,6), level=0.95, legend=FALSE,
             do.call(qqplot, ellipsis);
             qqline(ellipsis$y, distribution=function(p) qinvgauss(p, mean=1, dispersion=x$scale));
         }
+        else if(x$distribution=="dgamma"){
+            # Transform residuals for something meaningful
+            # This is not 100% accurate, because the scale should change together with mean...
+            if(!any(names(ellipsis)=="main")){
+                ellipsis$main <- "QQ-plot of Gamma distribution";
+            }
+            ellipsis$x <- qgamma(ppoints(500), shape=1/x$scale, scale=x$scale);
+
+            do.call(qqplot, ellipsis);
+            qqline(ellipsis$y, distribution=function(p) qgamma(p, shape=1/x$scale, scale=x$scale));
+        }
         else if(x$distribution=="dchisq"){
             message("Sorry, but we don't produce QQ plots for the Chi-Squared distribution");
         }
@@ -1434,7 +1447,7 @@ plot.greybox <- function(x, which=c(1,2,4,6), level=0.95, legend=FALSE,
         outliersID <- outliers$id;
         statistic <- outliers$statistic;
         # Analyse stuff in logarithms if the distribution is dinvgauss
-        if(x$distribution=="dinvgauss"){
+        if(any(x$distribution==c("dinvgauss","dgamma"))){
             ellipsis$x[] <- log(ellipsis$x);
             statistic[] <- log(statistic);
         }
@@ -1904,6 +1917,7 @@ print.summary.alm <- function(x, ...){
                       "dbcnorm" = paste0("Box-Cox Normal with lambda=",round(x$other$lambdaBC,digits)),
                       "dlogitnorm" = "Logit Normal",
                       "dinvgauss" = "Inverse Gaussian",
+                      "dgamma" = "Gamma",
                       "dchisq" = paste0("Chi-Squared with nu=",round(x$other$nu,digits)),
                       "dpois" = "Poisson",
                       "dnbinom" = paste0("Negative Binomial with size=",round(x$other$size,digits)),
@@ -1975,6 +1989,7 @@ print.summary.greybox <- function(x, ...){
                       "dbcnorm" = "Box-Cox Normal",
                       "dlogitnorm" = "Logit Normal",
                       "dinvgauss" = "Inverse Gaussian",
+                      "dgamma" = "Gamma",
                       "dchisq" = "Chi-Squared",
                       "dpois" = "Poisson",
                       "dnbinom" = "Negative Binomial",
@@ -2028,6 +2043,7 @@ print.summary.greyboxC <- function(x, ...){
                       "dbcnorm" = "Box-Cox Normal",
                       "dlogitnorm" = "Logit Normal",
                       "dinvgauss" = "Inverse Gaussian",
+                      "dgamma" = "Gamma",
                       "dchisq" = "Chi-Squared",
                       "dpois" = "Poisson",
                       "dnbinom" = "Negative Binomial",
@@ -2153,7 +2169,7 @@ rstandard.greybox <- function(model, ...){
         errors[residsToGo] <- ((errors[residsToGo] - mean(errors[residsToGo])) /
                          (model$scale^model$other$shape * obs / df)^{1/model$other$shape});
     }
-    else if(model$distribution=="dinvgauss"){
+    else if(any(model$distribution==c("dinvgauss","dgamma"))){
         errors[residsToGo] <- errors[residsToGo] / mean(errors[residsToGo]);
     }
     else{
@@ -2228,7 +2244,7 @@ rstudent.greybox <- function(model, ...){
             rstudentised[i] <- errors[i] / (sqrt(sum(errors[-i]^2) / df) * sqrt(3) / pi);
         }
     }
-    else if(model$distribution=="dinvgauss"){
+    else if(any(model$distribution==c("dinvgauss","dgamma"))){
         for(i in which(residsToGo)){
             rstudentised[i] <- errors[i] / mean(errors[-i]);
         }
@@ -2344,6 +2360,7 @@ outlierdummy.alm <- function(object, level=0.999, type=c("rstandard","rstudent")
                       "dinvgauss"=qinvgauss(c((1-level)/2, (1+level)/2), mean=1,
                                             dispersion=object$scale * nobs(object) /
                                                 (nobs(object)-nparam(object))),
+                      "dgamma"=qgamma(c((1-level)/2, (1+level)/2), shape=1/object$scale, scale=object$scale),
                       qnorm(c((1-level)/2, (1+level)/2), 0, 1));
     outliersID <- which(errors>statistic[2] | errors<statistic[1]);
     outliersNumber <- length(outliersID);
@@ -2597,7 +2614,7 @@ setMethod("extract", signature=className("summary.greyboxC","greybox"), definiti
 #### Predictions and forecasts ####
 
 #' @rdname predict.greybox
-#' @importFrom stats predict qchisq qlnorm qlogis qpois qnbinom qbeta
+#' @importFrom stats predict qchisq qlnorm qlogis qpois qnbinom qbeta qgamma
 #' @importFrom statmod qinvgauss
 #' @export
 predict.alm <- function(object, newdata=NULL, interval=c("none", "confidence", "prediction"),
@@ -2898,6 +2915,20 @@ predict.alm <- function(object, newdata=NULL, interval=c("none", "confidence", "
             greyboxForecast$upper[] <- exp(greyboxForecast$upper);
         }
     }
+    else if(object$distribution=="dgamma"){
+        greyboxForecast$mean <- exp(greyboxForecast$mean);
+        if(interval=="prediction"){
+            greyboxForecast$scale <- object$scale;
+            greyboxForecast$lower[] <- greyboxForecast$mean*qgamma(levelLow, shape=1/greyboxForecast$scale,
+                                                                   scale=greyboxForecast$scale);
+            greyboxForecast$upper[] <- greyboxForecast$mean*qgamma(levelUp, shape=1/greyboxForecast$scale,
+                                                                   scale=greyboxForecast$scale);
+        }
+        else if(interval=="confidence"){
+            greyboxForecast$lower[] <- exp(greyboxForecast$lower);
+            greyboxForecast$upper[] <- exp(greyboxForecast$upper);
+        }
+    }
     else if(object$distribution=="dlogis"){
         # Use the connection between the variance and scale in logistic distribution
         if(interval!="none"){

diff --git a/R/stepwise.R b/R/stepwise.R
@@ -66,7 +66,8 @@
 stepwise <- function(data, ic=c("AICc","AIC","BIC","BICc"), silent=TRUE, df=NULL,
                      method=c("pearson","kendall","spearman"),
                      distribution=c("dnorm","dlaplace","ds","dgnorm","dlogis","dt","dalaplace",
-                                    "dfnorm","dlnorm","dllaplace","dls","dbcnorm","dinvgauss",
+                                    "dfnorm","dlnorm","dllaplace","dls","dbcnorm",
+                                    "dinvgauss","dgamma",
                                     "dpois","dnbinom",
                                     "plogis","pnorm"),
                      occurrence=c("none","plogis","pnorm"), ...){

diff --git a/man/alm.Rd b/man/alm.Rd
diff --git a/man/stepwise.Rd b/man/stepwise.Rd