diff --git a/man/HarzTraffic_files/figure-html/unnamed-chunk-1-4.png b/man/HarzTraffic_files/figure-html/unnamed-chunk-1-4.png index c522f64..b2f0dda 100644 Binary files a/man/HarzTraffic_files/figure-html/unnamed-chunk-1-4.png and b/man/HarzTraffic_files/figure-html/unnamed-chunk-1-4.png differ diff --git a/man/gamlss2.html b/man/gamlss2.html index d2192fd..fed9a32 100644 --- a/man/gamlss2.html +++ b/man/gamlss2.html @@ -635,7 +635,7 @@

Examples

*-------- n = 610 df = 13.12 res.df = 596.88 Deviance = 4770.1554 Null Dev. Red. = 33.39% -AIC = 4796.3966 elapsed = 0.77sec +AIC = 4796.3966 elapsed = 0.80sec 
## plot estimated effects
 plot(b, which = "effects")
diff --git a/man/prodist.gamlss2.html b/man/prodist.gamlss2.html index 7c69c05..1c31f79 100644 --- a/man/prodist.gamlss2.html +++ b/man/prodist.gamlss2.html @@ -543,10 +543,10 @@

Examples

## simulate random numbers
 random(d, 5)
-
       r_1      r_2      r_3      r_4      r_5
-1 22.07090 29.92421 18.67442 29.23931 13.09939
-2 19.52497 41.95755 49.88034 28.46416 79.20217
-3 82.29180 72.53890 83.23728 76.62891 74.27616
+
        r_1       r_2       r_3       r_4      r_5
+1  31.99176  24.26454  29.24719  38.53560 22.14378
+2  32.35560  57.28431  41.84400  55.37877 53.11516
+3 100.10360 161.57604 152.17778 150.07732 94.09042
## density and distribution
 pdf(d, 50 * -2:2)
diff --git a/search.json b/search.json index 31dd305..331898a 100644 --- a/search.json +++ b/search.json @@ -69,7 +69,7 @@ "href": "vignettes/spatial.html#estimating-the-model", "title": "Spatial Effects", "section": "5 Estimating the Model", - "text": "5 Estimating the Model\nWe now estimate the spatial count model using the Negative Binomial distribution (NBI). The model includes smooth functions of altitude, year, and an interaction between altitude and year, as well as a (functional) random effect over regions for the time trend. Spatial effects are incorporated using the bs = \"mrf\" option.\n\n## estimate count model using the NBI family,\n## model formula is\nf <- ~ s(log(alt + 10)) + s(year) + s(id, bs = \"mrf\", xt = list(\"penalty\" = K), k = 100)\nf <- list(update(f, counts ~ .), f)\n\n## estimate model using BIC for shrinkage parameter selection\nb <- gamlss2(f, data = storms, family = NBI, criterion = \"BIC\")\n\nNote that estimation of this model takes some time, because of the quite high number of coefficients for the (functional) random effect ti(id, year, bs = c(\"re\", \"cr\"), k = c(157, 5)).\n\n## model summary\nsummary(b)\n\nCall:\ngamlss2(formula = counts ~ s(log(alt + 10)) + s(year) + s(id, \n bs = \"mrf\", xt = list(penalty = K), k = 100) | s(log(alt + \n 10)) + s(year) + s(id, bs = \"mrf\", xt = list(penalty = K), \n k = 100), data = storms, family = NBI, ... = pairlist(x = FALSE, criterion = \"BIC\"))\n---\nFamily: NBI \nLink function: mu = log, sigma = log\n*--------\nParameter: mu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 1.2999 0.0167 77.86 <2e-16 ***\n---\nSmooth terms:\n s(log(alt + 10)) s(year) s(id)\nedf 7.5077 2.2994 52.88\n*--------\nParameter: sigma \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.20330 0.02505 -8.117 6.6e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n---\nSmooth terms:\n s(log(alt + 10)) s(year) s(id)\nedf 7.4533 2.1207 24.712\n*--------\nn = 3494 df = 98.97 res.df = 3395.03\nDeviance = 16484.181 Null Dev. Red. = 14.46%\nAIC = 16682.1271 elapsed = 3.88sec\n\n\nModel calibration is checked using histogram, Q-Q plot, wormplot etc.\n\nplot(b, which = \"resid\")\n\n\n\n\n\n\n\n\nHere, the spatial effect is modeled as an MRF smooth (s(id, bs = \"mrf\")), where the penalty matrix K enforces spatial structure based on neighboring stations. We use the Bayesian Information Criterion (BIC) to select the optimal smoothing parameters.", + "text": "5 Estimating the Model\nWe now estimate the spatial count model using the Negative Binomial distribution (NBI). The model includes smooth functions of altitude, year, and an interaction between altitude and year, as well as a (functional) random effect over regions for the time trend. Spatial effects are incorporated using the bs = \"mrf\" option.\n\n## estimate count model using the NBI family,\n## model formula is\nf <- ~ s(log(alt + 10)) + s(year) + s(id, bs = \"mrf\", xt = list(\"penalty\" = K), k = 100)\nf <- list(update(f, counts ~ .), f)\n\n## estimate model using BIC for shrinkage parameter selection\nb <- gamlss2(f, data = storms, family = NBI, criterion = \"BIC\")\n\nNote that estimation of this model takes some time, because of the quite high number of coefficients for the (functional) random effect ti(id, year, bs = c(\"re\", \"cr\"), k = c(157, 5)).\n\n## model summary\nsummary(b)\n\nCall:\ngamlss2(formula = counts ~ s(log(alt + 10)) + s(year) + s(id, \n bs = \"mrf\", xt = list(penalty = K), k = 100) | s(log(alt + \n 10)) + s(year) + s(id, bs = \"mrf\", xt = list(penalty = K), \n k = 100), data = storms, family = NBI, ... = pairlist(x = FALSE, criterion = \"BIC\"))\n---\nFamily: NBI \nLink function: mu = log, sigma = log\n*--------\nParameter: mu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 1.2999 0.0167 77.86 <2e-16 ***\n---\nSmooth terms:\n s(log(alt + 10)) s(year) s(id)\nedf 7.5077 2.2994 52.88\n*--------\nParameter: sigma \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.20330 0.02505 -8.117 6.6e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n---\nSmooth terms:\n s(log(alt + 10)) s(year) s(id)\nedf 7.4533 2.1207 24.712\n*--------\nn = 3494 df = 98.97 res.df = 3395.03\nDeviance = 16484.181 Null Dev. Red. = 14.46%\nAIC = 16682.1271 elapsed = 3.92sec\n\n\nModel calibration is checked using histogram, Q-Q plot, wormplot etc.\n\nplot(b, which = \"resid\")\n\n\n\n\n\n\n\n\nHere, the spatial effect is modeled as an MRF smooth (s(id, bs = \"mrf\")), where the penalty matrix K enforces spatial structure based on neighboring stations. We use the Bayesian Information Criterion (BIC) to select the optimal smoothing parameters.", "crumbs": [ "Articles", "Spatial Effects" @@ -168,7 +168,29 @@ "href": "vignettes/survival.html", "title": "Survival Models", "section": "", - "text": "Here is the text for survival model introduction. Important is, how to estimate survival models with the GAMLSS model class.\n\npkg <- c(\"spBayesSurv\", \"gamlss.cens\", \"sf\")\nfor(p in pkg) {\n if(!(p %in% installed.packages())) install.packages(p)\n}\nlibrary(\"gamlss2\")\nlibrary(\"gamlss.cens\")\nlibrary(\"sf\")\ndata(\"LeukSurv\", package = \"spBayesSurv\")\n\nWe need to generate the censored family with\n\ngen.cens(WEI)\n\nA censored family of distributions from WEI has been generated \n and saved under the names: \n dWEIrc pWEIrc qWEIrc WEIrc \nThe type of censoring is right \n\nfam <- cens(WEI)\n\nFirst, we estimate time only models.\n\n## using gamlss2\nb1 <- gamlss2(Surv(time, cens) ~ 1, family = fam, data = LeukSurv)\n\nGAMLSS-RS iteration 1: Global Deviance = 12269.5961 eps = 0.009779 \nGAMLSS-RS iteration 2: Global Deviance = 12267.698 eps = 0.000154 \nGAMLSS-RS iteration 3: Global Deviance = 12267.6904 eps = 0.000000 \n\n## and now with survfit\ns1 <- survfit(Surv(time, cens) ~ 1, data = LeukSurv)\n\nPredict survival curves.\n\n## gamlss2 survival, predict parameters\npar <- predict(b1, newdata = LeukSurv)\n\n## compute survival probabilities\np1 <- 1 - family(b1)$p(Surv(LeukSurv$time, rep(1, nrow(LeukSurv))), par)\n\n## surfit() survival\nf <- stepfun(s1$time, c(1, s1$surv))\np2 <- f(LeukSurv$time)\n\nPlot estimated curves.\n\nmatplot(LeukSurv$time, cbind(p1, p2),\n type = \"l\", lty = 1, lwd = 2, col = 1:2,\n xlab = \"Time\", ylab = \"Survival Prob(t > Time)\",\n main = \"Estimated Survival Curves\")\nlegend(\"topright\", c(\"gamlss2\", \"survfit\"),\n lwd = 2, col = 1:2, bty = \"n\")\n\n\n\n\n\n\n\n\nNow, we estimate a full spatial survival model.\n\n## model formula\nf <- Surv(time, cens) ~ sex + s(age) + s(wbc) + s(tpi) + s(xcoord,ycoord,k=100) |\n sex + s(age) + s(wbc) + s(tpi) + s(xcoord,ycoord,k=100)\n\n## estimate model\nb2 <- gamlss2(f, family = fam, data = LeukSurv)\n\nGAMLSS-RS iteration 1: Global Deviance = 11863.604 eps = 0.042545 \nGAMLSS-RS iteration 2: Global Deviance = 11859.5156 eps = 0.000344 \nGAMLSS-RS iteration 3: Global Deviance = 11859.4914 eps = 0.000002 \n\n\nPlot estimated effects.\n\nplot(b2)\n\n\n\n\n\n\n\n\nResiduals diagnostic plots.\n\nplot(b2, which = \"resid\")\n\n\n\n\n\n\n\n\nVisualize the spatial effect.\n\n## read the map of new west england.\nfile <- system.file(\"otherdata/nwengland.bnd\", package = \"spBayesSurv\")\nd <- readLines(file)\n\n## transform the polygons to a list().\nid <- grep('\\\"', d, fixed = TRUE)\npolys <- list()\nfor(i in 1:length(id)) {\n j <- strsplit(d[id[i]], \",\")[[1]][2]\n if(i < length(id))\n polys[[j]] <- d[(id[i] + 1):(id[i + 1] - 1)]\n else\n polys[[j]] <- d[(id[i] + 1):length(d)]\n}\n\npolys <- lapply(polys, function(x) {\n tf <- tempfile()\n writeLines(x, tf)\n pol <- as.matrix(read.csv(tf, header = FALSE))\n unlink(tf)\n return(st_polygon(list(pol)))\n})\n\n## transform to sf object\npolys_sfc <- st_sfc(polys)\nmap <- st_sf(geometry = polys_sfc)\nmap$id <- names(polys)\nmap$district <- 1:nrow(map)\n\n## plot the map\npar(mar = rep(0, 4))\nplot(st_geometry(map))\n\n## sample coordinates for plotting\nco <- st_sample(map, size = 1000, type = \"regular\")\n\n## create new data for prediction\nnd <- as.data.frame(st_coordinates(co))\nnames(nd) <- c(\"xcoord\", \"ycoord\")\nnd$sex <- 1\nnd$wbc <- mean(LeukSurv$wbc)\nnd$tpi <- mean(LeukSurv$tpi)\nnd$age <- mean(LeukSurv$age)\n\n## predict parameters\npar <- predict(b2, newdata = nd)\n\n## compute survival probabilities\np2 <- 1 - family(b2)$p(Surv(rep(365, nrow(nd)), rep(1, nrow(nd))), par)\n\n## plot on map\nsp <- st_sf(geometry = co)\nsp$survprob <- p2\n\npar(mar = rep(0, 4))\nplot(st_geometry(map))\nplot(sp, pch = 15, cex = 1.3, add = TRUE)\nplot(st_geometry(map), add = TRUE)\n\n\n\n\n\n\n\n\n\n\n\n\nReferences\n\nRigby, R. A., and D. M. Stasinopoulos. 2005. “Generalized Additive Models for Location, Scale and Shape.” Journal of the Royal Statistical Society C 54 (3): 507–54. https://doi.org/10.1111/j.1467-9876.2005.00510.x.", + "text": "In survival analysis, we aim to model the time until an event occurs, which is typically represented as a random variable \\(T\\). The survival function \\(S(t)\\) indicates the probability that the event has not occurred by time \\(t\\), mathematically defined as \\[\nS(t) = P(T > t) = 1 - F(t),\n\\] where \\(F(t)\\) is the cumulative distribution function (CDF) of the random variable \\(T\\).\nThe hazard function \\(\\lambda(t)\\), which describes the instantaneous risk of the event occurring at time \\(t\\), is defined as \\[\n\\lambda(t) = \\frac{d(t)}{S(t)},\n\\] where \\(d(t)\\) is the probability density function (PDF) of \\(T\\).\nWithin the GAMLSS framework, we can specify a parametric form for the hazard function by selecting an appropriate distribution that fits the survival data, denoted as \\[\nd(t \\mid \\boldsymbol{\\theta}),\n\\] where \\(d( \\cdot )\\) represents a parametric density suitable for modeling survival times and \\(\\boldsymbol{\\theta} = (\\theta_1, \\ldots, \\theta_K)^\\top\\) are the parameters that need to be estimated.\nTo incorporate covariates into this model, we can express the parameters as functions of explanatory variables \\(\\theta_k(\\mathbf{x})\\) for \\(k = 1, \\ldots, K\\). We utilize GAM-type predictors represented by \\[\n\\eta_{k}(\\mathbf{x}) = f_1(\\mathbf{x}) + \\dots + f_{L_{k}}(\\mathbf{x}),\n\\] which are linked to the parameters through suitable link functions \\[\nh_{k}(\\theta_{k}(\\mathbf{x})) = \\eta_{k}(\\mathbf{x}).\n\\] The functions \\(f_l(\\cdot)\\) for \\(l = 1, \\ldots, L_{k}\\) can be nonlinear smooth functions (typically estimated using regression splines), linear effects, or random effects, among others. This flexible modeling approach allows for a richer representation of the relationship between covariates and the distribution parameters compared to simple linear effects.\nIn GAMLSS, model estimation is generally performed using (penalized) maximum likelihood estimation (MLE). The likelihood function for a survival model with right-censored data for \\(i = 1, \\ldots, n\\) observations can be expressed as \\[\nL(\\boldsymbol{\\beta}, \\boldsymbol{\\theta}) = \\prod_{i=1}^{n} \\left[ d(t_i \\mid \\boldsymbol{\\theta}(\\mathbf{x}_i)) \\right]^{\\delta_i} \\left[ S(t_i \\mid \\boldsymbol{\\theta}(\\mathbf{x}_i)) \\right]^{1 - \\delta_i},\n\\] where \\(\\delta_i\\) is the censoring indicator for the \\(i\\)-th observation. It takes the value of 1 if the event is observed (i.e., not censored) and 0 if it is censored. This likelihood function effectively accounts for both observed and censored data, facilitating the estimation of model parameters.", + "crumbs": [ + "Articles", + "Survival Models" + ] + }, + { + "objectID": "vignettes/survival.html#leuksurv-data-example", + "href": "vignettes/survival.html#leuksurv-data-example", + "title": "Survival Models", + "section": "1 LeukSurv data example", + "text": "1 LeukSurv data example\n\npkg <- c(\"spBayesSurv\", \"gamlss.cens\", \"sf\", \"stars\")\nfor(p in pkg) {\n if(!(p %in% installed.packages())) install.packages(p)\n}\nlibrary(\"gamlss2\")\nlibrary(\"gamlss.cens\")\nlibrary(\"sf\")\nlibrary(\"stars\")\ndata(\"LeukSurv\", package = \"spBayesSurv\")\n\nWe need to generate the censored family with\n\ngen.cens(WEI)\n\nA censored family of distributions from WEI has been generated \n and saved under the names: \n dWEIrc pWEIrc qWEIrc WEIrc \nThe type of censoring is right \n\nfam <- cens(WEI)\nfam <- fam(mu.link = \"log\")\n\nFirst, we estimate time only models.\n\n## using gamlss2\nb1 <- gamlss2(Surv(time, cens) ~ 1, family = fam, data = LeukSurv)\n\nGAMLSS-RS iteration 1: Global Deviance = 12269.5961 eps = 0.009779 \nGAMLSS-RS iteration 2: Global Deviance = 12267.698 eps = 0.000154 \nGAMLSS-RS iteration 3: Global Deviance = 12267.6904 eps = 0.000000 \n\n## and now with survfit\ns1 <- survfit(Surv(time, cens) ~ 1, data = LeukSurv)\n\nPredict survival curves.\n\n## gamlss2 survival, predict parameters\npar <- predict(b1, newdata = LeukSurv)\n\n## compute survival probabilities\np1 <- 1 - family(b1)$p(Surv(LeukSurv$time, rep(1, nrow(LeukSurv))), par)\n\n## surfit() survival\nf <- stepfun(s1$time, c(1, s1$surv))\np2 <- f(LeukSurv$time)\n\nPlot estimated curves.\n\nmatplot(LeukSurv$time, cbind(p1, p2),\n type = \"l\", lty = 1, lwd = 2, col = 1:2,\n xlab = \"Time\", ylab = \"Survival Prob(t > Time)\",\n main = \"Estimated Survival Curves\")\nlegend(\"topright\", c(\"gamlss2\", \"survfit\"),\n lwd = 2, col = 1:2, bty = \"n\")", + "crumbs": [ + "Articles", + "Survival Models" + ] + }, + { + "objectID": "vignettes/survival.html#spatial-survival-model", + "href": "vignettes/survival.html#spatial-survival-model", + "title": "Survival Models", + "section": "2 Spatial survival model", + "text": "2 Spatial survival model\nNow, we estimate a full spatial survival model.\n\n## model formula\nf <- Surv(time, cens) ~ sex + s(age) + s(wbc) + s(tpi) + s(xcoord, ycoord) |\n sex + s(age) + s(wbc) + s(tpi) + s(xcoord, ycoord)\n\n## estimate model\nb2 <- gamlss2(f, family = fam, data = LeukSurv)\n\nGAMLSS-RS iteration 1: Global Deviance = 11869.8787 eps = 0.042038 \nGAMLSS-RS iteration 2: Global Deviance = 11866.6901 eps = 0.000268 \nGAMLSS-RS iteration 3: Global Deviance = 11866.619 eps = 0.000005 \n\n\nPlot estimated effects.\n\nplot(b2)\n\n\n\n\n\n\n\n\nResiduals diagnostic plots.\n\nplot(b2, which = \"resid\")\n\n\n\n\n\n\n\n\nVisualize the spatial effect.\n\n## read the map of new west england.\nfile <- system.file(\"otherdata/nwengland.bnd\", package = \"spBayesSurv\")\nd <- readLines(file)\n\n## transform the polygons to a list().\nid <- grep('\\\"', d, fixed = TRUE)\npolys <- list()\nfor(i in 1:length(id)) {\n j <- strsplit(d[id[i]], \",\")[[1]][2]\n if(i < length(id))\n polys[[j]] <- d[(id[i] + 1):(id[i + 1] - 1)]\n else\n polys[[j]] <- d[(id[i] + 1):length(d)]\n}\n\npolys <- lapply(polys, function(x) {\n tf <- tempfile()\n writeLines(x, tf)\n pol <- as.matrix(read.csv(tf, header = FALSE))\n unlink(tf)\n return(st_polygon(list(pol)))\n})\n\n## transform to sf object\npolys_sfc <- st_sfc(polys)\nmap <- st_sf(geometry = polys_sfc)\nmap$id <- names(polys)\nmap$district <- 1:nrow(map)\n\n## plot the map\npar(mar = rep(0, 4))\nplot(st_geometry(map))\n\n\n\n\n\n\n\n\n\n## sample coordinates for plotting\nco <- st_sample(map, size = 10000, type = \"regular\")\n\n## create new data for prediction\nnd <- as.data.frame(st_coordinates(co))\nnames(nd) <- c(\"xcoord\", \"ycoord\")\nnd$sex <- 1\nnd$wbc <- mean(LeukSurv$wbc)\nnd$tpi <- mean(LeukSurv$tpi)\nnd$age <- mean(LeukSurv$age)\n\n## predict parameters\npar <- predict(b2, newdata = nd)\n\n## compute survival probabilities\np2 <- 1 - family(b2)$p(Surv(rep(365, nrow(nd)), rep(1, nrow(nd))), par)\n\n## plot on map\nsp <- st_sf(geometry = co)\nsp$survprob <- p2\n\n## reference system, only used for plotting\nsp <- st_set_crs(sp, 4326)\nmap <- st_set_crs(map, 4326)\n\n## plot as raster map\nspr <- st_as_stars(sp)\n\npar(mar = c(0, 0, 5, 0))\nplot(spr, main = \"Survival Probabilities Prob(t > 365)\",\n reset = FALSE, nbreaks = 100)\nplot(st_geometry(map), add = TRUE)\nbox()", "crumbs": [ "Articles", "Survival Models" @@ -409,7 +431,7 @@ "href": "man/prodist.gamlss2.html", "title": "gamlss2", "section": "", - "text": "Methods for gamlss2 model objects for extracting fitted (in-sample) or predicted (out-of-sample) probability distributions as distributions3 objects.\n\n\n\n## S3 method for class 'gamlss2'\nprodist(object, ...)\n\n\n\n\n\n\n\nobject\n\n\nA model object of class gamlss2.\n\n\n\n\n…\n\n\nArguments passed on to predict.gamlss2, e.g., newdata.\n\n\n\n\n\n\nTo facilitate making probabilistic forecasts based on gamlss2 model objects, the prodist method extracts fitted or predicted probability distribution objects. Internally, the predict.gamlss2 method is used first to obtain the distribution parameters (mu, sigma, tau, nu, or a subset thereof). Subsequently, the corresponding distribution object is set up using the GAMLSS class from the gamlss.dist package, enabling the workflow provided by the distributions3 package (see Zeileis et al. 2022).\nNote that these probability distributions only reflect the random variation in the dependent variable based on the model employed (and its associated distributional assumption for the dependent variable). This does not capture the uncertainty in the parameter estimates.\n\n\n\nAn object of class GAMLSS inheriting from distribution.\n\n\n\nZeileis A, Lang MN, Hayes A (2022). “distributions3: From Basic Probability to Probabilistic Regression.” Presented at useR! 2022 - The R User Conference. Slides, video, vignette, code at https://www.zeileis.org/news/user2022/.\n\n\n\nGAMLSS, predict.gamlss2\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## packages, code, and data\nlibrary(\"distributions3\")\ndata(\"cars\", package = \"datasets\")\n\n## fit heteroscedastic normal GAMLSS model\n## stopping distance (ft) explained by speed (mph)\nm <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO)\n\nGAMLSS-RS iteration 1: Global Deviance = 407.3541 eps = 0.125497 \nGAMLSS-RS iteration 2: Global Deviance = 405.7146 eps = 0.004024 \nGAMLSS-RS iteration 3: Global Deviance = 405.6978 eps = 0.000041 \nGAMLSS-RS iteration 4: Global Deviance = 405.6976 eps = 0.000000 \n\n## obtain predicted distributions for three levels of speed\nd <- prodist(m, newdata = data.frame(speed = c(10, 20, 30)))\nprint(d)\n\n 1 2 \n\"GAMLSS NO(mu = 23.04, sigma = 10.06)\" \"GAMLSS NO(mu = 59.04, sigma = 18.51)\" \n 3 \n\"GAMLSS NO(mu = 96.35, sigma = 33.95)\" \n\n## obtain quantiles (works the same for any distribution object 'd' !)\nquantile(d, 0.5)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = FALSE)\n\n q_0.05 q_0.5 q_0.95\n1 6.486962 23.03912 39.59128\n2 28.589641 59.03607 89.48250\n3 40.504887 96.34896 152.19303\n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = TRUE)\n\n 1 2 3 \n 6.486962 59.036073 152.193030 \n\n## visualization\nplot(dist ~ speed, data = cars)\nnd <- data.frame(speed = 0:240/4)\nnd$dist <- prodist(m, newdata = nd)\nnd$fit <- quantile(nd$dist, c(0.05, 0.5, 0.95))\nmatplot(nd$speed, nd$fit, type = \"l\", lty = 1, col = \"slategray\", add = TRUE)\n\n\n\n\n\n\n\n## moments\nmean(d)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nvariance(d)\n\n 1 2 3 \n 101.2639 342.6244 1152.6558 \n\n## simulate random numbers\nrandom(d, 5)\n\n r_1 r_2 r_3 r_4 r_5\n1 22.07090 29.92421 18.67442 29.23931 13.09939\n2 19.52497 41.95755 49.88034 28.46416 79.20217\n3 82.29180 72.53890 83.23728 76.62891 74.27616\n\n## density and distribution\npdf(d, 50 * -2:2)\n\n d_-100 d_-50 d_0 d_50 d_100\n1 1.365786e-34 1.440750e-13 0.0028836944 0.001095127 7.891037e-15\n2 2.012473e-18 6.289547e-10 0.0001332376 0.019131662 1.862073e-03\n3 6.414012e-10 1.084300e-06 0.0002095201 0.004627633 1.168286e-02\n\ncdf(d, 50 * -2:2)\n\n p_-100 p_-50 p_0 p_50 p_100\n1 1.116699e-34 1.961566e-13 0.0110254856 0.99631019 1.0000000\n2 4.279141e-18 1.923739e-09 0.0007128545 0.31271491 0.9865531\n3 3.661574e-09 8.139843e-06 0.0022705648 0.08609812 0.5428194\n\n## Poisson example\ndata(\"FIFA2018\", package = \"distributions3\")\nm2 <- gamlss2(goals ~ s(difference), data = FIFA2018, family = PO)\n\nGAMLSS-RS iteration 1: Global Deviance = 355.3922 eps = 0.045332 \nGAMLSS-RS iteration 2: Global Deviance = 355.3922 eps = 0.000000 \n\nd2 <- prodist(m2, newdata = data.frame(difference = 0))\nprint(d2)\n\n 1 \n\"GAMLSS PO(mu = 1.237)\" \n\nquantile(d2, c(0.05, 0.5, 0.95))\n\n[1] 0 1 3\n\n## note that log_pdf() can replicate logLik() value\nsum(log_pdf(prodist(m2), FIFA2018$goals))\n\n[1] -177.6961\n\nlogLik(m2)\n\n'log Lik.' -177.6961 (df=2.005144)", + "text": "Methods for gamlss2 model objects for extracting fitted (in-sample) or predicted (out-of-sample) probability distributions as distributions3 objects.\n\n\n\n## S3 method for class 'gamlss2'\nprodist(object, ...)\n\n\n\n\n\n\n\nobject\n\n\nA model object of class gamlss2.\n\n\n\n\n…\n\n\nArguments passed on to predict.gamlss2, e.g., newdata.\n\n\n\n\n\n\nTo facilitate making probabilistic forecasts based on gamlss2 model objects, the prodist method extracts fitted or predicted probability distribution objects. Internally, the predict.gamlss2 method is used first to obtain the distribution parameters (mu, sigma, tau, nu, or a subset thereof). Subsequently, the corresponding distribution object is set up using the GAMLSS class from the gamlss.dist package, enabling the workflow provided by the distributions3 package (see Zeileis et al. 2022).\nNote that these probability distributions only reflect the random variation in the dependent variable based on the model employed (and its associated distributional assumption for the dependent variable). This does not capture the uncertainty in the parameter estimates.\n\n\n\nAn object of class GAMLSS inheriting from distribution.\n\n\n\nZeileis A, Lang MN, Hayes A (2022). “distributions3: From Basic Probability to Probabilistic Regression.” Presented at useR! 2022 - The R User Conference. Slides, video, vignette, code at https://www.zeileis.org/news/user2022/.\n\n\n\nGAMLSS, predict.gamlss2\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## packages, code, and data\nlibrary(\"distributions3\")\ndata(\"cars\", package = \"datasets\")\n\n## fit heteroscedastic normal GAMLSS model\n## stopping distance (ft) explained by speed (mph)\nm <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO)\n\nGAMLSS-RS iteration 1: Global Deviance = 407.3541 eps = 0.125497 \nGAMLSS-RS iteration 2: Global Deviance = 405.7146 eps = 0.004024 \nGAMLSS-RS iteration 3: Global Deviance = 405.6978 eps = 0.000041 \nGAMLSS-RS iteration 4: Global Deviance = 405.6976 eps = 0.000000 \n\n## obtain predicted distributions for three levels of speed\nd <- prodist(m, newdata = data.frame(speed = c(10, 20, 30)))\nprint(d)\n\n 1 2 \n\"GAMLSS NO(mu = 23.04, sigma = 10.06)\" \"GAMLSS NO(mu = 59.04, sigma = 18.51)\" \n 3 \n\"GAMLSS NO(mu = 96.35, sigma = 33.95)\" \n\n## obtain quantiles (works the same for any distribution object 'd' !)\nquantile(d, 0.5)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = FALSE)\n\n q_0.05 q_0.5 q_0.95\n1 6.486962 23.03912 39.59128\n2 28.589641 59.03607 89.48250\n3 40.504887 96.34896 152.19303\n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = TRUE)\n\n 1 2 3 \n 6.486962 59.036073 152.193030 \n\n## visualization\nplot(dist ~ speed, data = cars)\nnd <- data.frame(speed = 0:240/4)\nnd$dist <- prodist(m, newdata = nd)\nnd$fit <- quantile(nd$dist, c(0.05, 0.5, 0.95))\nmatplot(nd$speed, nd$fit, type = \"l\", lty = 1, col = \"slategray\", add = TRUE)\n\n\n\n\n\n\n\n## moments\nmean(d)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nvariance(d)\n\n 1 2 3 \n 101.2639 342.6244 1152.6558 \n\n## simulate random numbers\nrandom(d, 5)\n\n r_1 r_2 r_3 r_4 r_5\n1 31.99176 24.26454 29.24719 38.53560 22.14378\n2 32.35560 57.28431 41.84400 55.37877 53.11516\n3 100.10360 161.57604 152.17778 150.07732 94.09042\n\n## density and distribution\npdf(d, 50 * -2:2)\n\n d_-100 d_-50 d_0 d_50 d_100\n1 1.365786e-34 1.440750e-13 0.0028836944 0.001095127 7.891037e-15\n2 2.012473e-18 6.289547e-10 0.0001332376 0.019131662 1.862073e-03\n3 6.414012e-10 1.084300e-06 0.0002095201 0.004627633 1.168286e-02\n\ncdf(d, 50 * -2:2)\n\n p_-100 p_-50 p_0 p_50 p_100\n1 1.116699e-34 1.961566e-13 0.0110254856 0.99631019 1.0000000\n2 4.279141e-18 1.923739e-09 0.0007128545 0.31271491 0.9865531\n3 3.661574e-09 8.139843e-06 0.0022705648 0.08609812 0.5428194\n\n## Poisson example\ndata(\"FIFA2018\", package = \"distributions3\")\nm2 <- gamlss2(goals ~ s(difference), data = FIFA2018, family = PO)\n\nGAMLSS-RS iteration 1: Global Deviance = 355.3922 eps = 0.045332 \nGAMLSS-RS iteration 2: Global Deviance = 355.3922 eps = 0.000000 \n\nd2 <- prodist(m2, newdata = data.frame(difference = 0))\nprint(d2)\n\n 1 \n\"GAMLSS PO(mu = 1.237)\" \n\nquantile(d2, c(0.05, 0.5, 0.95))\n\n[1] 0 1 3\n\n## note that log_pdf() can replicate logLik() value\nsum(log_pdf(prodist(m2), FIFA2018$goals))\n\n[1] -177.6961\n\nlogLik(m2)\n\n'log Lik.' -177.6961 (df=2.005144)", "crumbs": [ "Reference", "prodist.gamlss2" @@ -420,7 +442,7 @@ "href": "man/prodist.gamlss2.html#extracting-fitted-or-predicted-probability-distributions-from-gamlss2-models", "title": "gamlss2", "section": "", - "text": "Methods for gamlss2 model objects for extracting fitted (in-sample) or predicted (out-of-sample) probability distributions as distributions3 objects.\n\n\n\n## S3 method for class 'gamlss2'\nprodist(object, ...)\n\n\n\n\n\n\n\nobject\n\n\nA model object of class gamlss2.\n\n\n\n\n…\n\n\nArguments passed on to predict.gamlss2, e.g., newdata.\n\n\n\n\n\n\nTo facilitate making probabilistic forecasts based on gamlss2 model objects, the prodist method extracts fitted or predicted probability distribution objects. Internally, the predict.gamlss2 method is used first to obtain the distribution parameters (mu, sigma, tau, nu, or a subset thereof). Subsequently, the corresponding distribution object is set up using the GAMLSS class from the gamlss.dist package, enabling the workflow provided by the distributions3 package (see Zeileis et al. 2022).\nNote that these probability distributions only reflect the random variation in the dependent variable based on the model employed (and its associated distributional assumption for the dependent variable). This does not capture the uncertainty in the parameter estimates.\n\n\n\nAn object of class GAMLSS inheriting from distribution.\n\n\n\nZeileis A, Lang MN, Hayes A (2022). “distributions3: From Basic Probability to Probabilistic Regression.” Presented at useR! 2022 - The R User Conference. Slides, video, vignette, code at https://www.zeileis.org/news/user2022/.\n\n\n\nGAMLSS, predict.gamlss2\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## packages, code, and data\nlibrary(\"distributions3\")\ndata(\"cars\", package = \"datasets\")\n\n## fit heteroscedastic normal GAMLSS model\n## stopping distance (ft) explained by speed (mph)\nm <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO)\n\nGAMLSS-RS iteration 1: Global Deviance = 407.3541 eps = 0.125497 \nGAMLSS-RS iteration 2: Global Deviance = 405.7146 eps = 0.004024 \nGAMLSS-RS iteration 3: Global Deviance = 405.6978 eps = 0.000041 \nGAMLSS-RS iteration 4: Global Deviance = 405.6976 eps = 0.000000 \n\n## obtain predicted distributions for three levels of speed\nd <- prodist(m, newdata = data.frame(speed = c(10, 20, 30)))\nprint(d)\n\n 1 2 \n\"GAMLSS NO(mu = 23.04, sigma = 10.06)\" \"GAMLSS NO(mu = 59.04, sigma = 18.51)\" \n 3 \n\"GAMLSS NO(mu = 96.35, sigma = 33.95)\" \n\n## obtain quantiles (works the same for any distribution object 'd' !)\nquantile(d, 0.5)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = FALSE)\n\n q_0.05 q_0.5 q_0.95\n1 6.486962 23.03912 39.59128\n2 28.589641 59.03607 89.48250\n3 40.504887 96.34896 152.19303\n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = TRUE)\n\n 1 2 3 \n 6.486962 59.036073 152.193030 \n\n## visualization\nplot(dist ~ speed, data = cars)\nnd <- data.frame(speed = 0:240/4)\nnd$dist <- prodist(m, newdata = nd)\nnd$fit <- quantile(nd$dist, c(0.05, 0.5, 0.95))\nmatplot(nd$speed, nd$fit, type = \"l\", lty = 1, col = \"slategray\", add = TRUE)\n\n\n\n\n\n\n\n## moments\nmean(d)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nvariance(d)\n\n 1 2 3 \n 101.2639 342.6244 1152.6558 \n\n## simulate random numbers\nrandom(d, 5)\n\n r_1 r_2 r_3 r_4 r_5\n1 22.07090 29.92421 18.67442 29.23931 13.09939\n2 19.52497 41.95755 49.88034 28.46416 79.20217\n3 82.29180 72.53890 83.23728 76.62891 74.27616\n\n## density and distribution\npdf(d, 50 * -2:2)\n\n d_-100 d_-50 d_0 d_50 d_100\n1 1.365786e-34 1.440750e-13 0.0028836944 0.001095127 7.891037e-15\n2 2.012473e-18 6.289547e-10 0.0001332376 0.019131662 1.862073e-03\n3 6.414012e-10 1.084300e-06 0.0002095201 0.004627633 1.168286e-02\n\ncdf(d, 50 * -2:2)\n\n p_-100 p_-50 p_0 p_50 p_100\n1 1.116699e-34 1.961566e-13 0.0110254856 0.99631019 1.0000000\n2 4.279141e-18 1.923739e-09 0.0007128545 0.31271491 0.9865531\n3 3.661574e-09 8.139843e-06 0.0022705648 0.08609812 0.5428194\n\n## Poisson example\ndata(\"FIFA2018\", package = \"distributions3\")\nm2 <- gamlss2(goals ~ s(difference), data = FIFA2018, family = PO)\n\nGAMLSS-RS iteration 1: Global Deviance = 355.3922 eps = 0.045332 \nGAMLSS-RS iteration 2: Global Deviance = 355.3922 eps = 0.000000 \n\nd2 <- prodist(m2, newdata = data.frame(difference = 0))\nprint(d2)\n\n 1 \n\"GAMLSS PO(mu = 1.237)\" \n\nquantile(d2, c(0.05, 0.5, 0.95))\n\n[1] 0 1 3\n\n## note that log_pdf() can replicate logLik() value\nsum(log_pdf(prodist(m2), FIFA2018$goals))\n\n[1] -177.6961\n\nlogLik(m2)\n\n'log Lik.' -177.6961 (df=2.005144)", + "text": "Methods for gamlss2 model objects for extracting fitted (in-sample) or predicted (out-of-sample) probability distributions as distributions3 objects.\n\n\n\n## S3 method for class 'gamlss2'\nprodist(object, ...)\n\n\n\n\n\n\n\nobject\n\n\nA model object of class gamlss2.\n\n\n\n\n…\n\n\nArguments passed on to predict.gamlss2, e.g., newdata.\n\n\n\n\n\n\nTo facilitate making probabilistic forecasts based on gamlss2 model objects, the prodist method extracts fitted or predicted probability distribution objects. Internally, the predict.gamlss2 method is used first to obtain the distribution parameters (mu, sigma, tau, nu, or a subset thereof). Subsequently, the corresponding distribution object is set up using the GAMLSS class from the gamlss.dist package, enabling the workflow provided by the distributions3 package (see Zeileis et al. 2022).\nNote that these probability distributions only reflect the random variation in the dependent variable based on the model employed (and its associated distributional assumption for the dependent variable). This does not capture the uncertainty in the parameter estimates.\n\n\n\nAn object of class GAMLSS inheriting from distribution.\n\n\n\nZeileis A, Lang MN, Hayes A (2022). “distributions3: From Basic Probability to Probabilistic Regression.” Presented at useR! 2022 - The R User Conference. Slides, video, vignette, code at https://www.zeileis.org/news/user2022/.\n\n\n\nGAMLSS, predict.gamlss2\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## packages, code, and data\nlibrary(\"distributions3\")\ndata(\"cars\", package = \"datasets\")\n\n## fit heteroscedastic normal GAMLSS model\n## stopping distance (ft) explained by speed (mph)\nm <- gamlss2(dist ~ s(speed) | s(speed), data = cars, family = NO)\n\nGAMLSS-RS iteration 1: Global Deviance = 407.3541 eps = 0.125497 \nGAMLSS-RS iteration 2: Global Deviance = 405.7146 eps = 0.004024 \nGAMLSS-RS iteration 3: Global Deviance = 405.6978 eps = 0.000041 \nGAMLSS-RS iteration 4: Global Deviance = 405.6976 eps = 0.000000 \n\n## obtain predicted distributions for three levels of speed\nd <- prodist(m, newdata = data.frame(speed = c(10, 20, 30)))\nprint(d)\n\n 1 2 \n\"GAMLSS NO(mu = 23.04, sigma = 10.06)\" \"GAMLSS NO(mu = 59.04, sigma = 18.51)\" \n 3 \n\"GAMLSS NO(mu = 96.35, sigma = 33.95)\" \n\n## obtain quantiles (works the same for any distribution object 'd' !)\nquantile(d, 0.5)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = FALSE)\n\n q_0.05 q_0.5 q_0.95\n1 6.486962 23.03912 39.59128\n2 28.589641 59.03607 89.48250\n3 40.504887 96.34896 152.19303\n\nquantile(d, c(0.05, 0.5, 0.95), elementwise = TRUE)\n\n 1 2 3 \n 6.486962 59.036073 152.193030 \n\n## visualization\nplot(dist ~ speed, data = cars)\nnd <- data.frame(speed = 0:240/4)\nnd$dist <- prodist(m, newdata = nd)\nnd$fit <- quantile(nd$dist, c(0.05, 0.5, 0.95))\nmatplot(nd$speed, nd$fit, type = \"l\", lty = 1, col = \"slategray\", add = TRUE)\n\n\n\n\n\n\n\n## moments\nmean(d)\n\n 1 2 3 \n23.03912 59.03607 96.34896 \n\nvariance(d)\n\n 1 2 3 \n 101.2639 342.6244 1152.6558 \n\n## simulate random numbers\nrandom(d, 5)\n\n r_1 r_2 r_3 r_4 r_5\n1 31.99176 24.26454 29.24719 38.53560 22.14378\n2 32.35560 57.28431 41.84400 55.37877 53.11516\n3 100.10360 161.57604 152.17778 150.07732 94.09042\n\n## density and distribution\npdf(d, 50 * -2:2)\n\n d_-100 d_-50 d_0 d_50 d_100\n1 1.365786e-34 1.440750e-13 0.0028836944 0.001095127 7.891037e-15\n2 2.012473e-18 6.289547e-10 0.0001332376 0.019131662 1.862073e-03\n3 6.414012e-10 1.084300e-06 0.0002095201 0.004627633 1.168286e-02\n\ncdf(d, 50 * -2:2)\n\n p_-100 p_-50 p_0 p_50 p_100\n1 1.116699e-34 1.961566e-13 0.0110254856 0.99631019 1.0000000\n2 4.279141e-18 1.923739e-09 0.0007128545 0.31271491 0.9865531\n3 3.661574e-09 8.139843e-06 0.0022705648 0.08609812 0.5428194\n\n## Poisson example\ndata(\"FIFA2018\", package = \"distributions3\")\nm2 <- gamlss2(goals ~ s(difference), data = FIFA2018, family = PO)\n\nGAMLSS-RS iteration 1: Global Deviance = 355.3922 eps = 0.045332 \nGAMLSS-RS iteration 2: Global Deviance = 355.3922 eps = 0.000000 \n\nd2 <- prodist(m2, newdata = data.frame(difference = 0))\nprint(d2)\n\n 1 \n\"GAMLSS PO(mu = 1.237)\" \n\nquantile(d2, c(0.05, 0.5, 0.95))\n\n[1] 0 1 3\n\n## note that log_pdf() can replicate logLik() value\nsum(log_pdf(prodist(m2), FIFA2018$goals))\n\n[1] -177.6961\n\nlogLik(m2)\n\n'log Lik.' -177.6961 (df=2.005144)", "crumbs": [ "Reference", "prodist.gamlss2" @@ -651,7 +673,7 @@ "href": "man/gamlss2.html", "title": "gamlss2", "section": "", - "text": "Estimation of generalized additive models for location scale and shape (GAMLSS). The model fitting function gamlss2() provides flexible infrastructures to estimate the parameters of a response distribution. The number of distributional parameters is not fixed, see gamlss2.family. Moreover, gamlss2() supports all smooth term constructors from the mgcv package in addition to the classical model terms as provided by gamlss and gamlss.add.\n\n\n\ngamlss2(x, ...)\n\n## S3 method for class 'formula'\ngamlss2(formula, data, family = NO,\n subset, na.action, weights, offset, start = NULL,\n control = gamlss2_control(...), ...)\n\n## S3 method for class 'list'\ngamlss2(x, ...)\n\n\n\n\n\n\n\nformula\n\n\nA GAM-type formula or Formula. All smooth terms of the mgcv package are supported, see also formula.gam.\n\n\n\n\nx\n\n\nFor gamlss.list() x is a list of formulas.\n\n\n\n\ndata\n\n\nA data frame or list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which gamlss2 is called.\n\n\n\n\nfamily\n\n\nA gamlss.family or gamlss2.family object used to define distribution and the link functions of the parameters.\n\n\n\n\nsubset\n\n\nAn optional vector specifying a subset of observations to be used in the fitting process.\n\n\n\n\nna.action\n\n\nNA processing for setting up the model.frame.\n\n\n\n\nweights\n\n\nAn optional vector of prior weights to be used in the fitting process. Should be NULL or a numeric vector.\n\n\n\n\noffset\n\n\nThis can be used to specify an a priori known components to be included in the linear predictors during fitting. Please note that if only a single numeric vector is provided, the offset will be assigned to the first specified parameter of the distribution. In the case of multiple offsets, a data frame or list must be supplied. Each offset is assigned in the same order as the parameters of the distribution specified in the family object.\n\n\n\n\nstart\n\n\nStarting values for estimation algorithms.\n\n\n\n\ncontrol\n\n\nA list of control arguments, see gamlss2_control.\n\n\n\n\n…\n\n\nArguments passed to gamlss2_control.\n\n\n\n\n\n\nThe model fitting function gamlss2() provides flexible infrastructures for the estimation of GAMLSS.\n\n\nDistributional models are specified using family objects, either from the gamlss.dist package or using gamlss2.family objects.\n\n\nEstimation is carried out through a Newton-Raphson/Fisher scoring algorithm, see function RS. The estimation algorithms can also be exchanged using gamlss2_control. Additionally, if an optimizer is specified by the family object, this optimizer function will be employed for estimation.\n\n\nThe return value is determined by the object returned from the optimizer function, typically an object of class “gamlss2”. Default methods and extractor functions are available for this class. Nevertheless, users have the flexibility to supply their own optimizer function, along with user-specific methods tailored for the returned object.\n\n\n\n\n\nThe return value is determined by the object returned from the optimizer function. By default, the optimization is performed using the RS optimizer function (see gamlss2_control), yielding an object of class “gamlss2”. Default methods and extractor functions are available for this class.\n\n\n\nRigby RA, Stasinopoulos DM (2005). “Generalized Additive Models for Location, Scale and Shape (with Discussion).” Journal of the Royal Statistical Society, Series C (Applied Statistics), 54, 507–554. doi:10.1111/j.1467-9876.2005.00510.x\nRigby RA, Stasinopoulos DM, Heller GZ, De Bastiani F (2019). Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/9780429298547\nStasinopoulos DM, Rigby RA (2007). “Generalized Additive Models for Location Scale and Shape (GAMLSS) in R.” Journal of Statistical Software, 23(7), 1–46. doi:10.18637/jss.v023.i07\nStasinopoulos DM, Rigby RA, Heller GZ, Voudouris V, De Bastiani F (2017). Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/b21973\n\n\n\nRS, gamlss2_control, gamlss2.family\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## load the abdominal circumference data\ndata(\"abdom\", package = \"gamlss.data\")\n\n## specify the model Formula\nf <- y ~ s(x) | s(x) | s(x) | s(x)\n\n## estimate model\nb <- gamlss2(f, data = abdom, family = BCT)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002 \n\n## model summary\nsummary(b)\n\nCall:\ngamlss2(formula = f, data = abdom, family = BCT)\n---\nFamily: BCT \nLink function: mu = identity, sigma = log, nu = identity, tau = log\n*--------\nParameter: mu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 226.334 1.257 180 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 4.551\n*--------\nParameter: sigma \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -2.92264 0.01101 -265.5 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 2.5639\n*--------\nParameter: nu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.18021 0.04599 -3.918 9.95e-05 ***\n---\nSmooth terms:\n s(x)\nedf 1.0015\n*--------\nParameter: tau \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 2.6548 0.0144 184.4 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n---\nSmooth terms:\n s(x)\nedf 1.0042\n*--------\nn = 610 df = 13.12 res.df = 596.88\nDeviance = 4770.1554 Null Dev. Red. = 33.39%\nAIC = 4796.3966 elapsed = 0.77sec\n\n## plot estimated effects\nplot(b, which = \"effects\")\n\n\n\n\n\n\n\n## plot diagnostics\nplot(b, which = \"resid\")\n\n\n\n\n\n\n\n## predict parameters\npar <- predict(b)\n\n## predict quantiles\npq <- sapply(c(0.05, 0.5, 0.95), function(q) family(b)$q(q, par))\n\n## visualize\nplot(y ~ x, data = abdom, pch = 19,\n col = rgb(0.1, 0.1, 0.1, alpha = 0.3))\nmatplot(abdom$x, pq, type = \"l\", lwd = 2,\n lty = 1, col = 4, add = TRUE)\n\n\n\n\n\n\n\n## use of starting values\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10))\n\nGAMLSS-RS iteration 1: Global Deviance = 4789.7797 eps = 0.556012 \nGAMLSS-RS iteration 2: Global Deviance = 4774.5657 eps = 0.003176 \nGAMLSS-RS iteration 3: Global Deviance = 4771.4193 eps = 0.000658 \nGAMLSS-RS iteration 4: Global Deviance = 4769.9947 eps = 0.000298 \nGAMLSS-RS iteration 5: Global Deviance = 4769.9556 eps = 0.000008 \n\n## fix some parameters\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10),\n fixed = c(nu = TRUE, tau = TRUE))\n\nGAMLSS-RS iteration 1: Global Deviance = 4799.422 eps = 0.555118 \nGAMLSS-RS iteration 2: Global Deviance = 4795.2807 eps = 0.000862 \nGAMLSS-RS iteration 3: Global Deviance = 4795.2668 eps = 0.000002 \n\n## estimated coefficients (intercepts)\ncoef(m)\n\n mu.p.(Intercept) sigma.p.(Intercept) nu.p.(Intercept) tau.p.(Intercept) \n 226.347632 -2.922923 0.000000 2.302585 \n\n## starting values using full predictors\nm <- gamlss2(f, data = abdom, family = BCT,\n start = fitted(m))\n\nGAMLSS-RS iteration 1: Global Deviance = 4902.0767 eps = 0.372276 \nGAMLSS-RS iteration 2: Global Deviance = 4775.8465 eps = 0.025750 \nGAMLSS-RS iteration 3: Global Deviance = 4775.0302 eps = 0.000170 \nGAMLSS-RS iteration 4: Global Deviance = 4774.9582 eps = 0.000015 \nGAMLSS-RS iteration 5: Global Deviance = 4774.9519 eps = 0.000001 \n\n## same with\nm <- gamlss2(f, data = abdom, family = BCT,\n start = m)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002", + "text": "Estimation of generalized additive models for location scale and shape (GAMLSS). The model fitting function gamlss2() provides flexible infrastructures to estimate the parameters of a response distribution. The number of distributional parameters is not fixed, see gamlss2.family. Moreover, gamlss2() supports all smooth term constructors from the mgcv package in addition to the classical model terms as provided by gamlss and gamlss.add.\n\n\n\ngamlss2(x, ...)\n\n## S3 method for class 'formula'\ngamlss2(formula, data, family = NO,\n subset, na.action, weights, offset, start = NULL,\n control = gamlss2_control(...), ...)\n\n## S3 method for class 'list'\ngamlss2(x, ...)\n\n\n\n\n\n\n\nformula\n\n\nA GAM-type formula or Formula. All smooth terms of the mgcv package are supported, see also formula.gam.\n\n\n\n\nx\n\n\nFor gamlss.list() x is a list of formulas.\n\n\n\n\ndata\n\n\nA data frame or list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which gamlss2 is called.\n\n\n\n\nfamily\n\n\nA gamlss.family or gamlss2.family object used to define distribution and the link functions of the parameters.\n\n\n\n\nsubset\n\n\nAn optional vector specifying a subset of observations to be used in the fitting process.\n\n\n\n\nna.action\n\n\nNA processing for setting up the model.frame.\n\n\n\n\nweights\n\n\nAn optional vector of prior weights to be used in the fitting process. Should be NULL or a numeric vector.\n\n\n\n\noffset\n\n\nThis can be used to specify an a priori known components to be included in the linear predictors during fitting. Please note that if only a single numeric vector is provided, the offset will be assigned to the first specified parameter of the distribution. In the case of multiple offsets, a data frame or list must be supplied. Each offset is assigned in the same order as the parameters of the distribution specified in the family object.\n\n\n\n\nstart\n\n\nStarting values for estimation algorithms.\n\n\n\n\ncontrol\n\n\nA list of control arguments, see gamlss2_control.\n\n\n\n\n…\n\n\nArguments passed to gamlss2_control.\n\n\n\n\n\n\nThe model fitting function gamlss2() provides flexible infrastructures for the estimation of GAMLSS.\n\n\nDistributional models are specified using family objects, either from the gamlss.dist package or using gamlss2.family objects.\n\n\nEstimation is carried out through a Newton-Raphson/Fisher scoring algorithm, see function RS. The estimation algorithms can also be exchanged using gamlss2_control. Additionally, if an optimizer is specified by the family object, this optimizer function will be employed for estimation.\n\n\nThe return value is determined by the object returned from the optimizer function, typically an object of class “gamlss2”. Default methods and extractor functions are available for this class. Nevertheless, users have the flexibility to supply their own optimizer function, along with user-specific methods tailored for the returned object.\n\n\n\n\n\nThe return value is determined by the object returned from the optimizer function. By default, the optimization is performed using the RS optimizer function (see gamlss2_control), yielding an object of class “gamlss2”. Default methods and extractor functions are available for this class.\n\n\n\nRigby RA, Stasinopoulos DM (2005). “Generalized Additive Models for Location, Scale and Shape (with Discussion).” Journal of the Royal Statistical Society, Series C (Applied Statistics), 54, 507–554. doi:10.1111/j.1467-9876.2005.00510.x\nRigby RA, Stasinopoulos DM, Heller GZ, De Bastiani F (2019). Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/9780429298547\nStasinopoulos DM, Rigby RA (2007). “Generalized Additive Models for Location Scale and Shape (GAMLSS) in R.” Journal of Statistical Software, 23(7), 1–46. doi:10.18637/jss.v023.i07\nStasinopoulos DM, Rigby RA, Heller GZ, Voudouris V, De Bastiani F (2017). Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/b21973\n\n\n\nRS, gamlss2_control, gamlss2.family\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## load the abdominal circumference data\ndata(\"abdom\", package = \"gamlss.data\")\n\n## specify the model Formula\nf <- y ~ s(x) | s(x) | s(x) | s(x)\n\n## estimate model\nb <- gamlss2(f, data = abdom, family = BCT)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002 \n\n## model summary\nsummary(b)\n\nCall:\ngamlss2(formula = f, data = abdom, family = BCT)\n---\nFamily: BCT \nLink function: mu = identity, sigma = log, nu = identity, tau = log\n*--------\nParameter: mu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 226.334 1.257 180 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 4.551\n*--------\nParameter: sigma \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -2.92264 0.01101 -265.5 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 2.5639\n*--------\nParameter: nu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.18021 0.04599 -3.918 9.95e-05 ***\n---\nSmooth terms:\n s(x)\nedf 1.0015\n*--------\nParameter: tau \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 2.6548 0.0144 184.4 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n---\nSmooth terms:\n s(x)\nedf 1.0042\n*--------\nn = 610 df = 13.12 res.df = 596.88\nDeviance = 4770.1554 Null Dev. Red. = 33.39%\nAIC = 4796.3966 elapsed = 0.80sec\n\n## plot estimated effects\nplot(b, which = \"effects\")\n\n\n\n\n\n\n\n## plot diagnostics\nplot(b, which = \"resid\")\n\n\n\n\n\n\n\n## predict parameters\npar <- predict(b)\n\n## predict quantiles\npq <- sapply(c(0.05, 0.5, 0.95), function(q) family(b)$q(q, par))\n\n## visualize\nplot(y ~ x, data = abdom, pch = 19,\n col = rgb(0.1, 0.1, 0.1, alpha = 0.3))\nmatplot(abdom$x, pq, type = \"l\", lwd = 2,\n lty = 1, col = 4, add = TRUE)\n\n\n\n\n\n\n\n## use of starting values\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10))\n\nGAMLSS-RS iteration 1: Global Deviance = 4789.7797 eps = 0.556012 \nGAMLSS-RS iteration 2: Global Deviance = 4774.5657 eps = 0.003176 \nGAMLSS-RS iteration 3: Global Deviance = 4771.4193 eps = 0.000658 \nGAMLSS-RS iteration 4: Global Deviance = 4769.9947 eps = 0.000298 \nGAMLSS-RS iteration 5: Global Deviance = 4769.9556 eps = 0.000008 \n\n## fix some parameters\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10),\n fixed = c(nu = TRUE, tau = TRUE))\n\nGAMLSS-RS iteration 1: Global Deviance = 4799.422 eps = 0.555118 \nGAMLSS-RS iteration 2: Global Deviance = 4795.2807 eps = 0.000862 \nGAMLSS-RS iteration 3: Global Deviance = 4795.2668 eps = 0.000002 \n\n## estimated coefficients (intercepts)\ncoef(m)\n\n mu.p.(Intercept) sigma.p.(Intercept) nu.p.(Intercept) tau.p.(Intercept) \n 226.347632 -2.922923 0.000000 2.302585 \n\n## starting values using full predictors\nm <- gamlss2(f, data = abdom, family = BCT,\n start = fitted(m))\n\nGAMLSS-RS iteration 1: Global Deviance = 4902.0767 eps = 0.372276 \nGAMLSS-RS iteration 2: Global Deviance = 4775.8465 eps = 0.025750 \nGAMLSS-RS iteration 3: Global Deviance = 4775.0302 eps = 0.000170 \nGAMLSS-RS iteration 4: Global Deviance = 4774.9582 eps = 0.000015 \nGAMLSS-RS iteration 5: Global Deviance = 4774.9519 eps = 0.000001 \n\n## same with\nm <- gamlss2(f, data = abdom, family = BCT,\n start = m)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002", "crumbs": [ "Reference", "gamlss2" @@ -662,7 +684,7 @@ "href": "man/gamlss2.html#generalized-additive-models-for-location-scale-and-shape", "title": "gamlss2", "section": "", - "text": "Estimation of generalized additive models for location scale and shape (GAMLSS). The model fitting function gamlss2() provides flexible infrastructures to estimate the parameters of a response distribution. The number of distributional parameters is not fixed, see gamlss2.family. Moreover, gamlss2() supports all smooth term constructors from the mgcv package in addition to the classical model terms as provided by gamlss and gamlss.add.\n\n\n\ngamlss2(x, ...)\n\n## S3 method for class 'formula'\ngamlss2(formula, data, family = NO,\n subset, na.action, weights, offset, start = NULL,\n control = gamlss2_control(...), ...)\n\n## S3 method for class 'list'\ngamlss2(x, ...)\n\n\n\n\n\n\n\nformula\n\n\nA GAM-type formula or Formula. All smooth terms of the mgcv package are supported, see also formula.gam.\n\n\n\n\nx\n\n\nFor gamlss.list() x is a list of formulas.\n\n\n\n\ndata\n\n\nA data frame or list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which gamlss2 is called.\n\n\n\n\nfamily\n\n\nA gamlss.family or gamlss2.family object used to define distribution and the link functions of the parameters.\n\n\n\n\nsubset\n\n\nAn optional vector specifying a subset of observations to be used in the fitting process.\n\n\n\n\nna.action\n\n\nNA processing for setting up the model.frame.\n\n\n\n\nweights\n\n\nAn optional vector of prior weights to be used in the fitting process. Should be NULL or a numeric vector.\n\n\n\n\noffset\n\n\nThis can be used to specify an a priori known components to be included in the linear predictors during fitting. Please note that if only a single numeric vector is provided, the offset will be assigned to the first specified parameter of the distribution. In the case of multiple offsets, a data frame or list must be supplied. Each offset is assigned in the same order as the parameters of the distribution specified in the family object.\n\n\n\n\nstart\n\n\nStarting values for estimation algorithms.\n\n\n\n\ncontrol\n\n\nA list of control arguments, see gamlss2_control.\n\n\n\n\n…\n\n\nArguments passed to gamlss2_control.\n\n\n\n\n\n\nThe model fitting function gamlss2() provides flexible infrastructures for the estimation of GAMLSS.\n\n\nDistributional models are specified using family objects, either from the gamlss.dist package or using gamlss2.family objects.\n\n\nEstimation is carried out through a Newton-Raphson/Fisher scoring algorithm, see function RS. The estimation algorithms can also be exchanged using gamlss2_control. Additionally, if an optimizer is specified by the family object, this optimizer function will be employed for estimation.\n\n\nThe return value is determined by the object returned from the optimizer function, typically an object of class “gamlss2”. Default methods and extractor functions are available for this class. Nevertheless, users have the flexibility to supply their own optimizer function, along with user-specific methods tailored for the returned object.\n\n\n\n\n\nThe return value is determined by the object returned from the optimizer function. By default, the optimization is performed using the RS optimizer function (see gamlss2_control), yielding an object of class “gamlss2”. Default methods and extractor functions are available for this class.\n\n\n\nRigby RA, Stasinopoulos DM (2005). “Generalized Additive Models for Location, Scale and Shape (with Discussion).” Journal of the Royal Statistical Society, Series C (Applied Statistics), 54, 507–554. doi:10.1111/j.1467-9876.2005.00510.x\nRigby RA, Stasinopoulos DM, Heller GZ, De Bastiani F (2019). Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/9780429298547\nStasinopoulos DM, Rigby RA (2007). “Generalized Additive Models for Location Scale and Shape (GAMLSS) in R.” Journal of Statistical Software, 23(7), 1–46. doi:10.18637/jss.v023.i07\nStasinopoulos DM, Rigby RA, Heller GZ, Voudouris V, De Bastiani F (2017). Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/b21973\n\n\n\nRS, gamlss2_control, gamlss2.family\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## load the abdominal circumference data\ndata(\"abdom\", package = \"gamlss.data\")\n\n## specify the model Formula\nf <- y ~ s(x) | s(x) | s(x) | s(x)\n\n## estimate model\nb <- gamlss2(f, data = abdom, family = BCT)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002 \n\n## model summary\nsummary(b)\n\nCall:\ngamlss2(formula = f, data = abdom, family = BCT)\n---\nFamily: BCT \nLink function: mu = identity, sigma = log, nu = identity, tau = log\n*--------\nParameter: mu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 226.334 1.257 180 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 4.551\n*--------\nParameter: sigma \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -2.92264 0.01101 -265.5 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 2.5639\n*--------\nParameter: nu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.18021 0.04599 -3.918 9.95e-05 ***\n---\nSmooth terms:\n s(x)\nedf 1.0015\n*--------\nParameter: tau \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 2.6548 0.0144 184.4 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n---\nSmooth terms:\n s(x)\nedf 1.0042\n*--------\nn = 610 df = 13.12 res.df = 596.88\nDeviance = 4770.1554 Null Dev. Red. = 33.39%\nAIC = 4796.3966 elapsed = 0.77sec\n\n## plot estimated effects\nplot(b, which = \"effects\")\n\n\n\n\n\n\n\n## plot diagnostics\nplot(b, which = \"resid\")\n\n\n\n\n\n\n\n## predict parameters\npar <- predict(b)\n\n## predict quantiles\npq <- sapply(c(0.05, 0.5, 0.95), function(q) family(b)$q(q, par))\n\n## visualize\nplot(y ~ x, data = abdom, pch = 19,\n col = rgb(0.1, 0.1, 0.1, alpha = 0.3))\nmatplot(abdom$x, pq, type = \"l\", lwd = 2,\n lty = 1, col = 4, add = TRUE)\n\n\n\n\n\n\n\n## use of starting values\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10))\n\nGAMLSS-RS iteration 1: Global Deviance = 4789.7797 eps = 0.556012 \nGAMLSS-RS iteration 2: Global Deviance = 4774.5657 eps = 0.003176 \nGAMLSS-RS iteration 3: Global Deviance = 4771.4193 eps = 0.000658 \nGAMLSS-RS iteration 4: Global Deviance = 4769.9947 eps = 0.000298 \nGAMLSS-RS iteration 5: Global Deviance = 4769.9556 eps = 0.000008 \n\n## fix some parameters\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10),\n fixed = c(nu = TRUE, tau = TRUE))\n\nGAMLSS-RS iteration 1: Global Deviance = 4799.422 eps = 0.555118 \nGAMLSS-RS iteration 2: Global Deviance = 4795.2807 eps = 0.000862 \nGAMLSS-RS iteration 3: Global Deviance = 4795.2668 eps = 0.000002 \n\n## estimated coefficients (intercepts)\ncoef(m)\n\n mu.p.(Intercept) sigma.p.(Intercept) nu.p.(Intercept) tau.p.(Intercept) \n 226.347632 -2.922923 0.000000 2.302585 \n\n## starting values using full predictors\nm <- gamlss2(f, data = abdom, family = BCT,\n start = fitted(m))\n\nGAMLSS-RS iteration 1: Global Deviance = 4902.0767 eps = 0.372276 \nGAMLSS-RS iteration 2: Global Deviance = 4775.8465 eps = 0.025750 \nGAMLSS-RS iteration 3: Global Deviance = 4775.0302 eps = 0.000170 \nGAMLSS-RS iteration 4: Global Deviance = 4774.9582 eps = 0.000015 \nGAMLSS-RS iteration 5: Global Deviance = 4774.9519 eps = 0.000001 \n\n## same with\nm <- gamlss2(f, data = abdom, family = BCT,\n start = m)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002", + "text": "Estimation of generalized additive models for location scale and shape (GAMLSS). The model fitting function gamlss2() provides flexible infrastructures to estimate the parameters of a response distribution. The number of distributional parameters is not fixed, see gamlss2.family. Moreover, gamlss2() supports all smooth term constructors from the mgcv package in addition to the classical model terms as provided by gamlss and gamlss.add.\n\n\n\ngamlss2(x, ...)\n\n## S3 method for class 'formula'\ngamlss2(formula, data, family = NO,\n subset, na.action, weights, offset, start = NULL,\n control = gamlss2_control(...), ...)\n\n## S3 method for class 'list'\ngamlss2(x, ...)\n\n\n\n\n\n\n\nformula\n\n\nA GAM-type formula or Formula. All smooth terms of the mgcv package are supported, see also formula.gam.\n\n\n\n\nx\n\n\nFor gamlss.list() x is a list of formulas.\n\n\n\n\ndata\n\n\nA data frame or list or environment containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which gamlss2 is called.\n\n\n\n\nfamily\n\n\nA gamlss.family or gamlss2.family object used to define distribution and the link functions of the parameters.\n\n\n\n\nsubset\n\n\nAn optional vector specifying a subset of observations to be used in the fitting process.\n\n\n\n\nna.action\n\n\nNA processing for setting up the model.frame.\n\n\n\n\nweights\n\n\nAn optional vector of prior weights to be used in the fitting process. Should be NULL or a numeric vector.\n\n\n\n\noffset\n\n\nThis can be used to specify an a priori known components to be included in the linear predictors during fitting. Please note that if only a single numeric vector is provided, the offset will be assigned to the first specified parameter of the distribution. In the case of multiple offsets, a data frame or list must be supplied. Each offset is assigned in the same order as the parameters of the distribution specified in the family object.\n\n\n\n\nstart\n\n\nStarting values for estimation algorithms.\n\n\n\n\ncontrol\n\n\nA list of control arguments, see gamlss2_control.\n\n\n\n\n…\n\n\nArguments passed to gamlss2_control.\n\n\n\n\n\n\nThe model fitting function gamlss2() provides flexible infrastructures for the estimation of GAMLSS.\n\n\nDistributional models are specified using family objects, either from the gamlss.dist package or using gamlss2.family objects.\n\n\nEstimation is carried out through a Newton-Raphson/Fisher scoring algorithm, see function RS. The estimation algorithms can also be exchanged using gamlss2_control. Additionally, if an optimizer is specified by the family object, this optimizer function will be employed for estimation.\n\n\nThe return value is determined by the object returned from the optimizer function, typically an object of class “gamlss2”. Default methods and extractor functions are available for this class. Nevertheless, users have the flexibility to supply their own optimizer function, along with user-specific methods tailored for the returned object.\n\n\n\n\n\nThe return value is determined by the object returned from the optimizer function. By default, the optimization is performed using the RS optimizer function (see gamlss2_control), yielding an object of class “gamlss2”. Default methods and extractor functions are available for this class.\n\n\n\nRigby RA, Stasinopoulos DM (2005). “Generalized Additive Models for Location, Scale and Shape (with Discussion).” Journal of the Royal Statistical Society, Series C (Applied Statistics), 54, 507–554. doi:10.1111/j.1467-9876.2005.00510.x\nRigby RA, Stasinopoulos DM, Heller GZ, De Bastiani F (2019). Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/9780429298547\nStasinopoulos DM, Rigby RA (2007). “Generalized Additive Models for Location Scale and Shape (GAMLSS) in R.” Journal of Statistical Software, 23(7), 1–46. doi:10.18637/jss.v023.i07\nStasinopoulos DM, Rigby RA, Heller GZ, Voudouris V, De Bastiani F (2017). Flexible Regression and Smoothing: Using GAMLSS in R, Chapman and Hall/CRC. doi:10.1201/b21973\n\n\n\nRS, gamlss2_control, gamlss2.family\n\n\n\n\nlibrary(\"gamlss2\")\n\n\n## load the abdominal circumference data\ndata(\"abdom\", package = \"gamlss.data\")\n\n## specify the model Formula\nf <- y ~ s(x) | s(x) | s(x) | s(x)\n\n## estimate model\nb <- gamlss2(f, data = abdom, family = BCT)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002 \n\n## model summary\nsummary(b)\n\nCall:\ngamlss2(formula = f, data = abdom, family = BCT)\n---\nFamily: BCT \nLink function: mu = identity, sigma = log, nu = identity, tau = log\n*--------\nParameter: mu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 226.334 1.257 180 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 4.551\n*--------\nParameter: sigma \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -2.92264 0.01101 -265.5 <2e-16 ***\n---\nSmooth terms:\n s(x)\nedf 2.5639\n*--------\nParameter: nu \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -0.18021 0.04599 -3.918 9.95e-05 ***\n---\nSmooth terms:\n s(x)\nedf 1.0015\n*--------\nParameter: tau \n---\nCoefficients:\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) 2.6548 0.0144 184.4 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n---\nSmooth terms:\n s(x)\nedf 1.0042\n*--------\nn = 610 df = 13.12 res.df = 596.88\nDeviance = 4770.1554 Null Dev. Red. = 33.39%\nAIC = 4796.3966 elapsed = 0.80sec\n\n## plot estimated effects\nplot(b, which = \"effects\")\n\n\n\n\n\n\n\n## plot diagnostics\nplot(b, which = \"resid\")\n\n\n\n\n\n\n\n## predict parameters\npar <- predict(b)\n\n## predict quantiles\npq <- sapply(c(0.05, 0.5, 0.95), function(q) family(b)$q(q, par))\n\n## visualize\nplot(y ~ x, data = abdom, pch = 19,\n col = rgb(0.1, 0.1, 0.1, alpha = 0.3))\nmatplot(abdom$x, pq, type = \"l\", lwd = 2,\n lty = 1, col = 4, add = TRUE)\n\n\n\n\n\n\n\n## use of starting values\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10))\n\nGAMLSS-RS iteration 1: Global Deviance = 4789.7797 eps = 0.556012 \nGAMLSS-RS iteration 2: Global Deviance = 4774.5657 eps = 0.003176 \nGAMLSS-RS iteration 3: Global Deviance = 4771.4193 eps = 0.000658 \nGAMLSS-RS iteration 4: Global Deviance = 4769.9947 eps = 0.000298 \nGAMLSS-RS iteration 5: Global Deviance = 4769.9556 eps = 0.000008 \n\n## fix some parameters\nm <- gamlss2(f, data = abdom, family = BCT,\n start = c(mu = 200, sigma = 0.1, nu = 0, tau = 10),\n fixed = c(nu = TRUE, tau = TRUE))\n\nGAMLSS-RS iteration 1: Global Deviance = 4799.422 eps = 0.555118 \nGAMLSS-RS iteration 2: Global Deviance = 4795.2807 eps = 0.000862 \nGAMLSS-RS iteration 3: Global Deviance = 4795.2668 eps = 0.000002 \n\n## estimated coefficients (intercepts)\ncoef(m)\n\n mu.p.(Intercept) sigma.p.(Intercept) nu.p.(Intercept) tau.p.(Intercept) \n 226.347632 -2.922923 0.000000 2.302585 \n\n## starting values using full predictors\nm <- gamlss2(f, data = abdom, family = BCT,\n start = fitted(m))\n\nGAMLSS-RS iteration 1: Global Deviance = 4902.0767 eps = 0.372276 \nGAMLSS-RS iteration 2: Global Deviance = 4775.8465 eps = 0.025750 \nGAMLSS-RS iteration 3: Global Deviance = 4775.0302 eps = 0.000170 \nGAMLSS-RS iteration 4: Global Deviance = 4774.9582 eps = 0.000015 \nGAMLSS-RS iteration 5: Global Deviance = 4774.9519 eps = 0.000001 \n\n## same with\nm <- gamlss2(f, data = abdom, family = BCT,\n start = m)\n\nGAMLSS-RS iteration 1: Global Deviance = 4774.4683 eps = 0.534345 \nGAMLSS-RS iteration 2: Global Deviance = 4770.229 eps = 0.000887 \nGAMLSS-RS iteration 3: Global Deviance = 4770.1663 eps = 0.000013 \nGAMLSS-RS iteration 4: Global Deviance = 4770.1554 eps = 0.000002", "crumbs": [ "Reference", "gamlss2" diff --git a/vignettes/gamlss2_files/figure-html/unnamed-chunk-9-1.png b/vignettes/gamlss2_files/figure-html/unnamed-chunk-9-1.png index 1f8032b..252b702 100644 Binary files a/vignettes/gamlss2_files/figure-html/unnamed-chunk-9-1.png and b/vignettes/gamlss2_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/vignettes/spatial.html b/vignettes/spatial.html index 27623b2..5e5416f 100644 --- a/vignettes/spatial.html +++ b/vignettes/spatial.html @@ -601,7 +601,7 @@

+AIC = 16682.1271 elapsed = 3.92sec

Model calibration is checked using histogram, Q-Q plot, wormplot etc.

diff --git a/vignettes/spatial_files/figure-html/unnamed-chunk-7-1.png b/vignettes/spatial_files/figure-html/unnamed-chunk-7-1.png index c9ce49e..3130a08 100644 Binary files a/vignettes/spatial_files/figure-html/unnamed-chunk-7-1.png and b/vignettes/spatial_files/figure-html/unnamed-chunk-7-1.png differ diff --git a/vignettes/survival.html b/vignettes/survival.html index 9e7ab3e..cf734a2 100644 --- a/vignettes/survival.html +++ b/vignettes/survival.html @@ -121,6 +121,35 @@ } } + + + + @@ -415,8 +444,15 @@
-

First, we estimate time only models.

@@ -489,18 +545,21 @@

Survival Models

+ +
+

2 Spatial survival model

Now, we estimate a full spatial survival model.

## model formula
-f <- Surv(time, cens) ~ sex + s(age) + s(wbc) + s(tpi) + s(xcoord,ycoord,k=100) |
-  sex + s(age) + s(wbc) + s(tpi) + s(xcoord,ycoord,k=100)
+f <- Surv(time, cens) ~ sex + s(age) + s(wbc) + s(tpi) + s(xcoord, ycoord) |
+  sex + s(age) + s(wbc) + s(tpi) + s(xcoord, ycoord)
 
 ## estimate model
 b2 <- gamlss2(f, family = fam, data = LeukSurv)
-
GAMLSS-RS iteration  1: Global Deviance = 11863.604 eps = 0.042545     
-GAMLSS-RS iteration  2: Global Deviance = 11859.5156 eps = 0.000344     
-GAMLSS-RS iteration  3: Global Deviance = 11859.4914 eps = 0.000002
+
GAMLSS-RS iteration  1: Global Deviance = 11869.8787 eps = 0.042038     
+GAMLSS-RS iteration  2: Global Deviance = 11866.6901 eps = 0.000268     
+GAMLSS-RS iteration  3: Global Deviance = 11866.619 eps = 0.000005

Plot estimated effects.

@@ -558,33 +617,7 @@

Survival Models

## plot the map par(mar = rep(0, 4)) -plot(st_geometry(map)) - -## sample coordinates for plotting -co <- st_sample(map, size = 1000, type = "regular") - -## create new data for prediction -nd <- as.data.frame(st_coordinates(co)) -names(nd) <- c("xcoord", "ycoord") -nd$sex <- 1 -nd$wbc <- mean(LeukSurv$wbc) -nd$tpi <- mean(LeukSurv$tpi) -nd$age <- mean(LeukSurv$age) - -## predict parameters -par <- predict(b2, newdata = nd) - -## compute survival probabilities -p2 <- 1 - family(b2)$p(Surv(rep(365, nrow(nd)), rep(1, nrow(nd))), par) - -## plot on map -sp <- st_sf(geometry = co) -sp$survprob <- p2 - -par(mar = rep(0, 4)) -plot(st_geometry(map)) -plot(sp, pch = 15, cex = 1.3, add = TRUE) -plot(st_geometry(map), add = TRUE) +plot(st_geometry(map))
@@ -593,9 +626,52 @@

Survival Models

+
+
## sample coordinates for plotting
+co <- st_sample(map, size = 10000, type = "regular")
+
+## create new data for prediction
+nd <- as.data.frame(st_coordinates(co))
+names(nd) <- c("xcoord", "ycoord")
+nd$sex <- 1
+nd$wbc <- mean(LeukSurv$wbc)
+nd$tpi <- mean(LeukSurv$tpi)
+nd$age <- mean(LeukSurv$age)
+
+## predict parameters
+par <- predict(b2, newdata = nd)
+
+## compute survival probabilities
+p2 <- 1 - family(b2)$p(Surv(rep(365, nrow(nd)), rep(1, nrow(nd))), par)
+
+## plot on map
+sp <- st_sf(geometry = co)
+sp$survprob <- p2
+
+## reference system, only used for plotting
+sp <- st_set_crs(sp, 4326)
+map <- st_set_crs(map, 4326)
+
+## plot as raster map
+spr <- st_as_stars(sp)
+
+par(mar = c(0, 0, 5, 0))
+plot(spr, main = "Survival Probabilities Prob(t > 365)",
+  reset = FALSE, nbreaks = 100)
+plot(st_geometry(map), add = TRUE)
+box()
+
+
+
+

+
+
+
+
+

References

diff --git a/vignettes/survival_files/figure-html/unnamed-chunk-6-1.png b/vignettes/survival_files/figure-html/unnamed-chunk-6-1.png index 4e74226..5bed2c3 100644 Binary files a/vignettes/survival_files/figure-html/unnamed-chunk-6-1.png and b/vignettes/survival_files/figure-html/unnamed-chunk-6-1.png differ diff --git a/vignettes/survival_files/figure-html/unnamed-chunk-7-1.png b/vignettes/survival_files/figure-html/unnamed-chunk-7-1.png index 51596e5..decec37 100644 Binary files a/vignettes/survival_files/figure-html/unnamed-chunk-7-1.png and b/vignettes/survival_files/figure-html/unnamed-chunk-7-1.png differ diff --git a/vignettes/survival_files/figure-html/unnamed-chunk-8-1.png b/vignettes/survival_files/figure-html/unnamed-chunk-8-1.png index ae51aae..a2fa68f 100644 Binary files a/vignettes/survival_files/figure-html/unnamed-chunk-8-1.png and b/vignettes/survival_files/figure-html/unnamed-chunk-8-1.png differ diff --git a/vignettes/survival_files/figure-html/unnamed-chunk-9-1.png b/vignettes/survival_files/figure-html/unnamed-chunk-9-1.png new file mode 100644 index 0000000..63cbbbf Binary files /dev/null and b/vignettes/survival_files/figure-html/unnamed-chunk-9-1.png differ