## ## ## AUTHORS: JON WAKEFIELD, NICKY BEST AND SEBASTIEN HANEUSE, APRIL 6TH, 2006. ## PrettyPoly <- function(y, poly, nrepeats, ncut=1000, nlevels=10, lower=NULL, upper=NULL ) { # First replicate data so that no of polygons equals no of data items if (length(nrepeats)!=length(y)) cat("Different array lengths in RepeatPoly\n") x <- rep(y,nrepeats) ## Create id indicators for coloring scheme ## if (is.null(lower)) lower <- min(x) #- 0.0001 if (is.null(upper)) upper <- max(x) #+ 0.0001 id <- cut(x, breaks=seq(from=lower, to=upper, length=(ncut+1))) id <- as.numeric(id) id[is.na(id)] <- 0 id <- id + 1 ## Make the scale of the two axes the same ## xrnge <- range(poly[, 1], na.rm= T) yrnge <- range(poly[, 2], na.rm=T) xd <- xrnge[2] - xrnge[1] yd <- yrnge[2] - yrnge[1] if(xd > yd) { xplot <- xrnge yplot <- NULL yplot[1] <- ((yrnge[2] + yrnge[1])/2) - xd/2 yplot[2] <- ((yrnge[2] + yrnge[1])/2) + xd/2 } if(xd <= yd) { yplot <- yrnge xplot <- NULL xplot[1] <- ((xrnge[2] + xrnge[1])/2) - yd/2 xplot[2] <- ((xrnge[2] + xrnge[1])/2) + yd/2 } ## Set colours to grey scale ## palette(gray(seq(1,0, len=(ncut+1)))) ## Plot ## def.par <- par(no.readonly = TRUE) layout(matrix(c(1,2),ncol=2,nrow=1), heights=c(.3,.3), widths=c(.4,.1)) ## plot(poly, xlim=xplot, ylim=yplot, type="n", xlab="Eastings (km)",ylab="Northings (km)") polygon(poly, col=id, density=-1) ## Legend ## plot(c(0,1), c(0,1), type="n", axes=F, xlab="", ylab="") xlims <- rep(0, nlevels) atpts <- rep(0, nlevels) for(i in 1:nlevels) { # xlims[i] <- format(min(x)+(i-1)*(max(x)-min(x))/(nlevels-1),digits=2) xlims[i] <- format(lower+(i-1)*(upper-lower)/(nlevels-1),digits=2) atpts[i] <- (i-1)/(nlevels-1) } axis(2, at=c(atpts[1:nlevels]), labels=c(xlims[1:nlevels])) yb <- seq(0, (nlevels-2)/(nlevels-1), 1/(nlevels-1)) yt <- seq(1/(nlevels-1), 1, 1/(nlevels-1)) xl <- rep(0, nlevels-1) xr <- rep(1, nlevels-1) gr <- seq(0, 1, 1/nlevels) gr <- max(gr) - gr rect(xleft=xl, ybottom=yb, xright=xr, ytop=yt, col=gray(gr), border=T) ## Reset defaults ## palette("default") par(def.par) invisible() } # AUTHOR: JON WAKEFIELD, FEB 16TH, 2006. # # Function to result empirical Bayes estimates. The parameters of the # negative binomial marginal are estimated via MLE, MASS library is needed. # Inputs are: # Y = no of cases in each area # E = expected nos in each area # Xmat = n x p matrix of ecological covariates where n is the no of areas # and p is the number of covariates (not including the intercept). # If missing, intercept only model is fitted. # Outputs are: # RR = ecological EB posterior mean relative risk estimates # RRmed = ecological EB posterior median relative risk estimates # beta = MLEs of the regression coefficients # alpha = MLE of neg bin dispersion parameter, where the marginal mean # is given by mu(1+mu/alpha) # SMR = standardized moratility/morbifity ratio, Y/E. # eBayes <- function(Y,E,Xmat=NULL){ if (is.null(Xmat)) {mod <- glm.nb(Y ~ 1+offset(log(E)), link=log)} else {mod <- glm.nb(Y ~ Xmat+offset(log(E)), link=log)} alpha <- mod\$theta muhat <- mod\$fitted/E wgt <- E*muhat/(alpha+E*muhat); SMR <- Y/E RR <- as.numeric(wgt*SMR + (1-wgt)*muhat) RRmed <- qgamma(0.5,alpha+Y,(alpha+E*muhat)/muhat) list(RR=RR,RRmed=RRmed,beta=mod\$coeff,alpha=alpha,SMR=SMR) } # # AUTHOR: JON WAKEFIELD, FEB 16TH, 2006. # This function produces plots of emprical Bayes posterior densities which # are gamma distributions with parameters # alpha+Y # and # (alpha+E*mu)/mu # where mu = exp(x beta). The SMRs are drawn on for comparison. # EBpostdens <- function(Y,E,alpha,beta,Xrow=NULL,lower=NULL,upper=NULL,main=""){ if (is.null(Xrow)) Xrow <- matrix(c(1),nrow=1,ncol=1) xvals <- seq(lower,upper,.01) mu <- as.numeric(exp(Xrow %*% beta)) yvals <- dgamma(xvals,alpha+Y,(alpha+E*mu)/mu) plot(xvals,yvals,type="n", xlab=expression(theta),ylab="EB density",cex.lab=1.2) title(paste(main)) lines(xvals,yvals) lines(c(Y/E,Y/E),c(0,max(yvals)),lty=2) } # # AUTHOR: JON WAKEFIELD, FEB 16TH, 2006. # This function produces the probabilities of exceeding a threshold given a # gamma distributions with parameters # alpha+Y # and # (alpha+E*mu)/mu # where mu = exp(x beta). The SMRs are drawn on for comparison. # EBpostthresh <- function(Y,E,alpha,beta,Xrow=NULL,rrthresh){ if (is.null(Xrow)) Xrow <- matrix(rep(1,length(Y)),nrow=length(Y),ncol=1) mu <- as.numeric(exp(Xrow %*% beta)) thresh <- 1-pgamma(rrthresh,alpha+Y,(alpha+E*mu)/mu) return(thresh=thresh) } # # Function to produce maps for Ohio. If type="e" equal-spaced intervals, # if type="q" based on quantiles; if the former then lower and upper can # be specified. # April 10th, 2006. # OhioMap <-function(data, ncol=5, figmain="", digits=5, type="e", lower=NULL, upper=NULL) { if (is.null(lower)) lower <- min(data) if (is.null(upper)) upper <- max(data) if (type=="q"){p <- seq(0,1,length=ncol+1) br <- round(quantile(data,probs=p),2)} if (type=="e"){br <- round(seq(lower,upper,length=ncol+1),2)} shading <- gray((ncol-1):0/(ncol-1)) data.grp <- findInterval(data,vec=br,rightmost.closed=T,all.inside=T) data.shad <- shading[data.grp] map("county", "ohio", fill=TRUE, col=data.shad) leg.txt<-paste("[",br[ncol],",",br[ncol+1],"]",sep="") for(i in (ncol-1):1){ leg.txt<-append(leg.txt,paste("[",br[i],",",br[i+1],")",sep=""),) } leg.txt<-rev(leg.txt) legend(-81.4,39.4,legend=leg.txt,fill=shading,bty="n",ncol=1,cex=.8) title(main=figmain,cex=1.5) invisible() } # # AUTHOR: JON WAKEFIELD, FEB 16TH, 2006. # # Function to obtain parameters, mu and sigma, of a lognormal prior, given # quantiles theta1 and theta2, and distribution function values q1 and q2: # Pr( theta < theta1 ) = q1 and Pr( theta < theta2 ) = q2 # lnprior <- function(theta1,theta2,q1,q2){ zq1 <- qnorm(q1) zq2 <- qnorm(q2) mu <- log(theta1)*zq2/(zq2-zq1) - log(theta2)*zq1/(zq2-zq1) sigma <- (log(theta1)-log(theta2))/(zq1-zq2) list(mu=mu,sigma=sigma) }