"psp2dM"<- function(y, M, p1, p2, M.type = "stacked", M1.index = NULL, M2.index = NULL, Pars, ridge.adj = 0, M.pred = M, y.predicted = NULL, x.lab = "X1", y.lab = "X2", z.lab = "A.hat", coef.plot = T, image.plot = T, se.bands = T, family = "gaussian", link = "default", m.binomial = NULL, wts = NULL, r.gamma = NULL, int = F) { # P-spline regression with 2-D regressors # Input: # y: response (of length m) # M (stacked): (mp1) x p2 matrix stacked with m regressor matrices each of p1xp2 # M: (unfolded) m x (p1 x p2) regressor matrix with m regressor rows, each length p1xp2 # M.type: "stacked" or "unfolded" # p1, p2: dimensions of regressor matrices # M1.index (M2.index): index for rows (cols) of regressor matrix; default is simple seq # Pars: 2 rows with P-spline parameters: [min max nseg deg lambda pdeg] # M.pred (stacked): qp1 x p2 inputs for q new predicitons # M.pred (unfolded): q x (p1 x p2) inputs for q new predictions # y.predicted= a vector of responses from a cv data set (assoc. with cv M.pred) # ridge.adj: small ridge penalty to stabilize estimation, can be zero # x.lab, y.lab, z.lab: "character" labels for estimated alpha coefficient surface # coef.plot: T or F for perspective plot of coefficient surface # image.plot: T or F for image plot of coefficient surface # se.bands: T or F for twice standard error surface plots # family: "gaussian", "binomial", "poisson", "Gamma" # link: "logit", "probit", "log", "sqrt", "loglog", "cloglog", "identity", "inverse" # wts: non-negative weights (can be zero) # m.binomial: number of trials associated with binomial r.v. (can vary) # r.gamma: gamma scale parameter # int: T or F for intercept column of ones # # Output: list with elements # coef: tensor product P-spline coefficients # summary.predicted: predicted value at new 2-D regressor locations (with 2 se bands) # cv: cross-validation statistic # aic, dev, df.residual # eff.dim: effective df of estimated 2-D estimated coefficient surface # perspective plot and image plot of estimated alpha coefficient surface # # Support functions needed: bsplbase(), pspline2d.checker(), pspline.fitter() # # Paul Eilers (2000) and Brian Marx (2001, 2002) (c) # Prepare bases for estimation m <- length(y) if(missing(wts)) { wts <- rep(1, m) } if(missing(m.binomial)) { m.binomial <- rep(1, m) } if(missing(r.gamma)) { r.gamma <- rep(1, m) } parms <- pspline2d.checker(family, link, Pars[1, 4], Pars[2, 4], Pars[1, 6], Pars[2, 6], Pars[1, 3], Pars[2, 3], Pars[1, 5], Pars[2, 5], ridge.adj, wts) family <- parms$family link <- parms$link ridge.adj <- parms$ridge.adj wts <- parms$wts Pars[1, 3:6] <- c(parms$ps.intervals1, parms$degree1, parms$lambda1, parms$order1) Pars[2, 3:6] <- c(parms$ps.intervals2, parms$degree2, parms$lambda2, parms$order2) M <- X <- as.matrix(M) if(M.type == "stacked") { p1. <- nrow(M)/m p2. <- ncol(M) if(p1 != p1. | p2 != p2.) { warning(paste( "recheck input p1 or p2: at least one dimension does not match" )) } x <- as.vector(t(M)) X <- matrix(x, m, p1 * p2, byrow = T) } if(missing(M1.index)) { M1.index <- 1:p1 } if(missing(M2.index)) { M2.index <- 1:p2 } oM1 <- outer(rep(1, p2), M1.index) B1 <- bsplbase(as.vector(oM1), Pars[1, ]) oM2 <- outer(M2.index, rep(1, p1)) B2 <- bsplbase(as.vector(oM2), Pars[2, ]) n1 <- ncol(B1) n2 <- ncol(B2) # Compute tensor products for estimated alpha surface B1. <- kronecker(B1, t(rep(1, n2))) B2. <- kronecker(t(rep(1, n1)), B2) B. <- B1. * B2. # Construct penalty matrices #browser() #----- # B3 <- bsplbase(M1.index, Pars[1, ]) # B4 <- bsplbase(M2.index, Pars[2, ]) # V <- matrix(rep(0, m * n1 * n2), m, n1 * n2) # for(i in 1:m) { # V.in <- (t(B3) %*% M[((i - 1) * p1 + 1):(i * p1), # ] %*% B4) # V[i, ] <- as.vector((V.in)) # } d1 <- Pars[1, 6] D1 <- diag(n1) if(d1 != 0) { for(j in 1:d1) { D1 <- diff(D1) } } lambda1 <- Pars[1, 5] P1 <- sqrt(lambda1) * kronecker(D1, diag(n2)) d2 <- Pars[2, 6] D2 <- diag(n2) if(d2 != 0) { for(j in 1:d2) { D2 <- diff(D2) } } lambda2 <- Pars[2, 5] P2 <- sqrt(lambda2) * kronecker(diag(n1), D2) Pen <- rbind(P1, P2) p.ridge <- NULL if(ridge.adj > 0) { nix.ridge <- rep(0, n1 * n2) p.ridge <- sqrt(ridge.adj) * diag(n1 * n2) } # Data augmentation and regression z1 <- rep(0, n2 * (n1 - d1)) z2 <- rep(0, n1 * (n2 - d2)) Q <- X %*% B. n.col <- ncol(Q) if(int) { Q <- cbind(rep(1, nrow(Q)), Q) Pen <- cbind(rep(0, nrow(Pen)), Pen) if(ridge.adj > 0) { p.ridge <- cbind(rep(0, nrow(p.ridge)), p.ridge) } } ps.fit <- pspline.fitter(family, link, n.col, m.binomial, r.gamma, y, b = Q, Pen, p.ridge, nix = c(z1, z2), nix.ridge = rep(0, n1 * n2 ), ridge.adj, wts) mu <- ps.fit$mu pcoef <- ps.fit$coef bin.percent.correct <- bin.0percent <- bin.1percent <- NULL if(family == "binomial") { count1 <- count2 <- pcount1 <- pcount2 <- 0 p.hat <- mu/m.binomial for(ii in 1:m) { if(p.hat[ii] > 0.5) { count1 <- y[ii] count1 <- pcount1 + count1 pcount1 <- count1 } if(p.hat[ii] <= 0.5) { count2 <- m.binomial[ii] - y[ii] count2 <- pcount2 + count2 pcount2 <- count2 } } bin.percent.correct <- (count1 + count2)/sum(m.binomial) bin.1percent <- count1/sum(m.binomial[p.hat > 0.5]) bin.0percent <- count2/sum(m.binomial[p.hat <= 0.5]) } w <- ps.fit$w e <- 1e-009 h <- hat(ps.fit$f$qr, intercept = F)[1:m] trace <- eff.dim <- sum(h) if(family == "binomial") { dev <- 2 * sum((y + e) * log((y + e)/mu) + (m.binomial - y + e) * log((m.binomial - y + e)/(m.binomial - mu))) dispersion.parm <- 1 cv <- NULL } if(family == "poisson") { dev <- 2 * sum(y * log(y + e) - y - y * log(mu) + mu) dispersion.parm <- 1 cv <- NULL } if(family == "Gamma") { dev <- -2 * sum(r.gamma * (log((y + e)/mu) - ((y - mu)/mu))) ave.dev <- dev/m dispersion.parm <- (ave.dev * (6 + ave.dev))/(6 + 2 * ave.dev) cv <- NULL } cv <- press.mu <- press.e <- var.c <- NULL if(family == "gaussian") { dev <- sum(ps.fit$f$residuals[1:m]^2) dispersion.parm <- dev/(m - trace) press.e <- ps.fit$f$residuals[1:m]/(1 - h) cv <- sqrt(sum((press.e)^2)/(m)) press.mu <- y - press.e } aic <- dev + 2 * trace w.aug <- c(w, (c(z1, z2) + 1)) if(int) { A.hat <- B. %*% pcoef[2:(n.col + 1)] yint <- ps.fit$coef[1] } if(!int) { yint <- NULL A.hat <- B. %*% pcoef } A.hatm <- matrix(A.hat, p1, p2, byrow = T) if(coef.plot) { i.1 <- M1.index in1 <- 1:length(M1.index) i.2 <- M2.index in2 <- 1:length(M2.index) if(length(M1.index) > 100) { i.1 <- round(seq(from = min(M1.index), to = max( M1.index), length = 100)) in1 <- round(seq(from = 1, to = p1, length = 50)) } if(length(M2.index) > 100) { i.2 <- round(seq(from = min(M2.index), to = max( M2.index), length = 100)) in2 <- round(seq(from = 1, to = p2, length = 50)) } persp(M2.index[in2], (M1.index[in1]), t(A.hatm[in1, in2]), xlab = y.lab, ylab = x.lab, zlab = z.lab) if(image.plot) { image(M2.index[in2], M1.index[in1], t(A.hatm[in1, in2]), xlab = x.lab, ylab = y.lab, sub = "A hat") } if(se.bands) { half.meat <- sqrt(c(w)) * Q meat <- t(half.meat) %*% half.meat if(ridge.adj > 0) { bread <- solve(meat + t(Pen) %*% Pen + t( p.ridge) %*% p.ridge) } if(ridge.adj == 0) { bread <- solve(meat + t(Pen) %*% Pen) } half.sw <- half.meat %*% bread[, (1 + int):(n.col + int )] var.c <- t(half.sw) %*% half.sw half.lunch <- half.sw %*% t(B.) ones <- 0 * y + 1 var.Ahat <- ones %*% (half.lunch * half.lunch) stdev.Ahat <- sqrt(dispersion.parm) * t(sqrt(var.Ahat)) pivot <- 2 * stdev.Ahat upper <- A.hat + pivot lower <- A.hat - pivot L.hatm <- matrix(lower, p1, p2, byrow = T) U.hatm <- matrix(upper, p1, p2, byrow = T) persp(M2.index[in2], M1.index[in1], t(U.hatm[in1, in2]), xlab = x.lab, ylab = y.lab, zlab = "2 se Upper Surface") if(image.plot) { image(M2.index[in2], M1.index[in1], t(U.hatm[ in1, in2]), xlab = x.lab, ylab = y.lab, sub = "2 se Upper Surface") } persp(M2.index[in2], M1.index[in1], t(L.hatm[in1, in2]), xlab = x.lab, ylab = y.lab, zlab = "2 se Lower Surface") if(image.plot) { image(M2.index[in2], M1.index[in1], t(L.hatm[ in1, in2]), xlab = x.lab, ylab = y.lab, sub = "2 se Lower Surface") } } } summary.predicted <- NULL cv.predicted <- eta.predicted <- avediff.pred <- NULL if(!missing(M.pred)) { half.meat <- sqrt(c(w)) * Q meat <- t(half.meat) %*% half.meat if(ridge.adj > 0) { bread <- solve(meat + t(Pen) %*% Pen + t(p.ridge) %*% p.ridge) } if(ridge.adj == 0) { bread <- solve(meat + t(Pen) %*% Pen) } half.sw <- half.meat %*% bread[, (1 + int):(n.col + int)] var.c <- t(half.sw) %*% half.sw q <- nrow(M.pred) X.p <- as.matrix(M.pred) if(M.type == "stacked") { M.pred <- as.matrix(M.pred) q <- nrow(M.pred)/p1 X.p <- matrix(x, q, p1 * p2, byrow = T) } if(!int) { eta.predicted <- X.p %*% as.vector(A.hat) var.pred <- X.p %*% B. %*% var.c %*% t(X.p %*% B.) } if(int) { var.c <- t(bread) %*% t(half.meat) %*% half.meat %*% bread one.xpred.b <- cbind(rep(1, q), (X.p %*% B.)) eta.predicted <- X.p %*% A.hat + yint var.pred <- one.xpred.b %*% var.c %*% t(one.xpred.b) } stdev.pred <- as.vector(sqrt(diag(var.pred))) stdev.pred <- sqrt(dispersion.parm) * stdev.pred pivot <- as.vector(2 * stdev.pred) upper <- eta.predicted + pivot lower <- eta.predicted - pivot summary.predicted <- cbind(lower, eta.predicted, upper) if(!missing(y.predicted)) { if(family == "gaussian") { cv.predicted <- sqrt(sum((y.predicted - eta.predicted)^2)/(length(y.predicted))) avediff.pred <- (sum(y.predicted - eta.predicted))/length(y.predicted) } } bin.percent.correct <- bin.0percent <- bin.1percent <- NULL if(link == "logit") { summary.predicted <- 1/(1 + exp( - summary.predicted)) count1 <- count2 <- pcount1 <- pcount2 <- 0 p.hat <- exp(eta.predicted)/(1 + exp(eta.predicted)) if(!missing(y.predicted)) { y.predicted <- as.vector(y.predicted) for(ii in 1:length(eta.predicted)) { if(p.hat[ii] > 0.5) { count1 <- y.predicted[ii] count1 <- pcount1 + count1 pcount1 <- count1 } if(p.hat[ii] <= 0.5) { count2 <- 1 - y.predicted[ii] count2 <- pcount2 + count2 pcount2 <- count2 } } bin.percent.correct <- (count1 + count2)/length( y.predicted) bin.1percent <- count1/sum(y.predicted) bin.0percent <- count2/(length(y.predicted) - sum(y.predicted)) } } if(link == "probit") { summary.predicted <- apply(summary.predicted, c(1, 2), pnorm) } if(link == "cloglog") { summary.predicted <- (1 - exp( - exp(summary.predicted) )) } if(link == "loglog") { summary.predicted <- exp( - exp( - summary.predicted)) } if(link == "sqrt") { summary.predicted <- summary.predicted^2 } if(link == "log") { summary.predicted <- exp(summary.predicted) } if(link == "recipical") { summary.predd <- 1/(summary.predicted) summary.predicted[, 1] <- summary.predd[, 3] summary.predicted[, 3] <- summary.predd[, 1] summary.predd <- NULL } summary.predicted <- as.matrix(summary.predicted) dimnames(summary.predicted) <- list(NULL, c("-2std_Lower", "Predicted", "+2std_Upper")) } P <- list(coef = pcoef, Pars = Pars, yint = yint, int = int, family, link, dev = dev, aic = aic, bin.percent.correct = bin.percent.correct, bin.0 = bin.0percent, bin.1 = bin.1percent, df.resid = m - trace, dispersion.parm = dispersion.parm, mu = mu, press.mu = press.mu, summary.predicted = summary.predicted, eta.predicted = eta.predicted, avediff.pred = avediff.pred, ridge.adj = ridge.adj, cv = cv, cv.predicted = cv.predicted, eff.dim = eff.dim) P } "bsplbase"<- function(x, bpars) { # Compute a B-spline basis # Input: # x: abcissae # bpars: B-spline parameters: xmin, xmax, nseg, degree (= one less than "order") # Output: # base: matrix with nrow = length(x) and nseg + degree columns # # Paul Eilers, 2000 dx <- (bpars[2] - bpars[1])/bpars[3] knots <- seq(bpars[1] - bpars[4] * dx, bpars[2] + bpars[4] * dx, by = dx) base <- as.matrix(spline.des(knots, x, bpars[4] + 1, 0 * x)$design) base } "pspline2d.checker"<- function(family, link, degree1, degree2, order1, order2, ps.intervals1, ps.intervals2, lambda1, lambda2, ridge.adj, wts) { if(link == "default" && family == "gaussian") { link <- "identity" } if(link == "default" && family == "poisson") { link <- "log" } if(link == "default" && family == "binomial") { link <- "logit" } if(link == "default" && family == "Gamma") { link <- "log" } if(family != "binomial" && family != "gaussian" && family != "poisson" && family != "Gamma") { warning(paste("Improper FAMILY option. Choose: gaussian, poisson, binomial or Gamma" )) } if((family == "binomial") && (link != "logit" && link != "probit" && link != "cloglog" && link != "loglog")) { warning(paste("Improper LINK option with family=binomial. Choose: logit, probit, loglog, cloglog" )) } if((family == "Gamma") && (link != "log" && link != "recipical" && link != "identity")) { warning(paste("Improper LINK option with family=Gamma. Choose: recipical, log, identity" )) } if((family == "poisson") && (link != "log" && link != "sqrt" && link != "identity")) { warning(paste("Improper LINK option with family=poisson. Choose: log, sqrt, identity" )) } if((family == "gaussian") && (link != "identity")) { warning(paste("Improper LINK option with family=gaussian. Choose: identity" )) } if(degree1 < 0) { degree1 <- 1 warning(paste("degree1 must be non-neg integer: have used 1")) } if(order1 < 0) { order1 <- 0 warning(paste("order1 must be non-neg integer: have used 0")) } if(ps.intervals1 < 2) { ps.intervals1 <- 2 warning(paste("ps.intervals1 must be positive integer, > 1: have used 2" )) } if(lambda1 < 0) { lambda1 <- 0 warning(paste("lambda1 cannot be negative: have used 0")) } if(degree2 < 0) { degree2 <- 1 warning(paste("degree2 must be non-neg integer: have used 1")) } if(order2 < 0) { order2 <- 0 warning(paste("order2 must be non-neg integer: have used 0")) } if(ps.intervals2 < 2) { ps.intervals2 <- 2 warning(paste("ps.intervals2 must be positive integer, > 1: have used 2" )) } if(lambda2 < 0) { lambda2 <- 0 warning(paste("lambda2 cannot be negative: have used 0")) } if(ridge.adj < 0) { ridge.adj <- 0 warning(paste("ridge.adj cannot be negative: have used 0")) } if(min(wts) < 0) { warning(paste("At least one weight entry is negative")) } llist <- list(family = family, link = link, degree1 = degree1, order1 = order1, ps.intervals1 = ps.intervals1, lambda1 = lambda1, degree2 = degree2, order2 = order2, ps.intervals2 = ps.intervals2, lambda2 = lambda2, ridge.adj = ridge.adj, wts = wts) return(llist) } "pspline.fitter"<- function(family, link, n.col, m.binomial, r.gamma, y, b, p, p.ridge, nix, nix.ridge, ridge.adj, wts, ...) { coef.est <- rep(1, ncol(b)) if(family == "binomial") { mu <- (y + 0.5 * m.binomial)/2 } if(family == "Gamma" || family == "poisson") { mu <- (y + 3) } if(family == "gaussian") { mu <- rep(mean(y), length(y)) } it <- 0 repeat { if(it == 0) { if(link == "identity") { eta <- mu } if(link == "log") { eta <- log(mu) } if(link == "sqrt") { eta <- sqrt(mu) } if(link == "logit") { eta <- log(mu/(m.binomial - mu)) } if(link == "recipical") { eta <- 1/mu } if(link == "probit") { eta <- qnorm(mu/m.binomial) } if(link == "cloglog") { eta <- log( - log(1 - mu/m.binomial)) } if(link == "loglog") { eta <- - log( - log(mu/m.binomial)) } } it <- it + 1 if(it > 25) break if(link == "identity") { mu <- eta h.prime <- 1 } if(link == "log") { mu <- exp(eta) h.prime <- mu } if(link == "sqrt") { mu <- eta^2 h.prime <- 2 * eta } if(link == "logit") { mu <- m.binomial/(1 + exp( - eta)) h.prime <- mu * (1 - mu/m.binomial) } if(link == "recipical") { mu <- 1/eta h.prime <- - (mu^2) } if(link == "probit") { mu <- m.binomial * pnorm(eta) h.prime <- m.binomial * dnorm(eta) } if(link == "cloglog") { mu <- m.binomial * (1 - exp( - exp(eta))) h.prime <- (m.binomial) * exp(eta) * exp( - exp(eta)) } if(link == "loglog") { mu <- m.binomial * exp( - exp( - eta)) h.prime <- m.binomial * exp( - eta) * exp( - exp( - eta )) } if(family == "gaussian") { w <- rep(1, length(y)) } if(family == "poisson") { w <- h.prime^2/mu } if(family == "binomial") { w <- h.prime^2/(mu * (1 - mu/m.binomial)) } if(family == "Gamma") { w <- (r.gamma * h.prime^2)/mu^2 } u <- (y - mu)/h.prime + eta if(ridge.adj > 0) { f <- lsfit(rbind(b, p, p.ridge), c(u, nix, nix.ridge), wt = c(wts, nix + 1, nix.ridge + 1) * c(w, (nix + 1), (nix.ridge + 1)), intercept = F) } if(ridge.adj == 0) { f <- lsfit(rbind(b, p), c(u, nix), wt = c(wts, nix + 1) * c(w, (nix + 1)), intercept = F) } coef.old <- coef.est coef.est <- as.vector(f$coef) d.coef <- max(abs((coef.est - coef.old)/coef.old)) if(d.coef < 1e-008) break print(c(it, d.coef)) eta <- b %*% coef.est } if(it > 24) { warning(paste("parameter estimates did NOT converge in 25 iterations" )) } llist <- list(coef = coef.est, mu = mu, f = f, w = w * wts) return(llist) }