#### Plots ##########################################################################################
# Quantile-comparison plots (J. Fox)
# last modified 30 September 2009 by J. Fox
# November 2009 by S. Weisberg -- changed to use showLabels for point identification
# 14 April 2010: set id.n = 0. J. Fox
# 1 June 2010: set reps=100 in qqPlot.lm. J. Fox
# 28 June 2010: fixed labeling bug S. Weisberg
# 11 March 2011: moved up ... argument. J. Fox
qqPlot.default <- function(x, distribution="norm", ..., ylab=deparse(substitute(x)),
xlab=paste(distribution, "quantiles"), main=NULL, las=par("las"),
envelope=.95,
col=palette()[1], col.lines=palette()[2], lwd=2, pch=1, cex=par("cex"),
line=c("quartiles", "robust", "none"),
labels = if(!is.null(names(x))) names(x) else seq(along=x),
id.method = "y",
id.n = if(id.method[1]=="identify") Inf else 0,
id.cex=1, id.col=palette()[1], grid=TRUE, noShow=TRUE)
{
res=list()
line <- match.arg(line)
good <- !is.na(x)
ord <- order(x[good])
ord.x <- x[good][ord]
ord.lab <- labels[good][ord]
q.function <- eval(parse(text=paste("q", distribution, sep="")))
d.function <- eval(parse(text=paste("d", distribution, sep="")))
n <- length(ord.x)
P <- ppoints(n)
z <- q.function(P, ...)
if (!noShow) plot(z, ord.x, type="n", xlab=xlab, ylab=ylab, main=main, las=las)
if(grid & !noShow){
grid(lty=1, equilogs=FALSE)
box()}
if (!noShow) points(z, ord.x, col=col, pch=pch, cex=cex)
if (line == "quartiles" || line == "none"){
Q.x <- quantile(ord.x, c(.25,.75))
Q.z <- q.function(c(.25,.75), ...)
b <- (Q.x[2] - Q.x[1])/(Q.z[2] - Q.z[1])
a <- Q.x[1] - b*Q.z[1]
if (!noShow) abline(a, b, col=col.lines, lwd=lwd)
res$a=a
res$b=b
}
if (line=="robust") {
coef <- coef(rlm(ord.x ~ z))
a <- coef[1]
b <- coef[2]
if (!noShow) abline(a, b)
res$a=a
res$b=b
}
conf <- if (envelope == FALSE) .95 else envelope
zz <- qnorm(1 - (1 - conf)/2)
SE <- (b/d.function(z, ...))*sqrt(P*(1 - P)/n)
fit.value <- a + b*z
upper <- fit.value + zz*SE
lower <- fit.value - zz*SE
res$z=z
res$lower=lower
res$upper=upper
if (envelope != FALSE & !noShow) {
lines(z, upper, lty=2, lwd=lwd, col=col.lines)
lines(z, lower, lty=2, lwd=lwd, col=col.lines)
}
if (!noShow) showLabels(z, ord.x, labels=ord.lab,
id.method = id.method, id.n = id.n, id.cex=id.cex, id.col=id.col)
invisible(res)}
setMethod('plot', signature(x='FLsz',y="missing"),
function(x,y,...){
jk=any(unlist(qapply(x, function(x)
ifelse("jacknife" %in% names(attributes(x)), return(attributes(x)$jacknife), return(FALSE)))))
if (!jk) if (!is.null(attributes(x)$jacknife)) jk=attributes(x)$jacknife
if (jk) p.=plotJK(x,y,...) else
p.=ggplot(merge(getHat(x),as.data.frame(ln2z(x),drop=T)))+
geom_point( aes(year,data,group=iter)) +
#expand_limits(y = 0) +
theme_ms(size=12.5) +
scale_y_continuous(name="Z") +
scale_x_continuous(name="Year") +
geom_line(aes(year,z, group=bit+iter*10000),size=1.5) +
geom_line(aes(year,z+sd,group=bit+iter*10000),colour="red") +
geom_line(aes(year,z-sd,group=bit+iter*10000),colour="red") +
geom_smooth(aes(year,data,group=iter),se=FALSE)
print(p.)
invisible(p.)})
plotJK= function(x,y,...){
dat=merge(getHat(x),as.data.frame(ln2z(x),drop=T))
p.=ggplot(subset(dat,iter>1))+
geom_point( aes(year,data,colour=factor(year))) +
#expand_limits(y = 0) +
theme_ms(size=12.5) +
scale_y_continuous(name="Z") +
scale_x_continuous(name="Year") +
geom_line(aes(year,z, group=iter+bit*100000,colour=factor(iter)),size=1.5) +
geom_smooth(aes(year,data,group=factor(iter),colour=factor(iter)),span=0.75,lty=2,size=0.5,se=FALSE) +
geom_line(aes(year,z,group=bit),size=2,data=subset(dat,iter==1)) +
geom_smooth(aes(year,data),span=0.75,lty=1,size=1.0,se=FALSE,colour="black",data=subset(dat,iter==1)) +
opts(legend.position = "none")
}
setGeneric('diags', function(object,...)
standardGeneric('diags'))
setMethod("diags", signature(object="FLsz"),
function(object,plotIt=TRUE) {
rs1=getHat(object)
rs2=model.frame(FLQuants(y=object@obs,yHat=object@hat,residual=object@residuals))[,c("year","iter","y","yHat","residual")]
res=merge(rs1,rs2,by=c("year","iter"))
rsdls=function(x){
x$residualLag <-c(x$residual[-1],NA)
qq. = qqnorm(c(x$residual),plot.it=FALSE)
x$qqx= qq.$x
x$qqy= qq.$y
return(x)}
res=ddply(res, .(iter), rsdls)[,c("year","bit","iter","sd","z","y","yHat","residual","residualLag","qqx","qqy")]
if (!plotIt) return(res)
tmp1=cbind(res[,c("iter","residual","residualLag","bit")],title="AR(1) Residuals")
names(tmp1)[2:3]=c("x","y")
ref1=cbind(x=range(tmp1$x,na.rm=TRUE),y=predict(lm(y~x,data=tmp1,na.action=na.omit),data.frame(x=range(tmp1$x,na.rm=TRUE))))
qq=qqPlot.default(res$residual)
tmp2=cbind(res[,c("iter","qqx","qqy","bit")],title="Normal Q-Q Plot")
names(tmp2)[2:3]=c("x","y")
ref2=data.frame("x"=c(range(tmp2$x,na.rm=TRUE)))
ref2$y =ref2$x*qq$b + qq$a
tmp3=cbind(res[,c("iter","year","residual","bit")],title="Residuals by Year")
names(tmp3)[2:3]=c("x","y")
ref3=data.frame(x=range(tmp3$x,na.rm=TRUE),y=c(0,0))
tmp4=cbind(res[,c("iter","z","residual","bit")],title="Residuals by Z")
names(tmp4)[2:3]=c("x","y")
ref4=data.frame(x=range(tmp4$x,na.rm=TRUE),y=c(0,0))
tmp5=cbind(res[,c("iter","yHat","residual","bit")],title="Residuals by Length Hat")
names(tmp5)[2:3]=c("x","y")
ref5=data.frame(x=range(tmp5$x,na.rm=TRUE),y=c(0,0))
tmp6=cbind(res[,c("iter","year","z","bit")],title="Z by Year")
names(tmp6)[2:3]=c("x","y")
ttl=c("Z by Year","Normal Q-Q Plot","AR(1) Residuals","Residuals by Year","Residuals by Z","Residuals by Length Hat")
tmp=transform(rbind(tmp1,tmp2,tmp3,tmp4,tmp5,tmp6), bit=factor(bit), title=factor(title,levels=ttl))
ref=transform(cbind(rbind(ref1,ref2,ref3,ref4,ref5),title=rep(ac(unique(tmp$title)[-6]),each=2)), title=factor(title,levels=ttl))
tmp=transform(tmp,ymin=pmin(y,0),ymax=pmax(y,0))
tmp[ac(tmp$title) %in% ttl[1:2] ,c("ymin","ymax")]=NA
p.=ggplot(tmp)+geom_point(aes(x,y,colour=bit),position = position_dodge(width = 0.010))+
geom_smooth(aes(x,y),se=FALSE) +
geom_path(aes(x,y),data=ref, col="red") +
geom_linerange(aes(x,ymin=ymin,ymax=ymax),position = position_dodge(width = 0.010)) +
geom_line(aes(z,lower),data=data.frame(title="Normal Q-Q Plot",qq[3:4]),colour="red",lty=2)+
geom_line(aes(z,upper),data=data.frame(title="Normal Q-Q Plot",qq[c(3,5)]),colour="red",lty=2)+
facet_wrap(~title,scale="free")
print(p.)
invisible(p.)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.