# s - (n x p) matrix of locations # x - (n x 2) matrix of location-variable index # Y - (n x d) matrix of response values # Z - (n x q) matrix of additional covariates or a list of (n x q_i) matrices multiKrig <- function(s,Y,Z, cov.function="multi.cov",cov.args=NULL, wght.function="multi.wght",wght.args=NULL, verbose=FALSE){ # spatial locations s <- as.matrix(s) if(is.null(dimnames(s))){ nms.s <- paste('s[,',1:ncol(s),']',sep='') } else { nms.s <- dimnames(s)[[2]] } # covariates if(!is.null(Z)){ if(is.list(Z)){ Z <- lapply(Z,as.matrix) } else { Z <- as.matrix(Z) if(is.null(dimnames(Z))){ nms.Z <- paste('Z[,',1:ncol(Z),']',sep='') } else { nms.Z <- dimnames(Z)[[2]] } } } # response variables Y <- as.matrix(Y) if(is.null(dimnames(Y))){ nms.Y <- paste("V",1:ncol(Y),sep="") } else { nms.Y <- dimnames(Y)[[2]] } d <- ncol(Y) n <- nrow(Y) Y <- c(Y) # set-up x matrix x <- expand.grid(1:n,1:d) # set-up Z matrix if(is.null(Z)){ nZ <- kronecker(diag(d),s) } else { if(is.list(Z)){ qq <- c(0,cumsum(sapply(Z,ncol)+ncol(s))) nZ <- matrix(0,d*nrow(s),qq[length(qq)]) for(i in 1:length(Z)){ nZ[((i-1)*n+1):(i*n),(qq[i]+1):qq[i+1]] <- cbind(s,Z[[i]]) } } else { nZ <- kronecker(diag(d),cbind(s,Z)) } } # add spatial locations to covariance arguments cov.args$s <- s # call to Krig obj <- Krig(x=x,Y=c(Y),Z=nZ,method='REML', null.function="multi.null", cov.function="multi.cov",cov.args=cov.args, wght.function='multi.wght',wght.args=wght.args, verbose=verbose) # set up multiKrig object obj$dimY <- c(n,d) obj$nms.Y <- nms.Y obj$nms.s <- nms.s obj$nms.Z <- nms.Z class(obj) <- "multiKrig" # out return(obj) } predict.multiKrig <- function(obj,s,Z=NULL){ n <- obj$dimY[1] d <- obj$dimY[2] nx <- expand.grid((n+1):(n+nrow(s)),1:d) obj$args$s <- rbind(obj$args$s,s) if(is.null(Z)){ nZ <- kronecker(diag(d),s) } else { if(is.list(Z)){ qq <- c(0,cumsum(sapply(Z,ncol)+ncol(s))) nZ <- matrix(0,d*nrow(s),qq[length(qq)]) for(i in 1:length(Z)){ nZ[((i-1)*n+1):(i*n),(qq[i]+1):qq[i+1]] <- cbind(s,Z[[i]]) } } else { nZ <- kronecker(diag(d),as.matrix(cbind(s,Z))) } } out <- predict.Krig(obj,x=nx,Z=nZ) return(out) } summary.multiKrig <- function(obj){ n <- obj$dimY[1] d <- obj$dimY[2] if(!is.list(Z)){ if(obj$cov.function=='multi.cov' && obj$wght.function=='multi.wght'){ X <- do.call(obj$null.function,list(x=obj$x,Z=obj$Z)) O1 <- obj$rhohat*multi.cov(obj$x,obj$x,marginal=FALSE,C=NA, s=obj$args$s,theta=obj$args$theta,rho=obj$args$rho,smoothness=obj$args$smoothness) sp <- c(1,obj$wght.args$sp) S <- matrix(0,d,d) S[upper.tri(S,diag=TRUE)] <- sp S <- (S+t(S)) diag(S) <- diag(S)/2 S <- obj$shat.MLE^2*S O2 <- kronecker(S,diag(n)) O <- O1 + O2 V <- solve(t(X)%*%solve(O,X)) beta <- obj$d mbeta <- matrix(NA,nrow=length(beta)/d,ncol=d) mv <- matrix(NA,nrow=length(beta)/d,ncol=d) w <- cbind(1,contr.treatment(d)) mbeta[1,] <- w %*% beta[1:d] mv[1,] <- sqrt(diag(w %*% V[1:d,1:d] %*% t(w))) mbeta[-1,] <- matrix(beta[(d+1):length(beta)],ncol=d) mv[-1,] <- matrix(sqrt(diag(V)[(d+1):length(beta)]),ncol=d) betaout <- vector("list",d) for(i in 1:d){ betaout[[i]] <- cbind(mbeta[,i],mv[,i],mbeta[,i]/mv[,i]) dimnames(betaout[[i]]) <- list(c("Int",obj$nms.s,obj$nms.Z), c("Coef","SE","Z")) } names(betaout) <- obj$nms.Y } else { beta <- obj$d betaout <- matrix(NA,nrow=length(beta)/d,ncol=d) betaout[1,] <- cumsum(beta[1:d]) betaout[-1,] <- matrix(beta[(d+1):length(beta)],ncol=d) dimnames(betaout) <- list(c("Int",obj$nms.s,obj$nms.Z), obj$nms.Y) } } else { betaout <- obj$d } cat('\n') cat('regression coefficients:\n') print(betaout) cat('\n') cat('spatial covariance function:',obj$cov.function,'\n') if(obj$cov.function=='multi.cov'){ cat('rho:\n') print(diag(obj$rhohat*c(1,obj$args$rho))) } cat('\n') cat('error weight function:',obj$wght.function,'\n') if(obj$wght.function=='multi.wght'){ cat('S:\n') sp <- c(1,obj$wght.args$sp) S <- matrix(0,d,d) S[upper.tri(S,diag=TRUE)] <- sp S <- (S+t(S)) diag(S) <- diag(S)/2 print(obj$shat.MLE^2*S) } invisible() } print.multiKrig <- function(obj){ summary(obj) } plot.multiKrig <- function(obj,ask=FALSE){ Y <- matrix(obj$y,nrow=obj$dimY[1],ncol=obj$dimY[2]) p <- matrix(obj$fitted.values,nrow=obj$dimY[1],ncol=obj$dimY[2]) r <- matrix(obj$residuals,nrow=obj$dimY[1],ncol=obj$dimY[2]) op <- par(mar=c(5,4,4,2)+0.1) for(i in 1:obj$dimY[2]){ zl <- range(c(Y[,i],p[,i])) quilt.plot(obj$args$s[,1],obj$args$s[,2],Y[,i],zlim=zl, xlab=obj$nms.s[1],ylab=obj$nms.s[2], main=paste(obj$nms.Y[i],': observations',sep='')) quilt.plot(obj$args$s[,1],obj$args$s[,2],p[,i],zlim=zl, xlab=obj$nms.s[1],ylab=obj$nms.s[2], main=paste(obj$nms.Y[i],': fitted values',sep='')) m <- max(abs(r[,i])) zl <- c(-m,m) quilt.plot(obj$args$s[,1],obj$args$s[,2],r[,i],zlim=zl, xlab=obj$nms.s[1],ylab=obj$nms.s[2], main=paste(obj$nms.Y[i],': residuals',sep='')) } par(op) } multi.null <- function(x,Z=NULL,drop.Z=FALSE){ data <- data.frame(a=as.factor(x[,2])) X <- model.matrix(~a,data=data,contrasts=list(a="contr.treatment")) if (is.null(Z) | drop.Z) { return(X) } else { return(cbind(X,Z)) } } # No spatial cross-correlatons multi.cov <- function(x1,x2,marginal=FALSE,C=NA,s,theta,rho,smoothness){ if(marginal){ return(rep(1,d*nrow(x1))) } ind <- unique(x1[,2]) temp <- stationary.cov(s[x1[,1][x1[,2]==1],], s[x2[,1][x2[,2]==1],], Covariance='Matern',theta=theta,smoothness=smoothness) if(length(ind)>1){ for(i in 2:length(ind)){ temp2 <- rho[i-1]*stationary.cov(s[x1[,1][x1[,2]==i],], s[x2[,1][x2[,2]==i],], Covariance='Matern',theta=theta,smoothness=smoothness) d1 <- dim(temp) d2 <- dim(temp2) temp <- rbind(cbind(temp,matrix(0,d1[1],d2[2])), cbind(matrix(0,d1[1],d2[2]),temp2)) } } if(is.na(C[1])){ return(temp) } else { return(temp%*%C) } } # Kronecker form multi.wght <- function(x,sp){ n <- length(unique(x[,1])) d <- length(unique(x[,2])) sp <- c(1,sp) S <- matrix(0,d,d) S[upper.tri(S,diag=TRUE)] <- sp S <- (S+t(S)) diag(S) <- diag(S)/2 return(kronecker(solve(S),diag(n))) } # A simple 2 variable example multi.est <- function(pms,tol=0.0001, s=s,x=x,Y=Y,Z=NULL){ flag <- TRUE while(flag){ old <- c(pms$rho,pms$sp) out <- multiKrig(s=s,Y=Y[,1:2],Z=Z, cov.function='multi.cov', cov.args=list(theta=2.5,rho=pms$rho[2]/pms$rho[1],smoothness=1.5), wght.function='multi.wght', wght.args=list(sp=pms$sp[2:3]/pms$sp[1])) pms$rho[1] <- out$rhohat # update S u <- matrix(out$residuals,ncol=ncol(Y)) S <- t(u)%*%u/nrow(u) pms$sp <- c(S[1,1],S[1,2],S[2,2]) #cat('one:',pms$rho,pms$sp,'\n') out <- multiKrig(s=s,Y=Y[,2:1],Z=Z, cov.function='multi.cov', cov.args=list(theta=2.5,rho=pms$rho[1]/pms$rho[2],smoothness=1.5), wght.function='multi.wght', wght.args=list(sp=pms$sp[2:1]/pms$sp[3])) pms$rho[2] <- out$rhohat # update S u <- matrix(out$residuals,ncol=ncol(Y)) S <- t(u)%*%u/nrow(u) pms$sp <- c(S[2,2],S[1,2],S[1,1]) #cat('two:',pms$rho,pms$sp,'\n') new <- c(pms$rho,pms$sp) crit <- max(abs((new-old)/old)) if(crit < tol){ flag <- FALSE } cat(crit,pms$rho,pms$sp,'\n') } return(list(out=out,rho=pms$rho,sp=pms$sp)) }