Movatterモバイル変換

Generic function for Accumulated Local Effect (ALE)

Description

Generic function for producing ALE effects, to be plottied using theplot generic.

Usage

ALE(o, ...)

Arguments

References

Apley, D.W., and Zhu, J, 2016. Visualizing the effects of predictor variables in blackbox supervised learning models. arXiv preprint arXiv:1612.08468.

See Also

ALE.gam

Create Accumulated Local Effects (ALE) for GAMs

Description

Create Accumulated Local Effects (ALE) for GAMs

Usage

## S3 method for class 'gam'ALE(  o,  x,  newdata = NULL,  type = "link",  nbin = 40,  bins = NULL,  oind = 1,  center = 1,  ...)

Arguments

Value

An object of classALEXD, whereX is the number of dimensions, which can be plottedusingplot.ALEXD (onlyX=1 is provided at the moment).

Author(s)

Matteo Fasiolo and Christian Capezza, with some internal code having been adapted from the ALEPlotpackage of Dan Apley.

References

Apley, D.W., and Zhu, J, 2016. Visualizing the effects of predictor variables in blackbox supervised learning models. arXiv preprint arXiv:1612.08468.

See Also

plot.ALE1D

Examples

# Example using Tweedie distributionlibrary(mgcViz)set.seed(3)n<-400## Simulate data...dat <- gamSim(1,n=n,dist="poisson",scale=.2)dat$y <- rTweedie(exp(dat$f),p=1.3,phi=.5) ## Tweedie response## Fit a fixed p Tweedie, with wrong link ...b <- gam(list(y~s(x0)+s(x1)+s(x2)+s(x3),~1,~1), family=twlss(), data=dat)plot(ALE(b, "x2", type = "response", oind = 1))# With manually chosen binsplot(ALE(b, "x2", type = "response", oind = 1,          bins = c(0.1, 0.25, 0.5, 0.6, 0.9, 0.95, 0.99, 1)))

Fit a GAM model and get a gamViz object

Description

These are wrapper that fits a GAM model usingmgcv::gam ormgcv::bam andconverts it to agamViz object using thegetViz function.It is essentially a shortcut.

Usage

bamV(  formula,  family = gaussian(),  data = list(),  method = "fREML",  aGam = list(),  aViz = list())gamV(  formula,  family = gaussian(),  data = list(),  method = "REML",  aGam = list(),  aViz = list())

Arguments

Value

An object of class "gamViz" which can, for instance, be plotted usingplot.gamViz.

Examples

##### gam example# Simulate datalibrary(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(1,n=1000,dist="normal",scale=2)# Fit GAM and get gamViz objectb <- gamV(y~s(x0)+s(x1, x2)+s(x3), data = dat,           aGam = list(scale = 2), aViz = list("nsim" = 20))# This is equivalent to doing# 1. Fit GAM# b <- gam(y~s(x0)+s(x1, x2)+s(x3), data=dat, method="REML", scale = 2)# 2. Convert to gamViz object# b <- getViz(b, nsim = 20)# Either way, we plot first and third effects by doingprint(plot(b, select = c(1, 3)), pages = 1)##### bam example# Simulate datalibrary(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(1,n=2000,dist="normal",scale=2)# Fit using bam() and get gamViz objectb <- bamV(y~s(x0)+s(x1, x2)+s(x3), data = dat,           aGam = list(discrete = TRUE), aViz = list("nsim" = 0))          # Either way, we plot first and third effects by doingprint(plot(b, select = c(2)), pages = 1)

Some diagnostics for a fitted gam model

Description

Takes a fitted GAM model and produces some diagnostic information about the fittingprocedure and results. The default is to produce 4 residual plots, some information aboutthe convergence of the smoothness selection optimization, and to run diagnostic tests ofwhether the basis dimension choises are adequate.

Usage

## S3 method for class 'gamViz'check(  obj,  type = c("auto", "deviance", "pearson", "response", "tunif", "tnormal"),  k.sample = 5000,  k.rep = 200,  maxpo = 10000,  a.qq = list(),  a.hist = list(),  a.respoi = list(),  ...)

Arguments

Details

This is a essentially a re-write ofmgcv::gam.check usingggplot2. Seemgcv::gam.check for details.

Value

An object of classcheckGam, which is simply a list ofggplot objects.

Examples

library(mgcViz)set.seed(0)dat <- gamSim(1, n = 200)b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat)b <- getViz(b)# Checks using default optionscheck(b)# Change some algorithmic and graphical parameterscheck(b,      a.qq = list(method = "tnorm",                   a.cipoly = list(fill = "light blue")),       a.respoi = list(size = 0.2),       a.hist = list(bins = 10))

Checking GAM simulated residuals or responses

Description

This function extracts the residuals or responses of a fitted GAM model, then it comparestheir distribution with that of model-based simulations.

Usage

check0D(  o,  type = "auto",  maxpo = 10000,  na.rm = TRUE,  trans = NULL,  useSim = TRUE)

Arguments

Value

An object of classc("plotSmooth", "gg").

Examples

# The variance of the response distribution changes along x2 library(mgcViz)n  <- 400x1 <- runif(n, -1, 1)x2 <- runif(n, -1, 1)dat <- data.frame("x1" = x1, "x2" = x2,                  "y" = sin(3*x1) + 0.5 * x2^2 + pmax(x2, 0.2) * rnorm(n))# Fit model with constant variance and perform posterior simulations (post = TRUE) # which take into account smoothing parameter uncertainty (unconditional = TRUE)b <- gamV(y ~ s(x1)+s(x2), data = dat,           aViz = list(nsim = 50, post = TRUE, unconditional = TRUE))# Histogram of simulated vs observed residuals: the latter are fat tailedcheck0D(b) + l_hist() + l_rug()# Histogram of simulated 4th central moment (~ kurtosis) of simulated residuals.# The vertical line is the 4th moment of the observed residualscheck0D(b, trans = function(.y) mean((.y - mean(.y))^4)) + l_dens1D() + l_vline() + l_rug()# Residuals look very fat tails, but the real problem here is the heteroscedasticity # which can be diagnosted using check1D(b, "x2") + l_gridCheck1D(sd)

Checking GAM residuals or responses along one covariate

Description

This function extracts the residuals of a fitted GAM model, and ordersthem according to the value of a single covariate. Then several visual residuals diagnosticscan be plotted by adding layers.

Usage

check1D(  o,  x,  type = "auto",  maxpo = 10000,  na.rm = TRUE,  trans = NULL,  useSim = TRUE)

Arguments

Value

The function will return an object of classc("plotSmooth", "gg"), unless argumentx is alist. In that case the function will return an object of classc("plotGam", "gg") containinga checking plot for each variable.

Examples

### Example 1: diagnosing heteroscedasticitylibrary(mgcViz);set.seed(4124)n <- 1e4x <- rnorm(n); y <- rnorm(n);# Residuals are heteroscedastic w.r.t. xob <- (x)^2 + (y)^2 + (0.2*abs(x) + 1)  * rnorm(n)b <- bam(ob ~ s(x,k=30) + s(y, k=30), discrete = TRUE)b <- getViz(b)# Look at residuals along "x"ck <- check1D(b, "x", type = "tnormal")# Can't see that muchck + l_dens(type = "cond", alpha = 0.8) + l_points() + l_rug(alpha = 0.2)# Some evidence of heteroscedasticityck + l_densCheck()# Compare observed residuals std dev with that of simulated data,# heteroscedasticity is clearly visibleb <- getViz(b, nsim = 50)check1D(b, "x") + l_gridCheck1D(gridFun = sd, showReps = TRUE)# This also works with factor or logical datafac <- sample(letters, n, replace = TRUE)logi <- sample(c(TRUE, FALSE), n, replace = TRUE)b <- bam(ob ~ s(x,k=30) + s(y, k=30) + fac + logi, discrete = TRUE)b <- getViz(b, nsim = 50)# Look along "fac"ck <- check1D(b, "fac") ck + l_points() + l_rug() ck + l_gridCheck1D(gridFun = sd)# Look along "logi"ck <- check1D(b, "logi") ck + l_points() + l_rug() ck + l_gridCheck1D(gridFun = sd)

Checking GAM residuals along two covariates

Description

This function extracts the residuals of a fitted GAM model, and plots themaccording to the values of two covariates. Then several visual residuals diagnosticscan be plotted by adding layers.

Usage

check2D(  o,  x1,  x2,  type = "auto",  maxpo = 10000,  na.rm = TRUE,  trans = NULL,  useSim = TRUE)

Arguments

Value

The function will return an object of classc("plotSmooth", "gg"), unless argumentsx1 and/orx2 are lists. If they are lists of the same length, then the function will return an object of classc("plotGam", "gg") containing a checking plot for each pair of variables. Ifx1 is a listandx2 is not specified, the function will return an object of classc("plotGam", "gg") containinga plot for each unique combination of the variables inx1.

Examples

## Not run: library(mgcViz);#### Example 1: Rosenbrock function# Simulate datan <- 1e4X <- data.frame("x1"=rnorm(n, 0.5, 0.5), "x2"=rnorm(n, 1.5, 1))X$y <- (1-X$x1)^2 + 100*(X$x2 - X$x1^2)^2 + rnorm(n, 0, 2)b <- bam(y ~ te(x1, x2, k = 5), data = X, discrete = TRUE)b <- getViz(b, nsim = 50)# Plot joint density of observed covariate x1 and x2check2D(b, x1 = "x1", x2 = "x2") + l_rug() + l_dens(type="joint", alpha=0.6) + l_points()# Look at how mean of residuals varies across x1 and x2check2D(b, x1 = "x1", x2 = "x2") + l_gridCheck2D() + l_points()# Can't see much in previous plot, let's zoom in central area, where most# data is. Here we can clearly see that the mean model is mispecifiedcheck2D(b, x1 = "x1", x2 = "x2") + l_gridCheck2D(bw = c(0.05, 0.1)) +                                   xlim(-1, 1) + ylim(0, 3)# Fit can be improved by increasing k in the bam() call#### Example 2: checking along factor variables# Simulate data where variance changes along factor variable "fac"n <- 1e4X <- data.frame("x1"=rnorm(n, 0.5, 0.5), "x2"=rnorm(n, 1.5, 1))X$fac <- as.factor( sample(letters, n, replace = TRUE) )X$fac2 <- as.factor( sample(c("F1", "F2", "F3", "F4", "F5"), n, replace = TRUE) )X$y <- (1-X$x1)^2 + 5*(X$x2 - X$x1^2)^2 + 0.1*as.numeric(X$fac) * rnorm(n, 0, 2)b <- bam(y ~ te(x1, x2, k = 5) + fac + fac2, data = X, discrete = TRUE)b <- getViz(b, nsim = 50)# Check standard deviation of residuals along covariates "x1" and "fac"a <- check2D(b, x1 = "x2", x2 = "fac")a + l_gridCheck2D(gridFun = sd) + l_rug() + l_points() # Points and rug are jittered by default, but we can over-ride thisa + l_rug(position = position_jitter(width = 0, height = 0)) +   l_points(position = position_jitter(width = 0, height = 0)) # Check standard deviation of residuals along the two factor variablesa <- check2D(b, x1 = "fac", x2 = "fac2")a + l_gridCheck2D(gridFun = sd, bw = c(1, 4)) + l_rug() + l_points() ## End(Not run)

Getting the CDF of a gam family

Description

Some methods implemented inmgcViz require the c.d.f. of the response distribution.This function takes a family object as input and returns the same object, but with the cdffunction added to the$cdf slot. Mainly for internal use.

Usage

fix.family.cdf(fam)

Arguments

Fit a GAMM or GAMM4 model and get a gamViz object

Description

These are wrappers that fit GAM models usingmgcv::gamm orgamm4::gamm4 andconvert them to agamViz object using thegetViz function.It is essentially a shortcut.

Usage

gamm4V(  formula,  random = NULL,  family = gaussian(),  data = list(),  REML = TRUE,  aGam = list(),  aViz = list(),  keepGAMObj = FALSE)gammV(  formula,  random = NULL,  family = gaussian(),  data = list(),  method = "REML",  aGam = list(),  aViz = list(),  keepGAMObj = FALSE)

Arguments

Details

WARNING: Model comparisons (e.g. withanova) should only be done using the mixed model part as described ingamm4::gamm4.Formgcv::gamm please refer to the original help file.

Value

An object of class "gamViz" which can, for instance, be plotted usingplot.gamViz. Also the object has the following additional elements:

• lme mixed model as inmgcv::gamm

• mer mixed model as ingamm4::gamm4

• gam a copy of the gamViz Object if settingkeepGAMObj = TRUE.

Examples

##### gam examplelibrary(mgcViz)# Simulate datadat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth## Now add 20 level random effect `fac'...dat$fac <- fac <- as.factor(sample(1:20,400,replace=TRUE))dat$y <- dat$y + model.matrix(~fac-1) %*% rnorm(20) * 0.5br <- gammV(y~s(x0)+x1+s(x2), data=dat,random=list(fac=~1))summary(br)plot(br)summary(br$lme)## Not run: ## gamm4::gamm4 examplebr4 <- gamm4V(y~s(x0)+x1+s(x2),data=dat,random=~(1|fac))summary(br4)plot(br4)summary(br4$mer)## End(Not run)

Convert gamViz object to gamObject

Description

Function for converting agamViz object to agamObject.It is essentially the inverse of thegetViz function.

Usage

getGam(o)

Arguments

Examples

library(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(1,n=1000,dist="normal",scale=2)b <- gam(y~s(x0)+s(x1, x2)+s(x3), data=dat, method="REML")a <- getViz(b)identical(b, getGam(a)) # Must be TRUE

Converting gam objects to gamViz objects

Description

This function convertsgam objects intogamViz objects,for whichmgcViz provides several plotting methods.

Usage

getViz(o, nsim = 0, post = FALSE, newdata, ...)

Arguments

Value

An object of classgamViz.

Examples

library(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(1,n=1000,dist="normal",scale=2)b <- gam(y~s(x0)+s(x1, x2)+s(x3), data=dat, method="REML")b <- getViz(b, nsim = 20)str(b$store$sim) # Simulated responses now stored hereplot(sm(b,1)) + l_fitLine() + l_ciLine() + l_rug() + l_points()plot(sm(b,2)) + l_rug() + l_fitRaster() + l_fitContour()

Plotting plotSmooth objects on a grid

Description

This is a wrapper forgridExtra::grid.arrange, which allows toplot severalplotSmooth objects on a grid.

Usage

gridPrint(...)

Arguments

Details

This function simply extracts theggplot objects contained in anyobject of classplotSmooth and passes them togridExtra::grid.arrange.

Examples

library(mgcViz)n  <- 1e3x1 <- rnorm(n)x2 <- rnorm(n)dat <- data.frame("x1" = x1, "x2" = x2,                  "y" = sin(x1) + 0.5 * x2^2 + pmax(x2, 0.2) * rnorm(n))b <- bam(y ~ s(x1)+s(x2), data = dat, method = "fREML", discrete = TRUE)b <- getViz(b)o1 <- plot( sm(b, 1) ) + l_fitLine() + l_ciLine() o2 <- plot( sm(b, 2) ) + l_fitLine() + l_ciLine() qpl <- qq(b)# All on one page, method 1:gridPrint(o1, o2, qpl, ncol = 2)# All on one page, method 2:gridPrint(grobs = list(o1, o2, qpl), ncol = 2)# Works also when some ggplot objects are presentgridPrint(o1, o2, qpl, ggplot(), ncol = 2)

Add boundaries to smooth effect plot

Description

This layer adds boundaries to a smooth effect plot.

Usage

l_bound(n = 200, ...)

Arguments

Value

An object of classgamLayer.

See Also

plot.sos.smooth

Adding confidence intervals to barplots

Description

This layer adds confidence intervals to barplots, such as those produced by factor effects GAM.

Usage

l_ciBar(level = 0.95, mul = NULL, ...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.ptermFactor for examples.

Adding confidence intervals to effect plot

Description

This layer adds confidence interval lines to smooth, random or parametric effect plots.

Usage

l_ciLine(level = 0.95, mul = NULL, ...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.1D,plot.ptermNumeric orplot.random.effect for examples.

Adding confidence band to effect plots

Description

This layer adds a polygon representing the confidence band of asmooth, random or parametric effect plots.

Usage

l_ciPoly(level = 0.95, mul = NULL, ...)

Arguments

Value

An object of classgamLayer

See Also

Seeplot.mgcv.smooth.1D,plot.ptermNumeric orplot.random.effect for examples.

Cluster and plot smooth effects

Description

This layers clusters several smooth effects and plots the cluster centers.

Usage

l_clusterLine(centers, cluFun = kmeans, a.clu = list(), ...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.fs.interaction.1D for examples.

Adding coordinate lines

Description

This layers adds coordinate contours to smooth effect plots. It is mainlyuseful for smooth-on-the-sphere plots.

Usage

l_coordContour(brLO = c(-9:9 * 20), brLA = c(-8:8 * 10), ...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.sos.smooth for examples.

Adding density estimate to a plot

Description

This layer adds a density estimate to a plot. It is mainly awrapper aroundggplot2::geom_density.

Usage

l_dens1D(...)

Arguments

Value

An object of classgamLayer.

See Also

Seecheck0D for examples.

Adding density estimate heatmap

Description

This layer adds a 2D density estimate heat-map to a plot.For 1D effect plots, it adds either the conditional density of the partialresiduals,p(r|x), or the joint densityp(r, x). For 2Deffect plots it adds eitherp(x1|x2) orp(x1, x2), wherex1 andx2 are the relevant covariates.

Usage

l_dens2D(type, n = c(50, 50), bw = NULL, tol = 1e-06, trans = sqrt, ...)l_dens(type, n = c(50, 50), bw = NULL, tol = 1e-06, trans = sqrt, ...)

Arguments

Details

The density function is estimated using the fast binned kernel density estimationmethods provided by theKernSmooth package, hence this function should beable to handle relatively large datasets (~ 10^6 observations).

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.1D,plot.mgcv.smooth.2D andcheck1D for examples.

Checking residuals conditional density

Description

This layer calculates and plots how the empirical conditional density ofthe residuals, r, differs from its theoretical or model-based counterpart,along a covariate, x.

Usage

l_densCheck(n = c(80, 80), bw = NULL, tol = 1e-06, dFun = NULL, ...)

Arguments

Details

This layer is mainly meant to work together with thecheck1D function. If check1D() is called withresidual type == "tunif" or "tnormal", thenl_densCheck compares the conditional distributionof the residuals with Unif(0, 1) or N(0, 1). By changing the distance functiondFunone could of course change both the distance metric and the reference distribution (see Examples below).

WARNING: if check1D() is called with type != "tunif" or "tnormal", then the default distance used by l_densCheck is
dFun <- function(.ed, .gr, .y){
d <- dnorm(.gr, 0, sd=sd(.y)) # sd=sd(.y) !!!
d <- sqrt(.ed)-sqrt(d)
return(sign(d)*abs(d)^(1/3))
}
so the residuals are standardized using their own std devsd(.y).Hencel_densCheck might not detect that the mean estimated varianceunder the fitted model is different from the residuals variance.Hence it is safer to use residual types "tunif" or "tnormal", or acustomized distance function dFun (see below for an example on how to do this).

Value

An object of classgamLayer.

Examples

library(mgcViz);# Dataset where variance increases linearly with x2, for x2 > 0.2n <- 1e3x1 <- rnorm(1e3)x2 <- rnorm(1e3)dat <- data.frame("x1"=x1,                   "x2"=x2, "y"=sin(x1) + 0.5*x2^2 + pmax(x2, 0.2)*rnorm(n))b <- gam(y ~ s(x1)+s(x2), data=dat)b <- getViz(b)# (Red) Blue indicates area where the empirical density # of the residuals is (lower) higher than it should be under # the model (residuals should be N(0, sigma) here).# Here there are clear signs of heteroscedasticity: # the conditional variance is is increasing for x2 > 0.2. check1D(b, "x2", type = "tnormal") + l_densCheck() + l_rug()# Suppose we want to compare the conditional density of the standardized residuals# not with a Gaussian, but with a Student-t density with 3 degree of freedom.# We could achieve this as follows:myDistance <- function(.ed, .gr, .y){  d <- dt(.gr / sd(.y), df = 3)  d <- abs( sqrt(.ed) - sqrt(d) ) # We are using absolute difference between sqrt-densities }check1D(b, "x2", type = "response") + l_densCheck(dFun = myDistance) + l_rug()# NB comparing with a Student density is not useful for this example, but it illustrates# how both the distance function and the reference density can be customized.

Adding barplot to effect plots

Description

This layer adds a barplot to an effect plots. Mainly useful for factor orbinary effect plots.

Usage

l_fitBar(a.aes = list(), ...)

Arguments

Value

an object of classgamLayer.

See Also

Seeplot.ptermFactor for examples.

Adding fitted effect contour lines

Description

This layer adds the contour lines corresponding to a fitted multidimensional effect.

Usage

l_fitContour(...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.2D,plot.mgcv.smooth.MD,plot.sos.smooth andplotSlice for examples.

Adding density strip of fitted effect

Description

This layer adds a conditional posterior density strip to 1D smooth effects plots.With the default colour scale, the opacity is proportional to the conditional density of the fittedeffects, under the usual Gaussian approximation the posterior.

Usage

l_fitDens(n = 50, level = 0.95, trans = identity, ...)

Arguments

Details

See Bowman (2018) for explanations about the advantages of density strips, relativeto plots including the mean fit + confidence intervals.

Value

An object of class gamLayer.

References

Bowman, D. W (2018). Graphics for uncertainty. Journal of the Royal Statistical Society: Series A.

Examples

library(mgcViz)set.seed(44)dat <- gamSim(1,n=400,dist="normal",scale=2)b <- gamV(y~s(x0)+x1+s(x2)+s(x3),data=dat)plot(sm(b, 1)) + l_fitDens() + l_fitLine()plot(pterm(b, 1)) + l_fitDens(trans = function(x) x^0.25) + l_fitLine()

Add fitted smooth effect curve

Description

This layer add lines representing a single or a group ofparametric or smooth 1D effects.

Usage

l_fitLine(...)

Arguments

Details

When used in conjuction withplot.fs.interaction.1D, which plots smooth effectsof typebs="fs", this function uses transparency to avoid over-plotting.This can be avoided by settingalpha = 1 in the call tol_fitLine.

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.1D,plot.ptermNumeric, orplot.fs.interaction.1D for examples.

Adding points representing the fitted effect

Description

This function adds points representing the fitted effect. Mainlyuseful for plotting factor effects.

Usage

l_fitPoints(...)

Arguments

Value

an object of classgamLayer.

See Also

Seeplot.ptermFactor for examples.

Adding raster representing the fitted effect

Description

This layer adds a raster or heat-map representing a fitted multidimensional effect.

Usage

l_fitRaster(pTrans = function(.p) 1, noiseup = FALSE, mul = 1, ...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.2D,plot.sos.smooth orplotSlice for examples.

Adding glyphs to 2D plots

Description

This layer adds glyphs or subplots to 2D plots. It is mainly meant tobe used withcheck2D and to produce residuals checks.

Usage

l_glyphs2D(  glyFun,  ggLay = "geom_points",  n = c(4, 4),  mapping = NULL,  data = NULL,  polar = FALSE,  height = ggplot2::rel(0.95),  width = ggplot2::rel(0.95),  y_scale = identity,  x_scale = identity,  ...)

Arguments

Value

An object of classgamLayer.

See Also

check2D.

Examples

library(mgcViz);set.seed(4124)n <- 1e4dat <- data.frame("x1" = rnorm(n), "x2" = rnorm(n))# Residuals are heteroscedastic w.r.t. x1dat$y <- (dat$x1)^2 + (dat$x2)^2 + (1*abs(dat$x1) + 1)  * rnorm(n)b <- bam(y ~ s(x1,k=30) + s(x2, k=30), data = dat, discrete = TRUE)b <- getViz(b)pl <- check2D(b, x1 = "x1", x2 = "x2", type = "tnormal") +   l_points(colour = "blue", alpha = 0.5)# Look at distributions of residuals across x1 and x2# Approach 1: using binned kernel density estimate# Colour indicates whether we have more that 50 obs in that binglyFun <- function(.d){  .r <- .d$z  .qq <- as.data.frame( density(.r)[c("x", "y")], n = 100 )  .qq$colour <- rep(ifelse(length(.r)>50, "black", "red"), nrow(.qq))  return( .qq )}pl + l_glyphs2D(glyFun = glyFun, ggLay = "geom_path", n = c(8, 8),                 mapping = aes(x=gx, y=gy, group = gid, colour = I(colour)),                  height=1.5, width = 1) # Approach 2: using binned worm-plots. These are simply rotated QQplots.# An horizontal plot indicates well specified residual model. # Increasing (decreasing) worm indicates over (under) dispersionglyFun <- function(.d){  n <- nrow(.d)  px <- qnorm( (1:n - 0.5)/(n) )  py <- sort( .d$z )  clr <- if(n > 50) { "black" } else { "red" }  clr <- rep(clr, n)  return( data.frame("x" = px, "y" = py - px, "colour" = clr))}pl + l_glyphs2D(glyFun = glyFun, ggLay = "geom_point", n = c(10, 10),                mapping = aes(x=gx, y=gy, group = gid, colour = I(colour)),                height=2, width = 1, size = 0.2)

Binning and checking GAM residuals

Description

This layer bins the residuals, r, according to the value of the correspondingcovariate, x. Then the residuals in each bin are summarized using ascalar-valued statistic. Confidence intervals for the statistic correspondingto each bin can be obtained by simulating residuals from the fitted GAMmodel, binning and summarizing them. Mainly useful in conjuction withcheck1D.

Usage

l_gridCheck1D(  gridFun = NULL,  n = 20,  level = 0.8,  stand = "none",  showReps = TRUE,  showObs = TRUE,  ...)

Arguments

Value

An object of classgamLayer

Examples

library(mgcViz);set.seed(4124)n <- 1e4x <- rnorm(n); y <- rnorm(n);# Residuals are heteroscedastic w.r.t. xob <- (x)^2 + (y)^2 + (0.2*abs(x) + 1)  * rnorm(n)b <- bam(ob ~ s(x,k=30) + s(y, k=30), discrete = TRUE)b <- getViz(b, nsim = 50)# Don't see much by looking at meancheck1D(b, "x") + l_gridCheck1D()# Heteroscedasticity clearly visible herecheck1D(b, "x") + l_gridCheck1D(gridFun = sd, stand = "sc") # <- we are scaling and centering# Last point on the right of the rug seems to indicate that a bin is missing.# It is not an error, only on observation falls in that bin, hence the# standard deviation is not defined there.

Binning and checking GAM residuals

Description

This layer bins the residuals, r, according to the value of the correspondingcovariates, x1 and x2. Then the residuals in each bin are summarized using ascalar-valued statistic. Confidence intervals for the statistic correspondingto each bin can be obtained by simulating residuals from the fitted GAMmodel, which are then binned and summarized. Mainly useful in conjuction withcheck2D.

Usage

l_gridCheck2D(gridFun = mean, bw = c(NA, NA), stand = TRUE, binFun = NULL, ...)

Arguments

Value

An object of classgamLayer

Examples

library(mgcViz);set.seed(4124)n <- 1e4x <- rnorm(n); y <- rnorm(n);# Residuals are heteroscedastic w.r.t. xob <- (x)^2 + (y)^2 + (1*abs(x) + 1)  * rnorm(n)b <- bam(ob ~ s(x,k=30) + s(y, k=30), discrete = TRUE)b <- getViz(b, nsim = 50)# Don't see much by looking at meancheck2D(b, "x", "y") + l_gridCheck2D(gridFun = mean, bw = c(0.4, 0.4))# Variance pattern along x-axis clearer nowcheck2D(b, "x", "y") + l_gridCheck2D(gridFun = sd, bw = c(0.4, 0.4))

Checking sign of residuals along one covariate

Description

This layer is mainly useful when checking quantile GAMs fitted using theqgampackage. The residuals, r, are binned according to the corresponding value of acovariate, x. Then the proportions of negative residuals within each bin are calculated, andcompared with the theoretical value,qu. Confidence intervals for the proportionof negative residuals can be derived using binomial quantiles (under an independenceassumption). To be used in conjuction withcheck1D.

Usage

l_gridQCheck1D(qu = NULL, n = 20, level = 0.8, ...)

Arguments

Value

An object of classgamLayer

Examples

# Simulate some datalibrary(mgcViz)set.seed(3841)dat <- gamSim(1,n=400,dist="normal",scale=2)dat$fac <- as.factor( sample(letters[1:8], nrow(dat), replace = TRUE) ) fit <- qgam(y~s(x1)+s(x2)+s(x3)+fac, data=dat, err = 0.05, qu = 0.4)fit <- getViz(fit)# "x0" effect is missing, but should be there. l_gridQCheck1D shows# that fraction of negative residuals is quite different from the theoretical 0.4# in several places along "x0".check1D(fit, dat$x0) + l_gridQCheck1D(qu = 0.4, n = 20)# The problem gets better if s(x0) is added to the model.# Works also with factor variablescheck1D(fit, "fac") + l_gridQCheck1D(qu = 0.4)

Binning and checking QGAM residuals

Description

This layer bins the residuals, r, according to the value of the correspondingcovariates, x1 and x2. Then we calculate the proportion of negative residualsin each bin, which should not deviate too much from the theoretical proportion (eg 0.5 ifwe fit the median). Mainly useful in conjuction withcheck2D.

Usage

l_gridQCheck2D(qu = NULL, bw = c(NA, NA), stand = TRUE, binFun = NULL, ...)

Arguments

Value

An object of classgamLayer

Examples

library(mgcViz);set.seed(4124)n <- 4e2 dat <- data.frame(x = rnorm(n), y = rnorm(n))# Simulate some data, residuals are heteroscedastic w.r.t. xdat$ob <- (dat$x)^2 + (dat$y)^2 + (0.2*abs(dat$x) + 1)  * rnorm(n)b <- qgamV(ob ~ x + s(y), qu = 0.3, data = dat)# We have a residual pattern along x (increase n above to # see the problem more clearly) check2D(b, "x", "y") + l_gridQCheck2D(qu = 0.3, bw = c(0.4, 0.4))# We need a smooth wrt x to make the pattern disappear## Not run: b1 <- qgamV(ob ~ s(x) + s(y), qu = 0.3, data = dat) check2D(b1, "x", "y") + l_gridQCheck2D(qu = 0.3, bw = c(0.4, 0.4))## End(Not run)

Adding histogram to a plot

Description

This layer adds a histogram to a plot. It is mainly awrapper aroundggplot2::geom_histogram.

Usage

l_hist(...)

Arguments

Value

An object of classgamLayer.

See Also

Seecheck0D for examples.

Add points to plot

Description

This layers add points to smooth, parametric or random effect plots.It can also be used to add points to the output ofcheck1Dandcheck2D. The meaning of the added points, which could representresiduals or covariate values, should be clear from context.

Usage

l_points(...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.1D,plot.mgcv.smooth.2D,check1D orcheck2D for examples.

Add polygons to effect plots

Description

This layers adds polygons to plots and it is mainly usefuls forplotting Markov random field smooths.

Usage

l_poly(...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.mrf.smooth for examples.

Adding contour of p-values

Description

This function adds contour lines proportional to the p-valueof a multidimensional smooth effects. It is useful for checkingwhere (across covariates x1 and x2) the fitted smooth is significantlydifferent from zero.

Usage

l_pvContour(pTrans = identity, ...)

Arguments

Value

An object of classgamLayer.

See Also

SeeplotDiff.mgcv.smooth.2D andplotDiff.sos.smooth for examples.

Adding raster or heat-map of p-values

Description

This function adds a raster or heat-map proportional to the p-valueof a multidimensional smooth effects. It is useful for checkingwhere (across covariates x1 and x2) the fitted smooth is significantlydifferent from zero.

Usage

l_pvRaster(pTrans = identity, ...)

Arguments

Value

An object of classgamLayer.

See Also

SeeplotDiff.mgcv.smooth.2D andplotDiff.sos.smooth for examples.

Adding rug to margins of a plot

Description

This layer adds a rug plot to the margins of a plot. It is mainly awrapper aroundggplot2::geom_rug. Notice that for factor effectsplots the rug is jittered by default.

Usage

l_rug(...)

Arguments

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.1D,plot.mgcv.smooth.2D orcheck1D for examples.

Add simulated smooth effect curves

Description

This layer adds curves representing smooth effects simulated from theposterior distribution.

Usage

l_simLine(...)

Arguments

Details

This function uses transparency to avoid over-plotting.This can be avoided by settingalpha = 1 in the call tol_simLine.

Value

An object of classgamLayer.

See Also

Seeplot.mgcv.smooth.1D for examples.

Adding vertical line to a plot

Description

This layer adds a vertical to a plot. It is mainly awrapper aroundggplot2::geom_vline.

Usage

l_vline(...)

Arguments

Value

An object of classgamLayer.

See Also

Seecheck0D for examples.

Lists available layers for plotSmooth objects

Description

This function takes as input an object of classplotSmooth andreturns the names of all the possible visual layers thatcould be used with that object.

Usage

listLayers(o)

Arguments

Value

A vector containing the names of the available layers.

Examples

library(mgcViz)n  <- 400x1 <- rnorm(n)x2 <- rnorm(n)dat <- data.frame("x1" = x1, "x2" = x2,                  "y" = sin(x1) + 0.5 * x2^2 + rnorm(n))b <- gam(y ~ x1+s(x2), data = dat, method = "REML")b <- getViz(b)# List layers available for parametric effect ploto <- plot( pterm(b, 1) )listLayers(o)# List layers available for smooth effect ploto <- plot( sm(b, 1) )listLayers(o)# List layers available for checking ploto <- check1D(b, x1)listLayers(o)

Fit multiple QGAM models and get a mgamViz object

Description

These are wrapper that fits multple QGAM models usingqgam::mqgam andconverts it to amgamViz object using thegetViz function.It is essentially a shortcut.

Usage

mqgamV(form, data, qu, lsig = NULL, err = NULL, aQgam = list(), aViz = list())

Arguments

Value

An object of class "mgamViz" which can, for instance, be plotted usingplot.mgamViz.

Examples

library(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(2,n=500,dist="normal",scale=0.25)$data# Fit GAM and get gamViz objectb <- mqgamV(y~s(x) + s(z) + I(x*z), data = dat, qu = c(0.25, 0.5, 0.75),            aQgam = list(argGam = list(select = TRUE)), aViz = list("nsim" = 0))# This is equivalent to doing# 1. Fit QGAM# b <- mqgam(y~s(x) + s(z) + I(x*z), data=dat, #            qu = c(0.25, 0.5, 0.75), argGam = list(select = TRUE))# 2. Convert to gamViz object# b <- getViz(b, nsim = 0)# Either way, we all effects by doingprint(plot(b, allTerms = TRUE), pages = 1)

Plot 1D Accumulated Local Effects (ALE)

Description

Plot 1D Accumulated Local Effects (ALE)

Usage

## S3 method for class 'ALE1D'plot(x, trans = identity, maxpo = 10000, nsim = 0, ...)

Arguments

Value

An objects of classplotSmooth.

Author(s)

Matteo Fasiolo and Christian Capezza, with some internal code having been adapted from the ALEPlotpackage of Dan Apley.

References

Apley, D.W., and Zhu, J, 2016. Visualizing the effects of predictor variables in blackbox supervised learning models. arXiv preprint arXiv:1612.08468.

Examples

library(mgcViz)# Here x1 and x2 are very correlated, but only # x1 has influence of the responseset.seed(4141)n <- 1000X <- rmvn(n, c(0, 0), matrix(c(1, 0.9, 0.9, 1), 2, 2))y <- X[ , 1] + 0.2 * X[ , 1]^2 + rnorm(n, 0, 0.8)dat <- data.frame(y = y, x1 = X[ , 1], x2 = X[ , 2])fit <- gam(y ~ te(x1, x2), data = dat)# Marginal plot suggests that E(y) depends on x2, but# this is due to the correlation between x1 and x2...plot(dat$x2, fit$fitted.values)# ... in fact ALE effect of x2 is flat ...plot(ALE(fit, "x2")) + l_ciPoly() + l_fitLine() + l_rug()# ... while ALE effect of x1 is strong.plot(ALE(fit, "x1", center = 2), nsim = 20) +   l_simLine() + l_fitLine()

Plotting one dimensional smooth factor interactions

Description

This method should be used to plot smooth effectsof class"fs.interaction.1D", that is smooth constructedusing the basisbs="tp". Seemgcv::s.

Usage

## S3 method for class 'fs.interaction.1D'plot(x, n = 100, xlim = NULL, trans = identity, ...)

Arguments

Value

An object of classc("plotSmooth", "gg").

Examples

library(mgcViz)set.seed(0)## simulate data...f0 <- function(x) 2 * sin(pi * x)f1 <- function(x, a = 2, b = -1) exp(a * x) + bf2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 *  (10 * x)^3 * (1 - x)^10n <- 500; nf <- 25fac <- sample(1:nf, n, replace = TRUE)x0 <- runif(n); x1 <- runif(n); x2 <- runif(n)a <- rnorm(nf) * .2 + 2; b <- rnorm(nf) * .5f <- f0(x0) + f1(x1, a[fac], b[fac]) + f2(x2)fac <- factor(fac)y <- f + rnorm(n) * 2## so response depends on global smooths of x0 and## x2, and a smooth of x1 for each level of fac.## fit model (note p-values not available when fit## using gamm)...bm <- gamm(y ~ s(x0)+ s(x1, fac, bs = "fs", k = 5) + s(x2, k = 20))v <- getViz(bm$gam)# Plot with fitted effects and changing title plot(sm(v, 2)) + l_fitLine(alpha = 0.6) + labs(title = "Smooth factor interactions")# Changing plotting limitsplot(sm(v, 2)) + l_fitLine() + ylim(-0.5, 0.5) + xlim(0.25, 0.75)# Change line type and remove legendplot(sm(v, 2)) + l_fitLine(size = 1.3, linetype="dotted") +                  theme(legend.position="none")                 # Clustering smooth effects in 3 groupsplot(sm(v, 2)) + l_fitLine(colour = "grey") +                  l_clusterLine(centers = 3, a.clu = list(nstart = 100))

Basic GAM plotting

Description

This function is themgcViz equivalent ofplot.gam. It is the workhorse of themgcViz package, and allows plotting (almost) any type of smooth,parametric or random effects. It is basically a wrapper around plottingmethods that are specific to individual smooth effect classes (such asplot.mgcv.smooth.1D andplot.random.effect).

Usage

## S3 method for class 'gamViz'plot(x, n = 100, n2 = 40, select = NULL, allTerms = FALSE, ...)

Arguments

Value

An object of classc("plotGam", "gg").

Examples

library(mgcViz)######## Basic example# Simulate some data and fit modelset.seed(2)  dat <- gamSim(1,n=1e3,dist="normal",scale=2)b <- bam(y~s(x0)+s(x1, x2)+s(x3), data=dat)b <- getViz(b)# Default smooth effect plottingprint(plot(b), ask = FALSE)# Now on one page and with out title on the second plotprint(plot(b) + labs(title = NULL), pages = 1) # So far we used default layers, added in the printing phase, but# we might want to specify our own layers. Here we is how to do itpl <- plot(b) + l_points() + l_fitLine(linetype = 3) + l_fitContour() +   l_ciLine(colour = 2) + theme_get() + labs(title = NULL)print(pl, pages = 1)# We might want to plot only the first smoothplot(b, select = 1) + l_dens(type = "cond") + l_fitLine() + l_ciLine()## Not run: ######## Example with "by variable" smooth effect# Simulate data and fit modeldat <- gamSim(4)b <- gam(y ~ fac+s(x2,by=fac)+s(x0),data=dat)b <- getViz(b)# print() only needed because we want to plot on a single pageprint(plot(b), pages = 1) print(plot(b, allTerms = TRUE), pages = 1) # Including also parametric effect######## Example with 3D smooth effect which cannot be plotted# Simulate data and fit modeln <- 5e3x <- rnorm(n); y <- rnorm(n); z <- rnorm(n); z2 <- rnorm(n)ob <- (x-z)^2 + (y-z)^2 + z2^3 + rnorm(n)b1 <- bam(ob ~ s(x, y, z) + s(z2), discrete = TRUE)b1 <- getViz(b1)# Only second effect get plottedplot(b1)# In fact this does not plot anythingplot(b1, select = 1)# For plotting effects with more than 2D, one we need specific method. # See ?plot.mgcv.smooth.MD######## Examples about plotting parametric effects# 1 Gaussian GAMset.seed(3)dat <- gamSim(1,n=2500,dist="normal",scale=20)dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) ) dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) ) bs <- "cr"; k <- 12b <- bam(y ~ x0 + x1 + I(x1^2) + s(x2,bs=bs,k=k) + fac + x3:fac + I(x1*x2) + logi +            s(x3, bs=bs),data=dat, discrete = TRUE)b <- getViz(b)# All effects in one page. Notably 'x3:fac' is missing: we have no methods# for plotting second order effects.print(plot(b, allTerms = TRUE), pages = 1)# Plotting only parametric effectsprint(plot(b, select = 3:9), pages = 1)# 2 GAMLSS Gaussian modellibrary(mgcv);library(MASS)mcycle$fac <- as.factor( sample(c("z", "k", "a", "f"), nrow(mcycle), replace = TRUE) ) b <- gam(list(accel~times + I(times^2) + s(times,k=10), ~ times + fac + s(times)),          data=mcycle,family=gaulss())b <- getViz(b)# All effects on one page: effect of second linear predictor end with '.1'print(plot(b, allTerms = TRUE), pages = 1)## End(Not run)

Plotting multiple quantile GAMs

Description

This function is similar toplot.gamViz, but it is usedto plot multiple quantile GAM models fitted usingmqgamV ormqgam.It allows plotting standards 1D and 2D smooths, and parametric effects,It is basically a wrapper around plotting methods that are specific toindividual smooth effect classes (such asplot.multi.mgcv.smooth.1D).

Usage

## S3 method for class 'mgamViz'plot(x, n = 100, n2 = 40, select = NULL, allTerms = FALSE, ...)

Arguments

Value

An object of classc("plotGam", "gg").

Examples

library(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(1,n=500,dist="normal",scale=2)dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) )dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) )# Fit GAM and get gamViz objectfit <- mqgamV(y ~ fac + s(x0) + s(x1, x2) + x3 + logi, data = dat,               qu = c(0.2, 0.4, 0.6, 0.8))print(plot(fit, select = 1:4, allTerms = T), pages = 1)## Not run: # Example where we are fitting the same model to different datasets, but# plotting the estimate effects together dat <- list()for(ii in 1:4){  # Simulate 4 datasets, we are adding 2 factor variables "fac" and "ref" just  # for illustrating the plotting method (the two factors have no effect on y)  n <- 1000  dat[[ii]] <- gamSim(1,n=n,dist="normal",scale=2)  dat[[ii]]$fac <- as.factor( sample(c("A1", "A2", "A3"), n, replace = TRUE) )  dat[[ii]]$ref <- as.factor( sample(letters[1:10], n, replace = TRUE) )}# Estimating model on each datasetmods <- list()for(ii in 1:4){  mods[[ii]] <- gamV(y~s(x0)+s(x1, x2)+x3+fac+s(ref, bs = "re"), data = dat[[ii]])}# Names will be used to identify the four models we have fitted names(mods) <- c("M1", "M2", "M3", "M4")# Plotting on the same plotsprint(plot.mgamViz(mods, allTerms = TRUE), pages = 1)## End(Not run)

Plotting one dimensional smooth effects

Description

Plotting method for one dimensional smooth effects.

Usage

## S3 method for class 'mgcv.smooth.1D'plot(  x,  n = 100,  xlim = NULL,  maxpo = 10000,  trans = identity,  unconditional = FALSE,  seWithMean = FALSE,  nsim = 0,  ...)## S3 method for class 'multi.mgcv.smooth.1D'plot(  x,  n = 100,  xlim = NULL,  maxpo = 10000,  trans = identity,  unconditional = FALSE,  seWithMean = FALSE,  asFact = NULL,  ...)

Arguments

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

library(mgcViz)n  <- 1e3x1 <- rnorm(n)x2 <- rnorm(n)dat <- data.frame("x1" = x1, "x2" = x2,                  "y" = sin(x1) + 0.5 * x2^2 + pmax(x2, 0.2) * rnorm(n))b <- bamV(y ~ s(x1)+s(x2), data = dat, method = "fREML", aGam = list(discrete = TRUE))o <- plot(sm(b, 1), nsim = 50) # 50 posterior simulations ## Not run: # Plot with fitted effect + posterior simulations + rug on x axis( o <- o + l_simLine() + l_fitLine(colour = "red") +        l_rug(alpha = 0.8) )# Add CI lines at 1*sigma and 5*sigma( o <- o + l_ciLine(mul = 1) + l_ciLine(mul = 5, colour = "blue", linetype = 2) )# Add partial residuals and change theme( o + l_points(shape = 19, size = 1, alpha = 0.2) + theme_classic() )# Get second effect ploto2 <- plot( sm(b, 2) )# Plot it with polygon for partial residualso2 + l_ciPoly(mul = 5, fill = "light blue") +   l_fitLine(linetype = 2, colour = "red")# Plot is with conditional density of partial residualso2 + l_dens(type = "cond", alpha = 0.9)  +   l_fitLine(linetype = 2, colour = "red")  ######### Quantile GAM example######### Fit modelb <- mqgamV(y ~ s(x1) + s(x2), qu = c(0.2, 0.5, 0.8), data = dat)plot(sm(b, 1)) + l_fitLine(linetype = 2) + l_rug(colour = "blue")## End(Not run)

Plotting two dimensional smooth effects

Description

Plotting method for two dimensional smooth effects.

Usage

## S3 method for class 'mgcv.smooth.2D'plot(  x,  n = 40,  xlim = NULL,  ylim = NULL,  maxpo = 10000,  too.far = 0.1,  trans = identity,  seWithMean = FALSE,  unconditional = FALSE,  ...)## S3 method for class 'multi.mgcv.smooth.2D'plot(  x,  n = 30,  xlim = NULL,  ylim = NULL,  maxpo = 10000,  too.far = 0.1,  trans = identity,  seWithMean = FALSE,  unconditional = FALSE,  a.facet = list(),  ...)

Arguments

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

library(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(1, n = 700, dist = "normal", scale = 2)b <- gam(y ~ s(x0) + s(x1, x2) + s(x3), data = dat, method = "REML")b <- getViz(b)# Plot 2D effect with noised-up raster, contour and rug for design points # Opacity is proportional to the significance of the effectplot(sm(b, 2)) + l_fitRaster(pTrans = zto1(0.05, 2, 0.1), noiseup = TRUE) +   l_rug() + l_fitContour()  # Plot contour of effect joint density of design pointsplot(sm(b, 2)) + l_dens(type = "joint") + l_points() + l_fitContour() +   coord_cartesian(expand = FALSE) # Fill the plot  #### Quantile GAM example###b <- mqgamV(y ~ s(x0) + s(x1, x2) + s(x3), qu = c(0.3, 0.7), data = dat)plot(sm(b, 2)) + l_fitRaster(noiseup = TRUE) + l_fitContour(colour = 2)

Plotting slice of higher-dimensional smooth effects

Description

This function plots a 2D slice of a higher-dimensional smooth effects.

Usage

## S3 method for class 'mgcv.smooth.MD'plot(  x,  fix,  n = 40,  xlim = NULL,  ylim = NULL,  maxpo = 10000,  too.far = c(0.1, NA),  trans = identity,  seWithMean = FALSE,  unconditional = FALSE,  ...)

Arguments

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

## Not run: ## 3D examplelibrary(mgcViz)n <- 1e3x <- rnorm(n); y <- rnorm(n); z <- rnorm(n)ob <- (x-z)^2 + (y-z)^2 + rnorm(n)b <- gam(ob ~ s(x, y, z))b <- getViz(b)# Plot one 2D sliceplot( sm(b, 1), fix = c("z"=0) ) + l_fitRaster(noiseup = TRUE, mul = 3) +   l_fitContour(linetype = 2) + l_points(shape =  2)## 4Dn <- 5e3x <- rnorm(n); y <- rnorm(n); z <- rnorm(n); z2 <- rnorm(n)ob <- (x-z)^2 + (y-z)^2 + z2^3 + rnorm(n)b1 <- bam(ob ~ s(x, y, z, z2), discrete = TRUE)b1 <- getViz(b1)# Plot one 2D sliceplot(sm(b1, 1), fix = c("z"=0, "z2"=1)) + l_fitRaster() + l_fitContour()## End(Not run)

Plotting Markov random field smooths

Description

This is the plotting method for Markov random field smooths.

Usage

## S3 method for class 'mrf.smooth'plot(x, trans = identity, seWithMean = FALSE, unconditional = FALSE, ...)

Arguments

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

library(mgcViz)## Load Columbus Ohio crime data (see ?columbus for details and credits)data(columb)       ## data framedata(columb.polys) ## district shapes listxt <- list(polys=columb.polys) ## neighbourhood structure info for MRFpar(mfrow=c(2,2))## First a full rank MRF...b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")b <- getViz(b) # Manual plotplot(sm(b, 1)) + l_poly(colour = 2) +   scale_fill_gradientn(colours = heat.colors(50))# Default plot plot(b)

Plotting factor or logical parametric effects

Description

These are the plotting methods for parametric factor or logical effects.

Usage

## S3 method for class 'multi.ptermFactor'plot(x, a.facet = list(), asFact = TRUE, ...)## S3 method for class 'multi.ptermLogical'plot(x, ...)## S3 method for class 'ptermFactor'plot(x, maxpo = 10000, trans = identity, ...)## S3 method for class 'ptermLogical'plot(x, maxpo = 10000, trans = identity, ...)

Arguments

Value

An object of classplotSmooth.

Examples

# Simulate data and fit GAMset.seed(3)dat <- gamSim(1,n=2000,dist="normal",scale=20)dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) )dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) )bs <- "cr"; k <- 12b <- gam(y~fac + s(x0) + s(x1) + s(x2) + s(x3) + logi, data=dat)o <- getViz(b, nsim = 0)# Extract factor terms and plot itpt <- pterm(o, 1)plot(pt) + l_ciBar() + l_fitPoints(colour = 2) + l_rug(alpha = 0.2)# Use barplot instead of pointspt <- pterm(o, 1)plot(pt) + l_fitBar() + l_ciBar()# Same with binary variblept <- pterm(o, 2)plot(pt) + l_fitPoints() + l_ciBar()

Plotting random effects

Description

This is the plotting method for random effects (simple random intercepts).

Usage

## S3 method for class 'multi.random.effect'plot(x, trans = identity, ...)## S3 method for class 'random.effect'plot(x, trans = identity, ...)

Arguments

Value

An object of classplotSmooth.

Examples

library(mgcViz)b <- gam(travel~s(Rail,bs="re"), data=Rail, method="REML")b <- getViz(b)plot(sm(b, 1)) + l_fitLine(colour = 2, linetype = 2) + l_points() +   l_ciLine(colour = 4, linetype = 3)plot(sm(b, 1)) + l_ciPoly() + l_points()# Defaultplot(b)#### Quantile GAM version###b <- mqgamV(travel~s(Rail,bs="re"), data=as.data.frame(Rail), qu = c(0.2, 0.4, 0.6, 0.8))plot(sm(b, 1)) + l_ciPoly() + l_points()# Defaultplot(b)

Plotting one dimensional nested effects

Description

This method should be used to plot smooth effectsof class"si.smooth.1D".

Usage

## S3 method for class 'nested1D'plot(  x,  inner = FALSE,  n = 100,  xlim = NULL,  ylim = NULL,  maxpo = 10000,  trans = identity,  ...)

Arguments

Value

An object of classc("plotSmooth", "gg").

Plotting parametric interactions

Description

This function is here only to deal with parametric interactions (eg x0:fact), whichcannot be plotted at the moment.

Usage

## S3 method for class 'ptermInteraction'plot(x, ...)## S3 method for class 'multi.ptermInteraction'plot(x, ...)

Arguments

Value

Currently it returnsNULL.

Plotting numeric parametric effects

Description

This is the plotting method for parametric numerical effects.

Usage

## S3 method for class 'ptermMatrixNumeric'plot(x, n = 100, xlim = NULL, trans = identity, ...)## S3 method for class 'multi.ptermNumeric'plot(x, ...)## S3 method for class 'ptermNumeric'plot(x, n = 100, xlim = NULL, maxpo = 10000, trans = identity, ...)

Arguments

Value

An object of classplotSmooth.

Examples

# Simulate data and fit GAMset.seed(3)dat <- gamSim(1,n=2000,dist="normal",scale=20)bs <- "cr"; k <- 12b <- gam(y ~  x0 + x1 + I(x1^2) + s(x2,bs=bs,k=k) +               I(x1*x2) + s(x3, bs=bs), data=dat)o <- getViz(b, nsim = 0)# Extract first terms and plot itpt <- pterm(o, 1)plot(pt, n = 60) + l_ciPoly() + l_fitLine() + l_ciLine()# Extract I(x1^2) terms and plot it with partial residualspt <- pterm(o, 3)plot(pt, n = 60) + l_ciPoly() + l_fitLine() + l_ciLine() + l_points()

Plotting smooths on the sphere

Description

This is the plotting method for smooth effects on the sphere.

Usage

## S3 method for class 'sos.smooth'plot(  x,  n = 40,  xlim = NULL,  ylim = NULL,  maxpo = 10000,  too.far = 0.1,  phi = 30,  theta = 30,  trans = identity,  scheme = 0,  seWithMean = FALSE,  unconditional = FALSE,  ...)

Arguments

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

library(mgcViz)set.seed(0)n <- 400f <- function(la,lo) { ## a test function...  sin(lo)*cos(la-.3)}## generate with uniform density on sphere...  lo <- runif(n)*2*pi-pi ## longitudela <- runif(3*n)*pi-pi/2ind <- runif(3*n)<=cos(la)la <- la[ind];la <- la[1:n]ff <- f(la,lo)y <- ff + rnorm(n)*.2 ## test data## generate data for plotting truth...lam <- seq(-pi/2,pi/2,length=30)lom <- seq(-pi,pi,length=60)gr <- expand.grid(la=lam,lo=lom)fz <- f(gr$la,gr$lo)zm <- matrix(fz,30,60)require(mgcv)dat <- data.frame(la = la *180/pi,lo = lo *180/pi,y=y)## fit spline on sphere model...bp <- gam(y~s(la,lo,bs="sos",k=60),data=dat)bp <- getViz(bp)# Plot on sphereplot(sm(bp, 1), scheme=0) + l_fitRaster() + l_fitContour() +    l_points(shape = 19) + l_rug() + l_coordContour() + l_bound()  # Plotting as in standard 2D plotsplot(sm(bp, 1), scheme=1) + l_fitRaster() + l_fitContour() +            l_points(shape = 19) + l_rug()

Generic plotting of differences

Description

Generic function for plotting differences between objects. Usefulmainly for plotting the differences between two smooth effects.

Usage

plotDiff(...)

Arguments

See Also

plotDiff.mgcv.smooth.1D,plotDiff.mgcv.smooth.2D,plotDiff.sos.smooth

Plotting differences between two 1D smooth effects

Description

This method can be used to plot the difference between two 1Dsmooth effects. Mainly meant to be used with by-factor smooths.

Usage

## S3 method for class 'mgcv.smooth.1D'plotDiff(s1, s2, n = 100, trans = identity, unconditional = FALSE, ...)

Arguments

Details

Let sd be the difference between the fitted smooths, that is: sd = s1 - s2.sd is a vector of length n, and its covariance matrix isCov(sd) = X1\where: X1 (X2) and Sig11 (Sig22) are the design matrix and the covariance matrixof the coefficients of s1 (s2), while Sig12 is the cross-covariance matrix betweenthe coefficients of s1 and s2. To get the confidence intervals we need only diag(Cov(sd)),which here is calculated efficiently (without computing the whole of Cov(sd)).

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

# Simulate data and add factors uncorrelated to the responselibrary(mgcViz)set.seed(6898)dat <- gamSim(1,n=1500,dist="normal",scale=20)dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) ) dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) ) bs <- "cr"; k <- 12b <- gam(y ~ s(x2,bs=bs,by = fac), data=dat)o <- getViz(b, nsim = 0)# Extract the smooths correspoding to "A1" and "A2" and plot their differences# with credible intervalsplotDiff(s1 = sm(o, 1), s2 = sm(o, 2)) + l_ciPoly() +   l_fitLine() + geom_hline(yintercept = 0, linetype = 2)

Plotting differences between two 2D smooth effects

Description

This method can be used to plot the difference between two 2Dsmooth effects. Mainly meant to be used with by-factor smooths.

Usage

## S3 method for class 'mgcv.smooth.2D'plotDiff(  s1,  s2,  n = 40,  too.far = 0.1,  trans = identity,  unconditional = FALSE,  ...)

Arguments

Details

Let sd be the difference between the fitted smooths, that is: sd = s1 - s2.sd is a vector of length n, and its covariance matrix isCov(sd) = X1\where: X1 (X2) and Sig11 (Sig22) are the design matrix and the covariance matrixof the coefficients of s1 (s2), while Sig12 is the cross-covariance matrix betweenthe coefficients of s1 and s2. To get the confidence intervals we need only diag(Cov(sd)),which here is calculated efficiently (without computing the whole of Cov(sd)).

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

# Simulate data and add factors uncorrelated to the responselibrary(mgcViz)set.seed(235)dat <- gamSim(1,n=1500,dist="normal",scale=20)dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) ) dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) ) bs <- "cr"; k <- 12b <- gam(y ~ s(x2, x1, by = fac), data=dat)o <- getViz(b, nsim = 0)# Extract the smooths correspoding to "A1" and "A2" and plot their differencepl <- plotDiff(s1 = sm(o, 1), s2 = sm(o, 2))pl + l_fitRaster() + l_fitContour()# Plot p-values for differences between the two smoothspl + l_pvRaster() + l_pvContour()

Plotting differences between two smooths on the sphere

Description

This method can be used to plot the difference between two smootheffects on the sphere. Mainly meant to be used with by-factor smooths.

Usage

## S3 method for class 'sos.smooth'plotDiff(  s1,  s2,  n = 40,  too.far = 0.1,  phi = 30,  theta = 30,  scheme = 0,  trans = identity,  unconditional = FALSE,  ...)

Arguments

Details

Let sd be the difference between the fitted smooths, that is: sd = s1 - s2.sd is a vector of length n, and its covariance matrix isCov(sd) = X1\where: X1 (X2) and Sig11 (Sig22) are the design matrix and the covariance matrixof the coefficients of s1 (s2), while Sig12 is the cross-covariance matrix betweenthe coefficients of s1 and s2. To get the confidence intervals we need only diag(Cov(sd)),which here is calculated efficiently (without computing the whole of Cov(sd)).

Value

An objects of classplotSmooth.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

## Not run: #### 1) Simulate data and add factors uncorrelated to the responselibrary(mgcViz)set.seed(0)n <- 500f <- function(la,lo) { ## a test function...  sin(lo)*cos(la-.3)}## generate with uniform density on sphere...  lo <- runif(n)*2*pi-pi ## longitudela <- runif(3*n)*pi-pi/2ind <- runif(3*n)<=cos(la)la <- la[ind];la <- la[1:n]ff <- f(la,lo)y <- ff + rnorm(n)*.2 ## test data## generate data for plotting truth...lam <- seq(-pi/2,pi/2,length=30)lom <- seq(-pi,pi,length=60)gr <- expand.grid(la=lam,lo=lom)fz <- f(gr$la,gr$lo)zm <- matrix(fz,30,60)dat <- data.frame(la = la *180/pi,lo = lo *180/pi,y=y)dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) ) #### 2) fit spline on sphere model...bp <- gam(y~s(la,lo,bs="sos",k=60, by = fac),data=dat)bp <- getViz(bp)# Extract the smooths correspoding to "A1" and "A3" and plot their difference# Using scheme = 0 pl0 <- plotDiff(s1 = sm(bp, 1), s2 = sm(bp, 3))pl0 + l_fitRaster() + l_fitContour() + l_coordContour() + l_bound()# Plot p-values for significance of differencespl0 + l_pvRaster() + l_pvContour(breaks=c(0.05, 0.1, 0.2, 0.3, 0.5))# Using scheme = 1pl1 <- plotDiff(s1 = sm(bp, 1), s2 = sm(bp, 2), scheme = 1) pl1 + l_fitRaster() + l_fitContour()# Plot p-values for significance of differencespl1 + l_pvRaster() + l_pvContour(breaks=c(0.05, 0.1, 0.2, 0.3, 0.5))## End(Not run)

Generic RGL plotting function

Description

Generic function for producing an interactive RGL plot.

Usage

plotRGL(x, ...)

Arguments

See Also

plotRGL.mgcv.smooth.2D,plotRGL.mgcv.smooth.MD

Visualizing 2D smooth effects in 3D

Description

This method plots an interactive 3D representation of a 2D smooth effect, usingthe rgl package.

Usage

## S3 method for class 'mgcv.smooth.2D'plotRGL(  x,  se = TRUE,  n = 40,  residuals = FALSE,  type = "auto",  maxpo = 1000,  too.far = 0,  xlab = NULL,  ylab = NULL,  main = NULL,  xlim = NULL,  ylim = NULL,  se.mult = 1,  trans = identity,  seWithMean = FALSE,  unconditional = FALSE,  ...)

Arguments

Value

ReturnsNULL invisibly.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

# Example 1: taken from ?mgcv::te, shows how tensor pruduct deals nicely with # badly scaled covariates (range of x 5% of range of z )library(mgcViz)# Simulate some datatest1 <- function(x,z,sx=0.3,sz=0.4) {   x <- x*20  (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+                 0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))}n <- 500old.par <- par(mfrow=c(2,2))x <- runif(n)/20;z <- runif(n);xs <- seq(0,1,length=30)/20;zs <- seq(0,1,length=30)pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30)))truth <- matrix(test1(pr$x,pr$z),30,30)f <- test1(x,z)y <- f + rnorm(n)*0.2# Fit with t.p.r.s. basis and plotb1 <- gam(y~s(x,z))plotRGL(sm(getViz(b1), 1))# Need to load rgl at this point## Not run: library(rgl)rgl.close() # Close# Fit with tensor products basis and plot (with residuals)b2 <- gam(y~te(x,z))plotRGL(sm(getViz(b2), 1), residuals = TRUE)# We can still work on the plot, for instance change the aspect ratioaspect3d(1, 2, 1)rgl.close() # Close## End(Not run)

Visualizing a 2D slice of a smooth effects in 3D

Description

This method plots an interactive 3D representation of a 2-dimensionalslice of a multi-dimensional smooth effect, using the rgl package.

Usage

## S3 method for class 'mgcv.smooth.MD'plotRGL(  x,  fix,  se = TRUE,  n = 40,  residuals = FALSE,  type = "auto",  maxpo = 1000,  too.far = c(0, NA),  xlab = NULL,  ylab = NULL,  main = NULL,  xlim = NULL,  ylim = NULL,  se.mult = 1,  trans = identity,  seWithMean = FALSE,  unconditional = FALSE,  ...)

Arguments

Value

ReturnsNULL invisibly.

References

Marra, G and S.N. Wood (2012) Coverage Properties of Confidence Intervals forGeneralized Additive Model Components. Scandinavian Journal of Statistics.

Examples

# Example 1: taken from ?mgcv::te, shows how tensor pruduct deals nicely with # badly scaled covariates (range of x 5% of range of z )library(mgcViz)n <- 1e3x <- rnorm(n); y <- rnorm(n); z <- rnorm(n)ob <- (x-z)^2 + (y-z)^2 + rnorm(n)b <- gam(ob ~ s(x, y, z))v <- getViz(b)plotRGL(sm(v, 1), fix = c("z" = 0))# Need to load rgl at this point## Not run: library(rgl)rgl.close() # CloseplotRGL(sm(v, 1), fix = c("z" = 1), residuals = TRUE)# We can still work on the plot, for instance change the aspect ratioaspect3d(1, 2, 1)rgl.close() # Close## End(Not run)

Plotting sequence of slices of 2D smooth effect

Description

This function allows to slice a multi-dimensional (D > 2) smooth effect,and to plot the resulting sequence of 2D slices in an array of plots.

Usage

plotSlice(x, fix, a.facet = list(), ...)

Arguments

Value

An objects of classplotSmooth.

Examples

## Not run: ### Example 1: plotting slices of 3D smooth# Simulate data and fit GAMlibrary(mgcViz)n <- 1e3x <- rnorm(n); y <- rnorm(n); z <- rnorm(n)ob <- (x-z)^2 + (y-z)^2 + rnorm(n)b <- gam(ob ~ s(x, y, z))v <- getViz(b)# Get plot of slices and add layerspl <- plotSlice(x = sm(v, 1),                 fix = list("z" = seq(-2, 2, length.out = 9)))pl + l_fitRaster() + l_fitContour() + l_points() + l_rug()# Over-ride default layoutpl <- plotSlice(x = sm(v, 1),                 fix = list("z" = seq(-2, 2, length.out = 9)),                 a.facet = list(nrow = 2))pl + l_fitRaster() + l_fitContour() + theme(panel.spacing = unit(0.5, "lines"))### Example 2: plotting slices of 4D smooth# Simulate data and fit GAMn <- 5e3x <- rnorm(n); y <- rnorm(n); z <- rnorm(n); z2 <- rnorm(n)ob <- (x-z)^2 + (y-z)^2 + z2^3 + rnorm(n)b <- bam(ob ~ s(x, y, z, z2), discrete = TRUE)v <- getViz(b)# Plot slices across "z" and "x"pl <- plotSlice(x = sm(v, 1),                 fix = list("z" = seq(-2, 2, length.out = 3), "x" = c(-1, 0, 1)))pl + l_fitRaster() + l_fitContour() + l_points() + l_rug()# Plot slices across "x", keeping "z" fixedpl <- plotSlice(x = sm(v, 1),                 fix = list("z" = 0, "x" = seq(-3, 3, length.out = 9)))pl + l_fitRaster() + l_fitContour() + l_points() + l_rug()## End(Not run)

Posterior simulation from a GAM object

Description

This method can be used to simulate vectors of responses from theGaussian posterior approximation of a gamObject.

Usage

postSim(  o,  nsim,  newdata,  trans = NULL,  method = "auto",  w = NULL,  offset = NULL,  savePar = FALSE,  ...)

Arguments

Value

IfsavePar == FALSE the function will return a matrix where each column is a vector ofsimulated responses or a transformed version of it. IfsavePar == TRUE it will returna list where the$simY entry will contain the simulated responses and$simBetathe simulated parameters.

Examples

library(mgcViz)library(MASS)b <- gam(accel~s(times, k=20), data=mcycle)# Simulate list of 10 vectors of responses from posterior, taking into# account smoothing parameters uncertainty (see ?vcov.gam)n <- 10sim <- postSim(o = b, nsim = n, unconditional = TRUE)# Posterior simulations in grey and data in redplot(rep(mcycle$times, n), as.vector(sim), col = "grey",      ylab = "Acceleration", xlab = "Times")points(mcycle$times, mcycle$accel, col = 2)# There is clear disagreement between simulations' and data's # conditional variance, which can be solved using flexible GAMLSS model:b <- gam(list(accel~s(times, k=20), ~s(times)), data=mcycle, family = gaulss)sim <- postSim(o = b, nsim = n)plot(rep(mcycle$times, n), as.vector(sim), col = "grey",      ylab = "Acceleration", xlab = "Times")points(mcycle$times, mcycle$accel, col = 2)

Printing the output of check.gamViz

Description

This method prints the output ofcheck.gamViz.

Usage

## S3 method for class 'checkGam'print(x, lay = NULL, ...)

Arguments

Value

Returns the output ofgrid.arrange, invisibly.

Printing the output of plot.gamViz

Description

This method prints the output ofplot.gamViz.

Usage

## S3 method for class 'plotGam'print(x, ask = TRUE, pages = NULL, addLay = TRUE, ...)

Arguments

Value

Returns the output ofgridExtra::grid.arrange, invisibly.

Printing plots of smooth effects

Description

This method prints objects of classplotSmooth.

Usage

## S3 method for class 'plotSmooth'print(x, addLay = TRUE, ...)

Arguments

Value

ReturnsNULL, invisibly.

Printing the output of qq.gamViz

Description

This method prints the output ofqq.gamViz.

Usage

## S3 method for class 'qqGam'print(x, ...)

Arguments

Value

ReturnsNULL, invisibly.

Extracting parametric effects from a GAM model

Description

This function can be used to extract a parametric effect from an object ofclassgamViz.

Usage

pterm(o, select)

Arguments

Value

An object of class "pTermSomething" where "Something" is substituted withthe class of the variable of interest. For instance if this "numeric", theptermwill return an object of class "ptermNumeric".

Examples

####### 1. Gaussian GAM library(mgcViz)set.seed(3)dat <- gamSim(1,n=1500,dist="normal",scale=20)dat$fac <- as.factor( sample(c("A1", "A2", "A3"), nrow(dat), replace = TRUE) ) dat$logi <- as.logical( sample(c(TRUE, FALSE), nrow(dat), replace = TRUE) ) bs <- "cr"; k <- 12b <- gam(y ~ x0 + x1 + I(x1^2) + s(x2,bs=bs,k=k) + fac + x3:fac + I(x1*x2) + logi,data=dat)o <- getViz(b)# Plot effect of 'x0'pt <- pterm(o, 1)plot(pt, n = 60) + l_ciPoly() + l_fitLine() + l_ciLine() + l_points()# Plot effect of 'x3'pt <- pterm(o, 1)plot(pt, n = 60) + l_fitLine() + l_ciLine(colour = 2)# Plot effect of 'fac'pt <- pterm(o, 4)plot(pt) + l_ciBar(colour = "blue") + l_fitPoints(colour = "red") +            l_rug(alpha = 0.3)# Plot effect of 'logi'pt <- pterm(o, 6)plot(pt) + l_fitBar(a.aes = list(fill = I("light blue"))) + l_ciBar(colour = "blue")# Plot effect of 'x3:fac': no method available yet available for second order termspt <- pterm(o, 7)plot(pt)## Not run: ####### 1. Continued: Quantile GAMsb <- mqgamV(y ~ x0 + x1 + I(x1^2) + s(x2,bs=bs,k=k) + x3:fac +               I(x1*x2) + logi, data=dat, qu = c(0.3, 0.5, 0.8))plot(pterm(b, 3)) + l_ciBar(colour = 2) + l_fitPoints()plot(pterm(b, 4)) + l_fitBar(colour = "blue", fill = 3) + l_ciBar(colour = 2) # Don't know how to plot this interactionplot(pterm(b, 6)) ####### 2. Gaussian GAMLSS modellibrary(MASS)mcycle$fac <- as.factor( sample(c("z", "k", "a", "f"), nrow(mcycle), replace = TRUE) ) b <- gam(list(accel~times + I(times^2) + s(times,k=10), ~ times + fac + s(times)),          data=mcycle,family=gaulss(), optimizer = "efs")o <- getViz(b)# Plot effect of 'I(times^2)' on mean: notice that partial residuals# are unavailable for GAMLSS models, hence l_point does not do anything here.pt <- pterm(o, 2)plot(pt) + l_ciPoly() + l_fitLine() + l_ciLine() + l_points()# Plot effect of 'times' in second linear predictor.# Notice that partial residuals are unavailable.pt <- pterm(o, 3)plot(pt) + l_ciPoly() + l_fitLine() + l_ciLine(linetype = 3) + l_rug()# Plot effect of 'fac' in second linear predictor.pt <- pterm(o, 4)plot(pt) + l_ciBar(colour = "blue") + l_fitPoints(colour = "red") +            l_rug() ## End(Not run)

Fit a QGAM model and get a gamViz object

Description

These are wrapper that fits a QGAM model usingqgam::qgam andconverts it to agamViz object using thegetViz function.It is essentially a shortcut.

Usage

qgamV(form, data, qu, lsig = NULL, err = NULL, aQgam = list(), aViz = list())

Arguments

Value

An object of class "gamViz" which can, for instance, be plotted usingplot.gamViz.

Examples

library(mgcViz)set.seed(2) ## simulate some data...dat <- gamSim(2,n=1000,dist="normal",scale=0.25)$data# Fit GAM and get gamViz objectb <- qgamV(y~s(x) + s(z) + I(x*z), data = dat, qu = 0.2,           aQgam = list(argGam = list(select = TRUE)), aViz = list("nsim" = 0))# This is equivalent to doing# 1. Fit QGAM# b <- qgam(y~s(x) + s(z) + I(x*z), data=dat, qu = 0.2, argGam = list(select = TRUE))# 2. Convert to gamViz object# b <- getViz(b, nsim = 0)# Either way, we all effects by doingprint(plot(b, allTerms = TRUE), pages = 1)

Generic QQ plots

Description

Generic function for producing QQ-plots.

Usage

qq(...)

Arguments

See Also

qq.gamViz

QQ plots for gam model residuals

Description

Takes a fitted gam object, converted usinggetViz, and produces QQ plots of its residuals(conditional on the fitted model coefficients and scale parameter). If the model distributionalassumptions are met then usually these plots should be close to a straight line (althoughdiscrete data can yield marked random departures from this line).

Usage

## S3 method for class 'gamViz'qq(  o,  rep = 10,  level = 0.8,  method = "auto",  type = "auto",  CI = "none",  worm = FALSE,  showReps = FALSE,  sortFun = NULL,  discrete = NULL,  ngr = 1000,  xlim = NULL,  ylim = NULL,  a.qqpoi = list(),  a.ablin = list(),  a.cipoly = list(),  a.replin = list(),  ...)

Arguments

Details

Heremethod = "simul1" corresponds to the algorithm described in section 2.1 of Augustin et al. (2012), whichinvolves direct simulations of residuals from the models. This requireso$family$rd to be defined.Settingmethod = "simul2" results in a cheaper method, described in section 2.2 of Augustin et al. (2012),which requireso$family$qf to be defined.

Value

An object of classc("qqGam", "plotSmooth", "gg").

References

Augustin, N.H., Sauleau, E.A. and Wood, S.N., 2012. On quantile quantile plots for generalized linear models.Computational Statistics & Data Analysis, 56(8), pp.2404-2409.

Examples

######## Example: simulate binomial datalibrary(mgcViz)set.seed(0)n.samp <- 400dat <- gamSim(1,n = n.samp, dist = "binary", scale = .33)p <- binomial()$linkinv(dat$f) ## binomial pn <- sample(c(1, 3), n.samp, replace = TRUE) ## binomial ndat$y <- rbinom(n, n, p)dat$n <- nlr.fit <- gam(y/n ~ s(x0) + s(x1) + s(x2) + s(x3)              , family = binomial, data = dat,              weights = n, method = "REML")lr.fit <- getViz(lr.fit)# Quick QQ-plot of deviance residualsqq(lr.fit, method = "simul2")# Same, but changing points share and type of reference listqq(lr.fit, method = "simul2",        a.qqpoi = list("shape" = 1), a.ablin = list("linetype" = 2))# Simulation based QQ-plot with reference bands qq(lr.fit, rep = 100, level = .9, CI = "quantile")# Simulation based QQ-plot, Pearson resids, all simulations lines shown qq(lr.fit, rep = 100, CI = "none", showReps = TRUE, type = "pearson",        a.qqpoi = list(shape=19, size = 0.5))### Now fit the wrong model and checkpif <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3)           , family = poisson, data = dat, method = "REML")pif <- getViz(pif)qq(pif, method = "simul2")qq(pif, rep = 100, level = .9, CI = "quantile")qq(pif, rep = 100, type = "pearson", CI = "none", showReps = TRUE,                a.qqpoi = list(shape=19, size = 0.5))######## Example: binary data model violation so gross that you see a problem ######## on the QQ ploty <- c(rep(1, 10), rep(0, 20), rep(1, 40), rep(0, 10), rep(1, 40), rep(0, 40))x <- 1:160b <- glm(y ~ x, family = binomial)class(b) <- c("gamViz", class(b)) # Tricking qq.gamViz to use it on a glm# Note that the next two are not necessarily similar under gross # model violation...qq(b, method = "simul2")qq(b, rep = 50, CI = "none", showReps = TRUE)### alternative modelb <- gam(y ~ s(x, k = 5), family = binomial, method = "ML")b <- getViz(b)qq(b, method = "simul2")qq(b, rep = 50, showReps = TRUE, CI = "none", shape = 19)## Not run: ########  "Big Data" example: set.seed(0)n.samp <- 50000dat <- gamSim(1,n=n.samp,dist="binary",scale=.33)p <- binomial()$linkinv(dat$f) ## binomial pn <- sample(c(1,3),n.samp,replace=TRUE) ## binomial ndat$y <- rbinom(n,n,p)dat$n <- nlr.fit <- bam(y/n ~ s(x0) + s(x1) + s(x2) + s(x3)              , family = binomial, data = dat,              weights = n, method = "fREML", discrete = TRUE)lr.fit <- getViz(lr.fit)# Turning discretization off (on by default for large datasets).set.seed(414) # Setting the seed because qq.gamViz is doing simulationso <- qq(lr.fit, rep = 10, method = "simul1", CI = "normal", showReps = TRUE,             discrete = F, a.replin = list(alpha = 0.1))o # This might take some time!# Using default discretizationset.seed(414)o <- qq(lr.fit, rep = 10, method = "simul1", CI = "normal", showReps = TRUE,             a.replin = list(alpha = 0.1))o # Much faster plotting!# Very coarse discretizationset.seed(414)o <- qq(lr.fit, rep = 10, method = "simul1", CI = "normal", showReps = TRUE,            ngr = 1e2, a.replin = list(alpha = 0.1), a.qqpoi = list(shape = 19))o # We can also zoom in at no extra costs (most work already done by qq.gamViz)zoom(o, xlim = c(-0.25, 0.25), showReps = TRUE, discrete = TRUE, a.replin = list(alpha = 0.2))## End(Not run)

Quantile-Quantile Plots

Description

This is a re-write of the QQ-plotting functions provided bystats, using theggplot2 library.qqnorm is a generic function the default method of which produces a normalQQ plot of the values in y.qqline adds a line to a “theoretical”, by default normal,quantile-quantile plot which passes through theprobs quantiles, by default thefirst and third quartiles.qqplot produces a QQ plot of two datasets.

Usage

qqnorm(  y,  ylim,  main = "Normal Q-Q Plot",  xlab = "Theoretical Quantiles",  ylab = "Sample Quantiles",  datax = FALSE)qqplot(  x,  y,  xlab = deparse(substitute(x)),  ylab = deparse(substitute(y)),  main = "Q-Q Plot")qqline(  y,  datax = FALSE,  distribution = qnorm,  probs = c(0.25, 0.75),  qtype = 7,  ...)

Arguments

Note

Help file is mainly fromstats::qqnorm since this is a rewrite ofstats::qqplot,stats::qqline andstats::qqnorm using the ggplot2 library.

Examples

library(mgcViz)y <- rt(500, df = 5)# Compare new and old version of qqnormstats::qqnorm(y)qqnorm(y)# Compare new and old version of qqplotx <- rt(200, df = 5)y <- rt(300, df = 5)stats::qqplot(x, y)qqplot(x, y)# add a qqline()ggplot2::last_plot() + qqline(y = rt(500, df = 4.8), col = "green") ## "QQ-Chisquare" : --------------------------y <- rchisq(500, df = 3)## Q-Q plot for Chi^2 data against true theoretical distribution:x <- qchisq(ppoints(500), df = 3)stats::qqplot(qchisq(ppoints(500), df = 3), rchisq(500, df = 3),      main = expression("Q-Q plot for" ~~ {chi^2}[nu == 3]))qqplot(qchisq(ppoints(500), df = 3), rchisq(500, df = 3),      main = expression("Q-Q plot for" ~~ {chi^2}[nu == 3])) + theme_bw()

Generalized Additive Model residuals

Description

Extension ofmgcv::residuals.gam. Returns residuals for a fitted GAM model object.Pearson, deviance, working and response residuals are available as in the method frommgcv, but this version also provides residual types "tunif" and "tnormal". The former areobtained using the cdf of the response distribution (if available). The latter are obtainedby further transforming the uniform residuals using the quantile function (i.e. the inversecdf) of a standard normal variable.

Usage

## S3 method for class 'gamViz'residuals(object, type = "deviance", ...)

Arguments

See Also

See alsomgcv::residuals.gam for details.

Generic shine function

Description

Generic function for taking an object and transforming it into ashiny app.

Usage

shine(o, ...)

Arguments

See Also

shine.qqGam

Shiny QQ-plots for GAMs

Description

This function takes the output ofqq.gamViz and transforms itinto an interactive shiny app.

Usage

## S3 method for class 'qqGam'shine(o, ...)

Arguments

Details

In RStudio, this function returns a call toqq.gamViz that reproduces the last plotrendered in the interactive shiny window.

Examples

## Not run: ## simulate binomial data...library(mgcv)library(mgcViz)set.seed(0)n.samp <- 400dat <- gamSim(1,n = n.samp, dist = "binary", scale = .33)p <- binomial()$linkinv(dat$f) ## binomial pn <- sample(c(1, 3), n.samp, replace = TRUE) ## binomial ndat$y <- rbinom(n, n, p)dat$n <- nlr.fit <- gam(y/n ~ s(x0) + s(x1) + s(x2) + s(x3)              , family = binomial, data = dat,              weights = n, method = "REML")lr.fit <- getViz(lr.fit)# Need to load shiny and miniUI at this point# launch shiny gagdetshine(qq(lr.fit)) ## End(Not run)

Simulating responses from a GAM object

Description

This method can be used to simulate vectors of responses from a gamObject.

Usage

## S3 method for class 'gam'simulate(  object,  nsim = 1,  seed = NULL,  method = "auto",  newdata,  u = NULL,  w = NULL,  offset = NULL,  trans = NULL,  ...)

Arguments

Value

A matrix where each column is a vector of simulated responses. The number of rowsis equal to the number of responses in the fitted object.

Examples

library(mgcViz)set.seed(2) ## simulate some data... dat <- gamSim(1,n=400,dist="normal",scale=2)b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)# Simulate three vectors of responsesmatplot(simulate(b, nsim = 3), pch = 19, col = c(1, 3, 4))

Extracting a smooth effect from a GAM model

Description

This function can be used to extract a smooth or random effect from an object ofclassgamViz.

Usage

sm(o, select)

Arguments

Value

An object representing a smooth effect.

See Also

SeegetViz for examples.

Generic zooming function

Description

Generic function for zooming, mainly meant to work with graphical objects.

Usage

zoom(o, ...)

Arguments

See Also

zoom.qqGam

Efficiently zooming on GAM QQ-plots

Description

This function allows to zoom into a QQ-plot produced byqq.gamViz,in a computationally efficient manner.

Usage

## S3 method for class 'qqGam'zoom(  o,  xlim = NULL,  ylim = NULL,  discrete = NULL,  ngr = 1000,  adGrid = TRUE,  CI = FALSE,  worm = FALSE,  showReps = FALSE,  a.qqpoi = list(),  a.ablin = list(),  a.cipoly = list(),  a.replin = list(),  ...)

Arguments

Examples

library(mgcViz);set.seed(0)n.samp <- 500dat <- gamSim(1,n=n.samp,dist="binary",scale=.33)p <- binomial()$linkinv(dat$f) ## binomial pn <- sample(c(1,3),n.samp,replace=TRUE) ## binomial ndat$y <- rbinom(n,n,p)dat$n <- nlr.fit <- bam(y/n ~ s(x0) + s(x1) + s(x2) + s(x3)              , family = binomial, data = dat,              weights = n, method = "REML")lr.fit <- getViz(lr.fit)set.seed(414)o <- qq(lr.fit, rep = 50, method = "simul1", CI = "normal")o # This is the whole qqplot# We can zoom in along x at little extra costs (most computation already done by qq.gamViz)zoom(o, xlim = c(0, 1), showReps = TRUE,      a.replin = list(alpha = 0.1), a.qqpoi =  list(shape = 19))

Constructing a decreasing function from (0,1) to (0,1)

Description

This function returns a non-increasing function from (0, 1) to (0, 1).It takes inputs o, a and m, and it returns the functionf(p)={z=max(0, p-o); return(max((1-z)^a, m))}. The function f(p) can be used,for instance, for transforming p-values before plotting them.

Usage

zto1(o, a, m)

Arguments

Value

A function whose parameters o, a and m have been fixed.

Examples

library(mgcViz)x <- seq(0, 1, by = 0.01)plot(x, zto1(0.05, 2, 0.1)(x), ylim = c(0, 1), type = 'l')lines(x, zto1(0.05, 1, 0.2)(x), col = 2)lines(x, zto1(0.1, 3, 0)(x), col = 3)