# -------------------------------------------------------------------------------------------------- LL.to.ARBUTM<- function(cm,lat,lon) # This function converts from Lat-Lon (decimal degrees) to the Universal # Transverse Mercator Coordinates and returns the new coordinates in a # 2-column matrix with x- and y- as columns. In this program, the # coordinates are calculated from a user supplied central meridian. # Coordinates are returned in kilometers from the western-most # longitude and the southern-most latitude observed in the data set. { # initialize some variables e2 <- 0.00676865799729 a <- 6378206.4 ep2 <- e2 / (1-e2) drc <- pi / 180 sc <- 0.9996 fe <- 500000 ftm <- 0.30480371 #calculate some frequently used values lar <- lat * drc ls <- sin(lar) ls2 <- ls^2 els2 <- ep2 * ls2 lc <- cos(lar) lc2 <- lc^2 lc3 <- lc^3 lc5 <- lc^5 elc2 <- ep2 * lc2 lt2 <- tan(lar)^2 lt4 <- lt2^2 # do the transformation v <- a/sqrt(1 - e2*ls2) p <- drc*(cm - lon) temp <- 5104.57388 - (lc2*(21.73607 - 0.11422*lc2)) r1 <- 6367399.689*(lar - ls*lc*0.000001*temp) r2 <- (v*ls*lc*p^2)/2 temp <- 5 - lt2 + 9*elc2 + (2*elc2)^2 r3 <- (v*ls*lc3*p^4*temp)/24 r4 <- v*lc*p temp <- 1 - lt2 + elc2 r5 <- (v*lc3*p^3*temp)/6 temp <- 61 - 58*lt2 + lt4 + 270*elc2 - 330*els2 ra6 <- (v*ls*lc5*p^6*temp)/720 temp <- 5 - 18*lt2 + lt4 + 14*elc2 - 58*els2 rb5 <- (v*lc5*p^5*temp)/120 northing <- sc*(r1 + r2 + r3 + ra6) easting <- -sc*(r4 + r5 + rb5) y<-(northing-min(northing))/1000 x<-(easting-min(easting))/1000 cbind(x,y) } #------------------- GET COORDINATES FROM ESRI SHAPE FILE shape.prj.coords <- function(rshape){ tmp <- matrix(0, nrow = length(rshape$att.data[,1]), ncol = 2) for (i in 1:length(rshape$att.data[,1])){ tmp[i,] <- rshape$Shape[[i]]$verts } data.frame(Centroid.lon = tmp[,1], Centroid.lat = tmp[,2]) } #------------------- SOME SIMPLE FUNCTIONS, USED WITH SIZE-COLOR GRAPHS log1 <- function(x){log(x+1)} negf <- function(x){-x} invf <- function(x){1/x} #------------------- GENERALIZED INVERSE OF A MATRIX mginv <- function(X, tol = sqrt(.Machine$double.eps)) { dnx <- dimnames(X) if(is.null(dnx)) dnx <- vector("list", 2) s <- svd(X) nz <- s$d > tol * s$d[1] structure( if(any(nz)) s$v[, nz] %*% (t(s$u[, nz])/s$d[nz]) else X, dimnames = dnx[2:1]) } #------------------- BASIC SIZE-COLOR GRAPH size.color.gr <- function(data, xcol = "x", ycol = "y", sample.size.col, response.col, sample.size.fun = eval, response.fun = eval, cex.start.sample.size = 1, cex.increment.sample.size = 0, nclass.response = 10, nclass.sample = 10, color.palette = NULL, ...) { lower.breaks <- matrix(0, nrow = nclass.response, ncol = 1) upper.breaks <- matrix(0, nrow = nclass.response, ncol = 1) if(length(color.palette) == 0) color.palette <- rainbow(nclass.response, start = .66, end = .99) col.stepi <- (max(sample.size.fun(data[,sample.size.col]),na.rm = TRUE) - min(sample.size.fun(data[,sample.size.col]),na.rm = TRUE))/nclass.sample + 0.0001/(nclass.sample - 1) imax <- min(sample.size.fun(data[,sample.size.col]),na.rm = TRUE) - 0.0001 col.stepj <- (max(response.fun(data[,response.col]),na.rm = TRUE) - min(response.fun(data[,response.col]),na.rm = TRUE))/nclass.response + 0.0001/(nclass.response - 1) jmax1 <- min(response.fun(data[,response.col]),na.rm = TRUE) - 0.0001 for (i in 1:nclass.sample){ imin <- imax imax <- imax + col.stepi indi <- sample.size.fun(data[,sample.size.col]) >= imin & sample.size.fun(data[,sample.size.col]) <= imax for (j in 1:nclass.response){ if(j == 1) jmax <- jmax1 jmin <- jmax lower.breaks[j] <- jmin jmax <- jmax + col.stepj upper.breaks[j] <- jmax indj <- response.fun(data[,response.col]) >= jmin & response.fun(data[,response.col]) <= jmax points(data[indi & indj,xcol],data[indi & indj,ycol], col = color.palette[j], pch = 19, cex = cex.start.sample.size + (i-1)*cex.increment.sample.size) } } data.frame(lower.breaks = lower.breaks, upper.breaks = upper.breaks) } #------------------- EMPIRICAL SEMIVARIOGRAM empirical.semivariogram<- function(data, x, y, var, nlag = 20, directions = c(0,45,90,135), tolerance = 22.5, inc = 0, maxlag = 1e32, nlagcutoff = 1, EmpVarMeth = "MethMoment") # EMPIRICAL SEMIVARIOGRAM FUNCTION # var1 is a matrix or data frame with x-coord in the first column # y-coord in the second column # z (response) in the third column { n1 <- length(data[,1]) # distance matrix among locations distance <- sqrt( ( matrix(data[,x],nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1) - matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(data[,x],nrow=1,ncol=n1) )^2 + ( matrix(data[,y],nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1) - matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(data[,y],nrow=1,ncol=n1) )^2 ) difx <- -(matrix(data[,y],nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1) - matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(data[,y],nrow=1,ncol=n1)) signind <- -(matrix(data[,x],nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1) - matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(data[,x],nrow=1,ncol=n1)) < 0 distance <- distance*1.0000000001 theta.deg<-acos(difx/distance)*180/pi # matrix of degrees clockwise from north between locations theta.deg[signind] <- 360-theta.deg[signind] diff2 <- ( matrix(data[,var],nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1) - matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(data[,var],nrow=1,ncol=n1) )^2 sqrtdiff <- sqrt(abs( matrix(data[,var],nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1) - matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(data[,var],nrow=1,ncol=n1) ) ) if(EmpVarMeth == "CovMean") temp4cov <- data[,var] - mean(data[,var]) else temp4cov <- data[,var] covprod <- (matrix(temp4cov,nrow=n1,ncol=1) %*% matrix(rep(1,times=n1),nrow=1,ncol=n1)) * (matrix(rep(1,times=n1),nrow=n1,ncol=1) %*% matrix(temp4cov,ncol=n1,nrow=1)) # convert to vectors distance <- matrix(distance, ncol = 1) theta.deg <- matrix(theta.deg, ncol = 1) diff2 <- matrix(diff2, ncol = 1) sqrtdiff <- matrix(sqrtdiff, ncol = 1) covprod <- matrix(covprod, ncol = 1) # trim off values greater than maxlag indmax <- distance <= maxlag distance <- distance[indmax,] theta.deg <- theta.deg[indmax,] diff2 <- diff2[indmax,] sqrtdiff <- sqrtdiff[indmax,] covprod <- covprod[indmax,] maxd<-max(distance) if( inc <= 0) inc <- maxd/nlag ind <- distance==0 ndir <- length(directions) store.results <- matrix(data = NA, ncol = 6, dimnames = list(NULL, c("distance", "gamma", "np", "azimuth", "hx", "hy"))) for (j in 1:ndir) { for ( i in 1:nlag){ if( (directions[j]-tolerance)<0 && (directions[j]+tolerance)>0 ) ind1 <- theta.deg >= 360+directions[j]-tolerance | theta.deg < directions[j]+tolerance else if( (directions[j]+tolerance)>360 && (directions[j]-tolerance)<360 ) ind1 <- theta.deg < directions[j]+tolerance-360 | theta.deg >= directions[j]-tolerance else ind1 <- theta.deg >= directions[j]-tolerance & theta.deg < directions[j]+tolerance ind<-distance>(i-1)*inc & distance<=i*inc & !is.na(theta.deg) & ind1 nclass <- sum(ind) if(EmpVarMeth == "MethMoment") cv <- mean(diff2[ind]) if(EmpVarMeth == "RobustMean") cv <- ((mean(sqrtdiff[ind]))^4)/(.457 + .494/sum(ind)) if(EmpVarMeth == "RobustMedian") cv <- (median(sqrtdiff[ind]))^4/.457 if(EmpVarMeth == "Covariance" | EmpVarMeth == "CovMean") cv <- mean(covprod[ind]) mean.dis <- mean(distance[ind]) if(nclass > 0) store.results<-rbind(store.results, c(mean.dis,cv,nclass,directions[j],0,0)) } } store.results[,"hx"]<-store.results[,"distance"]*sin(store.results[,"azimuth"]*pi/180) store.results[,"hy"]<-store.results[,"distance"]*cos(store.results[,"azimuth"]*pi/180) store.results[,"gamma"]<-store.results[,"gamma"]/2 ind <- store.results[,"np"] >= nlagcutoff store.results <- store.results[ind,] ind <- !is.na(store.results[,"hx"]) store.results <- store.results[ind,] as.data.frame(store.results) } # ---------- TRELLIS GRAPH OF (DIRECTIONAL) EMPIRICAL SEMIVARIOGRAM VALUES graph.semivariogram.trellis <- function(data, response.column, xcol, ycol, nlag = 20, directions = c(0, 45, 90, 135), tolerance = 22.5, inc = 0, maxlag = 1e32, nlagcutoff = 1, EmpVarMeth = "MethMoment", cex = 1) { emp.var.out <- empirical.semivariogram(data, x = xcol, y = ycol, var = response.column, nlag = nlag, directions = directions, tolerance = tolerance, inc = inc, nlagcutoff = nlagcutoff, maxlag = maxlag, EmpVarMeth = EmpVarMeth) win.graph() plot(c(1,1),type = "n") t.background <- trellis.par.get("background") t.background$col <- "white" trellis.par.set("background", t.background) xyplot(gamma ~ distance | azimuth, data=emp.var.out, pch = 19, col = "black", cex = cex, as.table = TRUE) } # ---------- TRELLIS GRAPH OF (DIRECTIONAL) EMPIRICAL SEMIVARIOGRAM/COVARIANCE VALUES WITH FITTED MODEL graph.emp.var.fit <- function(data, slm.geo.out, resp.var, x = "x", y = "y", nlag = 15, tolerance = 22.5, directions = c(0, 45, 90, 135), inc = 0, maxlag = 1e32, lwd = 1, EmpVarMeth = "MethMoment", nlagcutoff = 5, cex = 1) { ev <- empirical.semivariogram(data = data, var = resp.var, x = x, y = y, nlag = nlag, tolerance = tolerance, directions = directions, inc = inc, maxlag = maxlag, EmpVarMeth = EmpVarMeth, nlagcutoff = nlagcutoff) nugget <- 0 if(any(names(slm.geo.out$theta) == "nugget")) nugget <- slm.geo.out$theta[names(slm.geo.out$theta) == "nugget"] parsil <- 0 if(any(names(slm.geo.out$theta) == "parsil")) parsil <- slm.geo.out$theta[names(slm.geo.out$theta) == "parsil"] range <- 0 if(any(names(slm.geo.out$theta) == "range")) range <- slm.geo.out$theta[names(slm.geo.out$theta) == "range"] rotate <- 90 if(any(names(slm.geo.out$theta) == "rotate")) rotate <- slm.geo.out$theta[names(slm.geo.out$theta) == "rotate"] minorp <- 1 if(any(names(slm.geo.out$theta) == "minorp")) minorp <- slm.geo.out$theta[names(slm.geo.out$theta) == "minorp"] extrap <- 0 if(any(names(slm.geo.out$theta) == "extrap")) extrap <- slm.geo.out$theta[names(slm.geo.out$theta) == "extrap"] CorModel <- slm.geo.out$CorModel # rotate coordinates newx <- cos(rotate*.0174533)*ev[,"hx"] - sin(rotate*.0174533)*ev[,"hy"] newy <- sin(rotate*.0174533)*ev[,"hx"] + cos(rotate*.0174533)*ev[,"hy"] # scale coordinates by minor and major axes */ newx <- newx/(range*minorp) newy <- newy/range dismat <- sqrt(newx^2 + newy^2) # compute distance for the scaled and rotated coordinates */ if(CorModel == "Exponential") cm <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") cm <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") cm <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") cm <- CorModel.Gaussian(dismat) if(CorModel == "Stable") cm <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") cm <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") cm <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") cm <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") cm <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") cm <- CorModel.Circular(dismat) if(CorModel == "Spherical") cm <- CorModel.Spherical(dismat) if(CorModel == "Cubic") cm <- CorModel.Cubic(dismat) if(CorModel == "Penta") cm <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") cm <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") cm <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") cm <- CorModel.BesselJ(dismat, extrap) if(EmpVarMeth == "Covariance" | EmpVarMeth == "CovMean") {fit <- parsil*cm ylabtrell <- "Covariance" } else { fit <- nugget + parsil*(1 - cm) ylabtrell <- "Semivariogram" } ev <- data.frame(ev, fit = fit) t.background <- trellis.par.get("background") t.background$col <- "white" trellis.par.set("background", t.background) ev[,"azimuth"] <- as.factor(ev[,"azimuth"]) win.graph() plot(c(1,1), type = "n") xyplot(gamma + fit ~ distance | azimuth, data = ev, type = c("p","l"), pch = c(19,19), col = c("red","black"), ylab = ylabtrell, cex = cex, lwd = lwd, main = paste("Empirical and Fitted",ylabtrell), as.table = TRUE, allow.multiple = TRUE, key = list(lines = list( col = c("red","black"), lty = c(1,1), lwd = c(3,3) ), text = list(lab = c("Empirical", "Fit")) ) ) } #------------------- SIMULATE GEOSTATISTICAL DATA geostat.sim <- function(loc.data = spp.syst, xcol = "x", ycol = "y", parsil = 1, range = 1, nugget = 0, minorp = 1, rotate = 90, extrap = NULL, CorModel = "Exponential") { xcoord <- loc.data[,xcol] ycoord <- loc.data[,ycol] n <- length(xcoord) dismat <- dist.geo.ani(xcoord, ycoord, rotate, range, minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CovMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CovMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CovMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CovMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CovMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CovMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CovMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CovMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CovMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CovMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CovMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CovMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CovMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CovMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CovMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CovMat <- CorModel.BesselJ(dismat, extrap) CovMat <- parsil*CovMat + diag(nugget, nrow = n, ncol = n) data.frame(loc.data, z = t(chol(CovMat))%*%rnorm(n)) } # -------------------- COMPLETELY SPATIAL RANDOM LOCATIONS point.sim.csr <- function(npair = 100, lower.x.lim = -1, upper.x.lim = 1, lower.y.lim = -1, upper.y.lim = 1) { x.range <- upper.x.lim - lower.x.lim y.range <- upper.y.lim - lower.y.lim data.frame(x=lower.x.lim + runif(npair)*x.range, y=lower.y.lim + runif(npair)*y.range) } # -------------------- CLUSTER LOCATIONS WITH CSR DATA AND SEEDS point.sim.clus <- function(npair.random = 100, nseed = 0, n.per.seed = 1, seed.range = 1, lower.x.lim = -1, upper.x.lim = 1, lower.y.lim = -1, upper.y.lim = 1) { store.results <- data.frame(x = NULL, y = NULL) x.range <- upper.x.lim - lower.x.lim y.range <- upper.y.lim - lower.y.lim if(npair.random > 0) { random.setx <- lower.x.lim + runif(npair.random)*x.range random.sety <- lower.y.lim + runif(npair.random)*y.range store.results <- rbind(store.results, data.frame(x = random.setx, y = random.sety)) } if(nseed > 0) { random.seedx <- lower.x.lim + seed.range/2 + runif(nseed)*(x.range - seed.range) random.seedy <- lower.y.lim + seed.range/2 + runif(nseed)*(y.range - seed.range) random.seedx <- kronecker(random.seedx,rep(1,times=n.per.seed)) random.seedy <- kronecker(random.seedy,rep(1,times=n.per.seed)) clusx <- random.seedx + (2*runif(nseed*n.per.seed) - 1) * seed.range clusy <- random.seedy + (2*runif(nseed*n.per.seed) - 1) * seed.range store.results <- rbind(store.results, data.frame(cbind(x = clusx, y = clusy))) } store.results } # -------------------- SYSTEMATIC LOCATIONS point.sim.syst <- function(nrow = 10, ncol = 10, lower.x.lim = -1, upper.x.lim = 1, lower.y.lim = -1, upper.y.lim = 1) { x.range <- upper.x.lim - lower.x.lim y.range <- upper.y.lim - lower.y.lim y.mat <- lower.y.lim + y.range* (nrow - matrix(rep(1:nrow, times = ncol), nrow = nrow, ncol = ncol))/ (nrow - 1) x.mat <- lower.x.lim + x.range* (t(matrix(rep(1:ncol, times = nrow), nrow = ncol, ncol = nrow))-1)/ (ncol - 1) data.frame(x = matrix(x.mat, ncol = 1), y = matrix(y.mat, ncol = 1)) } # -------------------- LOCATIONS FROM INHIBITION PROCESS point.sim.inhibit <- function(npair = 100, irange=0.05, lower.x.lim = -1, upper.x.lim = 1, lower.y.lim = -1, upper.y.lim = 1) { spatpts <- data.frame(matrix(0,npair,2)) colnames(spatpts) <- c("x","y") inpts <- 1 x.range <- upper.x.lim - lower.x.lim y.range <- upper.y.lim - lower.y.lim while(inpts <= npair){ok <- 1 xpt <- lower.x.lim + x.range*runif(1) ypt <- lower.y.lim + y.range*runif(1) chk <- inpts-1 ichk <- 1 while(ichk <= chk){ distpt <- sqrt((xpt-spatpts[ichk,1])^2+(ypt-spatpts[ichk,2])^2) if(distpt < irange){ ok <- 0 ichk <- chk} ichk <- ichk+1 } if(ok == 1){spatpts[inpts,1:2] <- cbind(xpt,ypt) inpts <- inpts+1 } } spatpts } # -------------------- SCALED DISTANCE USING ANISOTROPY dist.geo.ani <- function(xcoord, ycoord, rotate, range, minorp) { if(length(xcoord) != length(ycoord)){ return(list(errstate = 1, err.message = "number of x-coord not equal to number of y-coord", err.extra = data.frame(n.xcoord = length(xcoord), n.ycoord = length(ycoord)))) } if(range < 0){ return(list(errstate = 1, err.message = "range parameter less than 0", err.extra = data.frame(range.parm = range))) } if(rotate < 0 || rotate > 180){ return(list(errstate = 1, err.message = "rotation parameter beyond 0 - 180 range", err.extra = data.frame(rotate.parm = rotate))) } if(minorp < 0 || minorp > 1){ return(list(errstate = 1, err.message = "minor range proportion beyond 0 - 1 range", err.extra = dataframe(minorp.parm = minorp))) } # total number of observations n <- length(xcoord) # expand all x-coordinates sx <- matrix(xcoord,ncol = 1) %*% matrix(rep(1,times = n), nrow = 1) # find difference in x-coordinates between all pairwise locations sxdif <- t(sx) - sx # expand all y-coordinates sy <- matrix(ycoord,ncol = 1) %*% matrix(rep(1,times = n), nrow = 1) # find difference in x-coordinates between all pairwise locations sydif <- t(sy) - sy # rotate coordinates newx <- cos(rotate*.0174533)*sxdif - sin(rotate*.0174533)*sydif newy <- sin(rotate*.0174533)*sxdif + cos(rotate*.0174533)*sydif # scale coordinates by minor and major axes */ newx <- newx/(range*minorp) newy <- newy/range # compute distance for the scaled and rotated coordinates */ sqrt(newx^2 + newy^2) } # --------------- SCALED DISTANCE USING ISOTROPY dist.geo.iso <- function(xcoord, ycoord, range) { # total number of observations n <- length(xcoord) distance <- matrix(0, n, n) distance[lower.tri(distance)] <- dist(as.matrix(cbind(xcoord,ycoord))) (distance + t(distance))/range } # --------------- SCALED DISTANCE BETWEEN TWO SETS OF DATA USING ANISTROPY dist.geo.ani.pred <- function(xrow, yrow, xcol, ycol, rotate, range, minorp) { # total number of observations for each set of coordinates n.rows <- length(xrow) n.cols <- length(xcol) # expand all x-coordinates sxr <- matrix(xrow, ncol = 1) %*% matrix(rep(1,times = n.cols), nrow = 1) sxc <- matrix(rep(1,times = n.rows), ncol = 1) %*% matrix(xcol, nrow = 1) syr <- matrix(yrow,ncol = 1) %*% matrix(rep(1,times = n.cols), nrow = 1) syc <- matrix(rep(1,times = n.rows), ncol = 1) %*% matrix(ycol, nrow = 1) # find difference in coordinates between all pairwise locations sxdif <- sxr - sxc sydif <- syr - syc # rotate coordinates newx <- cos(rotate*.0174533)*sxdif - sin(rotate*.0174533)*sydif newy <- sin(rotate*.0174533)*sxdif + cos(rotate*.0174533)*sydif # scale coordinates by minor and major axes */ newx <- newx/(range*minorp) newy <- newy/range # compute distance for the scaled and rotated coordinates */ sqrt(newx^2 + newy^2) } # --------------- A BUNCH OF SPATIAL CORRELATION MODELS CorModel.Exponential <- function(distance.matrix) { exp(-distance.matrix) } CorModel.ExpRadon2 <- function(distance.matrix) { (1 + distance.matrix)*exp(-distance.matrix) } CorModel.ExpRadon4 <- function(distance.matrix) { (1 + distance.matrix + distance.matrix^2/3)*exp(-distance.matrix) } CorModel.Gaussian <- function(distance.matrix) { exp(-distance.matrix^2) } CorModel.Stable <- function(distance.matrix, extrap) { exp(-distance.matrix^extrap) } CorModel.RationalQuad <- function(distance.matrix) { 1/(1+distance.matrix^2) } CorModel.CauchyGrav <- function(distance.matrix) { 1/sqrt(1+distance.matrix^2) } CorModel.CauchyMag <- function(distance.matrix) { 1/(sqrt(1+distance.matrix^2))^3 } CorModel.Cauchy <- function(distance.matrix, extrap) { 1/(1+distance.matrix^2)^extrap } CorModel.Circular <- function(distance.matrix) { d <- distance.matrix d[distance.matrix > 1] <- 0 CovMat <- 2*(acos(d) - d*sqrt(1 - d^2))/pi CovMat[distance.matrix >= 1] <- 0 CovMat } CorModel.Spherical <- function(distance.matrix) { CovMat <- (1 - 1.5*distance.matrix + 0.5*distance.matrix^3) CovMat[distance.matrix > 1] <- 0 CovMat } CorModel.Cubic <- function(distance.matrix) { CovMat <- (1 - 7*distance.matrix^2 + 35*distance.matrix^3/4 - 7*distance.matrix^5/2 + 3*distance.matrix^7/4) CovMat[distance.matrix > 1] <- 0 CovMat } CorModel.Penta <- function(distance.matrix) { CovMat <- (1 - 22*distance.matrix^2/3 + 33*distance.matrix^4 - 77*distance.matrix^5/2 + 33*distance.matrix^7/2 - 11*distance.matrix^9/2 + 5*distance.matrix^11/6) CovMat[distance.matrix > 1] <- 0 CovMat } CorModel.CardinalSine <- function(distance.matrix) { d <- distance.matrix d[distance.matrix == 0] <- 1 CorMat <- sin(d)/d CorMat[distance.matrix == 0] <- 1 CorMat } CorModel.BesselJ <- function(distance.matrix, extrap) { d <- distance.matrix d[distance.matrix == 0] <- 1 extrap <- extrap + 2.0 CorMat <- d^(1-extrap/2)*besselJ(d, extrap/2 - 1)*2^(extrap/2 - 1)*gamma(extrap/2) CorMat[distance.matrix == 0] <- 1 CorMat } CorModel.BesselK <- function(distance.matrix, extrap) { d <- distance.matrix d[distance.matrix == 0] <- 1 CorMat <- d^extrap*besselK(d, extrap)/(2^(extrap - 1)*gamma(extrap)) CorMat[distance.matrix == 0] <- 1 CorMat } #------------ CREATE PREDICTION SET FROM BOUNDARY FILE create.pred.data.from.rshape <- function(rshape.outline, npoints, xcol, ycol) { test.poly <- rshape.outline$Shapes[[1]]$verts if(length(rshape.outline$Shapes[[1]]$Pstart) > 1) test.poly <- test.poly[1:rshape.outline$Shapes[[1]]$Pstart[2],] win.graph() polymap(test.poly, add = FALSE) grid.wi.poly <- gridpts(test.poly,npts = npoints) grid.wi.poly <- as.data.frame(grid.wi.poly) names(grid.wi.poly) <- c(xcol,ycol) pointmap(as.points(grid.wi.poly), add=TRUE) grid.wi.poly } #------------ BLOCK KRIGING USING PREDICTION SET block.krige.Vi <- function(formula, data, data.pred, rshape.outline, resp.var, CorModel = "Exponential", xcol, ycol, covb, beta.hat, resid, Vi, V, theta, case) { if(case == "12"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- 1 rotate <- 90 if(length(theta) == 4) extrap <- theta[4] } if(case == "123"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- theta[4] rotate <- theta[5] if(length(theta) == 6) extrap <- theta[6] } n.pred <- length(data.pred[,1]) pred.out <- data.pred[,c(xcol,ycol)] pred.out[,3] <- rep(-999.9, times = length(data.pred[,1])) # need the response variable in the prediction data set in order to create Xdesign.pred data.pred[,resp.var] <- rep(-999.9, times = length(data.pred[,1])) names(pred.out)[3] <- resp.var pred.out[,4] <- rep(-999.9, times = length(data.pred[,1])) names(pred.out)[4] <- "pred.se" z.pred <- model.frame(formula, data.pred) z.pred <- data.pred[,names(z.pred)[1]] m.pred <- model.frame(formula, data.pred) Xdesign.pred <- model.matrix(formula, m.pred) z <- model.frame(formula, data) z <- matrix(data[,names(z)[1]],ncol = 1) m <- model.frame(formula, data) Xdesign <- model.matrix(formula, m) xcoord <- data[,xcol] ycoord <- data[,ycol] n <- length(xcoord) x.pred <- data.pred[,xcol] y.pred <- data.pred[,ycol] dismat <- dist.geo.ani.pred(xrow = xcoord, yrow = ycoord, xcol = x.pred, ycol = y.pred, rotate = rotate, range = range, minorp = minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") c.pred <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") c.pred <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") c.pred <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") c.pred <- CorModel.Gaussian(dismat) if(CorModel == "Stable") c.pred <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") c.pred <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") c.pred <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") c.pred <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") c.pred <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") c.pred <- CorModel.Circular(dismat) if(CorModel == "Spherical") c.pred <- CorModel.Spherical(dismat) if(CorModel == "Cubic") c.pred <- CorModel.Cubic(dismat) if(CorModel == "Penta") c.pred <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") c.pred <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") c.pred <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") c.pred <- CorModel.BesselJ(dismat, extrap) if(case == "12") c.pred <- parsil*c.pred if(case == "123") c.pred <- parsil*c.pred dismat <- dist.geo.ani.pred(xrow = x.pred, yrow = y.pred, xcol = x.pred, ycol = y.pred, rotate = rotate, range = range, minorp = minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") sig.pp <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") sig.pp <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") sig.pp <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") sig.pp <- CorModel.Gaussian(dismat) if(CorModel == "Stable") sig.pp <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") sig.pp <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") sig.pp <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") sig.pp <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") sig.pp <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") sig.pp <- CorModel.Circular(dismat) if(CorModel == "Spherical") sig.pp <- CorModel.Spherical(dismat) if(CorModel == "Cubic") sig.pp <- CorModel.Cubic(dismat) if(CorModel == "Penta") sig.pp <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") sig.pp <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") sig.pp <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") sig.pp <- CorModel.BesselJ(dismat, extrap) if(case == "12") sig.pp <- parsil*sig.pp + diag(nugget, nrow = n.pred) if(case == "123") sig.pp <- parsil*sig.pp + diag(nugget, nrow = n.pred) sig.AA <- mean(sig.pp) c.A <- matrix(rowMeans(c.pred), ncol = 1) x.A <- matrix(colMeans(Xdesign.pred), ncol = 1) lambda <- matrix(t(c.A) %*% Vi + (t(x.A) - t(c.A) %*% Vi %*% Xdesign) %*% covb %*% t(Xdesign) %*% Vi, ncol = 1) block.krige.est <- t(lambda) %*% z block.krige.se <- sqrt(sig.AA - 2*t(c.A) %*% lambda + t(lambda) %*% V %*% lambda) data.frame(block.prediction.estimate = block.krige.est, block.prediction.stderror = block.krige.se) } #------------ MAKE MAP OF VARIOUS RESIDUALS graph.data.resid <- function(rshape.outline, data, xcol, ycol, sample.size.col, response.col, outfilen, sample.size.fun = eval, response.fun = eval, nclass.response = 10, nclass.sample = 10, cex.start.sample.size = .5, cex.increment.sample.size = 0.3, leg.deci.dig = 2, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red") ) { postscript(outfilen, onefile = TRUE, horizontal = FALSE) layout(matrix(c(1,2), ncol = 2), widths = c(2,1)) par(mar=c(0,0,0,0)) plot(rshape.outline, xaxt = "n", yaxt = "n", fg="white", xlab = "", ylab = "") breaks.j <- size.color.gr(data = data, xcol = xcol, ycol = ycol, sample.size.col = sample.size.col, response.col = response.col, sample.size.fun = invf, response.fun = eval, nclass.response = nclass.response, nclass.sample = nclass.sample, cex.start.sample.size = cex.start.sample.size, cex.increment.sample.size = cex.increment.sample.size, color.palette = color.palette) plot(data.frame(x = 0, y = 0), type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") nleg <- length(breaks.j[,1]) dec.dig <- leg.deci.dig lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("to",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha,dash.j, upperb.cha) legend(x = -1, y = .25, pch = 19, cex = 1.1, legend = leg.lab, col = color.palette) if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ MAKE MAP OF KRIGING PREDICTIONS AND STANDARD ERRORS graph.pred.data.krige <- function(rshape.outline, data, data.pred, xcol, ycol, sample.size.col, response.col, graphic.device = "windows", output.filename, sample.size.fun = eval, response.fun = eval, nclass.response = 10, nclass.sample = 10, cex.start.sample.size = .5, cex.increment.sample.size = 0.3, leg.deci.dig = 2, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red"), ... ) { if(graphic.device == "postscript") postscript(output.filename, onefile = TRUE, horizontal = FALSE) if(graphic.device == "windows") win.graph() layout(matrix(c(1,2), ncol = 2), widths = c(2,1)) par(mar=c(0,0,0,0)) plot(rshape.outline, xaxt = "n", yaxt = "n", fg="white", xlab = "", ylab = "") breaks.j <- size.color.gr(data = data.pred, xcol = xcol, ycol = ycol, sample.size.col = sample.size.col, response.col = response.col, sample.size.fun = sample.size.fun, response.fun = response.fun, nclass.response = nclass.response, nclass.sample = nclass.sample, cex.start.sample.size = cex.start.sample.size, cex.increment.sample.size = cex.increment.sample.size, color.palette = color.palette) points(data[,xcol], data[,ycol], pch = 19, ...) plot(data.frame(x = 0, y = 0), type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") nleg <- length(breaks.j[,1]) dec.dig <- leg.deci.dig lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("to",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha,dash.j, upperb.cha) legend(x = -1, y = .25, pch = 19, cex = 1.1, legend = leg.lab, col = color.palette) if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ GENERALIZED R2 FOR SPATIAL MODELS R2g <- function(z, Xdesign, bhat, Vi) { SS.bhat <- t(z - Xdesign %*% bhat) %*% Vi %*% (z - Xdesign %*% bhat) mu.hat <- sum(Vi %*% z)/sum(Vi) SS.mu <- matrix(z - mu.hat, nrow = 1) %*% Vi %*% matrix(z - mu.hat, ncol = 1) 1 - SS.bhat/SS.mu } #------------ CROSSVALIDATION DIAGNOSTICS crossval <- function(z, X, V, Vi, n, p, covb, covbi, bhat) { cdd.out <- matrix(-999.9, nrow = n, ncol = 2) for(i in 1:n) { Vi.i <- Vi[(1:n) != i,(1:n) != i] - matrix(Vi[(1:n) != i,i],ncol = 1) %*% matrix(Vi[i,(1:n) != i],nrow = 1)/Vi[i,i] c.i <- matrix(V[(1:n) != i,i],ncol = 1) xi <- matrix(X[i,], ncol = 1) X.i <- X[(1:n) != i,] z.i <- matrix(z[(1:n) != i], ncol = 1) xxi <- xi - t(X.i) %*% Vi.i %*% c.i hhi <- t(xxi) %*% covb %*% xxi zzi <- z[i] - t(z.i) %*% Vi.i %*% c.i si <- V[i,i] - t(c.i) %*% Vi.i %*% c.i covb.i <- mginv(t(X.i) %*% Vi.i %*% X.i) lam <- t(c.i + X.i %*% covb.i %*% xxi) %*% Vi.i cdd.out[i,1] <- lam %*% z.i cdd.out[i,2] <- sqrt(si + t(xxi) %*% covb.i %*% xxi) } cdd.out <- as.data.frame(cdd.out) names(cdd.out) <- c("cv.pred","cv.se") cdd.out } #------------ INFLUENCE DIAGNOSTICS infldiagn <- function(z, V, X, n, p) { svd.out <- svd(V) V.5 <- svd.out$u %*% diag(1/sqrt(svd.out$d)) %*% t(svd.out$u) z.5 <- V.5 %*% z X.5 <- V.5 %*% X # Hat matrix, Montgomery and Peck, pg. 170 Hat <- X.5 %*% mginv(t(X.5) %*% X.5) %*% t(X.5) # standardized residuals, Montgomery and Peck, pg. 170 e <- (diag(n) - Hat) %*% z.5 # studentized residuals Montgomery and Peck, pg. 172 r <- e/sqrt(diag((diag(n) - Hat))) # leverage - diagonal hat elements scaled by 2p/n, Montgomery and Peck, pg. 182 hii <- diag(Hat) lever <- hii/(2*p/n) # Cook's D, Montgomery and Peck, pg. 182 Di <- r^2*hii/((1-hii)*p) mat <- as.data.frame(cbind(e,r,lever,Di)) colnames(mat) <- c("resid.stand", "resid.student", "leverage", "CooksD") mat } # ------------ CREATE ONE OF SEVERAL TYPES OF DESIGN MATRICES X.design <- function(formula, data, Xmethod = "treatments") { if(Xmethod == "treat.first.0") { options(contrasts=c("contr.treatment", "contr.treatment")) X <- model.matrix(formula, data = data) options(contrasts=c("contr.treatment", "contr.poly")) } if(Xmethod == "sum.to.0") { options(contrasts=c("contr.sum", "contr.sum")) X <- model.matrix(formula, data = data) options(contrasts=c("contr.treatment", "contr.poly")) } if(Xmethod == "over.par") { terms.list <- attr(terms(formula),"term.labels") X <- NULL if(attr(terms(formula),"intercept") == 1) X <- rep(1, times = length(data[,1])) for (i in 1:length(terms.list)) { form1 <- formula(paste("~ ", terms.list[i], " - 1")) X <- cbind(X,model.matrix(form1, data = data)) } } X } #------------ CREATE A MAP WITH LOCATIONS COLORED ACCORDING TO THEIR VALUE geo.graph.sizecoloroutline <- function(data, rshape.outline = NULL, xcol = xcol, ycol = ycol, color.col, color.fun = eval, size.col, size.fun = eval, cex.start.size = 1.5, cex.increment.size = .001, dec.dig, graphic.device, output.filename) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() layout(matrix(c(1,2), ncol = 2), widths = c(2,1)) par(mar=c(0,0,0,0)) if(!is.null(rshape.outline)) plot(rshape.outline, xaxt = "n", yaxt = "n", fg="white", xlab = "", ylab = "") else plot(data[,c(xcol,ycol)], type = "n") breaks.j <- size.color.gr(data = data, xcol = xcol, ycol = ycol, sample.size.col = size.col, response.col = color.col, sample.size.fun = size.fun, response.fun = color.fun, nclass.response = 10, nclass.sample = 10, cex.start.sample.size = cex.start.size, cex.increment.sample.size = cex.increment.size, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) plot(data.frame(x = 0, y = 0), type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") nleg <- length(breaks.j[,1]) dec.dig <- dec.dig lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("to",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha, dash.j, upperb.cha) legend(x = -1, y = .5, pch = 19, cex = 1.1, legend = leg.lab, col = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ CREATE A (MULTIVARIATE) SCATTERPLOT graph.scatter <- function(data, output.filename, graphic.device, ...) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() plot(data, pch = 19, main = "Scatter Plots", ...) if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ CREATE A HISTOGRAM graph.histogram <- function(data, response.column, output.filename, graphic.device) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() hist(data[,response.column],lwd = 3, col = "blue", nclass = 15, xlim = c(min(data[,response.column]),max(data[,response.column])), main = "Histogram", xlab = "Value Classes") if(graphic.device == "postscript") dev.off(dev.cur()) } # ---------- CREATE A BOXPLOT graph.boxplot <- function(data, response.column, output.filename, graphic.device) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() boxplot(data[,response.column], col = "blue", pch = 19, cex.lab = 1.5, cex = 2, main = "Box Plot") if(graphic.device == "postscript") dev.off(dev.cur()) } # ---------- CREATE A QQPLOT graph.QQplot <- function(data, response.column, output.filename, graphic.device) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() qqnorm(data[,response.column], pch = 19, cex = 1.2, cex.lab = 1.5, main = "QQ Plot") if(graphic.device == "postscript") dev.off(dev.cur()) } # ---------- CREATE EMPIRICAL SEMIVARIOGRAM ROSE graph.semivariogram.rose <- function(data, response.column, xcol, ycol, legend.decimals = 2, maxlag.divisor = 4, EmpVarMeth = "MethMoment", nlagcutoff = 3, graphic.device, output.filename) { maxlag <- max(max(data[,xcol]) - min(data[,xcol]),max(data[,ycol]) - min(data[,ycol]))/maxlag.divisor emp.var.out <- empirical.semivariogram(data, x = xcol, y = ycol, var = response.column, nlag = 12, directions = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170), tolerance = 10, maxlag = maxlag, nlagcutoff = nlagcutoff, EmpVarMeth = EmpVarMeth) emp.var.out1 <- emp.var.out emp.var.out1[,"hx"] <- -emp.var.out1[,"hx"] emp.var.out1[,"hy"] <- -emp.var.out1[,"hy"] emp.var.out <- rbind(emp.var.out,emp.var.out1) if(graphic.device == "postscript") { postscript(output.filename, horizontal = FALSE) par(mar=c(2,4,1,4)) layout(matrix(c(1,2), nrow = 2), heights = c(2,1)) legx <- -.25 cex.inc <- .3} if(graphic.device == "windows") { win.graph() par(mar=c(2,9,2,9)) layout(matrix(c(1,2), nrow = 2), heights = c(1.2,1)) legx <- -.45 cex.inc <- .2} plot(x = emp.var.out[,"hx"], y = emp.var.out[,"hy"], type = "n", main = paste("Empirical Covariogram Rose - ",EmpVarMeth," Method"), xlab = "x", ylab = "y") breaks.j <- size.color.gr(data = emp.var.out, xcol = "hx", ycol = "hy", sample.size.col = "np", response.col = "gamma", cex.start.sample.size = 1, cex.increment.sample.size = cex.inc, nclass.response = 10, nclass.sample = 10, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) nleg <- length(breaks.j[,1]) dec.dig <- legend.decimals lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("-",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha,dash.j, upperb.cha) plot(data.frame(x = 0, y = 0), type = "n", xlim = c(-1,1), ylim = c(-1,1), xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") legend(x = legx, y = 1, pch = 19, cex = 1.1, legend = leg.lab, col = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) if(graphic.device == "postscript") dev.off(dev.cur()) } # ---------- CREATE EMPIRICAL SEMIVARIOGRAM/COVARIANCE ROSE WITH FITTED MODEL graph.covariogram.fitted.rose <- function(data, slm.geo.out, response.column, xcol, ycol, legend.decimals = 2, maxlag.divisor = 4, EmpVarMeth = "MethMoment", nlagcutoff = 3, graphic.device, output.filename) { maxlag <- max(max(data[,xcol]) - min(data[,xcol]),max(data[,ycol]) - min(data[,ycol]))/maxlag.divisor xy <- point.sim.syst(nrow = 41, ncol = 41, lower.x.lim = -maxlag, upper.x.lim = maxlag, lower.y.lim = -maxlag, upper.y.lim = maxlag) xy[,2] <- -xy[,2] emp.var.out <- empirical.semivariogram(data, x = xcol, y = ycol, var = response.column, nlag = 12, directions = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170), tolerance = 10, maxlag = maxlag, nlagcutoff = nlagcutoff, EmpVarMeth = EmpVarMeth) emp.var.out1 <- emp.var.out emp.var.out1[,"hx"] <- -emp.var.out1[,"hx"] emp.var.out1[,"hy"] <- -emp.var.out1[,"hy"] emp.var.out <- rbind(emp.var.out,emp.var.out1) if(graphic.device == "postscript") { postscript(output.filename, horizontal = FALSE) par(mar=c(2,4,1,4)) layout(matrix(c(1,2), nrow = 2), heights = c(2,1)) legx <- -.25 cex.inc <- .3} if(graphic.device == "windows") { win.graph() par(mar=c(2,9,2,9)) layout(matrix(c(1,2), nrow = 2), heights = c(1.2,1)) legx <- -.45 cex.inc <- .2} plot(x = emp.var.out[,"hx"], y = emp.var.out[,"hy"], type = "n", main = paste("Empirical/Fitted",EmpVarMeth,"Rose"), xlab = "x", ylab = "y") breaks.j <- size.color.gr(data = emp.var.out, xcol = "hx", ycol = "hy", sample.size.col = "np", response.col = "gamma", cex.start.sample.size = 1, cex.increment.sample.size = cex.inc, nclass.response = 10, nclass.sample = 10, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) nugget <- 0 if(any(names(slm.geo.out$theta) == "nugget")) nugget <- slm.geo.out$theta[names(slm.geo.out$theta) == "nugget"] parsil <- 0 if(any(names(slm.geo.out$theta) == "parsil")) parsil <- slm.geo.out$theta[names(slm.geo.out$theta) == "parsil"] range <- 0 if(any(names(slm.geo.out$theta) == "range")) range <- slm.geo.out$theta[names(slm.geo.out$theta) == "range"] rotate <- 90 if(any(names(slm.geo.out$theta) == "rotate")) rotate <- slm.geo.out$theta[names(slm.geo.out$theta) == "rotate"] minorp <- 1 if(any(names(slm.geo.out$theta) == "minorp")) minorp <- slm.geo.out$theta[names(slm.geo.out$theta) == "minorp"] extrap <- 0 if(any(names(slm.geo.out$theta) == "extrap")) extrap <- slm.geo.out$theta[names(slm.geo.out$theta) == "extrap"] CorModel <- slm.geo.out$CorModel # rotate coordinates newx <- cos(rotate*.0174533)*xy[,1] - sin(rotate*.0174533)*xy[,2] newy <- sin(rotate*.0174533)*xy[,1] + cos(rotate*.0174533)*xy[,2] # scale coordinates by minor and major axes */ newx <- newx/(range*minorp) newy <- newy/range dismat <- sqrt(newx^2 + newy^2) # compute distance for the scaled and rotated coordinates */ if(CorModel == "Exponential") cm <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") cm <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") cm <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") cm <- CorModel.Gaussian(dismat) if(CorModel == "Stable") cm <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") cm <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") cm <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") cm <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") cm <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") cm <- CorModel.Circular(dismat) if(CorModel == "Spherical") cm <- CorModel.Spherical(dismat) if(CorModel == "Cubic") cm <- CorModel.Cubic(dismat) if(CorModel == "Penta") cm <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") cm <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") cm <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") cm <- CorModel.BesselJ(dismat, extrap) if(EmpVarMeth == "Covariance" | EmpVarMeth == "CovMean") {fit <- parsil*cm ylabtrell <- "Covariance" } else { fit <- nugget + parsil*(1 - cm) ylabtrell <- "Semivariogram" } fit <- matrix(fit, ncol = 41, byrow = TRUE) x <- maxlag*(-20:20)/20 y <- maxlag*(-20:20)/20 contour(x,y,fit, add = TRUE, drawlabels = FALSE) nleg <- length(breaks.j[,1]) dec.dig <- legend.decimals lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("-",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha,dash.j, upperb.cha) plot(data.frame(x = 0, y = 0), type = "n", xlim = c(-1,1), ylim = c(-1,1), xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") legend(x = legx, y = 1, pch = 19, cex = 1.1, legend = leg.lab, col = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) if(graphic.device == "postscript") dev.off(dev.cur()) } # ---------- CREATE SCATTER PLOT OF CROSSVALIDATION PREDICTIONS VERSUS OBSERVED VALUES geo.graph.crossval.pred.vs.obs <- function(resid.data, response.column, output.filename, graphic.device) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() minop <- min(resid.data[,response.column], resid.data[,"cv.pred"]) maxop <- max(resid.data[,response.column], resid.data[,"cv.pred"]) plot(matrix(c(minop,minop,maxop,maxop), nrow = 2, byrow = TRUE), type = "l", xlab = "Measured Value", ylab = "Predicted Value", col = "blue", lwd = 2, main = "Cross-validation Scatter Plot") points(resid.data[,response.column], resid.data[,"cv.pred"], pch = 19, cex = 1.5) if(graphic.device == "postscript") dev.off(dev.cur()) } # ---------- FORMAT FIXED EFFECTS TABLE TO INCLUDE FACTOR LEVELS SET TO 0 fixed.effect.table <- function(formula, data, fixed.effects.estimates) { fee <- fixed.effects.estimates form1 <- formula fee.i <- 0 fee.xtra <- NULL levels.lab <- NULL if(attr(terms(form1),"intercept") == 1) { levels.lab <- data.frame( factor = "intercept", levels = "", fee[1,] ) fee.i <- 1 } terms.table <- attr(terms(form1),"factors") if(length(attr(terms(form1),"term.labels")) > 0) { if(sum(terms.table[,1]) == 1 & attr(terms(form1),"intercept") == 0) terms.table[terms.table[,1] == 1,1] <- 2 for (j in 1:length(terms.table[1,])) { ind.j <- terms.table[,j] > 0 n.interact <- sum(ind.j) lab.i <- NULL for (i in 1:n.interact) { if(!any(rownames(terms.table)[ind.j][i] == colnames(data))) {levels.i <- matrix("") if(i == 1) lab.i <- matrix(levels.i, ncol = 1) if(i > 1) { lab.t <- NULL for (k in 2:i) lab.t <- cbind(lab.t, matrix(rep(lab.i[,k-1], times = length(levels.i)), ncol = 1)) lab.tmp <- matrix(rep(levels.i, each = length(lab.i[,1])), ncol = 1) lab.i <- cbind(lab.t,lab.tmp) } } else if(is.factor(data[,rownames(terms.table)[ind.j][i]])) { levels.i <- levels(data[,rownames(terms.table)[ind.j][i]]) if(i == 1) lab.i <- matrix(levels.i, ncol = 1) if(i > 1) { lab.t <- NULL for (k in 2:i) lab.t <- cbind(lab.t, matrix(rep(lab.i[,k-1], times = length(levels.i)), ncol = 1)) lab.tmp <- matrix(rep(levels.i, each = length(lab.i[,1])), ncol = 1) lab.i <- cbind(lab.t,lab.tmp) } } else {levels.i <- matrix("") if(i == 1) lab.i <- matrix(levels.i, ncol = 1) if(i > 1) { lab.t <- NULL for (k in 2:i) lab.t <- cbind(lab.t, matrix(rep(lab.i[,k-1], times = length(levels.i)), ncol = 1)) lab.tmp <- matrix(rep(levels.i, each = length(lab.i[,1])), ncol = 1) lab.i <- cbind(lab.t,lab.tmp) } } } cum.ind <- rep(TRUE, times = length(lab.i[,1])) levels.lab.i <- NULL for (i in 1:n.interact) { if(!any(rownames(terms.table)[ind.j][i] == colnames(data))) cum.ind <- cum.ind & TRUE else if(is.factor(data[,rownames(terms.table)[ind.j][i]])){ if (terms.table[ind.j,j][i] == 1) { cum.ind <- cum.ind & lab.i[,i] != levels(data[,rownames(terms.table)[ind.j][i]])[1] } } else cum.ind <- cum.ind & TRUE if(i == 1) levels.lab.i <- lab.i[,1] if(i > 1) levels.lab.i <- paste(levels.lab.i, ",", lab.i[,i], sep = "") } n.fee <- sum(cum.ind) n.j <- length(cum.ind) fee.xtra <- matrix(c(0,NA,NA,NA,NA), nrow = n.j, ncol = 5, byrow = TRUE) fee.xtra[cum.ind,1:5] <- matrix(unlist(fee[(fee.i + 1):(fee.i + n.fee), 1:5]), ncol = 5) fee.xtra1 <- data.frame(factor = colnames(terms.table)[j], levels = levels.lab.i, parameters = fee.xtra[,1], std.err = fee.xtra[,2], df = fee.xtra[,3], t.value = fee.xtra[,4], prob.t = fee.xtra[,5]) if(is.null(levels.lab)) levels.lab <- fee.xtra1 else levels.lab <- rbind(levels.lab, fee.xtra1) fee.i <- fee.i + n.fee } } rownames(levels.lab) <- 1:length(levels.lab[,1]) colnames(levels.lab) <- c("effect", "levels", "estimate", "std.err", "df", "t.value", "prob.t") levels.lab } #------------ FUNCTION FOR ESTIMATES AND CONTRASTS, CHECKS FOR ESTIMABILITY estimates.and.contrasts <- function(formula = formula, data = data, est.con.mat = est.con.mat, est.con.lab = est.con.label, fixed.effects.estimates = fixed.effects.estimates, cov.fix.eff = cov.fix.eff, p = p, n = n) { Xop <- X.design(formula = formula, data = data, Xmethod = "over.par") estimability <- matrix(FALSE, nrow = length(est.con.mat[,1]), ncol = 1) for(i in 1:length(est.con.mat[,1])) { p1 <- sum(svd(rbind(Xop, est.con.mat[i,]))$d > 1e-10) estimability[i,] <- p == p1 } estimates.contrasts <- est.con.mat %*% matrix(fixed.effects.estimates[,3], ncol = 1) rownames(estimates.contrasts) <- 1:length(est.con.mat[,1]) estimates.contrasts[!estimability] <- NA ec.var <- est.con.mat %*% cov.fix.eff %*% t(est.con.mat) std.err <- sqrt(diag(ec.var)) std.err[!estimability] <- NA estimates.contrasts <- cbind(estimates.contrasts, std.err) np <- rep(n - p, times = length(est.con.mat[,1])) np[!estimability] <- NA estimates.contrasts <- cbind(estimates.contrasts, np) estimates.contrasts <- cbind(estimates.contrasts, estimates.contrasts[,1]/estimates.contrasts[,2]) estimates.contrasts <- cbind(estimates.contrasts, round(100000*(1- pt(abs(estimates.contrasts[,4]), df = n - p))*2)/100000) estimates.contrasts <- as.data.frame(estimates.contrasts) colnames(estimates.contrasts) <- c("estimate", "std.err", "df", "t.value", "prob.t") estimates.contrasts <- data.frame(label = est.con.lab, estimable = estimability, estimates.contrasts) estimates.contrasts } #------------ RETURN CHARACTER STRING OF CASE OF GEOSTATISTICAL MODEL geo.case <- function(use.spatial, use.nugget, use.anisotropy, use.locrep) { if(use.spatial == FALSE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == FALSE) case <- "2" if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == FALSE) case <- "12" if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == TRUE & use.locrep == FALSE) case <- "123" if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == TRUE) case <- "124" if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == TRUE & use.locrep == TRUE) case <- "1234" case } #------------ ESTIMATE COVARIANCE PARAMETERS est.covparm1 <- function(z, Xdesign, data, xcol = "x", ycol = "y", EstMeth = "REML", CorModel = "Exponential", group = NULL, use.spatial = TRUE, use.locrep = FALSE, use.nugget = TRUE, use.anisotropy = FALSE) { # -- CASE 2 PURE NUGGET if(use.spatial == FALSE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == FALSE) { n <- length(z) p <- sum(svd(Xdesign)$d>1e-10) XX <- t(Xdesign) %*% Xdesign covb <- mginv(XX) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(Xdesign) %*% z nmult <- n - p if(EstMeth == "ML") nmult <- n sigma2 <- as.double(t(z - Xdesign %*% b.hat) %*% (z - Xdesign %*% b.hat) / nmult) theta.opt <- sigma2 f1 <- n*log(sigma2) + nmult f2 <- 0 if(EstMeth == "REML") f2 <- sum(log(Re(eigen(XX/sigma2)$values[1:p]))) m2LL <- f1 + f2 + nmult*log(2*pi) CovMat <- sigma2*diag(n) Vi <- diag(n)/sigma2 names(theta.opt) <- "nugget" # Asycov <- Asy.cov(theta.opt, m2LL.geo12.noprof, z = z, Xdesign = Xdesign, group = group, # xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, # CorModel = CorModel) Asycov <- NULL outpt <- list(theta = theta.opt, m2LL = m2LL, Asycov = Asycov, V = CovMat, Vi = Vi) } # -- CASE 1,2 SPATIAL MODEL + NUGGET if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == FALSE) { if(length(group) > 0) group <- data[,"group"] xcoord <- data[,xcol] ycoord <- data[,ycol] names.theta <- c("parsil", "range", "nugget") if(CorModel == "Stable" || CorModel == "Cauchy" || CorModel == "BesselJ" || CorModel == "BesselK") names.theta <- c(names.theta, "extrap") sc1 <- var(z) sc2 <- 4*sqrt((var(xcoord) + var(ycoord))/2) n <- length(z) p <- sum(svd(Xdesign)$d>1e-10) lnt <- length(names.theta) if(lnt == 3) theta.mat <- cbind( matrix(c(.1,.5,.9) %x% c(1,1,1), ncol = 1), matrix(c(1,1,1) %x% c(.2,1,1.78), ncol = 1) ) if(lnt == 4) theta.mat <- cbind( matrix(c(.1,.89) %x% c(1,1,1,1), ncol = 1), matrix(c(1,1) %x% c(.2,1.78) %x% c(1,1), ncol = 1), matrix(c(1,1,1,1) %x% c(.5, 1.5), ncol = 1) ) best.m2LL <- 1e32 for (i in 1:length(theta.mat[,1])) { temp.m2LL <- m2LL.geo12(theta = theta.mat[i,], z = z, Xdesign = Xdesign, group = group, xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2) if(temp.m2LL < best.m2LL) { theta <- theta.mat[i,] best.m2LL <- temp.m2LL } } parmest <- optim(theta, m2LL.geo12, z = z, Xdesign = Xdesign, xcoord = xcoord, ycoord = ycoord, group = group, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2, hessian = TRUE) phi <- parmest$par[1] range <- parmest$par[2]*sc2 if(length(names.theta) == 4) extrap <- parmest$par[3] if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) CovMat <- matrix(0, nrow = n, ncol = n) Vi <- matrix(0, nrow = n, ncol = n) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.iso(xcoord = xi, ycoord = yi, range = range) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 if(sum(ind0 - ng) > 0) return(list(errstate = 1, errmessage = "Need model with replicated locations")) CovMat.temp <- phi*CorMat + (1 - phi)*diag(ni) CovMat[gind,gind] <- CovMat.temp Vi[gind,gind] <- solve(CovMat.temp) } covb <- mginv(t(Xdesign) %*% Vi %*% Xdesign) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(Xdesign) %*% Vi %*% z nmult <- n - p if(EstMeth == "ML") nmult <- n sigma2 <- as.double(t(z - Xdesign %*% b.hat) %*% Vi %*% (z - Xdesign %*% b.hat) / nmult) CovMat <- sigma2*CovMat Vi <- Vi/sigma2 parsil <- phi*sigma2 nugget <- (1 - phi)*sigma2 theta.opt <- c(parsil, range, nugget) if(length(names.theta) == 4) theta.opt <- c(theta.opt, extrap) names(theta.opt) <- names.theta # Asycov <- Asy.cov(theta.opt, m2LL.geo12.noprof, z = z, Xdesign = Xdesign, group = group, # xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, # CorModel = CorModel) Asycov <- NULL outpt <- list(theta = theta.opt, m2LL = parmest$value, Asycov = Asycov, V = CovMat, Vi = Vi) } # -- CASE 1,2,3 SPATIAL MODEL + NUGGET + ANISOTROPY if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == TRUE & use.locrep == FALSE) { if(length(group) > 0) group <- data[,"group"] xcoord <- data[,xcol] ycoord <- data[,ycol] names.theta <- c("parsil", "range", "nugget", "minorp", "rotate") if(CorModel == "Stable" || CorModel == "Cauchy" || CorModel == "BesselJ" || CorModel == "BesselK") names.theta <- c(names.theta, "extrap") sc1 <- var(z) sc2 <- 4*sqrt((var(xcoord) + var(ycoord))/2) n <- length(z) p <- sum(svd(Xdesign)$d>1e-10) lnt <- length(names.theta) if(lnt == 5) theta.mat <- cbind( matrix(c(.1,.89) %x% rep(1,times = 8), ncol = 1), matrix(c(1,1) %x% c(.2,1.78) %x% c(1,1,1,1), ncol = 1), matrix(c(1,1,1,1) %x% c(.1, .89) %x% c(1,1), ncol = 1), matrix(rep(1, times = 8) %x% c(.1, .89), ncol = 1) ) if(lnt == 6) theta.mat <- cbind( matrix(c(.1,.89) %x% rep(1,times = 16), ncol = 1), matrix(c(1,1) %x% c(.2,1.78) %x% rep(1, times = 8), ncol = 1), matrix(rep(1, times = 4) %x% c(.1, .89) %x% c(1,1,1,1), ncol = 1), matrix(rep(1, times = 8) %x% c(.1, .89) %x% c(1,1), ncol = 1), matrix(rep(1, times = 16) %x% c(.5, 1.5), ncol = 1) ) best.m2LL <- 1e32 for (i in 1:length(theta.mat[,1])) { temp.m2LL <- m2LL.geo123(theta = theta.mat[i,], z = z, Xdesign = Xdesign, group = group, xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2) if(temp.m2LL < best.m2LL) { theta <- theta.mat[i,] best.m2LL <- temp.m2LL } } parmest <- optim(theta, m2LL.geo123, z = z, Xdesign = Xdesign, xcoord = xcoord, ycoord = ycoord, group = group, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2, hessian = TRUE) phi <- parmest$par[1] range <- parmest$par[2]*sc2 minorp <- parmest$par[3] rotate <- parmest$par[4]*180 if(length(names.theta) == 6) extrap <- parmest$par[5] if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) CovMat <- matrix(0, nrow = n, ncol = n) Vi <- matrix(0, nrow = n, ncol = n) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.ani(xcoord = xi, ycoord = yi, range = range, minorp = minorp, rotate = rotate) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 if(sum(ind0 - ng) > 0) return(list(errstate = 1, errmessage = "Need model with replicated locations")) CovMat.temp <- phi*CorMat + (1 - phi)*diag(ni) CovMat[gind,gind] <- CovMat.temp Vi[gind,gind] <- solve(CovMat.temp) } covb <- mginv(t(Xdesign) %*% Vi %*% Xdesign) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(Xdesign) %*% Vi %*% z nmult <- n - p if(EstMeth == "ML") nmult <- n sigma2 <- as.double(t(z - Xdesign %*% b.hat) %*% Vi %*% (z - Xdesign %*% b.hat) / nmult) CovMat <- sigma2*CovMat Vi <- Vi/sigma2 parsil <- phi*sigma2 nugget <- (1 - phi)*sigma2 theta.opt <- c(parsil, range, nugget, minorp, rotate) if(length(names.theta) == 6) theta.opt <- c(theta.opt, extrap) names(theta.opt) <- names.theta # Asycov <- Asy.cov(theta.opt, m2LL.geo123.noprof, z = z, Xdesign = Xdesign, group = group, # xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, # CorModel = CorModel) Asycov <- NULL outpt <- list(theta = theta.opt, m2LL = parmest$value, Asycov = Asycov, V = CovMat, Vi = Vi) } # -- CASE 1,2,4 SPATIAL MODEL + NUGGET + REPLICATED LOCATIONS if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == TRUE) { if(length(group) > 0) group <- data[,"group"] xcoord <- data[,xcol] ycoord <- data[,ycol] names.theta <- c("parsil", "range", "nugget", "locrep") if(CorModel == "Stable" || CorModel == "Cauchy" || CorModel == "BesselJ" || CorModel == "BesselK") names.theta <- c(names.theta, "extrap") sc1 <- var(z) sc2 <- 4*sqrt((var(xcoord) + var(ycoord))/2) n <- length(z) p <- sum(svd(Xdesign)$d>1e-10) lnt <- length(names.theta) if(lnt == 4) theta.mat <- cbind( matrix(c(.1,.89) %x% rep(1,times = 16), ncol = 1), matrix(c(1,1) %x% c(.2,1.78) %x% rep(1, times = 8), ncol = 1), matrix(rep(1, times = 4) %x% c(.1, .89) %x% c(1,1,1,1), ncol = 1) ) if(lnt == 5) theta.mat <- cbind( matrix(c(.1,.89) %x% rep(1,times = 32), ncol = 1), matrix(c(1,1) %x% c(.2,1.78) %x% rep(1, times = 16), ncol = 1), matrix(rep(1, times = 4) %x% c(.1, .89) %x% rep(1, times = 8), ncol = 1), matrix(rep(1, times = 8) %x% c(.5, 1.5), ncol = 1) ) best.m2LL <- 1e32 for (i in 1:length(theta.mat[,1])) { temp.m2LL <- m2LL.geo124(theta = theta.mat[i,], z = z, Xdesign = Xdesign, group = group, xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2) if(temp.m2LL < best.m2LL) { theta <- theta.mat[i,] best.m2LL <- temp.m2LL } } parmest <- optim(theta, m2LL.geo124, z = z, Xdesign = Xdesign, xcoord = xcoord, ycoord = ycoord, group = group, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2, hessian = TRUE) phi1 <- parmest$par[1] range <- parmest$par[2]*sc2 phi2 <- parmest$par[3] if(length(names.theta) == 5) extrap <- parmest$par[4] if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) CovMat <- matrix(0, nrow = n, ncol = n) Vi <- matrix(0, nrow = n, ncol = n) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.iso(xcoord = xi, ycoord = yi, range = range) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 CovMat.temp <- phi1*CorMat + phi2*ind0 + (1 - phi1 - phi2)*diag(ni) CovMat[gind,gind] <- CovMat.temp Vi[gind,gind] <- solve(CovMat.temp) } covb <- mginv(t(Xdesign) %*% Vi %*% Xdesign) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(Xdesign) %*% Vi %*% z nmult <- n - p if(EstMeth == "ML") nmult <- n sigma2 <- as.double(t(z - Xdesign %*% b.hat) %*% Vi %*% (z - Xdesign %*% b.hat) / nmult) CovMat <- sigma2*CovMat Vi <- Vi/sigma2 parsil <- phi1*sigma2 locrep <- phi2*sigma2 nugget <- (1 - phi1 - phi2)*sigma2 theta.opt <- c(parsil, range, nugget, locrep) if(length(names.theta) == 5) theta.opt <- c(theta.opt, extrap) names(theta.opt) <- names.theta # Asycov <- Asy.cov(theta.opt, m2LL.geo123.noprof, z = z, Xdesign = Xdesign, group = group, # xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, # CorModel = CorModel) Asycov <- NULL outpt <- list(theta = theta.opt, m2LL = parmest$value, Asycov = Asycov, V = CovMat, Vi = Vi) } # -- CASE 1,2,3,4 SPATIAL MODEL + NUGGET + ANISOTROPY + REPLICATED LOCATIONS if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == TRUE & use.locrep == TRUE) { if(length(group) > 0) group <- data[,"group"] xcoord <- data[,xcol] ycoord <- data[,ycol] names.theta <- c("parsil", "range", "nugget", "minorp", "rotate", "locrep") if(CorModel == "Stable" || CorModel == "Cauchy" || CorModel == "BesselJ" || CorModel == "BesselK") names.theta <- c(names.theta, "extrap") sc1 <- var(z) sc2 <- 4*sqrt((var(xcoord) + var(ycoord))/2) n <- length(z) p <- sum(svd(Xdesign)$d>1e-10) lnt <- length(names.theta) if(lnt == 6) theta.mat <- cbind( matrix(c(.1,.89) %x% rep(1,times = 16), ncol = 1), matrix(c(1,1) %x% c(.2,1.78) %x% rep(1, times = 8), ncol = 1), matrix(rep(1, times = 4) %x% c(.1, .89) %x% c(1,1,1,1), ncol = 1), matrix(rep(1, times = 8) %x% c(.1, .89) %x% c(1,1), ncol = 1), matrix(rep(1, times = 16) %x% c(.1, .89), ncol = 1) ) if(lnt == 7) theta.mat <- cbind( matrix(c(.1, .89) %x% rep(1,times = 216), ncol = 1), matrix(c(1,1) %x% c(.2, 1, 1.78) %x% rep(1, times = 72), ncol = 1), matrix(rep(1, times = 6) %x% c(.1, .5, .89) %x% rep(1, times = 24), ncol = 1), matrix(rep(1, times = 18) %x% c(.1, .5, .89) %x% rep(1, times = 8), ncol = 1), matrix(rep(1, times = 54) %x% c(.1, .89) %x% rep(1, times = 4), ncol = 1), matrix(rep(1, times = 108) %x% c(.4, .8, 1.2, 1.6), ncol = 1) ) best.m2LL <- 1e32 for (i in 1:length(theta.mat[,1])) { temp.m2LL <- m2LL.geo1234(theta = theta.mat[i,], z = z, Xdesign = Xdesign, group = group, xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2) if(temp.m2LL < best.m2LL) { theta <- theta.mat[i,] best.m2LL <- temp.m2LL } } parmest <- optim(theta, m2LL.geo1234, z = z, Xdesign = Xdesign, xcoord = xcoord, ycoord = ycoord, group = group, EstMeth = EstMeth, CorModel = CorModel, sc1 = sc1, sc2 = sc2, hessian = TRUE) phi1 <- parmest$par[1] range <- parmest$par[2]*sc2 minorp <- parmest$par[3] rotate <- parmest$par[4]*180 phi2 <- parmest$par[5] if(length(names.theta) == 7) extrap <- parmest$par[6] if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) CovMat <- matrix(0, nrow = n, ncol = n) Vi <- matrix(0, nrow = n, ncol = n) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.ani(xcoord = xi, ycoord = yi, range = range, minorp = minorp, rotate = rotate) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 CovMat.temp <- phi1*CorMat + phi2*ind0 + (1 - phi1 - phi2)*diag(ni) CovMat[gind,gind] <- CovMat.temp Vi[gind,gind] <- solve(CovMat.temp) } covb <- mginv(t(Xdesign) %*% Vi %*% Xdesign) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(Xdesign) %*% Vi %*% z nmult <- n - p if(EstMeth == "ML") nmult <- n sigma2 <- as.double(t(z - Xdesign %*% b.hat) %*% Vi %*% (z - Xdesign %*% b.hat) / nmult) CovMat <- sigma2*CovMat Vi <- Vi/sigma2 parsil <- phi1*sigma2 locrep <- phi2*sigma2 nugget <- (1 - phi1 - phi2)*sigma2 theta.opt <- c(parsil, range, nugget, minorp, rotate, locrep) if(length(names.theta) == 7) theta.opt <- c(theta.opt, extrap) names(theta.opt) <- names.theta # Asycov <- Asy.cov(theta.opt, m2LL.geo123.noprof, z = z, Xdesign = Xdesign, group = group, # xcoord = xcoord, ycoord = ycoord, EstMeth = EstMeth, # CorModel = CorModel) Asycov <- NULL outpt <- list(theta = theta.opt, m2LL = parmest$value, Asycov = Asycov, V = CovMat, Vi = Vi) } outpt } # ---- PROFILED MINUS 2 TIME LOG-LIKELIHOOD, # --------- CASE 1,2,4 SPATIAL MODEL + NUGGET + REPLICATED LOCATIONS m2LL.geo124 <- function(theta, z, Xdesign, group = NULL, xcoord, ycoord, EstMeth = "REML", CorModel = "Exponential", sc1, sc2) { phi1 <- theta[1] range <- theta[2]*sc2 phi2 <- theta[3] if(length(theta) == 4) extrap <- theta[4] z <- z/sqrt(sc1) if(phi1 <= 0 || range <= 0 || phi2 < 0 || phi1 + phi2 > 1) return(1e32) if(length(theta) == 6) if(extrap <= 0.00 || extrap > 2.00) return(1e32) n <- length(z) if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) p <- sum(svd(Xdesign)$d>1e-10) nX <- length(Xdesign[1,]) f1 <- 0 XQX <- matrix(0, nrow = ng*nX, ncol = nX) XQz <- matrix(0, nrow = nX, ncol = ng) zQzsum <- 0 XQXsum <- matrix(0, nrow = nX, ncol = nX) XQzsum <- matrix(0, nrow = nX, ncol = 1) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.iso(xcoord = xi, ycoord = yi, range = range) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 CovMat <- phi1*CorMat + phi2*ind0 + (1 - phi1 - phi2)*diag(ni) # inverse of covariance matrix Vinv <- solve(CovMat) # store some temporary value f1 <- f1 + sum(log(Re(eigen(CovMat)$values))) zQzsum <- zQzsum + t(zi) %*% Vinv %*% zi XQXi <- t(Xi) %*% Vinv %*% Xi XQXsum <- XQXsum + XQXi XQX[((i-1)*nX + 1):(i*nX),] <- XQXi XQzi <- t(Xi) %*% Vinv %*% zi XQzsum <- XQzsum + XQzi XQz[,i] <- XQzi } # estimate beta hat b.hat <- mginv(XQXsum) %*% XQzsum t1 <- 0 t2 <- 0 for(i in 1:ng) { t1 <- t1 + t(XQz[,i]) %*% b.hat t2 <- t2 + t(b.hat) %*% XQX[((i-1)*nX + 1):(i*nX),] %*% b.hat } # parts of the (restricted) likelihood f2 <- zQzsum - 2*t1 + t2 f3 <- 0 if(EstMeth == "REML") f3 <- sum(log(Re(eigen(XQXsum)$values[1:p]))) # if maximum likelihood, a couple of the values above need to be changed nmult <- n - p if(EstMeth == "ML") nmult <- n # the final (restricted) profiled likelihood; this is minus 2 time log likelihood return( f1 + nmult*log(f2*sc1) + f3 + nmult*log(2*pi/nmult) + nmult ) } # ---- PROFILED MINUS 2 TIME LOG-LIKELIHOOD, # --------- CASE 1,2,3, 4 SPATIAL MODEL + NUGGET + ANISTROPY + REPLICATED LOCATIONS m2LL.geo1234 <- function(theta, z, Xdesign, group = NULL, xcoord, ycoord, EstMeth = "REML", CorModel = "Exponential", sc1, sc2) { phi1 <- theta[1] range <- theta[2]*sc2 minorp <- theta[3] rotate <- theta[4]*180 phi2 <- theta[5] if(length(theta) == 6) extrap <- theta[6] z <- z/sqrt(sc1) if(phi1 <= 0 || range <= 0 || minorp <= 0 || minorp > 1 || rotate < 0 || rotate > 180 || phi2 <= 0 || phi1 + phi2 > 1) return(1e32) if(length(theta) == 6) if(extrap <= 0.00 || extrap > 2.00) return(1e32) n <- length(z) if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) p <- sum(svd(Xdesign)$d>1e-10) nX <- length(Xdesign[1,]) f1 <- 0 XQX <- matrix(0, nrow = ng*nX, ncol = nX) XQz <- matrix(0, nrow = nX, ncol = ng) zQzsum <- 0 XQXsum <- matrix(0, nrow = nX, ncol = nX) XQzsum <- matrix(0, nrow = nX, ncol = 1) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.ani(xcoord = xi, ycoord = yi, range = range, minorp = minorp, rotate = rotate) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 CovMat <- phi1*CorMat + phi2*ind0 + (1 - phi1 - phi2)*diag(ni) # inverse of covariance matrix Vinv <- solve(CovMat) # store some temporary value f1 <- f1 + sum(log(Re(eigen(CovMat)$values))) zQzsum <- zQzsum + t(zi) %*% Vinv %*% zi XQXi <- t(Xi) %*% Vinv %*% Xi XQXsum <- XQXsum + XQXi XQX[((i-1)*nX + 1):(i*nX),] <- XQXi XQzi <- t(Xi) %*% Vinv %*% zi XQzsum <- XQzsum + XQzi XQz[,i] <- XQzi } # estimate beta hat b.hat <- mginv(XQXsum) %*% XQzsum t1 <- 0 t2 <- 0 for(i in 1:ng) { t1 <- t1 + t(XQz[,i]) %*% b.hat t2 <- t2 + t(b.hat) %*% XQX[((i-1)*nX + 1):(i*nX),] %*% b.hat } # parts of the (restricted) likelihood f2 <- zQzsum - 2*t1 + t2 f3 <- 0 if(EstMeth == "REML") f3 <- sum(log(Re(eigen(XQXsum)$values[1:p]))) # if maximum likelihood, a couple of the values above need to be changed nmult <- n - p if(EstMeth == "ML") nmult <- n # the final (restricted) profiled likelihood; this is minus 2 time log likelihood return( f1 + nmult*log(f2*sc1) + f3 + nmult*log(2*pi/nmult) + nmult ) } # ---- PROFILED MINUS 2 TIME LOG-LIKELIHOOD # ---------- CASE 1,2,3 SPATIAL MODEL + NUGGET + ANISTROPY m2LL.geo123 <- function(theta, z, Xdesign, group = NULL, xcoord, ycoord, EstMeth = "REML", CorModel = "Exponential", sc1, sc2) { phi <- theta[1] range <- theta[2]*sc2 minorp <- theta[3] rotate <- theta[4]*180 if(length(theta) == 5) extrap <- theta[5] z <- z/sqrt(sc1) if(phi <= 0 || phi > 1 || range <= 0 || minorp <= 0 || minorp > 1 || rotate < 0 || rotate > 180) return(1e32) if(length(theta) == 5) if(extrap <= 0.00 || extrap > 2.00) return(1e32) n <- length(z) if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) p <- sum(svd(Xdesign)$d>1e-10) nX <- length(Xdesign[1,]) f1 <- 0 XQX <- matrix(0, nrow = ng*nX, ncol = nX) XQz <- matrix(0, nrow = nX, ncol = ng) zQzsum <- 0 XQXsum <- matrix(0, nrow = nX, ncol = nX) XQzsum <- matrix(0, nrow = nX, ncol = 1) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.ani(xcoord = xi, ycoord = yi, range = range, minorp = minorp, rotate = rotate) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 CovMat <- phi*CorMat + (1 - phi)*diag(ni) # inverse of covariance matrix Vinv <- solve(CovMat) # store some temporary value f1 <- f1 + sum(log(Re(eigen(CovMat)$values))) zQzsum <- zQzsum + t(zi) %*% Vinv %*% zi XQXi <- t(Xi) %*% Vinv %*% Xi XQXsum <- XQXsum + XQXi XQX[((i-1)*nX + 1):(i*nX),] <- XQXi XQzi <- t(Xi) %*% Vinv %*% zi XQzsum <- XQzsum + XQzi XQz[,i] <- XQzi } # estimate beta hat b.hat <- mginv(XQXsum) %*% XQzsum t1 <- 0 t2 <- 0 for(i in 1:ng) { t1 <- t1 + t(XQz[,i]) %*% b.hat t2 <- t2 + t(b.hat) %*% XQX[((i-1)*nX + 1):(i*nX),] %*% b.hat } # parts of the (restricted) likelihood f2 <- zQzsum - 2*t1 + t2 f3 <- 0 if(EstMeth == "REML") f3 <- sum(log(Re(eigen(XQXsum)$values[1:p]))) # if maximum likelihood, a couple of the values above need to be changed nmult <- n - p if(EstMeth == "ML") nmult <- n # the final (restricted) profiled likelihood; this is minus 2 time log likelihood return( f1 + nmult*log(f2*sc1) + f3 + nmult*log(2*pi/nmult) + nmult ) } # ---- PROFILED MINUS 2 TIME LOG-LIKELIHOOD, # --------- CASE 1,2 SPATIAL MODEL + NUGGET m2LL.geo12 <- function(theta, z, Xdesign, group = NULL, xcoord, ycoord, EstMeth = "REML", CorModel = "Exponential", sc1, sc2) { phi <- theta[1] range <- theta[2]*sc2 if(length(theta) == 3) extrap <- theta[3] z <- z/sqrt(sc1) if(phi <= 0 || phi > 1 || range <= 0) return(1e32) if(length(theta) == 3) if(extrap <= 0.00 || extrap > 2.00) return(1e32) n <- length(z) if(length(group) < 1) group <- rep(1, times = n) glevels <- unique(group) ng <- length(glevels) p <- sum(svd(Xdesign)$d>1e-10) nX <- length(Xdesign[1,]) f1 <- 0 XQX <- matrix(0, nrow = ng*nX, ncol = nX) XQz <- matrix(0, nrow = nX, ncol = ng) zQzsum <- 0 XQXsum <- matrix(0, nrow = nX, ncol = nX) XQzsum <- matrix(0, nrow = nX, ncol = 1) for(i in 1:ng) { gind <- group == glevels[i] zi <- z[gind] xi <- xcoord[gind] yi <- ycoord[gind] Xi <- Xdesign[gind,] ni <- length(zi) dismat <- dist.geo.iso(xcoord = xi, ycoord = yi, range = range) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") CorMat <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") CorMat <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") CorMat <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") CorMat <- CorModel.Gaussian(dismat) if(CorModel == "Stable") CorMat <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") CorMat <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") CorMat <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") CorMat <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") CorMat <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") CorMat <- CorModel.Circular(dismat) if(CorModel == "Spherical") CorMat <- CorModel.Spherical(dismat) if(CorModel == "Cubic") CorMat <- CorModel.Cubic(dismat) if(CorModel == "Penta") CorMat <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") CorMat <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") CorMat <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") CorMat <- CorModel.BesselJ(dismat, extrap) # create the full covariance matrix with random effects for location and nugget ind0 <- dismat == 0 CovMat <- phi*CorMat + (1 - phi)*diag(ni) # inverse of covariance matrix Vinv <- solve(CovMat) # store some temporary value f1 <- f1 + sum(log(Re(eigen(CovMat)$values))) zQzsum <- zQzsum + t(zi) %*% Vinv %*% zi XQXi <- t(Xi) %*% Vinv %*% Xi XQXsum <- XQXsum + XQXi XQX[((i-1)*nX + 1):(i*nX),] <- XQXi XQzi <- t(Xi) %*% Vinv %*% zi XQzsum <- XQzsum + XQzi XQz[,i] <- XQzi } # estimate beta hat b.hat <- mginv(XQXsum) %*% XQzsum t1 <- 0 t2 <- 0 for(i in 1:ng) { t1 <- t1 + t(XQz[,i]) %*% b.hat t2 <- t2 + t(b.hat) %*% XQX[((i-1)*nX + 1):(i*nX),] %*% b.hat } # parts of the (restricted) likelihood f2 <- zQzsum - 2*t1 + t2 f3 <- 0 if(EstMeth == "REML") f3 <- sum(log(Re(eigen(XQXsum)$values[1:p]))) # if maximum likelihood, a couple of the values above need to be changed nmult <- n - p if(EstMeth == "ML") nmult <- n # the final (restricted) profiled likelihood; this is minus 2 time log likelihood return( f1 + nmult*log(f2*sc1) + f3 + nmult*log(2*pi/nmult) + nmult ) } #------------ FUNCTION FOR COMPUTING TYPE III GENERAL LINEAR HYPOTHESES hypoth <- function(formula, data, Vi) { termlabs <- attr(terms(formula),"term.labels") if(length(termlabs) == 0) return(NULL) z <- model.frame(formula, data) z <- data[,names(z)[1]] z.col <- as.character(attr(terms(formula, data = data),"variables"))[2] form1 <- termlabs[1] if(length(termlabs) > 1) for(i in 2:length(termlabs)) form1 <- paste(form1, "+", termlabs[i]) formula <- as.formula(paste(z.col, "~", form1)) X <- X.design(formula = formula, data = data, Xmethod = "sum.to.0") terms.table <- attr(terms(formula),"factors") n <- length(z) p <- sum(svd(X)$d>1e-10) covb.sum0 <- mginv(t(X) %*% Vi %*% X) b.sum0 <- covb.sum0 %*% t(X) %*% Vi %*% z nterms <- length(terms.table[1,]) cumcol <- 1 storage <- matrix(0, ncol = 4, nrow = nterms) nfactor <- 0 for(i in 1:nterms) { facts <- rownames(terms.table)[terms.table[,i] > 0] nfacts <- length(facts) nxcol <- 1 for(j in 1:nfacts) if(is.factor(data[,facts[j]])) nxcol <- nxcol*(length(levels(data[,facts[j]])) - 1) K <- matrix(0, nrow = p, ncol = nxcol) K[(cumcol + 1):(cumcol + nxcol), 1:nxcol] <- diag(nxcol) cumcol <- cumcol + nxcol storage[i,1] <- nxcol storage[i,2] <- n - p storage[i,3] <- t(t(K) %*% b.sum0) %*% mginv(t(K) %*% covb.sum0 %*% K) %*% (t(K) %*% b.sum0)/sum(svd(K)$d>1e-10) storage[i,4] <- round(100000*(1 - pf(storage[i,3],storage[i,1], storage[i,2])))/100000 } storage <- as.data.frame(storage) colnames(storage) <- c("Num.df", "Den.df", "F.value", "Prob.F") outpt <- data.frame(Effect = termlabs, storage) outpt } #------------ MAKE PREDICTIONS (KRIGE) USING PREDICTION SET pred.krige.Vi <- function(formula, data.pred, data, resp.var, xcol, ycol, CorModel = "Exponential", covb, Vi, theta, case) { if(case == "12"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- 1 rotate <- 90 if(length(theta) == 4) extrap <- theta[4] } if(case == "123"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- theta[4] rotate <- theta[5] if(length(theta) == 6) extrap <- theta[6] } n.pred <- length(data.pred[,1]) pred.out <- data.pred[,c(xcol,ycol)] pred.out[,3] <- rep(-999.9, times = length(data.pred[,1])) # need the response variable in the prediction data set in order to create Xdesign.pred data.pred[,resp.var] <- rep(-999.9, times = length(data.pred[,1])) names(pred.out)[3] <- resp.var pred.out[,4] <- rep(-999.9, times = length(data.pred[,1])) names(pred.out)[4] <- "pred.se" z.pred <- model.frame(formula, data.pred) z.pred <- data.pred[,names(z.pred)[1]] m.pred <- model.frame(formula, data.pred) Xdesign.pred <- model.matrix(formula, m.pred) z <- model.frame(formula, data) z <- matrix(data[,names(z)[1]],ncol = 1) m <- model.frame(formula, data) Xdesign <- model.matrix(formula, m) xcoord <- data[,xcol] ycoord <- data[,ycol] n <- length(xcoord) x.pred <- data.pred[,xcol] y.pred <- data.pred[,ycol] dismat <- dist.geo.ani.pred(xrow = xcoord, yrow = ycoord, xcol = x.pred, ycol = y.pred, rotate = rotate, range = range, minorp = minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") c.pred <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") c.pred <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") c.pred <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") c.pred <- CorModel.Gaussian(dismat) if(CorModel == "Stable") c.pred <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") c.pred <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") c.pred <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") c.pred <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") c.pred <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") c.pred <- CorModel.Circular(dismat) if(CorModel == "Spherical") c.pred <- CorModel.Spherical(dismat) if(CorModel == "Cubic") c.pred <- CorModel.Cubic(dismat) if(CorModel == "Penta") c.pred <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") c.pred <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") c.pred <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") c.pred <- CorModel.BesselJ(dismat, extrap) if(case == "12") c.pred <- parsil*c.pred if(case == "123") c.pred <- parsil*c.pred for(i in 1:n.pred) { mat1 <- matrix(Xdesign.pred[i,], ncol = 1) - t(Xdesign) %*% Vi %*% c.pred[,i] lam <- t(c.pred[,i] + Xdesign %*% covb %*% mat1) %*% Vi pred.var <- parsil + nugget - t(c.pred[,i]) %*% Vi %*% c.pred[,i] + t(mat1) %*% covb %*% mat1 pred.out[i,3] <- lam %*% z pred.out[i,4] <- sqrt(pred.var) } pred.out } #------------ KRIGING WEIGHTS FOR A GIVEN LOCATION graph.krige.wts <- function(formula, data.pred, data, rshape.outline = NULL, resp.var, xcol, ycol, slm.geo.out, graphic.device = "windows") { covb <- slm.geo.out$cov.fix.eff Vi <- slm.geo.out$Vi theta <- slm.geo.out$theta case <- slm.geo.out$case CorModel <- slm.geo.out$CorModel if(case == "12"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- 1 rotate <- 90 if(length(theta) == 4) extrap <- theta[4] } if(case == "123"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- theta[4] rotate <- theta[5] if(length(theta) == 6) extrap <- theta[6] } pred.out <- data.pred[,c(xcol,ycol)] pred.out[,3] <- -999.9 # need the response variable in the prediction data set in order to create Xdesign.pred data.pred[,resp.var] <- -999.9 names(pred.out)[3] <- resp.var pred.out[,4] <- -999.9 names(pred.out)[4] <- "pred.se" z.pred <- model.frame(formula, data.pred) z.pred <- data.pred[,names(z.pred)[1]] m.pred <- model.frame(formula, data.pred) Xdesign.pred <- model.matrix(formula, m.pred) z <- model.frame(formula, data) z <- matrix(data[,names(z)[1]],ncol = 1) m <- model.frame(formula, data) Xdesign <- model.matrix(formula, m) xcoord <- data[,xcol] ycoord <- data[,ycol] n <- length(xcoord) x.pred <- data.pred[,xcol] y.pred <- data.pred[,ycol] dismat <- dist.geo.ani.pred(xrow = xcoord, yrow = ycoord, xcol = x.pred, ycol = y.pred, rotate = rotate, range = range, minorp = minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") c.pred <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") c.pred <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") c.pred <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") c.pred <- CorModel.Gaussian(dismat) if(CorModel == "Stable") c.pred <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") c.pred <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") c.pred <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") c.pred <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") c.pred <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") c.pred <- CorModel.Circular(dismat) if(CorModel == "Spherical") c.pred <- CorModel.Spherical(dismat) if(CorModel == "Cubic") c.pred <- CorModel.Cubic(dismat) if(CorModel == "Penta") c.pred <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") c.pred <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") c.pred <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") c.pred <- CorModel.BesselJ(dismat, extrap) if(case == "12") c.pred <- parsil*c.pred if(case == "123") c.pred <- parsil*c.pred mat1 <- matrix(Xdesign.pred[1,], ncol = 1) - t(Xdesign) %*% Vi %*% c.pred[,1] lam <- t(c.pred[,1] + Xdesign %*% covb %*% mat1) %*% Vi pred.var <- parsil + nugget - t(c.pred[,1]) %*% Vi %*% c.pred[,1] + t(mat1) %*% covb %*% mat1 pred.out[1,3] <- lam %*% z pred.out[1,4] <- sqrt(pred.var) data[,"abslam"] <- matrix(abs(lam), ncol = 1) if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() layout(matrix(c(1,2), ncol = 2), widths = c(2,1)) par(mar=c(0,0,2,0)) if(!is.null(rshape.outline)) plot(rshape.outline, xaxt = "n", yaxt = "n", fg="white", xlab = "", ylab = "", main = "ABSOLUTE VALUE OF KRIGING WEIGHTS") else plot(data[,c(xcol,ycol)], type = "n", main = "ABSOLUTE VALUE OF KRIGING WEIGHTS") breaks.j <- size.color.gr(data = data, xcol = xcol, ycol = ycol, sample.size.col = "abslam", response.col = "abslam", sample.size.fun = eval, response.fun = eval, nclass.response = 10, nclass.sample = 10, cex.start.sample.size = 2, cex.increment.sample.size = .001, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) points(pred.out[1],pred.out[2], pch = "+", cex = 2) plot(data.frame(x = 0, y = 0), type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") nleg <- length(breaks.j[,1]) dec.dig <- 3 lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("to",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha, dash.j, upperb.cha) legend(x = -1, y = .5, pch = 19, cex = 1.1, legend = leg.lab, col = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) if(graphic.device == "postscript") dev.off(dev.cur()) } #------------------------- FINITE POPULATION BLOCK KRIGING FPBK <- function(formula, data, xcol, ycol, FPBK.col, theta, CorModel = "Exponential", case, V, Vi) { if(case == "2") { parsil <- 0 range <- 1 nugget <- theta[1] minorp <- 1 rotate <- 90 } if(case == "12"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- 1 rotate <- 90 if(length(theta) == 4) extrap <- theta[4] } if(case == "123"){ parsil <- theta[1] range <- theta[2] nugget <- theta[3] minorp <- theta[4] rotate <- theta[5] if(length(theta) == 6) extrap <- theta[6] } response.col <- as.character(attr(terms(formula, data = data),"variables"))[2] ind.sa <- !is.na(data[, response.col]) ind.un <- is.na(data[, response.col]) data.sa <- data[ind.sa,] if(is.factor(data.sa[,response.col])) data.sa[,response.col] <- as.numeric(as.character(data.sa[,response.col])) if(is.character(data.sa[,response.col])) data.sa[,response.col] <- as.numeric(data.sa[,response.col]) data.un <- data[ind.un,] data.un[,response.col] <- 1 m.un <- model.frame(formula, data.un) Xu <- model.matrix(formula, m.un) z.sa <- model.frame(formula, data.sa) z.sa <- matrix(data.sa[,names(z.sa)[1]],ncol = 1) m.sa <- model.frame(formula, data.sa) Xs <- model.matrix(formula, m.sa) x.sa <- data.sa[,xcol] y.sa <- data.sa[,ycol] n.sa <- sum(ind.sa) x.un <- data.un[,xcol] y.un <- data.un[,ycol] n.un <- sum(ind.un) B <- data[, FPBK.col] B[is.na(B)] <- 0 Bs <- B[ind.sa] Bu <- B[ind.un] dismat <- dist.geo.ani.pred(xrow = x.sa, yrow = y.sa, xcol = x.un, ycol = y.un, rotate = rotate, range = range, minorp = minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") c.pred <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") c.pred <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") c.pred <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") c.pred <- CorModel.Gaussian(dismat) if(CorModel == "Stable") c.pred <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") c.pred <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") c.pred <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") c.pred <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") c.pred <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") c.pred <- CorModel.Circular(dismat) if(CorModel == "Spherical") c.pred <- CorModel.Spherical(dismat) if(CorModel == "Cubic") c.pred <- CorModel.Cubic(dismat) if(CorModel == "Penta") c.pred <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") c.pred <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") c.pred <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") c.pred <- CorModel.BesselJ(dismat, extrap) if(case == "2") SU <- 0*c.pred if(case == "12") SU <- parsil*c.pred if(case == "123") SU <- parsil*c.pred dismat <- dist.geo.ani.pred(xrow = x.un, yrow = y.un, xcol = x.un, ycol = y.un, rotate = rotate, range = range, minorp = minorp) # compute correlation matrix for scaled distance matrix if(CorModel == "Exponential") c.pred <- CorModel.Exponential(dismat) if(CorModel == "ExpRadon2") c.pred <- CorModel.ExpRadon2(dismat) if(CorModel == "ExpRadon4") c.pred <- CorModel.ExpRadon4(dismat) if(CorModel == "Gaussian") c.pred <- CorModel.Gaussian(dismat) if(CorModel == "Stable") c.pred <- CorModel.Stable(dismat, extrap) if(CorModel == "RationalQuad") c.pred <- CorModel.RationalQuad(dismat) if(CorModel == "CauchyGrav") c.pred <- CorModel.CauchyGrav(dismat) if(CorModel == "CauchyMag") c.pred <- CorModel.CauchyMag(dismat) if(CorModel == "Cauchy") c.pred <- CorModel.Cauchy(dismat, extrap) if(CorModel == "Circular") c.pred <- CorModel.Circular(dismat) if(CorModel == "Spherical") c.pred <- CorModel.Spherical(dismat) if(CorModel == "Cubic") c.pred <- CorModel.Cubic(dismat) if(CorModel == "Penta") c.pred <- CorModel.Penta(dismat) if(CorModel == "CardinalSine") c.pred <- CorModel.CardinalSine(dismat) if(CorModel == "BesselK") c.pred <- CorModel.BesselK(dismat, extrap) if(CorModel == "BesselJ") c.pred <- CorModel.BesselJ(dismat, extrap) if(case == "2") UU <- nugget*diag(n.un) if(case == "12") UU <- nugget*diag(n.un) + parsil*c.pred if(case == "123") UU <- nugget*diag(n.un) + parsil*c.pred # predictions part1 <- Xs %*% mginv(t(Xs) %*% Vi %*% Xs) part2 <- t(Xu) - t(Xs) %*% Vi %*% SU D <- SU + part1 %*% part2 FF <- Vi %*% D Ao <- Bs + FF %*% Bu y.est <- t(Ao) %*% z.sa # variance part1 <- t(FF) %*% V %*% FF part2 <- t(FF) %*% SU y.var <- t(Bu) %*% ( part1 - part2 - t(part2) + UU ) %*% Bu y.se <- sqrt(y.var) data.frame(FPBK.est = y.est, FPBK.se = y.se) } # -------------- A FUNCTION TO DEFINE NEIGHBORS BASED ON A SHARED BORDER ndef.commonborder <- function(rshape) { npolys <- length(rshape$att.data[,1]) neigh.list <- matrix(NA, nrow = npolys, ncol = 100) n.list <- matrix(0, nrow = npolys, ncol = 1) # compare each line segment for all polygons to see if any of them overlap # if the returned value is 4, they share a line segment # note that the segments for one polygon may be going clockwise, while the other # may be going counterclockwise, so we need to switch the vertices on the second # polygon when checking for (i in 1:(npolys-1)) { ni <- 0 for (j in (i+1):npolys) { tot <- 0 for (k in (1:(length(rshape$Shapes[[i]]$verts[,1])-1))) { for (l in 1:length(rshape$Shapes[[j]]$verts[,1])) { if( l < length(rshape$Shapes[[j]]$verts[,1]) ) tot <- tot + sum(sum(rshape$Shapes[[i]]$verts[k:(k+1),] == rshape$Shapes[[j]]$verts[l:(l+1),]) == 4) if( l > 1 ) tot <- tot + sum(sum(rshape$Shapes[[i]]$verts[k:(k+1),] == rshape$Shapes[[j]]$verts[l:(l-1),]) == 4) } } if(tot > 0) { n.list[i,1] <- n.list[i,1] + 1 n.list[j,1] <- n.list[j,1] + 1 neigh.list[i,n.list[i,1]] <- j neigh.list[j,n.list[j,1]] <- i } } } for (i in 1:npolys) { if(i == 1) adj <- matrix(neigh.list[i,1:n.list[i]],ncol = 1) else adj <- rbind(adj, matrix(neigh.list[i,1:n.list[i]], ncol = 1) ) } num <- as.vector(n.list) adj <- as.vector(adj) list(adj = adj, num = num) } # -------------- A FUNCTION TO COMPUTE DISTANCE MATRIX dist.mat <- function(xcoord, ycoord) { # total number of observations n <- length(xcoord) distance <- matrix(0, n, n) distance[lower.tri(distance)] <- dist(as.matrix(cbind(xcoord,ycoord))) (distance + t(distance)) } # -------------- A FUNCTION TO DEFINE NEIGHBORS BASED ON A DISTANCE THRESHOLD ndef.disthresh <- function(data, xcol, ycol, threshold) { thresh <- threshold distmat <- dist.mat(data[,xcol], data[,ycol]) dist.ind <- (distmat < thresh)*1.0 diag(dist.ind) <- 0 n <- length(distmat[,1]) adj <- NULL num <- NULL for(i in 1:n) { neigh <- (1:n)*dist.ind[i,] neigh <- neigh[neigh > 0] adj <- rbind(adj,matrix(neigh, ncol = 1)) num <- rbind(num, length(neigh)) } list(adj = adj, num = num) } # ---------- A FUNCTION TO LINK NEIGHBORS WITH LINES ON A GRAPH lattice.link.gr <- function(data, neighbor.ind, num.neighbor, xcol = "x", ycol = "y", link.lwd = 3) { k <- 1 for (i in 1:length(num.neighbor)) { if(num.neighbor[i] > 0) { for (j in 1:num.neighbor[i]){ line.tmp <- rbind(data[i,c(xcol, ycol)], data[neighbor.ind[k],c(xcol, ycol)]) lines(line.tmp, lwd = link.lwd) k <- k + 1 } } } } #------------------- BASIC SIZE-COLOR GRAPH size.color.gr <- function(data, xcol = "x", ycol = "y", sample.size.col, response.col, sample.size.fun = eval, response.fun = eval, cex.start.sample.size = 1, cex.increment.sample.size = 0, nclass.response = 10, nclass.sample = 10, color.palette = NULL, ...) { lower.breaks <- matrix(0, nrow = nclass.response, ncol = 1) upper.breaks <- matrix(0, nrow = nclass.response, ncol = 1) if(length(color.palette) == 0) color.palette <- rainbow(nclass.response, start = .66, end = .99) col.stepi <- (max(sample.size.fun(data[,sample.size.col]),na.rm = TRUE) - min(sample.size.fun(data[,sample.size.col]),na.rm = TRUE))/nclass.sample + 0.0001/(nclass.sample - 1) imax <- min(sample.size.fun(data[,sample.size.col]),na.rm = TRUE) - 0.0001 col.stepj <- (max(response.fun(data[,response.col]),na.rm = TRUE) - min(response.fun(data[,response.col]),na.rm = TRUE))/nclass.response + 0.0001/(nclass.response - 1) jmax1 <- min(response.fun(data[,response.col]),na.rm = TRUE) - 0.0001 for (i in 1:nclass.sample){ imin <- imax imax <- imax + col.stepi indi <- sample.size.fun(data[,sample.size.col]) >= imin & sample.size.fun(data[,sample.size.col]) <= imax for (j in 1:nclass.response){ if(j == 1) jmax <- jmax1 jmin <- jmax lower.breaks[j] <- jmin jmax <- jmax + col.stepj upper.breaks[j] <- jmax indj <- response.fun(data[,response.col]) >= jmin & response.fun(data[,response.col]) <= jmax points(data[indi & indj,xcol],data[indi & indj,ycol], col = color.palette[j], pch = 19, cex = cex.start.sample.size + (i-1)*cex.increment.sample.size) } } data.frame(lower.breaks = lower.breaks, upper.breaks = upper.breaks) } # ---------- A SIZE-COLOR-LINK GRAPH FOR USE WITH LATTICE MODELS lat.graph.sizecolorlink <- function(data, linkdata = NULL, rshape.poly = NULL, xcol = xcol, ycol = ycol, neighbor.ind, num.neighbor, color.col, color.fun = eval, size.col, size.fun = eval, cex.start.size = 1.5, cex.increment.size = .001, dec.dig, graphic.device, output.filename) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() layout(matrix(c(1,2), ncol = 2), widths = c(2,1)) par(mar=c(0,0,0,0)) if(!is.null(rshape.poly)) { rpolys <- Map2poly(rshape.poly) plot(rpolys) } else plot(data[,c(xcol,ycol)], type = "n") if(!is.null(linkdata)) ldata <- linkdata else ldata <- data lattice.link.gr(data = ldata, neighbor.ind = neighbor.ind, num.neighbor = num.neighbor, xcol = xcol, ycol = ycol, link.lwd = 3) breaks.j <- size.color.gr(data = data, xcol = xcol, ycol = ycol, sample.size.col = size.col, response.col = color.col, sample.size.fun = size.fun, response.fun = color.fun, nclass.response = 10, nclass.sample = 10, cex.start.sample.size = cex.start.size, cex.increment.sample.size = cex.increment.size, color.palette = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) plot(data.frame(x = 0, y = 0), type = "n", xaxt = "n", yaxt = "n", bty = "n", xlab = "", ylab = "") nleg <- length(breaks.j[,1]) lowerb.cha <- as.character(as.numeric(as.integer(breaks.j[,"lower.breaks"]*10^dec.dig))/10^dec.dig) dash.j <- rep("to",times = nleg) upperb.cha <- as.character(as.numeric(as.integer(breaks.j[,"upper.breaks"]*10^dec.dig))/10^dec.dig) leg.lab <- paste(lowerb.cha, dash.j, upperb.cha) legend(x = -1, y = .5, pch = 19, cex = 1.1, legend = leg.lab, col = c("darkblue", "blue", "cyan3", "cyan", "lightgreen", "greenyellow", "yellow", "orange", "tomato", "red")) if(graphic.device == "postscript") dev.off(dev.cur()) } # ------------ CENTROIDS OF POLYGONS IN LATTICE GRAPH shape.prj.coords <- function(rshape){ tmp <- matrix(0, nrow = length(rshape$att.data[,1]), ncol = 2) for (i in 1:length(rshape$att.data[,1])){ tmp[i,] <- rshape$Shape[[i]]$verts } data.frame(Centroid.lon = tmp[,1], Centroid.lat = tmp[,2]) } # ------------ FREEMAN TUKEY TRANSFORMATION Free.Tukey <- function(data, numer, denom) { sqrt(1000*data[,numer]/data[,denom]) + sqrt(1000*(data[,numer]+1)/data[,denom]) } # ------------ CREATES A NEIGHBORHOOD MATRIX FROM ADJACENCY/NUMBER VECTORS Neighmat <- function(adj, num, n) { N.mat <- matrix( 0, nrow = n, ncol = n ) k <- 1 for (i in 1:n){ if(num[i] > 0) { for (j in 1:num[i]){ N.mat[i,adj[k]] <- 1 k <- k + 1 } } } N.mat } # ------------ PERFORMS ROW STANDARDIZATION ON A NEIGHBORHOOD MATRIX row.stand.fun <- function(Wmat) { ndata <- length(Wmat[,1]) onevec <- matrix(1, nrow = ndata, ncol = 1) ci.star <- rowSums(Wmat) ci.star[ci.star == 0] <-1 ci.star <- 1/ci.star Dvec <- ci.star Wmat <- Wmat * (ci.star %*% t(onevec)) list(Wmat = Wmat, Dvec = Dvec) } # ------------ ESTIMATE SIGMA AFTER ESTIMATING RHO sig2hat <- function(rho, z = z, X = X, n, p, ind, CorModel = "CAR", EstMeth = "REML", Wmat, Dvec) { nind <- length(ind) if (CorModel == "CAR"){ IrhoC <- as.matrix(diag(nind) - as.numeric(rho)*Wmat) Qi <- diag(as.vector(1/Dvec)) %*% IrhoC Q <- solve(Qi) } if (CorModel == "SAR"){ IrhoC <- as.matrix(diag(nind) - as.numeric(rho)*Wmat) Qi <- IrhoC %*% diag(as.vector(1/Dvec)) %*% t(IrhoC) Q <- solve(Qi) } if (CorModel == "MA"){ IrhoC1 <- as.matrix(diag(nind) + as.numeric(rho)*Wmat) Q <- IrhoC1 %*% diag(as.vector(Dvec)) %*% t(IrhoC1) Qi <- solve(Q) } if(nind > n) {Q <- Q[ind,ind] Qi <- solve(Q) } r <- z - X %*% mginv(t(X) %*% Qi %*% X) %*% t(X) %*% Qi %*% z if (EstMeth == "ML") return (t(r) %*% Qi %*% r/n) t(r) %*% Qi %*% r/(n - p) } # ------------ MINUS 2*LOG LIKELIHOOD FOR LATTICE MODELS WITH SIGMA PROFILED OUT m2LL.lat <- function(rho, z = z, X = X, n, p, ind, CorModel = "CAR", EstMeth = "REML", Wmat = Wmat, Dvec = Dvec) { nind <- length(ind) if (CorModel == "CAR"){ IrhoC <- as.matrix(diag(nind) - as.numeric(rho)*Wmat) Qi <- diag(as.vector(1/Dvec)) %*% IrhoC Q <- solve(Qi) } if (CorModel == "SAR"){ IrhoC <- as.matrix(diag(nind) - as.numeric(rho)*Wmat) Qi <- IrhoC %*% diag(as.vector(1/Dvec)) %*% t(IrhoC) Q <- solve(Qi) } if (CorModel == "MA"){ IrhoC1 <- as.matrix(diag(nind) + as.numeric(rho)*Wmat) Q <- IrhoC1 %*% diag(as.vector(Dvec)) %*% t(IrhoC1) Qi <- solve(Q) } if(nind > n) {Q <- Q[ind,ind] Qi <- solve(Q) } xpqix <- t(X) %*% Qi %*% X r <- z - X %*% mginv(xpqix) %*% t(X) %*% Qi %*% z rQir <- t(r) %*% Qi %*% r lndv <- sum(log(Re(eigen(Q)$values))) lndx <- sum(log(Re(eigen(xpqix)$values))) # if maximum likelihood, a couple of the values above need to be changed nmult <- n - p if(EstMeth == "ML") {nmult <- n lndx <- 0 } nmult*log(rQir) + lndv + lndx + nmult*log(2*pi/nmult) + nmult } #------------ MAKE PREDICTIONS (KRIGE) FOR ALL MISSING DATA pred.krige.lat <- function(z, X, X.pred, ind, V.all, Vi, covb) { n.pred <- sum(!ind) pred.out <- matrix(NA, nrow = n.pred, ncol = 2) pred.out[,2] <- rep(-999.9, times = n.pred) c.pred <- as.matrix(V.all[ind,!ind]) UU <- as.matrix(V.all[!ind,!ind]) for(i in 1:n.pred) { mat1 <- matrix(X.pred[i,], ncol = 1) - t(X) %*% Vi %*% c.pred[,i] lam <- t(c.pred[,i] + X %*% covb %*% mat1) %*% Vi pred.var <- UU[i,i] - t(c.pred[,i]) %*% Vi %*% c.pred[,i] + t(mat1) %*% covb %*% mat1 pred.out[i,1] <- lam %*% z pred.out[i,2] <- sqrt(pred.var) } pred.out } # ---------- CREATE SCATTER PLOT OF CROSSVALIDATION PREDICTIONS VERSUS OBSERVED VALUES lat.graph.crossval.pred.vs.obs <- function(data, obs.column, pred.col, output.filename = NULL, graphic.device = NULL) { if(graphic.device == "postscript") postscript(output.filename) if(graphic.device == "windows") win.graph() minop <- min(data[,obs.column], data[,pred.col]) maxop <- max(data[,obs.column], data[,pred.col]) plot(matrix(c(minop,minop,maxop,maxop), nrow = 2, byrow = TRUE), type = "l", xlab = "Measured Value", ylab = "Predicted Value", col = "blue", lwd = 2, main = "Cross-validation Scatter Plot") points(data[,obs.column], data[,pred.col], pch = 19, cex = 1.5) if(graphic.device == "postscript") dev.off(dev.cur()) } # ------------ ESTIMATE THE COVARIANCE PARAMETERS FOR LATTICE MODEL est.covparm.lat <- function(z, X, data, data.all, ind, n, p, ndef.adj, ndef.num, CorModel = "CAR", EstMeth = "REML", use.row.standard = TRUE, use.distance = FALSE, weight.col = NULL, xcol = NULL, ycol = NULL, alpha = 0) { nind <- length(ind) Wmat <- Neighmat(ndef.adj, ndef.num, length(ind)) Dvec <- matrix(rep(1, times = nind), ncol = 1) if(use.distance == TRUE) { dismat <- dist.mat(data.all[,xcol], data.all[,ycol]) D <- 1/dismat diag(D) <- 0 D <- D^alpha Wmat <- Wmat*D maxD <- max(Wmat) Wmat <- Wmat/maxD Dvec <- rep(1,times = length(ind)) } if(use.row.standard == TRUE) { temp <- row.stand.fun(Wmat) Wmat <- temp$Wmat Dvec <- temp$Dvec } if(!is.null(weight.col)) { ni <- matrix(data.all[,weight.col], ncol = 1) onevec <- matrix(rep(1, times = length(ni)), ncol = 1) Wmat <- Wmat * sqrt((onevec %*% t(ni))/(ni %*% t(onevec))) Dvec <- Dvec/ni } evals <- Re(eigen(Wmat)$values) rhomin <- 1/min(evals) rhomax <- 1/max(evals) opt.out <- optimize(m2LL.lat, z = z, X = X, n = n, p = p, ind = ind, CorModel = CorModel, EstMeth = EstMeth, Wmat = Wmat, Dvec = Dvec, lower = rhomin*.999999999, upper = rhomax*.9999999999) theta <- data.frame(rho = opt.out$minimum) m2LL <- opt.out$objective rho.restrict <- data.frame(rhomin = rhomin, rhomax = rhomax) sig2est <- as.double(sig2hat(rho = theta["rho"], z = z, X = X, n = n, p = p, ind, CorModel = CorModel, EstMeth = EstMeth, Wmat = Wmat, Dvec = Dvec)) theta <- data.frame(theta, sigma = sig2est) if (CorModel == "CAR"){ IrhoC <- as.matrix(diag(nind) - as.numeric(theta["rho"])*Wmat) Vi.all <- (diag(as.vector(1/Dvec)) %*% IrhoC)/sig2est V.all <- solve(Vi.all) } if (CorModel == "SAR"){ IrhoC <- as.matrix(diag(nind) - as.numeric(theta["rho"])*Wmat) Vi.all <- (IrhoC %*% diag(as.vector(1/Dvec)) %*% t(IrhoC))/sig2est V.all <- solve(Vi.all) } if (CorModel == "MA"){ IrhoC1 <- as.matrix(diag(nind) + as.numeric(theta["rho"])*Wmat) V.all <- (IrhoC1 %*% diag(as.vector(Dvec)) %*% t(IrhoC1))*sig2est Vi.all <- solve(V.all) } Vi <- Vi.all V <- V.all if(sum(!ind) > 0) { V <- V[ind,ind] Vi <- solve(V) } list(theta = theta, rho.restrict = rho.restrict, m2LL = m2LL, V = V, Vi = Vi, V.all = V.all, Vi.all = Vi.all) } #------------ FINITE POPULATION BLOCK KRIGING FOR LATTICEL DATA FPBK.lat <- function(z.sa, Xs, Xu, ind.sa, data, FPBK.col, V.all, Vi) { ind.un <- !ind.sa B <- data[, FPBK.col] B[is.na(B)] <- 0 Bs <- B[ind.sa] Bu <- B[ind.un] SU <- as.matrix(V.all[ind.sa,ind.un]) UU <- as.matrix(V.all[ind.un,ind.un]) V <- as.matrix(V.all[ind.sa,ind.sa]) # predictions part1 <- Xs %*% mginv(t(Xs) %*% Vi %*% Xs) part2 <- t(Xu) - t(Xs) %*% Vi %*% SU D <- SU + part1 %*% part2 FF <- Vi %*% D Ao <- Bs + FF %*% Bu y.est <- t(Ao) %*% z.sa # variance part1 <- t(FF) %*% V %*% FF part2 <- t(FF) %*% SU y.var <- t(Bu) %*% ( part1 - part2 - t(part2) + UU ) %*% Bu y.se <- sqrt(y.var) data.frame(FPBK.est = y.est, FPBK.se = y.se) } #------------ JOIN-COUNT STATISTIC WITH RANDOMIZATION TEST BW.join <- function(datavec.01, ndef.adj, ndef.num, nsim = 999) { ind <- !is.na(datavec.01) datavec.01 <- datavec.01[ind] n <- length(ind) Nmat <- Neighmat(ndef.adj, ndef.num, n) Nmat <- Nmat[ind,ind] BW.obs <- 0 n <- length(datavec.01) for (i in 1:n){ BW.obs <- BW.obs + sum(datavec.01[i] != datavec.01[Nmat[i,] == 1]) } storage.join <- matrix(NA, nrow = nsim, ncol = 1) for (j in 1:nsim) { rand.ord <- cbind(datavec.01,runif(n)) rand.ord <- rand.ord[order(rand.ord[,2]),][,1] BW <- 0 for (i in 1:n){ BW <- BW + sum(rand.ord[i] != rand.ord[Nmat[i,] == 1]) } storage.join[j] <- BW } storage.join <- storage.join[order(storage.join[,1]),] P.pos.join <- sum(BW.obs < storage.join)/nsim P.neg.join <- 1 - P.pos.join list(BW.obs = BW.obs, Pvalue.neg = P.neg.join, Pvalue.pos = P.pos.join, BW.null = storage.join) } #------------ GRAPH FOR JOIN-COUNT STATISTIC WITH RANDOMIZATION TEST lat.graph.BWjoin <- function(BW.join.out, graphic.device = "windows", output.filename) { if(graphic.device == "postscript") postscript(output.filename, onefile = TRUE, horizontal = FALSE) if(graphic.device == "windows") win.graph() hist(BW.join.out$BW.null, xlim= c(min(BW.join.out$BW.null,BW.join.out$BW.obs), max(BW.join.out$BW.null,BW.join.out$BW.obs)), nclass = 20, col = "blue", xlab = "BW Joins", cex.axis = 1.3, cex.lab = 1.5, main = "JOIN-COUNT STATISTIC") points(BW.join.out$BW.obs,1, pch = 19, cex = 2, col = "darkgreen") if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ GEARY'S C STATISTIC WITH RANDOMIZATION TEST Geary.c <- function(datavec, ndef.adj, ndef.num, nsim = 100) { ind <- !is.na(datavec) datavec <- datavec[ind] n <- length(ind) Wmat <- Neighmat(ndef.adj, ndef.num, n) Wmat <- Wmat[ind,ind] n <- length(datavec) zzt <- matrix(datavec ,nrow= length(datavec), ncol = length(datavec)) denom <- 2*sum(Wmat)*var(datavec) geary.obs <- 1 - sum((zzt - t(zzt))^2*Wmat)/denom storage.geary <- matrix(NA, nrow = nsim, ncol = 1) for (j in 1:nsim) { rand.ord <- cbind(datavec,runif(n)) rand.ord <- rand.ord[order(rand.ord[,2]),][,1] zzt <- matrix(rand.ord ,nrow= length(datavec), ncol = length(datavec)) storage.geary[j] <- 1- sum((zzt - t(zzt))^2*Wmat)/denom } storage.geary <- storage.geary[order(storage.geary[,1]),] P.pos.geary <- sum(geary.obs < storage.geary)/nsim P.neg.geary <- 1 - P.pos.geary list(geary.obs = geary.obs, Pvalue.neg = P.neg.geary, Pvalue.pos = P.pos.geary, geary.null = storage.geary) } #------------ GRAPH FOR GEARY'S C WITH RANDOMIZATION TEST lat.graph.Geary.c <- function(Geary.c.out, graphic.device = "windows", output.filename) { if(graphic.device == "postscript") postscript(output.filename, onefile = TRUE, horizontal = FALSE) if(graphic.device == "windows") win.graph() hist(Geary.c.out$geary.null, xlim= c(min(Geary.c.out$geary.null,Geary.c.out$geary.obs), max(Geary.c.out$geary.null,Geary.c.out$geary.obs)), nclass = 20, col = "blue", xlab = "Geary's c", cex.axis = 1.3, cex.lab = 1.5, main = "GEARY'S C") points(Geary.c.out$geary.obs,1, pch = 19, cex = 2, col = "darkgreen") if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ MORAN'S I STATISTIC WITH RANDOMIZATION TEST Moran.I <- function(datavec, ndef.adj, ndef.num, nsim = 100){ ind <- !is.na(datavec) datavec <- datavec[ind] Wmat <- Neighmat(ndef.adj, ndef.num, length(ind)) Wmat <- Wmat[ind,ind] x.bar <- mean(datavec) zzt <- matrix(datavec - x.bar ,nrow= length(datavec), ncol = length(datavec)) denom <- sum(Wmat)*var(datavec) moran.obs <- sum(zzt*t(zzt)*Wmat)/denom storage.moran <- matrix(NA, nrow = nsim, ncol = 1) for (j in 1:nsim) { rand.ord <- cbind(datavec,runif(length(datavec))) rand.ord <- rand.ord[order(rand.ord[,2]),][,1] zzt <- matrix(rand.ord - x.bar ,nrow= length(datavec), ncol = length(datavec)) storage.moran[j] <- sum(zzt*t(zzt)*Wmat)/denom } storage.moran <- storage.moran [order(storage.moran [,1]),] P.pos.moran <- sum(moran.obs < storage.moran)/nsim P.neg.moran <- 1 - P.pos.moran list(moran.obs = moran.obs, Pvalue.neg = P.neg.moran, Pvalue.pos = P.pos.moran, moran.null = storage.moran) } #------------ SCATTER PLOT GRAPH BASED ON MORAN'S I lat.graph.MoranI.scatter <- function(datavec, ndef.adj, ndef.num, graphic.device = "windows", output.filename) { if(graphic.device == "postscript") postscript(output.filename, onefile = TRUE, horizontal = FALSE) if(graphic.device == "windows") win.graph() ind <- !is.na(datavec) datavec1 <- datavec[ind] Nmat <- Neighmat(ndef.adj, ndef.num, length(ind)) Nmat <- Nmat[ind,ind] Wmat <- row.stand.fun(Nmat)$Wmat ind1 <- ndef.num[ind] > 0 datavec1 <- datavec1[ind1] Wmat <- Wmat[ind1, ind1] WZ <- Wmat %*% datavec1 minop <- min(datavec1, WZ) maxop <- max(datavec1, WZ) plot(matrix(c(minop,minop,maxop,maxop), nrow = 2, byrow = TRUE), type = "l", xlab = "Average of Neighbors", ylab = "Observed Value", col = "blue", lwd = 2, main = "MORAN'S I SCATTER PLOT", cex.axis = 1.3, cex.lab = 1.5) points(WZ, datavec1, pch = 19, col = "black", cex = 1.3) if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ MORAN'S I SCATTERGRAPH lat.graph.Moran.I <- function(Moran.I.out, graphic.device = "windows", output.filename) { if(graphic.device == "postscript") postscript(output.filename, onefile = TRUE, horizontal = FALSE) if(graphic.device == "windows") win.graph() hist(Moran.I.out$moran.null, xlim= c(min(Moran.I.out$moran.null,Moran.I.out$moran.obs), max(Moran.I.out$moran.null,Moran.I.out$moran.obs)), nclass = 20, col = "blue", xlab = "Moran's I", cex.axis = 1.3, cex.lab = 1.5, main = "MORAN'S I") points(Moran.I.out$moran.obs,1, pch = 19, cex = 2, col = "darkgreen") if(graphic.device == "postscript") dev.off(dev.cur()) } #------------ FUNCTION FOR FITTING A SPATIAL LINEAR MODEL TO GEOSTATISTICAL DATA slm.geo <- function(formula, data = data, xcol = xcol, ycol = xcol, group = NULL, EstMeth = "REML", CorModel = "Exponential", use.nugget = TRUE, use.spatial = TRUE, use.anisotropy = FALSE, use.locrep = FALSE, prediction.set = NULL, block.prediction.polygon = NULL, est.con.mat = NULL, est.con.label = NULL, use.FPBK = FALSE, FPBK.col = NULL) { if(use.anisotropy == TRUE & use.spatial == FALSE) return(list(errstate = 1, errmessage = "Cannot have anisotropy without spatial model")) if(use.locrep == TRUE & use.spatial == FALSE & use.nugget == FALSE) return(list(errstate = 1, errmessage = "Cannot have pure measurement error model")) if(use.spatial == FALSE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == FALSE) case <- "2" if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == FALSE & use.locrep == FALSE) case <- "12" if(use.spatial == TRUE & use.nugget == TRUE & use.anisotropy == TRUE & use.locrep == FALSE) case <- "123" case <- geo.case(use.spatial, use.nugget, use.anisotropy, use.locrep) response.col <- as.character(attr(terms(formula, data = data),"variables"))[2] data1 <- data[!is.na(data[,response.col]),] n.na <- 0 if(any(is.na(data[,response.col]))) { ind.na <- is.na(data[,response.col]) data.na <- data[ind.na,] n.na <- sum(ind.na) data.na[,response.col] <- 1 m.na <- model.frame(formula, data.na) X.na <- model.matrix(formula, m.na) } if(is.factor(data1[,response.col])) data1[,response.col] <- as.numeric(as.character(data1[,response.col])) if(is.character(data1[,response.col])) data1[,response.col] <- as.numeric(data1[,response.col]) z <- model.frame(formula, data1) z <- data1[,names(z)[1]] Xdesign <- X.design(formula = formula, data = data1, Xmethod = "treat.first.0") estcovparm.out <- est.covparm1(z, Xdesign, data = data1, xcol = xcol, ycol = ycol, EstMeth = EstMeth, CorModel = CorModel, group = group, use.spatial = use.spatial, use.nugget = use.nugget, use.anisotropy = use.anisotropy, use.locrep = use.locrep) n <- length(z) p <- sum(svd(Xdesign)$d>1e-10) Vi <- estcovparm.out$Vi V <- estcovparm.out$V covbi <- t(Xdesign) %*% Vi %*% Xdesign covb <- mginv(covbi) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(Xdesign) %*% Vi %*% z cv.out <- crossval(z = z, X = Xdesign, V = estcovparm.out$V, Vi = Vi, n = n, p = p, covb = covb, covbi = covbi, bhat = b.hat) cv.stats <- data.frame(bias = sum(z - cv.out[,1])/n, std.bias = sum((z - cv.out[,1])/sqrt(cv.out[,2]))/n, RMSPE = sqrt(sum((z - cv.out[,1])^2)/n), RAV = sqrt(sum(cv.out[,2]^2)/n), std.MSPE = sum((z - cv.out[,1])^2/cv.out[,2]^2)/n, cov.80 = sum(abs((z - cv.out[,1])/cv.out[,2]) < 1.281552)/n, cov.90 = sum(abs((z - cv.out[,1])/cv.out[,2]) < 1.644854)/n, cov.95 = sum(abs((z - cv.out[,1])/cv.out[,2]) < 1.959964)/n) infld <- infldiagn(z = z, V = V, X = Xdesign, n = n, p = p) fit <- data.frame(data1, fit = Xdesign %*% b.hat, resid = z - Xdesign %*% b.hat, infld, cv.out, cv.resid = z - cv.out[,1], cv.stndr = (z - cv.out[,1])/cv.out[,2]) m2LL <- estcovparm.out$m2LL if(EstMeth == "REML") { k <- length(estcovparm.out$theta) n1 <- n - p } if(EstMeth == "ML") { k <- p + length(estcovparm.out$theta) n1 <- n } AIC <- m2LL + 2*k AICc <- m2LL + 2*k*n1/(n1 - k - 1) BIC <- m2LL + log(n1)*k fixed.eff.est <- data.frame(parameters = b.hat, std.err = bhat.se, df = n - p, t.value = b.hat/bhat.se, prob.t = round(100000*(1- pt(abs(b.hat/bhat.se), df = n - p))*2)/100000 ) fixed.effects.estimates <- fixed.effect.table(formula = formula, data = data1, fixed.effects.estimates = fixed.eff.est) typeIII.hypoth <- hypoth(formula, data1, Vi) ind <- !is.na(fixed.effects.estimates[,4]) cov.fix.eff <- matrix(0, nrow = length(fixed.effects.estimates[,4]), ncol = length(fixed.effects.estimates[,4])) cov.fix.eff[ind,ind] <- covb pred.data.krige <- NULL if(!is.null(prediction.set)) pred.data.krige <- pred.krige.Vi(formula, data.pred = prediction.set, data = data1, resp.var = response.col, xcol = xcol, ycol = ycol, covb = covb, Vi = Vi, CorModel = CorModel, theta = estcovparm.out$theta, case = case) if(n.na > 0) {na.data.krige <- pred.krige.Vi(formula, data.pred = data.na, data = data1, resp.var = response.col, xcol = xcol, ycol = ycol, covb = covb, Vi = Vi, CorModel = CorModel, theta = estcovparm.out$theta, case = case) data <- data.frame(data, pred = rep(NA, times = length(data[,1])), pred.se = rep(NA, times = length(data[,1]))) data[ind.na, "pred"] <- na.data.krige[,3] data[ind.na, "pred.se"] <- na.data.krige[,4] colnames(data)[colnames(data) == "pred"] <- paste("pred.",response.col, sep = "") } block.krige <- NULL if(!is.null(block.prediction.polygon)){ if(is.null(prediction.set)) block.krige <- "Need prediction set" else block.krige <- block.krige.Vi( formula = formula, data = data1, data.pred = prediction.set, rshape.outline = block.prediction.polygon, resp.var = resp.var, CorModel = CorModel, xcol = xcol, ycol = ycol, covb = covb, Vi = Vi, V = estcovparm.out$V, resid = fit["resid"], beta.hat = b.hat, theta = estcovparm.out$theta, case = case) } estimates.contrasts <- NULL if(!is.null(est.con.mat)){ estimates.contrasts <- estimates.and.contrasts(formula = formula, data = data1, est.con.mat = est.con.mat, est.con.lab = est.con.label, fixed.effects.estimates = fixed.effects.estimates, cov.fix.eff = cov.fix.eff, p = p, n = n) } FPBK.out <- NULL if(use.FPBK == TRUE){ FPBK.out <- FPBK(formula = formula, data = data, xcol = xcol, ycol = ycol, FPBK.col = FPBK.col, theta = estcovparm.out$theta, CorModel = CorModel, case = case, V = V, Vi = Vi) } list( sample.size = n + n.na, obs.sample.size = n, missing.sample.size = n.na, CorModel = CorModel, EstMeth = EstMeth, fixed.effects.estimates = fixed.effects.estimates, typeIII.hypoth = typeIII.hypoth, R2g = R2g(z = z, Xdesign = Xdesign, bhat = b.hat, Vi = Vi), cov.fix.eff = cov.fix.eff, theta = estcovparm.out$theta, m2LL = m2LL, AIC = AIC, AICc = AICc, BIC = BIC, cv.stats = cv.stats, data = data, fit = fit, Predictions = pred.data.krige, Block.Prediction = block.krige, FPBK = FPBK.out, Estimates.Contrasts = estimates.contrasts, V = estcovparm.out$V, Vi = Vi, X = Xdesign, case = case ) } #------------ FUNCTION FOR FITTING A SPATIAL LINEAR MODEL TO LATTICE DATA slm.lat <- function(formula, data, ndef.adj, ndef.num, CorModel = "CAR", EstMeth = "REML", use.row.standard = TRUE, use.distance = FALSE, weight.col = NULL, xcol = NULL, ycol = NULL, create.smooth = FALSE, est.con.mat = NULL, est.con.label = NULL, alpha = 0, use.FPBK = FALSE, FPBK.col = NULL) { response.col <- as.character(attr(terms(formula, data = data),"variables"))[2] ind <- !is.na(data[,response.col]) data1 <- data[ind,] n.na <- 0 if(any(is.na(data[,response.col]))) { ind.na <- is.na(data[,response.col]) data.na <- data[ind.na,] n.na <- sum(ind.na) data.na[,response.col] <- 1 m.na <- model.frame(formula, data.na) X.na <- model.matrix(formula, m.na) } if(is.factor(data1[,response.col])) data1[,response.col] <- as.numeric(as.character(data1[,response.col])) if(is.character(data1[,response.col])) data1[,response.col] <- as.numeric(data1[,response.col]) z <- model.frame(formula, data1) z <- data1[,names(z)[1]] X <- X.design(formula = formula, data = data1, Xmethod = "treat.first.0") n <- length(z) p <- sum(svd(X)$d>1e-10) parmest <- est.covparm.lat(z, X, data = data1, data.all = data, ind, n, p, ndef.adj = ndef.adj, ndef.num = ndef.num, CorModel = CorModel, EstMeth = EstMeth, use.row.standard = use.row.standard, use.distance = use.distance, weight.col = weight.col, xcol = xcol, ycol = ycol, alpha = alpha) V <- parmest$V Vi <- parmest$Vi V.all <- parmest$V.all Vi.all <- parmest$Vi.all covbi <- t(X) %*% Vi %*% X covb <- mginv(covbi) bhat.se <- sqrt(diag(covb)) b.hat <- covb %*% t(X) %*% Vi %*% z cv.out <- crossval(z = z, X = X, V = V, Vi = Vi, n = n, p = p, covb = covb, covbi = covbi, bhat = b.hat) data <- data.frame(data, smooth = rep(NA, times = length(data[,1])), smooth.se =rep(NA, times = length(data[,1]))) data[ind,"smooth"] <- cv.out[,1] data[ind,"smooth.se"] <- cv.out[,2] cv.stats <- data.frame(bias = sum(z - cv.out[,1])/n, std.bias = sum((z - cv.out[,1])/sqrt(cv.out[,2]))/n, RMSPE = sqrt(sum((z - cv.out[,1])^2)/n), RAV = sqrt(sum(cv.out[,2]^2)/n), std.MSPE = sum((z - cv.out[,1])^2/cv.out[,2]^2)/n, cov.80 = sum(abs((z - cv.out[,1])/cv.out[,2]) < 1.281552)/n, cov.90 = sum(abs((z - cv.out[,1])/cv.out[,2]) < 1.644854)/n, cov.95 = sum(abs((z - cv.out[,1])/cv.out[,2]) < 1.959964)/n) infld <- infldiagn(z = z, V = V, X = X, n = n, p = p) fit <- data.frame(data1, fit = X %*% b.hat, resid = z - X %*% b.hat, infld, cv.out, cv.resid = z - cv.out[,1], cv.stndr = (z - cv.out[,1])/cv.out[,2]) m2LL <- parmest$m2LL if(EstMeth == "REML") { k <- length(parmest$theta) n1 <- n - p } if(EstMeth == "ML") { k <- p + length(parmest$theta) n1 <- n } AIC <- m2LL + 2*k AICc <- m2LL + 2*k*n1/(n1 - k - 1) BIC <- m2LL + log(n1)*k fixed.eff.est <- data.frame(parameters = b.hat, std.err = bhat.se, df = n - p, t.value = b.hat/bhat.se, prob.t = round(100000*(1- pt(abs(b.hat/bhat.se), df = n - p))*2)/100000 ) fixed.effects.estimates <- fixed.effect.table(formula = formula, data = data1, fixed.effects.estimates = fixed.eff.est) typeIII.hypoth <- hypoth(formula, data1, Vi) ind1 <- !is.na(fixed.effects.estimates[,4]) cov.fix.eff <- matrix(0, nrow = length(fixed.effects.estimates[,4]), ncol = length(fixed.effects.estimates[,4])) cov.fix.eff[ind1,ind1] <- covb if(n.na > 0) {na.data.krige <- pred.krige.lat(z = z, X = X, X.pred = X.na, ind = ind, V.all = V.all, Vi = Vi, covb = covb) data <- data.frame(data, pred = rep(NA, times = length(data[,1])), pred.se = rep(NA, times = length(data[,1]))) data[ind.na, "pred"] <- na.data.krige[,1] data[ind.na, "pred.se"] <- na.data.krige[,2] data[ind.na, "smooth"] <- na.data.krige[,1] data[ind.na, "smooth.se"] <- na.data.krige[,2] colnames(data)[colnames(data) == "pred"] <- paste("pred.",response.col, sep = "") } estimates.contrasts <- NULL if(!is.null(est.con.mat)){ estimates.contrasts <- estimates.and.contrasts(formula = formula, data = data1, est.con.mat = est.con.mat, est.con.lab = est.con.label, fixed.effects.estimates = fixed.effects.estimates, cov.fix.eff = cov.fix.eff, p = p, n = n) } FPBK.out <- NULL if(use.FPBK == TRUE){ FPBK.out <- FPBK.lat(z.sa = z, Xs = X, Xu = X.na, ind.sa = ind, data = data, FPBK.col = FPBK.col, V.all = V.all, Vi = Vi) } fit.4.data <- matrix(NA, nrow = n + n.na, ncol = 10) fit.4.data[ind,1:10] <- as.matrix(fit[,c("fit","resid","resid.stand","resid.student","leverage", "CooksD", "cv.pred", "cv.se", "cv.resid", "cv.stndr")]) fit.4.data <- as.data.frame(fit.4.data) colnames(fit.4.data) <- c("fit","resid","resid.stand","resid.student","leverage", "CooksD", "cv.pred", "cv.se", "cv.resid", "cv.stndr") data <- data.frame(data,fit.4.data) list( sample.size = n + n.na, obs.sample.size = n, missing.sample.size = n.na, CorModel = CorModel, EstMeth = EstMeth, fixed.effects.estimates = fixed.effects.estimates, typeIII.hypoth = typeIII.hypoth, R2g = R2g(z = z, Xdesign = X, bhat = b.hat, Vi = Vi), cov.fix.eff = cov.fix.eff, theta = parmest$theta, rho.restrict = parmest$rho.restrict, m2LL = m2LL, AIC = AIC, AICc = AICc, BIC = BIC, cv.stats = cv.stats, data = data, fit = fit, FPBK = FPBK.out, Estimates.Contrasts = estimates.contrasts, V = V, Vi = Vi, X = X ) }