# -------------------------------------------------------------------------- # PSS.r: Population Smoothing Splines method (Russu et al., IEEE TBME, 2011) # Copyright (C) 2011 Alberto Russu # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # -------------------------------------------------------------------------- PSS <- function(study.data, initial.values) { DOSE <- study.data$DOSE RESPONSE <- study.data$RESPONSE n.tot <- study.data$n.tot n.subjects <- study.data$n.subjects n.j <- study.data$n.j DOSE <- log(DOSE) RESPONSE <- log(RESPONSE) d.min <- min(DOSE) - log(2) DOSE <- DOSE - d.min min.dose <- min(DOSE) max.dose <- max(DOSE) ### HYPERPARAMETER ESTIMATION ### f <- function(x) { lambda2.bar <- 10^x[1] lambda2.tilde <- 10^x[2] sigma2 <- 10^x[3] sigma2.0 <- 10^x[4] t.grid <- DOSE R.bar <- matrix(NA, length(t.grid), length(t.grid)) for (i in 1 : length(t.grid)) for (j in i : length(t.grid)) { R.bar[i, j] <- r.bar.fun(t.grid[i], t.grid[j], lambda2.bar) R.bar[j, i] <- R.bar[i, j] } R.tilde <- matrix(0, length(DOSE), length(DOSE)) ptr <- 1 for (i in 1 : n.subjects) { t.grid <- DOSE[ptr:(ptr + n.j[i] - 1)] R <- matrix(NA, length(t.grid), length(t.grid)) for (i in 1 : length(t.grid)) for (j in i : length(t.grid)) { R[i, j] <- r.tilde.fun(t.grid[i], t.grid[j], lambda2.tilde, sigma2.0) R[j, i] <- R[i, j] } R.tilde[(ptr):(ptr+length(t.grid)-1), (ptr):(ptr+length(t.grid)-1)] <- R ptr <- ptr + length(t.grid) } phi <- cbind(matrix(1, length(DOSE), 1), DOSE) sigma.v <- sigma2 * diag(length(RESPONSE)) M <- R.bar + R.tilde + sigma.v inv.M <- solve(M) N <- solve(t(phi) %*% inv.M %*% phi) gamma.t.gamma <- RESPONSE %*% (inv.M - inv.M %*% phi %*% N %*% t(phi) %*% inv.M) %*% RESPONSE J <- log(det(M)) - log(det(N)) + gamma.t.gamma return(J) } ret <- optim(log10(initial.values), f) print(str(ret)) lambda2.bar <- 10^ret$par[1] lambda2.tilde <- 10^ret$par[2] sigma2 <- 10^ret$par[3] sigma2.0 <- 10^ret$par[4] t.grid <- DOSE R.bar <- matrix(NA, length(t.grid), length(t.grid)) for (i in 1 : length(t.grid)) for (j in i : length(t.grid)) { R.bar[i, j] <- r.bar.fun(t.grid[i], t.grid[j], lambda2.bar) R.bar[j, i] <- R.bar[i, j] } R.tilde <- matrix(0, length(DOSE), length(DOSE)) ptr <- 1 for (i in 1 : n.subjects) { t.grid <- DOSE[ptr:(ptr + n.j[i] - 1)] R <- matrix(NA, length(t.grid), length(t.grid)) for (i in 1 : length(t.grid)) for (j in i : length(t.grid)) { R[i, j] <- r.tilde.fun(t.grid[i], t.grid[j], lambda2.tilde, sigma2.0) R[j, i] <- R[i, j] } R.tilde[(ptr):(ptr+length(t.grid)-1), (ptr):(ptr+length(t.grid)-1)] <- R ptr <- ptr + length(t.grid) } phi <- cbind(matrix(1, length(DOSE), 1), DOSE) sigma.v <- sigma2 * diag(length(RESPONSE)) M <- R.bar + R.tilde + sigma.v inv.M <- solve(M) N <- solve(t(phi) %*% inv.M %*% phi) d <- N %*% t(phi) %*% inv.M %*% RESPONSE w <- inv.M %*% (RESPONSE - phi %*% d) ### PREDICTION (reconstruction of individual curves) ### pred.grid <- rep(seq(min.dose, max.dose, by = (max.dose - min.dose)/20), n.subjects) n.j.pred <- rep(length(pred.grid)/n.subjects, n.subjects) t.grid <- DOSE R.bar.pred <- matrix(0, length(pred.grid), length(t.grid)) u <- 1 for (k in 1 : n.subjects) { sub.grid <- pred.grid[u:(u + n.j.pred[k] - 1)] R <- matrix(NA, length(sub.grid), length(t.grid)) for (j in 1 : length(t.grid)) for (i in 1 : length(sub.grid)) R[i, j] <- r.bar.fun(sub.grid[i], t.grid[j], lambda2.bar) R.bar.pred[(u):(u+length(sub.grid)-1), ] <- R u <- u + n.j.pred[k] } R.tilde.pred <- matrix(0, length(pred.grid), length(w)) ptr <- 1 u <- 1 for (k in 1 : n.subjects) { t.grid <- DOSE[ptr:(ptr + n.j[k] - 1)] sub.grid <- pred.grid[u:(u + n.j.pred[k] - 1)] R <- matrix(NA, length(sub.grid), length(t.grid)) for (j in 1 : length(t.grid)) for (i in 1 : length(sub.grid)) R[i, j] <- r.tilde.fun(sub.grid[i], t.grid[j], lambda2.tilde, sigma2.0) R.tilde.pred[(u):(u+length(sub.grid)-1), (ptr):(ptr+length(t.grid)-1)] <- R ptr <- ptr + length(t.grid) u <- u + n.j.pred[k] } phi.pred <- cbind(matrix(1, length(pred.grid), 1), pred.grid) z.bar <- R.bar.pred %*% w + phi.pred %*% d z.j <- z.bar + R.tilde.pred %*% w pred.grid <- pred.grid + d.min pred.grid <- exp(pred.grid) z.j <- exp(z.j) ret <- list() ret$pred.grid <- pred.grid ret$z.j <- z.j ret$n.j.pred <- n.j.pred return(ret) } r.bar.fun <- function(t1, t2, lambda2.bar) { if (t1 <= t2) { res <- lambda2.bar * (t1^2/2) * (t2 - t1/3) } else { res <- lambda2.bar * (t2^2/2) * (t1 - t2/3) } return(res) } r.tilde.fun <- function(t1, t2, lambda2.tilde, sigma2.0) { if (t1 <= t2) { res <- sigma2.0 + lambda2.tilde * (t1^2/2) * (t2 - t1/3) } else { res <- sigma2.0 + lambda2.tilde * (t2^2/2) * (t1 - t2/3) } return(res) } tab <- read.csv2("example.csv") study.data <- list() study.data$DOSE <- tab[,1] study.data$RESPONSE <- tab[,3] study.data$n.tot <- dim(tab)[1] study.data$subjects <- unique(tab[,2]) study.data$n.subjects <- length(study.data$subjects) study.data$n.j <- array(NA, study.data$n.subjects) for (i in 1 : study.data$n.subjects) study.data$n.j[i] <- length(tab[tab[, 2] == study.data$subjects[i], 1]) initial.values <- c(10, 10, 1, 1) ret <- PSS(study.data, initial.values) windows() par(mfrow = c(3, 4)) u <- 1 for (i in 1 : study.data$n.subjects) { plot(tab[tab[, 2] == study.data$subjects[i],1], tab[tab[, 2] == study.data$subjects[i],3], type = "p", xlim = c(0, 1.1*max(tab[,1])), ylim = c(0, 1.1*max(tab[,3])), col = "black", pch = 19, xlab = "Dose", ylab = "Response") lines(ret$pred.grid[u : (u + ret$n.j.pred[i] - 1)], ret$z.j[u : (u + ret$n.j.pred[i] - 1)], col = "red", lwd = 2) u <- u + ret$n.j.pred[i] }