# function that calculates a model for predicting the outcome
# depending from a Dvolume, used by fit.Dvolume
glm.Dvolume <- function(dvh.matrix, InputVolume, formula) {
DVolume<<-Dvolume(dvh.matrix=dvh.matrix, Volume=InputVolume)
f <- update.formula(formula, ~ . + DVolume)
if (type.binary=="binomial") out.model <<- glm(formula=f, family=binomial(link = "logit"))
if (type.binary=="gaussian") out.model <<- glm(formula=f, family=gaussian(link = "identity"))
if (type.binary=="poisson") out.model <<- glm(formula=f, family=poisson(link = "log"))
output<-out.model$aic
return(output)
}
# function that finds Dvolume and other covariates affecting a given outcome
fit.Dvolume <- function(dvh.matrix, formula, num.steps.volume=100, type.binary=c("binomial", "gaussian", "poisson"),
range.min=NULL, range.max=NULL) {
require(data.table)
type.binary<<-match.arg(type.binary)
formula <- as.formula(formula)
# define the volume range for Dvolume calculation
# maximum Volume is the minimum volume of the maxima in the dvh matrix
if (is.null(range.max)) range.max<-min(apply(X=dvh.matrix[,2:ncol(dvh.matrix)], MARGIN=2, FUN=max))
# fix the number of the steps for calculating the dose range
if (is.null(range.min)) range.min<-min(dvh.matrix[,2:ncol(dvh.matrix)])
# create sequence of volumes, find the minimum of the maxima in the volume values
seq.volume <- seq(from=range.min, to=range.max, length.out=num.steps.volume)
# delete the zero if present in the sequence and steps forward to the following Dvolume
if (range.min==0) seq.volume<-seq.volume[2:length(seq.volume)]
# create matrix of Dvolumes
seq.Dvolume <- c()
for (n in 1:length(seq.volume)) {
seq.Dvolume <- cbind(seq.Dvolume, Dvolume(dvh.matrix=dvh.matrix, Volume=seq.volume[n]))
}
# create sequence of glm models for getting the minimum aic and starting optimization
temp.formula <- update.formula(formula, ~ . + D.volume)
result.matrix <- c()
for (n in 1:length(seq.volume)) {
D.volume <<- seq.Dvolume[,n]
if (type.binary=="binomial") family=binomial(link = "logit")
if (type.binary=="gaussian") family=gaussian(link = "identity")
if (type.binary=="poisson") family=poisson(link = "log")
temp.model <- glm(temp.formula, family)
temp.p <- c()
# create the temp vector of p values
if (type.binary=="gaussian") for (m in 1:length(summary(temp.model)$coefficients[, "Pr(>|t|)"])) {
temp.p <- c(temp.p, summary(temp.model)$coefficients[, "Pr(>|t|)"][m])
} else for (m in 1:length(summary(temp.model)$coefficients[, "Pr(>|z|)"])) {
temp.p <- c(temp.p, summary(temp.model)$coefficients[, "Pr(>|z|)"][m])
}
# add row in the matrix: Dvolume, aic, list of P values
temp.result<-c(seq.volume[n], temp.model$aic, temp.p)
result.matrix<-rbind(result.matrix, temp.result)
}
# typecast result.matrix
result.matrix<-as.matrix(result.matrix)
rm(D.volume, envir=.GlobalEnv)
# set the start Dvolume at minimum aic found from series of aic
start.Dvolume <- result.matrix [(which.min(x=result.matrix[,2])),1]
# finds the minimum value of aic for best modeling of Dvolume
out<-nlminb(start=start.Dvolume, objective=glm.Dvolume, dvh.matrix=dvh.matrix, formula=formula,
lower=range.min, upper=range.max)
# unname the matrix
result.matrix<-unname(result.matrix)
cnames<-c("Volume", "AIC")
temp.p <- c()
# create the temp vector of p values
if (type.binary=="gaussian") for (m in 1:length(summary(out.model)$coefficients[, "Pr(>|t|)"])) {
temp.p <- c(temp.p, summary(out.model)$coefficients[, "Pr(>|t|)"][m])
# add the column names
cnames <- c(cnames, paste("P Value ", names(out.model$coefficients[m])))
} else for (m in 1:length(summary(out.model)$coefficients[, "Pr(>|z|)"])) {
temp.p <- c(temp.p, summary(out.model)$coefficients[, "Pr(>|z|)"][m])
# add the column names
cnames <- c(cnames, paste("P Value ", names(out.model$coefficients[m])))
}
colnames(result.matrix) <- cnames
result.matrix <- rbind(result.matrix, c(out$par, out.model$aic, temp.p))
result.matrix <- data.table(result.matrix, key="Volume")
# new typecast after ordering table
result.matrix <- as.matrix(result.matrix)
output<-(list("model"=out.model, "Dvolume"=out$par, "dvhmatrix"=dvh.matrix,
"resultmatrix"=result.matrix))
class(output)<-"Dvolume.fit"
# clean environment from global variables
rm(type.binary, out.model, envir=.GlobalEnv)
return(output)
}
# output of fit.Dvolume
print.Dvolume.fit <- function (fittedmodel) {
cat("Dvolume = ", fittedmodel$Dvolume, " cc", "\n")
print(fittedmodel$model)
}
# summary of fit.Dvolume
summary.Dvolume.fit <- function (fittedmodel) {
cat("Dvolume = ", fittedmodel$Dvolume, " cc", "\n")
summary(fittedmodel$model)
}
# function for plotting the series of P values of covariates
# and AIC related to Dvolumes
plot.Dvolume.fit <- function (fittedmodel, plot.legend=FALSE) {
# draw the frame
par(mar=c(3, 14, 0, 0) + 0.1)
# plot the AIC curve
plot(fittedmodel$resultmatrix[,1], fittedmodel$resultmatrix[,2],
axes=F, ylim=c(min(fittedmodel$resultmatrix[,2]),max(fittedmodel$resultmatrix[,2])),
xlab="", ylab="", type="l",col="black", main="",xlim=c(min(fittedmodel$resultmatrix[,1]),max(fittedmodel$resultmatrix[,1])))
axis(2, ylim=c(min(fittedmodel$resultmatrix[,2],max(fittedmodel$resultmatrix[,2]))), col="black", lwd=1)
mtext(2,text=colnames(fittedmodel$resultmatrix)[2] ,line=2)
abline(v=fittedmodel$Dvolume, lty=1, col="red")
text(x=fittedmodel$Dvolume, y=max(fittedmodel$resultmatrix[,2]),
labels=paste(round(fittedmodel$Dvolume, digits=2), " cc"), pos=4)
legend.labels<-"AIC"
legend.lty <- 1
# draw the remaining covariates P values
linepos<-3.5
for (n in 3:ncol(fittedmodel$resultmatrix)){
par(new=T)
plot(fittedmodel$resultmatrix[,1], fittedmodel$resultmatrix[,n], lty=n-1, log="y",
axes=F, ylim=c(min(fittedmodel$resultmatrix[,n]),max(fittedmodel$resultmatrix[,n])),
xlab="", ylab="", type="l",col="black", main="",xlim=c(min(fittedmodel$resultmatrix[,1]),max(fittedmodel$resultmatrix[,1])))
axis(2, ylim=c(min(fittedmodel$resultmatrix[,n]),max(fittedmodel$resultmatrix[,n])),lwd=1,line=linepos)
mtext(2,text=colnames(fittedmodel$resultmatrix)[n] ,line=2+linepos)
legend.labels<-c(legend.labels, colnames(fittedmodel$resultmatrix)[n])
legend.lty <- c(legend.lty, n - 1)
linepos<-linepos + 3.5
}
# plot the bottom axis with text
axis(1, pretty(range(fittedmodel$resultmatrix[,1]),10), lwd=1)
mtext("D-Volume Value (cc)", side=1, col="black", line=2)
# plot the legend if required
if (plot.legend==TRUE) {
legend(x=min(fittedmodel$resultmatrix[,1]),y=median(fittedmodel$resultmatrix[,n]),
legend=legend.labels, lty=legend.lty)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.