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 f6531cd..6dbe50d 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 c8a4afd..d5030c6 100644 --- a/man/gamlss2.html +++ b/man/gamlss2.html @@ -629,7 +629,7 @@

Examples

*-------- n = 610 df = 13.12 res.df = 596.88 Deviance = 4770.1554 Null Dev. Red. = 33.39% -AIC = 4796.3966 elapsed = 0.79sec +AIC = 4796.3966 elapsed = 0.78sec 
## plot estimated effects
 plot(b, which = "effects")
diff --git a/man/prodist.gamlss2.html b/man/prodist.gamlss2.html index 77ed4bb..74e19a8 100644 --- a/man/prodist.gamlss2.html +++ b/man/prodist.gamlss2.html @@ -537,10 +537,10 @@

Examples

## simulate random numbers
 random(d, 5)
-
       r_1       r_2       r_3      r_4       r_5
-1 20.05192  12.01664  28.09310 17.86293  32.08695
-2 39.93538  54.87385  58.66783 84.24932  77.59599
-3 75.54687 126.27259 128.96634 74.50611 173.73709
+
        r_1      r_2       r_3      r_4      r_5
+1  21.37530 30.88462  25.30309 11.07565 21.63578
+2  36.46126 34.15373  73.37629 75.47346 42.17611
+3 102.67860 55.63671 100.03750 99.18436 65.02560
## density and distribution
 pdf(d, 50 * -2:2)
diff --git a/search.json b/search.json index 2041558..6ee1376 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 = 4.45sec\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 = 4.40sec\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" @@ -124,7 +124,7 @@ "href": "vignettes/mixture.html", "title": "Mixture Models", "section": "", - "text": "Mixture models are used for modeling data generated from multiple distinct processes, where each process can be described by a separate probability distribution. The challenge is that for each observation, we do not know which process generated it, and therefore we must infer the underlying latent classes (or components).\nThis vignette demonstrates how to fit a mixture of normal distributions where the mixing probabilities depend on a covariate.", + "text": "Mixture models are used for modeling data generated from multiple distinct processes, where each process can be described by a separate probability distribution. The challenge is that for each observation, we do not know which process generated it, and therefore we must infer the underlying latent components.\nThis vignette demonstrates how to fit a mixture of normal distributions where the mixing probabilities depend on a covariate.", "crumbs": [ "Articles", "Mixture Models" @@ -135,7 +135,7 @@ "href": "vignettes/mixture.html#simulating-data", "title": "Mixture Models", "section": "1 Simulating Data", - "text": "1 Simulating Data\nWe begin by simulating data where the response y is generated by one of two normal distributions, depending on a covariate x. The latent class membership is determined by a logistic model, where the mixing probability is a function of x.\n\nset.seed(123)\n\n## simulate covariate\nn <- 1000\nx <- sort(runif(n, -pi, pi))\n\n## logistic model for mixing probability (dependent on x)\n## this creates a probability that varies with x using a logistic function\nmix_prob <- 1 / (1 + exp(-1 - 0.5 * x))\n\n## simulate latent component assignment based on covariate-dependent probability\nz <- rbinom(n, 1, mix_prob)\n\n## generate response based on the latent component\ny <- ifelse(z == 1,\n sin(x) + rnorm(n, sd = 0.2),\n cos(x) + rnorm(n, sd = 0.2))\n\n## combine into a data frame\nd <- data.frame(\"x\" = x, \"y\" = y)\n\n## plot the data, color by latent component z\npar(mar = c(4, 4, 4, 1))\nplot(d, col = z + 1, main = \"Simulated Data by Latent Class\", \n xlab = \"x\", ylab = \"y\")", + "text": "1 Simulating Data\nWe begin by simulating data where the response y is generated by one of two normal distributions, depending on a covariate x. The latent component membership is determined by a logistic model, where the mixing probability is a function of x.\n\nset.seed(123)\n\n## simulate covariate\nn <- 1000\nx <- sort(runif(n, -pi, pi))\n\n## logistic model for mixing probability (dependent on x)\n## this creates a probability that varies with x using a logistic function\nmix_prob <- 1 / (1 + exp(-1 - 0.5 * x))\n\n## simulate latent component assignment based on covariate-dependent probability\nz <- rbinom(n, 1, mix_prob)\n\n## generate response based on the latent component\ny <- ifelse(z == 1,\n sin(x) + rnorm(n, sd = 0.2),\n cos(x) + rnorm(n, sd = 0.2))\n\n## combine into a data frame\nd <- data.frame(\"x\" = x, \"y\" = y)\n\n## plot the data, color by latent component z\npar(mar = c(4, 4, 4, 1))\nplot(d, col = z + 1, main = \"Simulated Data by Latent Component\", \n xlab = \"x\", ylab = \"y\")", "crumbs": [ "Articles", "Mixture Models" @@ -168,7 +168,7 @@ "href": "vignettes/mixture.html#predicting-and-visualizing-results", "title": "Mixture Models", "section": "4 Predicting and Visualizing Results", - "text": "4 Predicting and Visualizing Results\nWe predict the component-specific means (mu1, mu2) and the mixing probability (pi) from the fitted model. Then, we compare the predicted and true class probabilities and plot the fitted means along with the data.\nFirst, compute predicted parameters\n\np <- predict(b, type = \"parameter\")\n\nPlot estimated vs. true component probabilities\n\nplot(p$pi, mix_prob, main = \"Estimated vs. True Class Probabilities\", \n xlab = \"Estimated Probability\", ylab = \"True Probability\")\nabline(0, 1)\n\n\n\n\n\n\n\n## compute posterior class probabilities\nd1 <- dnorm(d$y, p$mu1, p$sigma1)\nd2 <- dnorm(d$y, p$mu2, p$sigma2)\ntotprob <- rowSums(cbind(p$pi * d1, (1 - p$pi) * d2))\nc1 <- p$pi * d1 / totprob\nc2 <- (1 - p$pi) * d2 / totprob\n\n## get classes\nclasses <- apply(cbind(c2, c1), 1, which.max)\n\n## plot the data and fitted component means\npar(mar = c(4, 4, 4, 1))\nplot(y ~ x, data = d,\n main = \"Fitted Component Means\",\n xlab = \"x\", ylab = \"y\", col = classes)\nlines(p$mu1 ~ x, col = 2, lwd = 4)\nlines(p$mu2 ~ x, col = 1, lwd = 4)", + "text": "4 Predicting and Visualizing Results\nWe predict the component-specific means (mu1, mu2) and the mixing probability (pi) from the fitted model. Then, we compare the predicted and true component probabilities and plot the fitted means along with the data.\nFirst, compute predicted parameters.\n\np <- predict(b, type = \"parameter\")\n\nPlot estimated vs. true component probabilities.\n\nplot(p$pi, mix_prob, main = \"Estimated vs. True Component Probabilities\", \n xlab = \"Estimated Probability\", ylab = \"True Probability\")\nabline(0, 1)\n\n\n\n\n\n\n\n## compute posterior component probabilities\nd1 <- dnorm(d$y, p$mu1, p$sigma1)\nd2 <- dnorm(d$y, p$mu2, p$sigma2)\ntotprob <- rowSums(cbind(p$pi * d1, (1 - p$pi) * d2))\nc1 <- p$pi * d1 / totprob\nc2 <- (1 - p$pi) * d2 / totprob\n\n## get components\ncomp <- apply(cbind(c2, c1), 1, which.max)\n\n## plot the data and fitted component means\npar(mar = c(4, 4, 4, 1))\nplot(y ~ x, data = d,\n main = \"Fitted Component Means\",\n xlab = \"x\", ylab = \"y\", col = comp)\nlines(p$mu1 ~ x, col = 2, lwd = 4)\nlines(p$mu2 ~ x, col = 1, lwd = 4)", "crumbs": [ "Articles", "Mixture Models" @@ -179,7 +179,7 @@ "href": "vignettes/mixture.html#summary", "title": "Mixture Models", "section": "5 Summary", - "text": "5 Summary\nIn this vignette, we have demonstrated how to fit a mixture model where the mixing probabilities depend on a covariate. We simulated data, defined a custom mixture family, and fitted the model to estimate the component means and mixing probabilities. Finally, we compared the predicted class probabilities with the true values and visualized the results.", + "text": "5 Summary\nIn this vignette, we have demonstrated how to fit a mixture model where the mixing probabilities depend on a covariate. We simulated data, defined a custom mixture family, and fitted the model to estimate the component means and mixing probabilities. Finally, we compared the predicted component probabilities with the true values and visualized the results.", "crumbs": [ "Articles", "Mixture Models" @@ -299,7 +299,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.79sec\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.78sec\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" @@ -310,7 +310,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.79sec\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.78sec\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" @@ -541,7 +541,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 20.05192 12.01664 28.09310 17.86293 32.08695\n2 39.93538 54.87385 58.66783 84.24932 77.59599\n3 75.54687 126.27259 128.96634 74.50611 173.73709\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 21.37530 30.88462 25.30309 11.07565 21.63578\n2 36.46126 34.15373 73.37629 75.47346 42.17611\n3 102.67860 55.63671 100.03750 99.18436 65.02560\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" @@ -552,7 +552,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 20.05192 12.01664 28.09310 17.86293 32.08695\n2 39.93538 54.87385 58.66783 84.24932 77.59599\n3 75.54687 126.27259 128.96634 74.50611 173.73709\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 21.37530 30.88462 25.30309 11.07565 21.63578\n2 36.46126 34.15373 73.37629 75.47346 42.17611\n3 102.67860 55.63671 100.03750 99.18436 65.02560\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" 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 67f2e8a..044674e 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/mixture.html b/vignettes/mixture.html index a70e2a5..ef398d0 100644 --- a/vignettes/mixture.html +++ b/vignettes/mixture.html @@ -431,11 +431,11 @@

Mixture Models

-

Mixture models are used for modeling data generated from multiple distinct processes, where each process can be described by a separate probability distribution. The challenge is that for each observation, we do not know which process generated it, and therefore we must infer the underlying latent classes (or components).

+

Mixture models are used for modeling data generated from multiple distinct processes, where each process can be described by a separate probability distribution. The challenge is that for each observation, we do not know which process generated it, and therefore we must infer the underlying latent components.

This vignette demonstrates how to fit a mixture of normal distributions where the mixing probabilities depend on a covariate.

1 Simulating Data

-

We begin by simulating data where the response y is generated by one of two normal distributions, depending on a covariate x. The latent class membership is determined by a logistic model, where the mixing probability is a function of x.

+

We begin by simulating data where the response y is generated by one of two normal distributions, depending on a covariate x. The latent component membership is determined by a logistic model, where the mixing probability is a function of x.

set.seed(123)
 
@@ -460,7 +460,7 @@ 

 ## plot the data, color by latent component z
 par(mar = c(4, 4, 4, 1))
-plot(d, col = z + 1, main = "Simulated Data by Latent Class", 
+plot(d, col = z + 1, main = "Simulated Data by Latent Component", 
      xlab = "x", ylab = "y")
@@ -529,14 +529,14 @@

4 Predicting and Visualizing Results

-

We predict the component-specific means (mu1, mu2) and the mixing probability (pi) from the fitted model. Then, we compare the predicted and true class probabilities and plot the fitted means along with the data.

-

First, compute predicted parameters

+

We predict the component-specific means (mu1, mu2) and the mixing probability (pi) from the fitted model. Then, we compare the predicted and true component probabilities and plot the fitted means along with the data.

+

First, compute predicted parameters.

p <- predict(b, type = "parameter")
-

Plot estimated vs. true component probabilities

+

Plot estimated vs. true component probabilities.

-
plot(p$pi, mix_prob, main = "Estimated vs. True Class Probabilities", 
+
plot(p$pi, mix_prob, main = "Estimated vs. True Component Probabilities", 
      xlab = "Estimated Probability", ylab = "True Probability")
 abline(0, 1)

 

@@ -546,21 +546,21 @@ 
## compute posterior class probabilities
+
## compute posterior component probabilities
 d1 <- dnorm(d$y, p$mu1, p$sigma1)
 d2 <- dnorm(d$y, p$mu2, p$sigma2)
 totprob <- rowSums(cbind(p$pi * d1, (1 - p$pi) * d2))
 c1 <- p$pi * d1  / totprob
 c2 <- (1 - p$pi) * d2  / totprob
 
-## get classes
-classes <- apply(cbind(c2, c1), 1, which.max)
+## get components
+comp <- apply(cbind(c2, c1), 1, which.max)
 
 ## plot the data and fitted component means
 par(mar = c(4, 4, 4, 1))
 plot(y ~ x, data = d,
   main = "Fitted Component Means",
-  xlab = "x", ylab = "y", col = classes)
+  xlab = "x", ylab = "y", col = comp)
 lines(p$mu1 ~ x, col = 2, lwd = 4)
 lines(p$mu2 ~ x, col = 1, lwd = 4)

 

@@ -574,7 +574,7 @@ 

 
5 Summary

-
In this vignette, we have demonstrated how to fit a mixture model where the mixing probabilities depend on a covariate. We simulated data, defined a custom mixture family, and fitted the model to estimate the component means and mixing probabilities. Finally, we compared the predicted class probabilities with the true values and visualized the results.

+
In this vignette, we have demonstrated how to fit a mixture model where the mixing probabilities depend on a covariate. We simulated data, defined a custom mixture family, and fitted the model to estimate the component means and mixing probabilities. Finally, we compared the predicted component probabilities with the true values and visualized the results.

 
 
 
diff --git a/vignettes/mixture_files/figure-html/unnamed-chunk-1-1.png b/vignettes/mixture_files/figure-html/unnamed-chunk-1-1.png
index d766115..0107f52 100644
Binary files a/vignettes/mixture_files/figure-html/unnamed-chunk-1-1.png and b/vignettes/mixture_files/figure-html/unnamed-chunk-1-1.png differ
diff --git a/vignettes/mixture_files/figure-html/unnamed-chunk-6-1.png b/vignettes/mixture_files/figure-html/unnamed-chunk-6-1.png
index 1cff4b9..9117023 100644
Binary files a/vignettes/mixture_files/figure-html/unnamed-chunk-6-1.png and b/vignettes/mixture_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/vignettes/spatial.html b/vignettes/spatial.html
index 94b2854..720a1c4 100644
--- a/vignettes/spatial.html
+++ b/vignettes/spatial.html
@@ -595,7 +595,7 @@ 

+AIC = 16682.1271 elapsed =  4.40sec

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 8ea89f6..326d8b7 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