Nothing
# R plotting routines for the package mgcv (c) Simon Wood 2000-2017
## With contributions from Henric Nilsson
in.out <- function(bnd,x) {
## tests whether point defined by each row of x is inside
## or outside boundary defined by bnd. bnd my be made up of multiple
## nested loops.
if (!is.matrix(x)) x <- matrix(x,1,2)
if (is.list(bnd)) { ## convert list of lists to matrix form
b1 <- bnd[[1]][[1]]
b2 <- bnd[[1]][[2]]
if (length(bnd)>1) for (i in 2:length(bnd)) {
b1 <- c(b1,NA,bnd[[i]][[1]])
b2 <- c(b2,NA,bnd[[i]][[2]])
}
bnd <- cbind(b1,b2)
}
## replace NA segment separators with a numeric code
lowLim <- min(bnd,na.rm=TRUE) - mean(abs(bnd),na.rm=TRUE)
ind <- is.na(rowSums(bnd))
bnd[ind,] <- lowLim
n <- nrow(bnd)
um <-.C(C_in_out,bx=as.double(bnd[,1]),by=as.double(bnd[,2]),break.code=as.double(lowLim),
x=as.double(x[,1]),y=as.double(x[,2]),inside=as.integer(x[,2]*0),nb=as.integer(n),
n=as.integer(nrow(x)))
as.logical(um$inside)
}
fix.family.qf <- function(fam) {
## add quantile function to family object
if (!inherits(fam, "family"))
stop("fam not a family object")
if (!is.null(fam$qf)) return(fam) ## already exists
family <- fam$family
if (family=="poisson") {
fam$qf <- function(p,mu,wt,scale) {
qpois(p,mu)
}
} else if (family=="binomial") {
fam$qf <- function(p,mu,wt,scale) {
if (all.equal(wt,ceiling(wt))!=TRUE) {
wt <- ceiling(wt)
warning("non-integer binomial denominator: quantiles incorrect")
}
qbinom(p,wt,mu)/(wt + as.numeric(wt==0))
}
} else if (family=="Gamma") {
fam$qf <- function(p,mu,wt,scale) {
qgamma(p,shape=1/scale,scale=mu*scale)
}
} else if (family=="gaussian") {
fam$qf <- function(p,mu,wt,scale) {
qnorm(p,mean=mu,sd=sqrt(scale/wt))
}
}
fam
}
fix.family.rd <- function(fam) {
## add random deviate function to family objet
if (!inherits(fam, "family"))
stop("fam not a family object")
if (!is.null(fam$rd)) return(fam) ## already exists
family <- fam$family
if (family=="poisson") {
fam$rd <- function(mu,wt,scale) {
rpois(length(mu),mu)
}
} else if (family=="binomial") {
fam$rd <- function(mu,wt,scale) {
rbinom(length(mu),wt,mu)/(wt + as.numeric(wt==0))
}
} else if (family=="Gamma") {
fam$rd <- function(mu,wt,scale) {
rgamma(length(mu),shape=1/scale,scale=mu*scale)
}
} else if (family=="gaussian") {
fam$rd <- function(mu,wt,scale) {
rnorm(length(mu),mean=mu,sd=sqrt(scale/wt))
}
} else if (family=="inverse.gaussian") {
fam$rd <- function(mu,wt,scale) {
rig(length(mu),mu,scale)
}
}
fam
}
qq.gam <- function(object, rep=0, level=.9,s.rep=10,
type=c("deviance","pearson","response"),
pch=".", rl.col=2, rep.col="gray80",...) {
## get deviance residual quantiles under good fit
type <- match.arg(type)
ylab <- paste(type,"residuals")
if (inherits(object,c("glm","gam"))) {
if (is.null(object$sig2)) object$sig2 <- summary(object)$dispersion
} else stop("object is not a glm or gam")
## in case of NA & na.action="na.exclude", we need the "short" residuals:
object$na.action <- NULL
D <- residuals(object,type=type)
if (object$method %in% c("PQL","lme.ML","lme.REML","lmer.REML","lmer.ML","glmer.ML")) {
## then it's come out of a gamm fitter and qq.gam can't see the random effects
## that would be necessary to get quantiles. Fall back to normal QQ plot.
qqnorm(D,ylab=ylab,pch=pch,...)
return()
}
lim <- Dq <- NULL
if (rep==0) {
fam <- fix.family.qf(object$family)
if (is.null(fam$qf))
rep <- 50 ## try simulation if quantile function not available
level <- 0
}
n <- length(D)
if (rep > 0) { ## simulate quantiles
fam <- fix.family.rd(object$family)
if (!is.null(fam$rd)) {
##d <- rep(0,0)
## simulate deviates...
dm <- matrix(0,n,rep)
for (i in 1:rep) {
yr <- fam$rd(object$fitted.values, object$prior.weights, object$sig2)
#di <- fam$dev.resids(yr,object$fitted.values,object$prior.weights)^.5*
# sign(yr-object$fitted.values)
object$y <- yr
dm[,i] <- sort(residuals(object,type=type))
#d <- c(d,sort(di))
}
# n <- length(D)
Dq <- quantile(as.numeric(dm),(1:n - .5)/n)
## now get simulation limits on QQ plot
#dm <- matrix(d,length(Dq),rep)
alpha <- (1-level)/2
if (alpha>.5||alpha<0) alpha <- .05
if (level>0&&level<1) lim <- apply(dm,1,FUN=quantile,p=c(alpha,1-alpha)) else
if (level >= 1) lim <- level
}
} else {
## ix <- sort.int(D,index.return=TRUE)$ix ## messes up under multiple ties!
#ix <- rank(D)
#U <- (ix-.5)/length(D) ## code used pre-randomization - not needed
U <- (1:n-.5)/n
if (!is.null(fam$qf)) {
dm <- matrix(0,n,s.rep)
for (i in 1:s.rep) {
U <- sample(U,n) ## randomize uniform quantiles w.r.t. obs
q0 <- fam$qf(U,object$fitted.values,object$prior.weights,object$sig2)
object$y <- q0
dm[,i] <- sort(residuals(object,type=type)) ## original proposal
}
Dq <- sort(rowMeans(dm))
}
}
if (!is.null(Dq))
{ qqplot(Dq,D,ylab=ylab,xlab="theoretical quantiles",ylim=range(c(lim,D)),
pch=pch,...)
abline(0,1,col=rl.col)
if (!is.null(lim)) {
if (level>=1) for (i in 1:rep) lines(Dq,dm[,i],col=rep.col) else {
n <- length(Dq)
polygon(c(Dq,Dq[n:1],Dq[1]),c(lim[1,],lim[2,n:1],lim[1,1]),col=rep.col,border=NA)
}
abline(0,1,col=rl.col)
}
points(Dq,sort(D),pch=pch,...)
return(invisible(Dq))
} else qqnorm(D,ylab=ylab,pch=pch,...)
} ## qq.gam
k.check <- function(b,subsample=5000,n.rep=400) {
## function to check k in a gam fit...
## does a randomization test looking for evidence of residual
## pattern attributable to covariates of each smooth.
m <- length(b$smooth)
if (m==0) return(NULL)
rsd <- residuals(b)
ve <- rep(0,n.rep)
p.val<-v.obs <- kc <- edf<- rep(0,m)
snames <- rep("",m)
n <- nrow(b$model)
if (n>subsample) { ## subsample to avoid excessive cost
ind <- sample(1:n,subsample)
modf <- b$model[ind,]
rsd <- rsd[ind]
} else modf <- b$model
nr <- length(rsd)
for (k in 1:m) { ## work through smooths
ok <- TRUE
b$smooth[[k]]$by <- "NA" ## can't deal with by variables
dat <- ExtractData(b$smooth[[k]],modf,NULL)$data
if (!is.null(attr(dat,"index"))||!is.null(attr(dat[[1]],"matrix"))||is.matrix(dat[[1]])) ok <- FALSE
if (ok) dat <- as.data.frame(dat)
snames[k] <- b$smooth[[k]]$label
ind <- b$smooth[[k]]$first.para:b$smooth[[k]]$last.para
kc[k] <- length(ind)
edf[k] <- sum(b$edf[ind])
nc <- b$smooth[[k]]$dim
if (ok && ncol(dat)>nc) dat <- dat[,1:nc,drop=FALSE] ## drop any by variables
for (j in 1:nc) if (is.factor(dat[[j]])) ok <- FALSE
if (!ok) {
p.val[k] <- v.obs[k] <- NA ## can't do this test with summation convention/factors
} else { ## normal term
if (nc==1) { ## 1-D term
e <- diff(rsd[order(dat[,1])])
v.obs[k] <- mean(e^2)/2
for (i in 1:n.rep) {
e <- diff(rsd[sample(1:nr,nr)]) ## shuffle
ve[i] <- mean(e^2)/2
}
p.val[k] <- mean(ve<v.obs[k])
v.obs[k] <- v.obs[k]/mean(rsd^2)
} else { ## multi-D
if (!is.null(b$smooth[[k]]$margin)) { ## tensor product (have to consider scaling)
## get the scale factors...
beta <- coef(b)[ind]
f0 <- PredictMat(b$smooth[[k]],dat)%*%beta
gr.f <- rep(0,ncol(dat))
for (i in 1:nc) {
datp <- dat;dx <- diff(range(dat[,i]))/1000
datp[,i] <- datp[,i] + dx
fp <- PredictMat(b$smooth[[k]],datp)%*%beta
gr.f[i] <- mean(abs(fp-f0))/dx
}
for (i in 1:nc) { ## rescale distances
dat[,i] <- dat[,i] - min(dat[,i])
dat[,i] <- gr.f[i]*dat[,i]/max(dat[,i])
}
}
nn <- 3
ni <- nearest(nn,as.matrix(dat))$ni
e <- rsd - rsd[ni[,1]]
for (j in 2:nn) e <- c(e,rsd-rsd[ni[,j]])
v.obs[k] <- mean(e^2)/2
for (i in 1:n.rep) {
rsdr <- rsd[sample(1:nr,nr)] ## shuffle
e <- rsdr - rsdr[ni[,1]]
for (j in 2:nn) e <- c(e,rsdr-rsdr[ni[,j]])
ve[i] <- mean(e^2)/2
}
p.val[k] <- mean(ve<v.obs[k])
v.obs[k] <- v.obs[k]/mean(rsd^2)
}
}
}
k.table <- cbind(kc,edf,v.obs, p.val)
dimnames(k.table) <- list(snames, c("k\'","edf","k-index", "p-value"))
k.table
} ## end of k.check
gam.check <- function(b, old.style=FALSE,
type=c("deviance","pearson","response"),
k.sample=5000,k.rep=200,
## arguments passed to qq.gam() {w/o warnings !}:
rep=0, level=.9, rl.col=2, rep.col="gray80", ...)
## takes a fitted gam object and produces some standard diagnostic plots
{
type <- match.arg(type)
resid <- residuals(b, type=type)
linpred <- if (is.matrix(b$linear.predictors)&&!is.matrix(resid))
napredict(b$na.action, b$linear.predictors[,1]) else
napredict(b$na.action, b$linear.predictors)
## if (b$method%in%c("GCV","GACV","UBRE","REML","ML","P-ML","P-REML","mle.REML","mle.ML","PQL")) {
if (is.null(.Platform$GUI) || .Platform$GUI != "RStudio") old.par <- par(mfrow=c(2,2))
if (old.style)
qqnorm(resid,...)
else
qq.gam(b, rep=rep, level=level, type=type, rl.col=rl.col, rep.col=rep.col, ...)
plot(linpred, resid,main="Resids vs. linear pred.",
xlab="linear predictor",ylab="residuals",...)
hist(resid,xlab="Residuals",main="Histogram of residuals",...)
fv <- if (inherits(b$family,"extended.family")) predict(b,type="response") else fitted(b)
if (is.matrix(fv)&&!is.matrix(b$y)) fv <- fv[,1]
plot(fv, napredict(b$na.action, b$y),
xlab="Fitted Values",ylab="Response",main="Response vs. Fitted Values",...)
gamm <- !(b$method%in%c("GCV","GACV","UBRE","REML","ML","P-ML","P-REML","fREML","NCV")) ## gamm `gam' object
#if (is.null(.Platform$GUI) || .Platform$GUI != "RStudio") par(old.par)
# return(invisible())
#}
## now summarize convergence information
if (gamm) {
cat("\n\'gamm\' based fit - care required with interpretation.")
cat("\nChecks based on working residuals may be misleading.")
} else {
cat("\nMethod:",b$method," Optimizer:",b$optimizer)
if (!is.null(b$outer.info)) { ## summarize convergence information
if (b$optimizer[2]%in%c("newton","bfgs"))
{ boi <- b$outer.info
cat("\n",boi$conv," after ",boi$iter," iteration",sep="")
if (boi$iter==1) cat(".") else cat("s.")
cat("\nGradient range [",min(boi$grad),",",max(boi$grad),"]",sep="")
cat("\n(score ",b$gcv.ubre," & scale ",b$sig2,").",sep="")
ev <- eigen(boi$hess)$values
if (min(ev)>0) cat("\nHessian positive definite, ") else cat("\n")
cat("eigenvalue range [",min(ev),",",max(ev),"].\n",sep="")
} else { ## just default print of information ..
cat("\n");print(b$outer.info)
}
} else { ## no sp, perf iter or AM case
if (length(b$sp)==0) ## no sp's estimated
cat("\nModel required no smoothing parameter selection")
else {
cat("\nSmoothing parameter selection converged after",b$mgcv.conv$iter,"iteration")
if (b$mgcv.conv$iter>1) cat("s")
if (!b$mgcv.conv$fully.converged)
cat(" by steepest\ndescent step failure.\n") else cat(".\n")
cat("The RMS",b$method,"score gradient at convergence was",b$mgcv.conv$rms.grad,".\n")
if (b$mgcv.conv$hess.pos.def)
cat("The Hessian was positive definite.\n") else cat("The Hessian was not positive definite.\n")
#cat("The estimated model rank was ",b$mgcv.conv$rank,
# " (maximum possible: ",b$mgcv.conv$full.rank,")\n",sep="")
}
}
if (!is.null(b$rank)) {
cat("Model rank = ",b$rank,"/",length(b$coefficients),"\n")
}
} ## if gamm
cat("\n")
## now check k
kchck <- k.check(b,subsample=k.sample,n.rep=k.rep)
if (!is.null(kchck)) {
cat("Basis dimension (k) checking results. Low p-value (k-index<1) may\n")
cat("indicate that k is too low, especially if edf is close to k\'.\n\n")
printCoefmat(kchck,digits=3);
}
if (is.null(.Platform$GUI) ||.Platform$GUI != "RStudio") par(old.par)
## } else plot(linpred,resid,xlab="linear predictor",ylab="residuals",...)
} ## end of gam.check
#############################################
## Start of plot method functions for smooths
#############################################
plot.random.effect <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
## plot method for a "random.effect" smooth class
if (is.null(P)) { ## get plotting information...
if (!x$plot.me) return(NULL) else { ## shouldn't or can't plot
raw <- data[x$term][[1]]
p <- x$last.para - x$first.para + 1
X <- diag(p) # prediction matrix for this term
if (is.null(xlab)) xlabel<- "Gaussian quantiles" else xlabel <- xlab
if (is.null(ylab)) ylabel <- "effects" else ylabel <- ylab
if (!is.null(main)) label <- main
return(list(X=X,scale=FALSE,se=FALSE,raw=raw,xlab=xlabel,ylab=ylabel,
main=label))
} ## end of basic plot data production
} else { ## produce plot
b <- as.numeric(trans(P$fit+shift))
qqnorm(b,main=P$main,xlab=P$xlab,ylab=P$ylab,...)
qqline(b)
} ## end of plot production
} ## end of plot.random.effect
repole <- function(lo,la,lop,lap) {
## painfully plodding function to get new lo, la relative to pole at
## lap,lop...
## x,y,z location of pole...
yp <- sin(lap)
xp <- cos(lap)*sin(lop)
zp <- cos(lap)*cos(lop)
## x,y,z location of meridian point for pole - i.e. point lat pi/2
## from pole on pole's lon.
ym <- sin(lap-pi/2)
xm <- cos(lap-pi/2)*sin(lop)
zm <- cos(lap-pi/2)*cos(lop)
## x,y,z locations of points in la, lo
y <- sin(la)
x <- cos(la)*sin(lo)
z <- cos(la)*cos(lo)
## get angle between points and new equatorial plane (i.e. plane orthogonal to pole)
d <- sqrt((x-xp)^2+(y-yp)^2+(z-zp)^2) ## distance from points to to pole
phi <- pi/2-2*asin(d/2)
## location of images of la,lo on (new) equatorial plane
## sin(phi) gives distance to plane, -(xp, yp, zp) is
## direction...
x <- x - xp*sin(phi)
y <- y - yp*sin(phi)
z <- z - zp*sin(phi)
## get distances to meridian point
d <- sqrt((x-xm)^2+(y-ym)^2+(z-zm)^2)
## angles to meridian plane (i.e. plane containing origin, meridian point and pole)...
theta <- (1+cos(phi)^2-d^2)/(2*cos(phi))
theta[theta < -1] <- -1; theta[theta > 1] <- 1
theta <- acos(theta)
## now decide which side of meridian plane...
## get points at extremes of hemispheres on either side
## of meridian plane....
y1 <- 0
x1 <- sin(lop+pi/2)
z1 <- cos(lop+pi/2)
y0 <- 0
x0 <- sin(lop-pi/2)
z0 <- cos(lop-pi/2)
d1 <- sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2)
d0 <- sqrt((x-x0)^2+(y-y0)^2+(z-z0)^2)
ii <- d0 < d1 ## index -ve lon hemisphere
theta[ii] <- -theta[ii]
list(lo=theta,la=phi)
} ## end of repole
lolaxy <- function(lo,la,theta,phi) {
## takes locations lo,la, relative to a pole at lo=theta, la=phi.
## theta, phi are expressed relative to plotting co-ordinate system
## with pole at top. Convert to x,y in plotting co-ordinates.
## all in radians!
er <- repole(-lo,la,-pi,phi)
er$lo <- er$lo - theta
y <- sin(er$la)
x <- cos(er$la)*sin(er$lo)
z <- cos(er$la)*cos(er$lo)
ind <- z<0
list(x=x[ind],y=y[ind])
} ## end of lolaxy
plot.sos.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(100),
contour.col=4,...) {
## plot method function for sos.smooth terms
if (scheme>1) return(plot.mgcv.smooth(x,P=P,data=data,label=label,se1.mult=se1.mult,se2.mult=se2.mult,
partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,
theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme-2,
hcolors=hcolors,contour.col=contour.col,...))
## convert location of pole in plotting grid to radians
phi <- phi*pi/180
theta <- theta*pi/180
## re-map to sensible values...
theta <- theta%%(2*pi)
if (theta>pi) theta <- theta - 2*pi
phi <- phi%%(2*pi)
if (phi > pi) phi <- phi - 2*pi
if (phi > pi/2) phi <- pi - phi
if (phi < -pi/2 ) phi <- -(phi+pi)
if (is.null(P)) { ## get plotting information...
if (!x$plot.me) return(NULL) ## shouldn't or can't plot
## get basic plot data
raw <- data[x$term]
if (rug) { ## need to project data onto plotting grid...
raw <- lolaxy(lo=raw[[2]]*pi/180,la=raw[[1]]*pi/180,theta,phi)
}
m <- round(n2*1.5)
ym <- xm <- seq(-1,1,length=m)
gr <- expand.grid(x=xm,y=ym)
r <- z <- gr$x^2+gr$y^2
z[z>1] <- NA
z <- sqrt(1-z)
## generate la, lo in plotting grid co-ordinates...
ind <- !is.na(z)
r <- r[ind]
la <- asin(gr$y[ind])
lo <- cos(la)
lo <- asin(gr$x[ind]/lo)
um <- repole(lo,la,theta,phi)
dat <- data.frame(la=um$la*180/pi,lo=um$lo*180/pi)
names(dat) <- x$term
if (x$by!="NA") dat[[x$by]] <- la*0+1
X <- PredictMat(x,dat) # prediction matrix for this term
## fix lo for smooth contouring
lo <- dat[[2]]
ii <- lo <= -177
lo[ii] <- lo[ii] <- 360 + lo[ii]
ii <- lo < -165 & lo > -177
ii <- ii | (abs(dat[[1]])>80)
lo[ii] <- NA
return(list(X=X,scale=FALSE,se=FALSE,raw=raw,xlab="",ylab="",main="",
ind=ind,xm=xm,ym=ym,lo=lo,la=dat[[1]]))
} else { ## do plot
op <- par(pty="s",mar=c(0,0,0,0))
m <- length(P$xm); zz <- rep(NA,m*m)
if (scheme == 0) {
col <- 1# "lightgrey
zz[P$ind] <- trans(P$fit+shift)
image(P$xm,P$ym,matrix(zz,m,m),col=hcolors,axes=FALSE,xlab="",ylab="",...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
}
zz[P$ind] <- P$la
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,levels=c(-8:8*10),col=col,...)
zz[P$ind] <- P$lo
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,levels=c(-8:9*20),col=col,...)
zz[P$ind] <- P$fit
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,col=contour.col,...)
} else if (scheme == 1) {
col <- 1
zz[P$ind] <- trans(P$fit+shift)
contour(P$xm,P$ym,matrix(zz,m,m),col=1,axes=FALSE,xlab="",ylab="",...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
}
zz[P$ind] <- P$la
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,levels=c(-8:8*10),col=col,lty=2,...)
zz[P$ind] <- P$lo
contour(P$xm,P$ym,matrix(zz,m,m),add=TRUE,levels=c(-8:9*20),col=col,lty=2,...)
theta <- seq(-pi/2,pi/2,length=200)
x <- sin(theta);y <- cos(theta)
x <- c(x,x[200:1]);y <- c(y,-y[200:1])
lines(x,y)
}
par(op)
}
} ## end plot.sos.smooth
poly2 <- function(x,col) {
## let x be a 2 col matrix defining some polygons.
## Different closed loop sections are separated by
## NA rows. This routine assumes that loops nested within
## other loops are holes (further nesting gives and island
## in hole, etc). Holes are left unfilled.
## The first polygon should not be a hole.
ind <- (1:nrow(x))[is.na(rowSums(x))] ## where are the splits?
if (length(ind)==0|| ind[1]==nrow(x)) polygon(x,col=col,border="black") else {
base <- x[1,]
xf <- x
xf[ind,1] <- base[1]
xf[ind,2] <- base[2]
if (!is.na(col)) polygon(xf,col=col,border=NA,fillOddEven=TRUE)
polygon(x,border="black")
}
} ## poly2
polys.plot <- function(pc,z=NULL,scheme="heat",lab="",...) {
## pc is a list of polygons defining area boundaries
## pc[[i]] is the 2 col matrix of vertex co-ords for polygons defining
## boundary of area i
## z gives the value associated with the area
## first find the axes ranges...
for (i in 1:length(pc)) {
yr <- range(pc[[i]][,2],na.rm=TRUE)
xr <- range(pc[[i]][,1],na.rm=TRUE)
if (i==1) {
ylim <- yr
xlim <- xr
} else {
if (yr[1]<ylim[1]) ylim[1] <- yr[1]
if (yr[2]>ylim[2]) ylim[2] <- yr[2]
if (xr[1]<xlim[1]) xlim[1] <- xr[1]
if (xr[2]>xlim[2]) xlim[2] <- xr[2]
}
} ## end of axes range loop
mar <- par("mar");
oldpar <- par(mar=c(2,mar[2],2,1))
if (is.null(z)) { ## no z value, no shading, no scale, just outlines...
plot(0,0,ylim=ylim,xlim=xlim,xaxt="n",yaxt="n",type="n",bty="n",ylab=lab,xlab="",...)
for (i in 1:length(pc)) {
poly2(pc[[i]],col=NA)
}
} else {
nz <- names(z)
npc <- names(pc)
if (!is.null(nz)&&!is.null(npc)) { ## may have to re-order z into pc order.
if (all.equal(sort(nz),sort(npc))!=TRUE) stop("names of z and pc must match")
z <- z[npc]
}
xmin <- xlim[1]
xlim[1] <- xlim[1] - .1 * (xlim[2]-xlim[1]) ## allow space for scale
n.col <- 100
if (scheme=="heat") scheme <- heat.colors(n.col+1) else
scheme <- gray(0:n.col/n.col)
zlim <- range(pretty(z))
## Now want a grey or color scale up the lhs of plot
## first scale the y range into the z range for plotting
for (i in 1:length(pc)) pc[[i]][,2] <- zlim[1] +
(zlim[2]-zlim[1])*(pc[[i]][,2]-ylim[1])/(ylim[2]-ylim[1])
ylim <- zlim
plot(0,0,ylim=ylim,xlim=xlim,type="n",xaxt="n",bty="n",xlab="",ylab=lab,...)
for (i in 1:length(pc)) {
coli <- round((z[i] - zlim[1])/(zlim[2]-zlim[1])*n.col)+1
poly2(pc[[i]],col=scheme[coli])
}
## now plot the scale bar...
xmin <- min(c(axTicks(1),xlim[1]))
dx <- (xlim[2]-xlim[1])*.05
x0 <- xmin-2*dx
x1 <- xmin+dx
dy <- (ylim[2]-ylim[1])/n.col
poly <- matrix(c(x0,x0,x1,x1,ylim[1],ylim[1]+dy,ylim[1]+dy,ylim[1]),4,2)
for (i in 1:n.col) {
polygon(poly,col=scheme[i],border=NA)
poly[,2] <- poly[,2] + dy
}
poly <- matrix(c(x0,x0,x1,x1,ylim[1],ylim[2],ylim[2],ylim[1]),4,2)
polygon(poly,border="black")
}
par(oldpar)
} ## polys.plot
plot.mrf.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
## plot method function for mrf.smooth terms, depends heavily on polys.plot, above
if (is.null(P)) { ## get plotting information...
if (!x$plot.me||is.null(x$xt$polys)) return(NULL) ## shouldn't or can't plot
## get basic plot data
raw <- data[x$term][[1]]
dat <- data.frame(x=factor(names(x$xt$polys),levels=levels(x$knots)))
names(dat) <- x$term
x$by <- "NA"
X <- PredictMat(x,dat) # prediction matrix for this term
if (is.null(xlab)) xlabel<- "" else xlabel <- xlab
if (is.null(ylab)) ylabel <- "" else ylabel <- ylab
return(list(X=X,scale=FALSE,se=FALSE,raw=raw,xlab=xlabel,ylab=ylabel,
main=label))
} else { ## do plot
if (scheme==0) scheme <- "heat" else scheme <- "grey"
polys.plot(x$xt$polys,trans(P$fit+shift),scheme=scheme,lab=P$main,...)
}
} ## end plot.mrf.smooth
plot.sz.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
## plot method for factor smooth interactions designed for models such as s(x) + s(fac,x) where
## the factor level dependent smooths are strictly deviations from the main effect smooth.
nf2i <- function(nf,i) {
k <- length(nf)
kk <- rep(0,k)
i <- i-1
for (j in k:1) {
kk[j] <- i %% nf[j] + 1
i <- i %/% nf[j]
}
kk
} ## nf2i
if (is.null(P)) { ## get plotting info
if (x$base$dim!=1) return(NULL) ## no method for base smooth dim > 1
raw <- data[x$base$term][[1]]
xx <- seq(min(raw),max(raw),length=n) # generate x sequence for prediction
nf <- unlist(lapply(x$flev,length)) ## levels per grouping factor
dat <- data.frame(rep(xx,prod(nf)))
k <- length(x$flev) ## number of factors
for (i in 1:k) {
re <- if (i<k) prod(nf[(i+1):k]) else 1
rs <- if (i>1) prod(nf[1:(i-1)]) else 1
dat[,i+1] <- factor(rep(rep(x$flev[[i]],each=re*n),rs),levels=x$flev[[i]])
}
names(dat) <- c(x$base$term,x$fterm)
if (x$by!="NA") { # deal with any by variables
dat[[x$by]] <- rep(1,n)
}
X <- PredictMat(x,dat)
if (is.null(xlab)) xlabel <- x$base$term else xlabel <- xlab
if (is.null(ylab)) ylabel <- label else ylabel <- ylab
return(list(X=X,scale=TRUE,se=TRUE,se.mult=se1.mult,raw=raw,xlab=xlabel,ylab=ylabel,
main="",x=xx,n=n,nf=nf))
} else { ## produce the plot
nft <- prod(P$nf) ## total number of curves
if (scheme!=1) {
kol <- hcl.colors(nft,palette = "viridis", alpha = .33) ## CI
lkol <- hcl.colors(nft,palette = "viridis", alpha = .66) ## mode
tkol <- hcl.colors(nft,palette = "viridis", alpha = 1) ## label
}
xlim <- range(P$x);dx <- xlim[2]-xlim[1]
xt <- xlim[1] + (1:nft-.5)*dx/nft ## text locations
ind <- 1:P$n; mind <- P$n:1
if(is.null(ylim)) ylim <- trans(range(c(P$fit+P$se,P$fit-P$se))+shift)
plot(P$x[ind],trans(P$fit[ind]+shift),ylim=ylim,xlab=P$xlab,ylab=P$ylab,type="n",...)
nfac <- length(P$nf) ## number of factors
kk <- rep(0,nfac) ## factor level index vector
if (scheme==1) {
for (i in 1:nft) {
ul <- trans(P$fit[ind] + P$se[ind]+shift)
ll <- trans(P$fit[ind] - P$se[ind]+shift)
lines(P$x,ul,col="grey",lty=i);lines(P$x,ll,col="grey",lty=i)
ii <- P$x < xt[i] - dx/30
yt <- approx(P$x,P$fit[ind],xt[i])$y
lines(P$x[ii],(P$fit[ind])[ii],lty=i,lwd=2)
text(xt[i],yt,paste(nf2i(P$nf,i),collapse="."))
ii <- P$x > xt[i] + dx/30
lines(P$x[ii],(P$fit[ind])[ii],lty=i,lwd=2)
ind <- ind + P$n; mind <- mind + P$n
}
} else {
for (i in 1:nft) {
ul <- trans(P$fit[ind] + P$se[ind]+shift)
ll <- trans(P$fit[mind] - P$se[mind]+shift)
polygon(c(P$x,P$x[P$n:1]),c(ul,ll),col=kol[i],border=kol[i])
yt <- approx(P$x,P$fit[ind],xt[i])$y
ii <- P$x < xt[i] - dx/30
lines(P$x[ii],(P$fit[ind])[ii],col=lkol[i])
text(xt[i],yt,paste(nf2i(P$nf,i),collapse="."),col=tkol[i])
ii <- P$x > xt[i] + dx/30
lines(P$x[ii],(P$fit[ind])[ii],col=lkol[i])
ind <- ind + P$n; mind <- mind + P$n
}
}
}
} ## end plot.sz.interaction
plot.fs.interaction <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,...) {
## plot method for simple smooth factor interactions...
if (is.null(P)) { ## get plotting info
if (x$dim!=1) return(NULL) ## no method for base smooth dim > 1
raw <- data[x$base$term][[1]]
xx <- seq(min(raw),max(raw),length=n) # generate x sequence for prediction
nf <- length(x$flev)
fac <- rep(x$flev,rep(n,nf))
dat <- data.frame(fac,xx,stringsAsFactors=TRUE)
names(dat) <- c(x$fterm,x$base$term)
if (x$by!="NA") { # deal with any by variables
dat[[x$by]] <- rep(1,n)
}
X <- PredictMat(x,dat)
if (is.null(xlab)) xlabel <- x$base$term else xlabel <- xlab
if (is.null(ylab)) ylabel <- label else ylabel <- ylab
return(list(X=X,scale=TRUE,se=FALSE,raw=raw,xlab=xlabel,ylab=ylabel,
main="",x=xx,n=n,nf=nf))
} else { ## produce the plot
ind <- 1:P$n
if(is.null(ylim)) ylim <- trans(range(P$fit)+shift)
plot(P$x[ind],trans(P$fit[ind]+shift),ylim=ylim,xlab=P$xlab,ylab=P$ylab,type="l",...)
if (P$nf>1) for (i in 2:P$nf) {
ind <- ind + P$n
if (scheme==0) lines(P$x,trans(P$fit[ind]+shift),lty=i,col=i) else
lines(P$x,trans(P$fit[ind]+shift),lty=i)
}
}
} ## end plot.fs.interaction
plot.mgcv.smooth <- function(x,P=NULL,data=NULL,label="",se1.mult=2,se2.mult=1,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(50),
contour.col=4,...) {
## default plot method for smooth objects `x' inheriting from "mgcv.smooth"
## `x' is a smooth object, usually part of a `gam' fit. It has an attribute
## 'coefficients' containing the coefs for the smooth, but usually these
## are not needed.
## Usually this function is called twice. First to set up P, then to compute the
## actual plot information including standard error bands, and then to actually
## plot...
## `P' is a list of plot data.
## If `P' is NULL (first call) then the routine should compute some of this plot data
## and return without plotting...
## * X the matrix mapping the smooth's coefficients to the values at
## which the smooth must be computed for plotting.
## * The values against which to plot.
## * `exclude' indicates rows of X%*%p to set to NA for plotting -- NULL for none.
## * se TRUE if plotting of the term can use standard error information.
## * se.mult - the multiplier of the standard error used to plot confidence bands
## * scale TRUE if the term should be considered by plot.gam if a common
## y scale is required.
## * any raw data information.
## * axis labels and plot titles
## As an alternative, P may contain a 'fit' field directly, in which case the
## very little processing is done outside the routine, except for partial residual
## computations.
## Alternatively return P as NULL if x should not be plotted.
## If P is not NULL (second call) it will contain the following...
## * fit - the values for plotting
## * se - standard errors of fit multiplied by se.mult
## * the values against which to plot
## * any raw data information
## * any partial.residuals
## `data' is a data frame containing the raw data for the smooth, usually the
## model.frame of the fitted gam. Can be NULL if P is not NULL.
## `label' is the term label, usually something like e.g. `s(x,12.34)'.
## Note that if ylim is supplied it should not be transformed using trans and shift.
#############################
sp.contour <- function(x,y,z,zse,xlab="",ylab="",zlab="",titleOnly=FALSE,
se.plot=TRUE,se.mult=1,trans=I,shift=0,...)
## function for contouring 2-d smooths with 1 s.e. limits
{ gap<-median(zse,na.rm=TRUE)
zr<-max(trans(z+zse+shift),na.rm=TRUE)-min(trans(z-zse+shift),na.rm=TRUE) # plotting range
n<-10
while (n>1 && zr/n<2.5*gap) n<-n-1
zrange<-c(min(trans(z-zse+shift),na.rm=TRUE),max(trans(z+zse+shift),na.rm=TRUE))
zlev<-pretty(zrange,n) ## ignore codetools on this one
yrange<-range(y);yr<-yrange[2]-yrange[1]
xrange<-range(x);xr<-xrange[2]-xrange[1]
ypos<-yrange[2]+yr/10
args <- as.list(substitute(list(...)))[-1]
args$x <- substitute(x);args$y <- substitute(y)
args$type="n";args$xlab<-args$ylab<-"";args$axes<-FALSE
do.call("plot",args)
cs<-(yr/10)/strheight(zlab);if (cs>1) cs<-1 # text scaling based on height
tl<-strwidth(zlab);
if (tl*cs>3*xr/10) cs<-(3*xr/10)/tl
args <- as.list(substitute(list(...)))[-1]
n.args <- names(args)
zz <- trans(z+shift) ## ignore codetools for this
args$x<-substitute(x);args$y<-substitute(y);args$z<-substitute(zz)
if (!"levels"%in%n.args) args$levels<-substitute(zlev)
if (!"lwd"%in%n.args) args$lwd<-2
if (!"labcex"%in%n.args) args$labcex<-cs*.65
if (!"axes"%in%n.args) args$axes <- FALSE
if (!"add"%in%n.args) args$add <- TRUE
do.call("contour",args)
if (is.null(args$cex.main)) cm <- 1 else cm <- args$cex.main
if (titleOnly) title(zlab,cex.main=cm) else
{ xpos<-xrange[1]+3*xr/10
xl<-c(xpos,xpos+xr/10); yl<-c(ypos,ypos)
lines(xl,yl,xpd=TRUE,lwd=args$lwd)
text(xpos+xr/10,ypos,zlab,xpd=TRUE,pos=4,cex=cs*cm,off=0.5*cs*cm)
}
if (is.null(args$cex.axis)) cma <- 1 else cma <- args$cex.axis
axis(1,cex.axis=cs*cma);axis(2,cex.axis=cs*cma);box();
if (is.null(args$cex.lab)) cma <- 1 else cma <- args$cex.lab
mtext(xlab,1,2.5,cex=cs*cma);mtext(ylab,2,2.5,cex=cs*cma)
if (!"lwd"%in%n.args) args$lwd<-1
if (!"lty"%in%n.args) args$lty<-2
if (!"col"%in%n.args) args$col<-2
if (!"labcex"%in%n.args) args$labcex<-cs*.5
zz <- trans(z+zse+shift)
args$z<-substitute(zz)
do.call("contour",args)
if (!titleOnly) {
xpos<-xrange[1]
xl<-c(xpos,xpos+xr/10)#;yl<-c(ypos,ypos)
lines(xl,yl,xpd=TRUE,lty=args$lty,col=args$col)
text(xpos+xr/10,ypos,paste("-",round(se.mult),"se",sep=""),xpd=TRUE,pos=4,cex=cs*cm,off=0.5*cs*cm)
}
if (!"lty"%in%n.args) args$lty<-3
if (!"col"%in%n.args) args$col<-3
zz <- trans(z - zse+shift)
args$z<-substitute(zz)
do.call("contour",args)
if (!titleOnly) {
xpos<-xrange[2]-xr/5
xl<-c(xpos,xpos+xr/10);
lines(xl,yl,xpd=TRUE,lty=args$lty,col=args$col)
text(xpos+xr/10,ypos,paste("+",round(se.mult),"se",sep=""),xpd=TRUE,pos=4,cex=cs*cm,off=0.5*cs*cm)
}
} ## end of sp.contour
if (is.null(P)) { ## get plotting information...
if (!x$plot.me||x$dim>4) return(NULL) ## shouldn't or can't plot
if (x$dim==1) { ## get basic plotting data for 1D terms
raw <- data[x$term][[1]]
if (is.null(xlim)) xx <- seq(min(raw),max(raw),length=n) else # generate x sequence for prediction
xx <- seq(xlim[1],xlim[2],length=n)
if (x$by!="NA") # deal with any by variables
{ by<-rep(1,n);dat<-data.frame(x=xx,by=by)
names(dat)<-c(x$term,x$by)
} else {
dat<-data.frame(x=xx);names(dat) <- x$term
} ## prediction data.frame finished
X <- PredictMat(x,dat) # prediction matrix for this term
if (is.null(xlab)) xlabel<- x$term else xlabel <- xlab
if (is.null(ylab)) ylabel <- label else ylabel <- ylab
if (is.null(xlim)) xlim <- range(xx)
return(list(X=X,x=xx,scale=TRUE,se=TRUE,raw=raw,xlab=xlabel,ylab=ylabel,
main=main,se.mult=se1.mult,xlim=xlim))
} else if (x$dim==2) { ## basic plot data for 2D terms
xterm <- x$term[1]
if (is.null(xlab)) xlabel <- xterm else xlabel <- xlab
yterm <- x$term[2]
if (is.null(ylab)) ylabel <- yterm else ylabel <- ylab
raw <- data.frame(x=as.numeric(data[xterm][[1]]),
y=as.numeric(data[yterm][[1]]))
n2 <- max(10,n2)
if (is.null(xlim)) xm <- seq(min(raw$x),max(raw$x),length=n2) else
xm <- seq(xlim[1],xlim[2],length=n2)
if (is.null(ylim)) ym <- seq(min(raw$y),max(raw$y),length=n2) else
ym <- seq(ylim[1],ylim[2],length=n2)
xx <- rep(xm,n2)
yy <- rep(ym,rep(n2,n2))
if (too.far>0)
exclude <- exclude.too.far(xx,yy,raw$x,raw$y,dist=too.far) else
exclude <- rep(FALSE,n2*n2)
if (x$by!="NA") { # deal with any by variables
by <- rep(1,n2^2);dat <- data.frame(x=xx,y=yy,by=by)
names(dat) <- c(xterm,yterm,x$by)
} else {
dat<-data.frame(x=xx,y=yy);names(dat)<-c(xterm,yterm)
} ## prediction data.frame complete
X <- PredictMat(x,dat) ## prediction matrix for this term
if (is.null(main)) {
main <- label
}
if (is.null(ylim)) ylim <- range(ym)
if (is.null(xlim)) xlim <- range(xm)
return(list(X=X,x=xm,y=ym,scale=FALSE,se=TRUE,raw=raw,xlab=xlabel,ylab=ylabel,
main=main,se.mult=se2.mult,ylim=ylim,xlim=xlim,exclude=exclude))
} else { ## basic plot data for 3 or 4 d terms
vname <- x$term
## if the smooth has margins and one is 2D then set that as the
## term for 2D plotting, rather than conditioning....
if (!is.null(x$margin)) {
for (i in 1:length(x$margin)) if (x$margin[[i]]$dim==2) {
vname <- x$margin[[i]]$term ## these are the variables to 2d plot
vname <- c(vname,x$term[!x$term%in%vname])
break;
}
}
## ... so first 2 terms in vname are the vars to plot in 2D.
## Now get the limits for plotting...
nv <- length(vname)
lo <- hi <- rep(0,nv)
for (i in 1:length(vname)) {
xx <- data[vname[i]][[1]]
lo[i] <- min(xx);hi[i] <- max(xx)
}
nc <- nr <- n3 ## set number cols and rows of plot
m <- n2 ## 2d plotting grid side
x1 <- seq(lo[1],hi[1],length=m)
x2 <- seq(lo[2],hi[2],length=m)
if (nv==3) {
x3 <- seq(lo[3],hi[3],length=nr*nc)
dat <- cbind(rep(x1,m*nr*nc),
rep(rep(x2,each=m*nr),nc),
x3[rep(rep((1:nr-1)*nc,each=m),m*nc) + rep(1:nc,each=m*m*nr)])
} else {
x3 <- seq(lo[3],hi[3],length=nr)
x4 <- seq(lo[4],hi[4],length=nc)
dat <- cbind(rep(x1,m*nr*nc),
rep(rep(x2,each=m*nr),nc),
rep(rep(x3,each=m),m*nc),
rep(x4,each=m*m*nr))
} ## 4D term end
if (x$by!="NA") {
dat <- data.frame(cbind(dat,1))
names(dat) <- c(vname,x$by)
} else {
dat <- data.frame(dat)
names(dat) <- vname
}
X <- PredictMat(x,dat) ## prediction matrix for this term
exclude <- if (too.far<=0) rep(FALSE,nrow(X)) else
exclude.too.far(dat[,1],dat[,2],data[vname[1]][[1]],data[vname[2]][[1]],dist=too.far)
if (is.null(main)) {
main <- label
}
return(list(X=X,scale=FALSE,se=FALSE,m=m,nc=nc,nr=nr,lo=lo,hi=hi,vname=vname,
main=main,exclude=exclude))
} ## end of 3/4 D case
} else { ## produce plot
if (se) { ## produce CI's
if (x$dim==1) {
if (scheme == 1) shade <- TRUE
ul <- P$fit + P$se ## upper CL
ll <- P$fit - P$se ## lower CL
if (scale==0&&is.null(ylim)) { ## get scale
ylimit<-c(min(ll),max(ul))
if (partial.resids) {
max.r <- max(P$p.resid,na.rm=TRUE)
if ( max.r> ylimit[2]) ylimit[2] <- max.r
min.r <- min(P$p.resid,na.rm=TRUE)
if (min.r < ylimit[1]) ylimit[1] <- min.r
}
}
ylimit <- if (is.null(ylim)) ylimit <- trans(ylimit + shift) else ylim
## plot the smooth...
if (shade) {
plot(P$x,trans(P$fit+shift),type="n",xlab=P$xlab,ylim=ylimit,
xlim=P$xlim,ylab=P$ylab,main=P$main,...)
polygon(c(P$x,P$x[n:1],P$x[1]),
trans(c(ul,ll[n:1],ul[1])+shift),col = shade.col,border = NA)
lines(P$x,trans(P$fit+shift),...)
} else { ## ordinary plot
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,ylim=ylimit,xlim=P$xlim,
ylab=P$ylab,main=P$main,...)
if (is.null(list(...)[["lty"]])) {
lines(P$x,trans(ul+shift),lty=2,...)
lines(P$x,trans(ll+shift),lty=2,...)
} else {
lines(P$x,trans(ul+shift),...)
lines(P$x,trans(ll+shift),...)
}
} ## ... smooth plotted
if (partial.resids&&(by.resids||x$by=="NA")) { ## add any partial residuals
if (length(P$raw)==length(P$p.resid)) {
if (is.null(list(...)[["pch"]]))
points(P$raw,trans(P$p.resid+shift),pch=".",...) else
points(P$raw,trans(P$p.resid+shift),...)
} else {
warning("Partial residuals do not have a natural x-axis location for linear functional terms")
}
} ## partial residuals finished
if (rug) {
if (jit) rug(jitter(as.numeric(P$raw)),...)
else rug(as.numeric(P$raw),...)
} ## rug plot done
} else if (x$dim==2) {
P$fit[P$exclude] <- NA
if (scheme == 1) { ## perspective plot
persp(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
zlab=P$main,ylim=P$ylim,xlim=P$xlim,theta=theta,phi=phi,...)
} else if (scheme==2||scheme==3) {
if (scheme==3) hcolors <- grey(0:50/50)
image(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,col=hcolors,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=contour.col,...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
}
} else { ## contour plot with error contours
sp.contour(P$x,P$y,matrix(P$fit,n2,n2),matrix(P$se,n2,n2),
xlab=P$xlab,ylab=P$ylab,zlab=P$main,titleOnly=!is.null(main),
se.mult=1,trans=trans,shift=shift,...)
if (rug) {
if (is.null(list(...)[["pch"]]))
points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
}
} ## contour plot done
} else if (x$dim<5) {
if (scheme==1) hcolors <- grey(0:50/50)
md.plot(P$fit,P$nr,P$nc,P$m,P$vname,P$lo,P$hi,hcolors=hcolors,scheme=scheme,P$main,...)
} else {
warning("no automatic plotting for smooths of more than two variables")
}
} else { ## no CI's
if (x$dim==1) {
if (scale==0&&is.null(ylim)) {
if (partial.resids) ylimit <- range(P$p.resid,na.rm=TRUE) else ylimit <-range(P$fit)
}
ylimit <- if (is.null(ylim)) ylimit <- trans(ylimit + shift) else ylim
plot(P$x,trans(P$fit+shift),type="l",xlab=P$xlab,
ylab=P$ylab,ylim=ylimit,xlim=P$xlim,main=P$main,...)
if (rug) {
if (jit) rug(jitter(as.numeric(P$raw)),...)
else rug(as.numeric(P$raw),...)
}
if (partial.resids&&(by.resids||x$by=="NA")) {
if (is.null(list(...)[["pch"]]))
points(P$raw,trans(P$p.resid+shift),pch=".",...) else
points(P$raw,trans(P$p.resid+shift),...)
}
} else if (x$dim==2) {
P$fit[P$exclude] <- NA
if (!is.null(main)) P$title <- main
if (scheme==1) {
persp(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
zlab=P$main,theta=theta,phi=phi,xlim=P$xlim,ylim=P$ylim,...)
} else if (scheme==2||scheme==3) {
if (scheme==3) hcolors <- grey(0:50/50)
image(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,col=hcolors,...)
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),add=TRUE,col=contour.col,...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
}
} else {
contour(P$x,P$y,matrix(trans(P$fit+shift),n2,n2),xlab=P$xlab,ylab=P$ylab,
main=P$main,xlim=P$xlim,ylim=P$ylim,...)
if (rug) {
if (is.null(list(...)[["pch"]])) points(P$raw$x,P$raw$y,pch=".",...) else
points(P$raw$x,P$raw$y,...)
}
}
} else if (x$dim<5) {
if (scheme==1) hcolors <- grey(0:50/50)
md.plot(P$fit,P$nr,P$nc,P$m,P$vname,P$lo,P$hi,hcolors=hcolors,scheme=scheme,P$main,...)
} else {
warning("no automatic plotting for smooths of more than four variables")
}
} ## end of no CI code
} ## end of plot production
} ## plot.mgcv.smooth
md.plot <- function(f,nr,nc,m,vname,lo,hi,hcolors,scheme,main,...) {
## multi-dimensional term plotter, called from plot.mgcv.smooth for
## 3 and 4 dimensional terms.
## *f is the plot data. See `basic plot data for 3 or 4 d terms'
## in plot.mgcv.smooth for details of the packing conventions
## (f = X %*% coefs).
## *nr and nc the number of rows and columns of plot panels
## *m each panel is m by m
## *vname contains the variable names
## *lo and hi are the arrays of axis limits
## *hcolors is the color palette for the image plot.
## *scheme indicates b/w or color
## *main is a title.
concol <- if (scheme==1) "white" else "black"
nv <- length(vname)
## insert NA breaks to separate the panels within a plot...
f1 <- matrix(NA,nr*m+nr-1,nc*m)
ii <- rep(1:m,nr) + rep(0:(nr-1)*(m+1),each=m)
f1[ii,] <- f
f <- matrix(NA,nr*m+nr-1,nc*m+nc-1)
ii <- rep(1:m,nc) + rep(0:(nc-1)*(m+1),each=m)
f[,ii] <- f1
xx <- seq(0,1,length=ncol(f))
yy <- seq(0,1,length=nrow(f))
image(xx,yy,t(f),axes=FALSE,xlab="",ylab="",col=hcolors)
contour(xx,yy,t(f),add=TRUE,col=concol)
dl <- list(...)
c1 <- if (is.null(dl[["cex"]])) 1 else dl[["cex"]]
c2 <- if (is.null(dl[["cex.axis"]])) .6 else dl[["cex.axis"]]
c3 <- if (is.null(dl[["cex.lab"]])) .9 else dl[["cex.lab"]]
if (nv==4) {
x3 <- seq(lo[3],hi[3],length=nr)
x4 <- seq(lo[4],hi[4],length=nc)
mtext(vname[4],1,1.7,cex=c1*c3) ## x label
mtext(vname[3],2,1.7,cex=c1*c3) ## y label
at=(1:nc-.5)/nc
lab <- format(x4,digits=2)
for (i in 1:nc) mtext(lab[i],1,at=at[i],line=.5,cex=c1*c3)
at=(1:nr-.5)/nr
lab <- format(x4,digits=2)
for (i in 1:nr) mtext(lab[i],2,at=at[i],line=.5,cex=c1*c3)
## now the 2d panel axes...
xr <- axisTicks(c(lo[2],hi[2]),log=FALSE,nint=4)
x0 <- ((nc-1)*(m+1)+1)/(nc*m+nc-1)
xt <- (xr-lo[2])/(hi[2]-lo[2])*(1-x0)+x0
axis(3,at=xt,labels=as.character(xr),cex.axis=c2,cex=c1)
xr <- axisTicks(c(lo[1],hi[1]),log=FALSE,nint=4)
x0 <- ((nr-1)*(m+1)+1)/(nr*m+nr-1)
xt <- (xr-lo[1])/(hi[1]-lo[1])*(1-x0)+x0
axis(4,at=xt,labels=as.character(xr),cex.axis=c2,cex=c1)
at <- (2*nc-3)/(2*nc)
mtext(vname[2],3,at=at,line=.5,cex=c1*c2)
at <- (2*nr-3)/(2*nr)
mtext(vname[1],4,at=at,line=.5,cex=c1*c2)
mtext(main,3,at=0,adj=0,line=1,cex=c1*c3)
} else {
x3 <- seq(lo[3],hi[3],length=nr*nc)
## get pretty ticks
xr <- axisTicks(c(lo[2],hi[2]),log=FALSE,nint=4)
x0 <- (m-1)/(nc*m+nc-1)
xt <- (xr-lo[2])/(hi[2]-lo[2])*x0
axis(1,at=xt,labels=as.character(xr),cex.axis=c2,cex=c1)
mtext(vname[2],1,at=x0/2,line=2,cex=c1*c2)
xr <- axisTicks(c(lo[1],hi[1]),log=FALSE,nint=4)
x0 <- (m-1)/(nr*m+nr-1)
xt <- (xr-lo[1])/(hi[1]-lo[1])*x0
axis(2,at=xt,labels=as.character(xr),cex.axis=c2,cex=c1)
mtext(vname[1],2,at=x0/2,line=2,cex=c1*c2)
lab <- c("",format(x3[-1],digits=2))
at=(1:nc-.5)/nc
for (i in 2:nc) mtext(lab[i],1,at=at[i],line=.5,cex=c1*c3)
mtext(parse(text=paste(vname[3],"%->% \" \"")),1,at=mean(at[2:nc]),line=2,cex=c1*c3)
ii <- ((nr-1)*nr+1):(nc*nr)
for (i in 1:nc) mtext(lab[ii[i]],3,at=at[i],line=.5,cex=c1*c3)
mtext(parse(text=paste(vname[3],"%->% \" \"")),3,at=mean(at),line=2,cex=c1*c3)
mtext(main,2,at=1/nr+0.5*(nr-1)/nr,line=1,cex=c1*c3)
}
} ## md.plot
plot.gam <- function(x,residuals=FALSE,rug=NULL,se=TRUE,pages=0,select=NULL,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,all.terms=FALSE,shade=FALSE,shade.col="gray80",
shift=0,trans=I,seWithMean=FALSE,unconditional=FALSE,by.resids=FALSE,scheme=0,...)
# Create an appropriate plot for each smooth term of a GAM.....
# x is a gam object
# rug determines whether a rug plot should be added to each plot
# se determines whether twice standard error bars are to be added
# pages is the number of pages over which to split output - 0 implies that
# graphic settings should not be changed for plotting
# scale -1 for same y scale for each plot
# 0 for different y scales for each plot
# n - number of x axis points to use for plotting each term
# n2 is the square root of the number of grid points to use for contouring
# 2-d terms.
{ ######################################
## Local function for producing labels
######################################
sub.edf <- function(lab,edf) {
## local function to substitute edf into brackets of label
## labels are e.g. smooth[[1]]$label
pos <- regexpr(":",lab)[1]
if (pos<0) { ## there is no by variable stuff
pos <- nchar(lab) - 1
lab <- paste(substr(lab,start=1,stop=pos),",",round(edf,digits=2),")",sep="")
} else {
lab1 <- substr(lab,start=1,stop=pos-2)
lab2 <- substr(lab,start=pos-1,stop=nchar(lab))
lab <- paste(lab1,",",round(edf,digits=2),lab2,sep="")
}
lab
} ## end of sub.edf
#########################
## start of main function
#########################
if (is.null(rug)) rug <- if (nrow(x$model)>10000) FALSE else TRUE
if (unconditional) {
if (is.null(x$Vc)) warning("Smoothness uncertainty corrected covariance not available") else
x$Vp <- x$Vc ## cov matrix reset to full Bayesian
}
w.resid <- NULL
if (length(residuals)>1) { # residuals supplied
if (length(residuals)==length(x$residuals))
w.resid <- residuals else
warning("residuals argument to plot.gam is wrong length: ignored")
partial.resids <- TRUE
} else partial.resids <- residuals # use working residuals or none
m <- length(x$smooth) ## number of smooth terms
if (length(scheme)==1) scheme <- rep(scheme,m)
if (length(scheme)!=m) {
warn <- paste("scheme should be a single number, or a vector with",m,"elements")
warning(warn)
scheme <- rep(scheme[1],m)
}
## array giving order of each parametric term...
order <- if (is.list(x$pterms)) unlist(lapply(x$pterms,attr,"order")) else attr(x$pterms,"order")
if (all.terms) # plot parametric terms as well
n.para <- sum(order==1) # plotable parametric terms
else n.para <- 0
if (se) ## sort out CI widths for 1 and 2D
{ if (is.numeric(se)) se2.mult <- se1.mult <- se else { se1.mult <- 2;se2.mult <- 1}
if (se1.mult<0) se1.mult<-0;if (se2.mult < 0) se2.mult <- 0
} else se1.mult <- se2.mult <-1
if (se && x$Vp[1,1] < 0) ## check that variances are actually available
{ se <- FALSE
warning("No variance estimates available")
}
if (partial.resids) { ## getting information needed for partial residuals...
if (is.null(w.resid)) { ## produce working resids if info available
if (is.null(x$residuals)||is.null(x$weights)) partial.resids <- FALSE else {
wr <- sqrt(x$weights)
w.resid <- x$residuals*wr#/mean(wr) # weighted working residuals
}
}
if (partial.resids) fv.terms <- predict(x,type="terms") ## get individual smooth effects
}
pd <- list(); ## plot data list
i <- 1 # needs a value if no smooths, but parametric terms ...
##################################################
## First the loop to get the data for the plots...
##################################################
if (m>0) for (i in 1:m) { ## work through smooth terms
first <- x$smooth[[i]]$first.para
last <- x$smooth[[i]]$last.para
edf <- sum(x$edf[first:last]) ## Effective DoF for this term
term.lab <- sub.edf(x$smooth[[i]]$label,edf)
#P <- plot(x$smooth[[i]],P=NULL,data=x$model,n=n,n2=n2,xlab=xlab,ylab=ylab,too.far=too.far,label=term.lab,
# se1.mult=se1.mult,se2.mult=se2.mult,xlim=xlim,ylim=ylim,main=main,scheme=scheme[i],...)
attr(x$smooth[[i]],"coefficients") <- x$coefficients[first:last] ## relevent coefficients
P <- plot(x$smooth[[i]],P=NULL,data=x$model,partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,n3=n3,
theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,label=term.lab,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
se1.mult=se1.mult,se2.mult=se2.mult,shift=shift,trans=trans,
by.resids=by.resids,scheme=scheme[i],...)
if (is.null(P)) pd[[i]] <- list(plot.me=FALSE) else if (is.null(P$fit)) {
p <- x$coefficients[first:last] ## relevent coefficients
offset <- attr(P$X,"offset") ## any term specific offset
## get fitted values ....
if (is.null(offset)) P$fit <- P$X%*%p else P$fit <- P$X%*%p + offset
if (!is.null(P$exclude)) P$fit[P$exclude] <- NA
if (se && P$se) { ## get standard errors for fit
## test whether mean variability to be added to variability (only for centred terms)
if (seWithMean && attr(x$smooth[[i]],"nCons")>0) {
if (length(x$cmX) < ncol(x$Vp)) x$cmX <- c(x$cmX,rep(0,ncol(x$Vp)-length(x$cmX)))
if (seWithMean==2) x$cmX[-(1:x$nsdf)] <- 0 ## variability of fixed effects mean only
X1 <- matrix(x$cmX,nrow(P$X),ncol(x$Vp),byrow=TRUE)
meanL1 <- x$smooth[[i]]$meanL1
if (!is.null(meanL1)) X1 <- X1 / meanL1
X1[,first:last] <- P$X
lpi <- attr(x$formula,"lpi")
if (is.null(lpi)) se.fit <- sqrt(pmax(0,rowSums(as(X1%*%x$Vp,"matrix")*X1))) else {
ii <- rep(0,0) ## only include constant uncertainty from relevant linear predictors
for (q in 1:length(lpi)) if (any(first:last%in%lpi[[q]])) ii <- c(ii,lpi[[q]])
se.fit <- sqrt(pmax(0,rowSums(as(X1[,ii]%*%x$Vp[ii,ii],"matrix")*X1[,ii])))
}
} else se.fit <- ## se in centred (or anyway unconstained) space only
sqrt(pmax(0,rowSums(as(P$X%*%x$Vp[first:last,first:last,drop=FALSE],"matrix")*P$X)))
if (!is.null(P$exclude)) se.fit[P$exclude] <- NA
} ## standard errors for fit completed
if (partial.resids) { P$p.resid <- fv.terms[,length(order)+i] + w.resid }
if (se && P$se) P$se <- se.fit*P$se.mult # Note multiplier
P$X <- NULL
P$plot.me <- TRUE
pd[[i]] <- P;rm(P)
} else { ## P$fit created directly
if (partial.resids) { P$p.resid <- fv.terms[,length(order)+i] + w.resid }
P$plot.me <- TRUE
pd[[i]] <- P;rm(P)
}
} ## end of data setup loop through smooths
##############################################
## sort out number of pages and plots per page
##############################################
n.plots <- n.para
if (m>0) for (i in 1:m) n.plots <- n.plots + as.numeric(pd[[i]]$plot.me)
if (n.plots==0) stop("No terms to plot - nothing for plot.gam() to do.")
if (pages>n.plots) pages<-n.plots
if (pages<0) pages<-0
if (pages!=0) # figure out how to display things
{ ppp<-n.plots%/%pages
if (n.plots%%pages!=0)
{ ppp<-ppp+1
while (ppp*(pages-1)>=n.plots) pages<-pages-1
}
# now figure out number of rows and columns
c <- r <- trunc(sqrt(ppp))
if (c<1) r <- c <- 1
if (c*r < ppp) c <- c + 1
if (c*r < ppp) r <- r + 1
oldpar<-par(mfrow=c(r,c))
} else
{ ppp<-1;oldpar<-par()}
#####################################
## get a common scale, if required...
#####################################
if (scale==-1&&is.null(ylim)) {
k <- 0
if (m>0) for (i in 1:m) if (pd[[i]]$plot.me&&pd[[i]]$scale) { ## loop through plot data
if (se&&length(pd[[i]]$se)>1) { ## require CIs on plots
ul<-pd[[i]]$fit+pd[[i]]$se
ll<-pd[[i]]$fit-pd[[i]]$se
if (k==0) {
ylim <- c(min(ll,na.rm=TRUE),max(ul,na.rm=TRUE));k <- 1
} else {
if (min(ll,na.rm=TRUE)<ylim[1]) ylim[1] <- min(ll,na.rm=TRUE)
if (max(ul,na.rm=TRUE)>ylim[2]) ylim[2] <- max(ul,na.rm=TRUE)
}
} else { ## no standard errors
if (k==0) {
ylim <- range(pd[[i]]$fit,na.rm=TRUE);k <- 1
} else {
if (min(pd[[i]]$fit,na.rm=TRUE)<ylim[1]) ylim[1] <- min(pd[[i]]$fit,na.rm=TRUE)
if (max(pd[[i]]$fit,na.rm=TRUE)>ylim[2]) ylim[2] <- max(pd[[i]]$fit,na.rm=TRUE)
}
}
if (partial.resids) {
ul <- max(pd[[i]]$p.resid,na.rm=TRUE)
if (ul > ylim[2]) ylim[2] <- ul
ll <- min(pd[[i]]$p.resid,na.rm=TRUE)
if (ll < ylim[1]) ylim[1] <- ll
} ## partial resids done
} ## loop end
ylim <- trans(ylim+shift)
} ## end of common scale computation
##############################################################
## now plot smooths, by calling plot methods with plot data...
##############################################################
if ((pages==0&&prod(par("mfcol"))<n.plots&&dev.interactive())||
pages>1&&dev.interactive()) ask <- TRUE else ask <- FALSE
if (!is.null(select)) {
ask <- FALSE
}
# if (ask) { ## asks before plotting
# oask <- devAskNewPage(TRUE)
# on.exit(devAskNewPage(oask))
# }
if (m>0) for (i in 1:m) if (pd[[i]]$plot.me&&(is.null(select)||i==select)) {
plot(x$smooth[[i]],P=pd[[i]],partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,n3=n3,
theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
shift=shift,trans=trans,by.resids=by.resids,scheme=scheme[i],...)
if (ask) { ## this is within loop so we don't get asked before it's necessary
oask <- devAskNewPage(TRUE)
on.exit(devAskNewPage(oask))
ask <- FALSE ## only need to do this once
}
} ## end of smooth plotting loop
####################################################
## Finally deal with any parametric term plotting...
####################################################
if (n.para>0) # plot parameteric terms
{ class(x) <- c("gam","glm","lm") # needed to get termplot to call model.frame.glm
if (is.null(select)) {
attr(x,"para.only") <- TRUE
termplot(x,se=se,rug=rug,col.se=1,col.term=1,main=attr(x$pterms,"term.labels"),...)
} else { # figure out which plot is required
if (select > m) {
## can't figure out how to get this to work with more than first linear predictor
## as termplots relies on matching terms to names in original data...
select <- select - m # i.e. which parametric term
term.labels <- attr(x$pterms,"term.labels")
term.labels <- term.labels[order==1]
if (select <= length(term.labels)) {
# if (interactive() && m &&i%%ppp==0)
termplot(x,terms=term.labels[select],se=se,rug=rug,col.se=1,col.term=1,...)
}
}
}
}
if (pages>0) par(oldpar)
invisible(pd)
} ## end plot.gam
exclude.too.far <- function(g1,g2,d1,d2,dist)
# if g1 and g2 are the co-ordinates of grid modes and d1,d2 are co-ordinates of data
# then this routine returns a vector with TRUE if the grid node is too far from
# any data and FALSE otherwise. Too far is judged using dist: a positive number indicating
# distance on the unit square into which the grid is scaled prior to calculation
{ mig<-min(g1)
d1<-d1-mig;g1<-g1-mig
mag<-max(g1)
d1<-d1/mag;g1<-g1/mag
mig<-min(g2)
d2<-d2-mig;g2<-g2-mig
mag<-max(g2)
d2<-d2/mag;g2<-g2/mag
# all now in unit square
n<-length(g1)
m<-length(d1)
if (length(g2)!=n) stop("grid vectors are different lengths")
if (m!=length(d2)) stop("data vectors are of different lengths")
if (dist<0) stop("supplied dist negative")
distance<-array(0,n)
o<-.C(C_MinimumSeparation,x=as.double(cbind(g1,g2)),n=as.integer(n), d=as.integer(2),
t=as.double(cbind(d1,d2)),m=as.integer(m),distance=as.double(distance))
res <- rep(FALSE,n)
res[o$distance > dist] <-TRUE
res
} ## exclude.too.far
vis.gam <- function(x,view=NULL,cond=list(),n.grid=30,too.far=0,col=NA,color="heat",
contour.col=NULL,se=-1,type="link",plot.type="persp",zlim=NULL,nCol=50,lp=1,...)
# takes a gam object and plots 2D views of it, supply ticktype="detailed" to get proper axis anotation
# (c) Simon N. Wood 23/2/03
{ fac.seq<-function(fac,n.grid)
# generates a sequence of factor variables of length n.grid
{ fn<-length(levels(fac));gn<-n.grid;
if (fn>gn) mf<-factor(levels(fac))[1:gn]
else
{ ln<-floor(gn/fn) # length of runs
mf<-rep(levels(fac)[fn],gn)
mf[1:(ln*fn)]<-rep(levels(fac),rep(ln,fn))
mf<-factor(mf,levels=levels(fac))
}
mf
}
# end of local functions
dnm <- names(list(...))
## basic issues in the following are that not all objects will have a useful `data'
## component, but they all have a `model' frame. Furthermore, `predict.gam' recognises
## when a model frame has been supplied
v.names <- names(x$var.summary) ## names of all variables
## Note that in what follows matrices in the parametric part of the model
## require special handling. Matrices arguments to smooths are different
## as they follow the summation convention.
if (is.null(view)) # get default view if none supplied
{ ## need to find first terms that can be plotted against
k <- 0;view <- rep("",2)
for (i in 1:length(v.names)) {
ok <- TRUE
if (is.matrix(x$var.summary[[i]])) ok <- FALSE else
if (is.factor(x$var.summary[[i]])) {
if (length(levels(x$var.summary[[i]]))<=1) ok <- FALSE
} else {
if (length(unique(x$var.summary[[i]]))==1) ok <- FALSE
}
if (ok) {
k <- k + 1;view[k] <- v.names[i]
}
if (k==2) break;
}
if (k<2) stop("Model does not seem to have enough terms to do anything useful")
} else {
if (sum(view%in%v.names)!=2) stop(gettextf("view variables must be one of %s",
paste(v.names, collapse = ", ")))
for (i in 1:2)
if (!inherits(x$var.summary[[view[i]]],c("numeric","factor")))
stop("Don't know what to do with parametric terms that are not simple numeric or factor variables")
}
ok <- TRUE
for (i in 1:2) if (is.factor(x$var.summary[[view[i]]])) {
if (length(levels(x$var.summary[[view[i]]]))<=1) ok <- FALSE
} else {
if (length(unique(x$var.summary[[view[i]]]))<=1) ok <- FALSE
}
if (!ok) stop(gettextf("View variables must contain more than one value. view = c(%s,%s).",
view[1], view[2]))
# now get the values of the variables which are not the arguments of the plotted surface
# Make dataframe....
if (is.factor(x$var.summary[[view[1]]]))
m1<-fac.seq(x$var.summary[[view[1]]],n.grid)
else { r1<-range(x$var.summary[[view[1]]]);m1<-seq(r1[1],r1[2],length=n.grid)}
if (is.factor(x$var.summary[[view[2]]]))
m2<-fac.seq(x$var.summary[[view[2]]],n.grid)
else { r2<-range(x$var.summary[[view[2]]]);m2<-seq(r2[1],r2[2],length=n.grid)}
v1<-rep(m1,n.grid);v2<-rep(m2,rep(n.grid,n.grid))
newd <- data.frame(matrix(0,n.grid*n.grid,0)) ## creating prediction data frame full of conditioning values
for (i in 1:length(x$var.summary)) {
ma <- cond[[v.names[i]]]
if (is.null(ma)) {
ma <- x$var.summary[[i]]
if (is.numeric(ma)) ma <- ma[2] ## extract median
}
if (is.matrix(x$var.summary[[i]])) newd[[i]] <- matrix(ma,n.grid*n.grid,ncol(x$var.summary[[i]]),byrow=TRUE)
else newd[[i]]<-rep(ma,n.grid*n.grid)
}
names(newd) <- v.names
newd[[view[1]]]<-v1
newd[[view[2]]]<-v2
# call predict.gam to get predictions.....
if (type=="link") zlab <- paste("linear predictor") ## ignore codetools
else if (type=="response") zlab <- type
else stop("type must be \"link\" or \"response\"")
fv <- predict.gam(x,newdata=newd,se.fit=TRUE,type=type)
z <- fv$fit # store NA free copy now
if (is.matrix(z)) {
lp <- min(ncol(z),max(1,round(lp)))
z <- z[,lp] ## retain selected linear predictor
fv$fit <- fv$fit[,lp]
fv$se.fit <- fv$se.fit[,lp]
}
if (too.far>0) { # exclude predictions too far from data
ex.tf <- exclude.too.far(v1,v2,x$model[,view[1]],x$model[,view[2]],dist=too.far)
fv$se.fit[ex.tf] <- fv$fit[ex.tf] <- NA
}
# produce a continuous scale in place of any factors
if (is.factor(m1)) {
m1<-as.numeric(m1);m1<-seq(min(m1)-0.5,max(m1)+0.5,length=n.grid)
}
if (is.factor(m2)) {
m2<-as.numeric(m2);m2<-seq(min(m1)-0.5,max(m2)+0.5,length=n.grid)
}
if (se<=0) {
old.warn<-options(warn=-1)
av <- matrix(c(0.5,0.5,rep(0,n.grid-1)),n.grid,n.grid-1)
options(old.warn)
# z is without any exclusion of gridpoints, so that averaging works nicely
max.z <- max(z,na.rm=TRUE)
z[is.na(z)] <- max.z*10000 # make sure NA's don't mess it up
z <- matrix(z,n.grid,n.grid) # convert to matrix
surf.col <- t(av)%*%z%*%av # average over tiles
surf.col[surf.col>max.z*2] <- NA # restore NA's
# use only non-NA data to set colour limits
if (!is.null(zlim)) {
if (length(zlim)!=2||zlim[1]>=zlim[2]) stop("Something wrong with zlim")
min.z <- zlim[1]
max.z <- zlim[2]
} else {
min.z <- min(fv$fit,na.rm=TRUE)
max.z <- max(fv$fit,na.rm=TRUE)
}
if (min.z==max.z) {min.z <- min.z-1;max.z <- max.z + 1}
surf.col <- surf.col-min.z
surf.col <- surf.col/(max.z-min.z)
surf.col <- round(surf.col*nCol)
con.col <- 1
if (color=="heat") { pal<-heat.colors(nCol);con.col<-4;}
else if (color=="topo") { pal<-topo.colors(nCol);con.col<-2;}
else if (color=="cm") { pal<-cm.colors(nCol);con.col<-1;}
else if (color=="terrain") { pal<-terrain.colors(nCol);con.col<-2;}
else if (color=="gray"||color=="bw") {pal <- gray(seq(0.1,0.9,length=nCol));con.col<-1}
else stop("color scheme not recognised")
if (is.null(contour.col)) contour.col <- con.col # default colour scheme
surf.col[surf.col<1] <- 1; surf.col[surf.col>nCol] <- nCol # otherwise NA tiles can get e.g. -ve index
if (is.na(col)) col <- pal[as.array(surf.col)]
z <- matrix(fv$fit,n.grid,n.grid)
if (plot.type=="contour")
{ stub <- paste(ifelse("xlab" %in% dnm, "" , ",xlab=view[1]"),
ifelse("ylab" %in% dnm, "" , ",ylab=view[2]"),
ifelse("main" %in% dnm, "" , ",main=zlab"),",...)",sep="")
if (color!="bw")
{ txt <- paste("image(m1,m2,z,col=pal,zlim=c(min.z,max.z)",stub,sep="") # assemble image() call
eval(parse(text=txt))
txt <- paste("contour(m1,m2,z,col=contour.col,zlim=c(min.z,max.z)",
ifelse("add" %in% dnm, "" , ",add=TRUE"),",...)" , sep="") # assemble contour() call
eval(parse(text=txt))
} else
{ txt <- paste("contour(m1,m2,z,col=1,zlim=c(min.z,max.z)",stub,sep="") # assemble contour() call
eval(parse(text=txt))
}
} else
{ stub <- paste(ifelse("xlab" %in% dnm, "" , ",xlab=view[1]"),
ifelse("ylab" %in% dnm, "" , ",ylab=view[2]"),
ifelse("zlab" %in% dnm, "" , ",zlab=zlab"),",...)",sep="")
if (color=="bw")
{ op <- par(bg="white")
txt <- paste("persp(m1,m2,z,col=\"white\",zlim=c(min.z,max.z) ",stub,sep="") # assemble persp() call
eval(parse(text=txt))
par(op)
} else
{ txt <- paste("persp(m1,m2,z,col=col,zlim=c(min.z,max.z)",stub,sep="") # assemble persp() call
eval(parse(text=txt))
}
}
} else # add standard error surfaces
{ if (color=="bw"||color=="gray")
{ subs <- paste("grey are +/-",se,"s.e.") ## ignore codetools
lo.col <- "gray" ## ignore codetools claims about this
hi.col <- "gray" ## ignore codetools
} else
{ subs <- paste("red/green are +/-",se,"s.e.")
lo.col <- "green"
hi.col <- "red"
}
if (!is.null(zlim)) {
if (length(zlim)!=2||zlim[1]>=zlim[2]) stop("Something wrong with zlim")
min.z <- zlim[1]
max.z <- zlim[2]
} else {
max.z <- max(fv$fit+fv$se.fit*se,na.rm=TRUE)
min.z <- min(fv$fit-fv$se.fit*se,na.rm=TRUE)
zlim<-c(min.z,max.z)
}
z <- fv$fit - fv$se.fit*se; z <- matrix(z,n.grid,n.grid)
if (plot.type=="contour") warning("sorry no option for contouring with errors: try plot.gam")
stub <- paste(ifelse("xlab" %in% dnm, "" , ",xlab=view[1]"),
ifelse("ylab" %in% dnm, "" , ",ylab=view[2]"),
ifelse("zlab" %in% dnm, "" , ",zlab=zlab"),
ifelse("sub" %in% dnm, "" , ",sub=subs"),
",...)",sep="")
txt <- paste("persp(m1,m2,z,col=col,zlim=zlim",
ifelse("border" %in% dnm, "" ,",border=lo.col"),
stub,sep="") # assemble persp() call
eval(parse(text=txt))
par(new=TRUE) # don't clean device
z <- fv$fit; z <- matrix(z,n.grid,n.grid)
txt <- paste("persp(m1,m2,z,col=col,zlim=zlim",
ifelse("border" %in% dnm, "" ,",border=\"black\""),
stub,sep="")
eval(parse(text=txt))
par(new=TRUE) # don't clean device
z <- fv$fit+se*fv$se.fit; z <- matrix(z,n.grid,n.grid)
txt <- paste("persp(m1,m2,z,col=col,zlim=zlim",
ifelse("border" %in% dnm, "" ,",border=hi.col"),
stub,sep="")
eval(parse(text=txt))
}
} ## vis.gam
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.