# Moved this outside the document for easy of reading # I often have those sections in here source("Setup_and_munge.R")
info <- sessionInfo() r_ver <- paste(info$R.version$major, info$R.version$minor, sep=".")
All analyses were performed using R (ver. r r_ver
)[R Core Team, 2013] and packages rms (ver. r info$otherPkgs$rms$Version
) [F. Harrell, 2014] for analysis, Gmisc for plot and table output (ver. r info$otherPkgs$Gmisc$Version
), and knitr (ver r info$otherPkgs$knitr$Version
) [Xie, 2013] for reproducible research.
We found r nrow(melanoma)
patients with malignant melanoma between the years r paste(range(melanoma$year), collapse=" and ")
. Patients were followed until the end of 1977, the median follow-up time was r sprintf("%.1f", median(melanoma$time_years))
years (range r paste(sprintf("%.1f", range(melanoma$time_years)), collapse=" to ")
years). Males were more common than females and had also a higher mortality rate.
table_data <- list() getT1Stat <- function(varname, digits=0){ getDescriptionStatsBy(melanoma[, varname], melanoma$status, add_total_col=TRUE, show_all_values=TRUE, hrzl_prop=TRUE, statistics=FALSE, html=TRUE, digits=digits) } # Get the basic stats table_data[["Sex"]] <- getT1Stat("sex") table_data[["Age<sup>†</sup>"]] <- getT1Stat("age") table_data[["Ulceration"]] <- getT1Stat("ulcer") table_data[["Thickness<sup>‡</sup>"]] <- getT1Stat("thickness", digits=1) # Now merge everything into a matrix # and create the rgroup & n.rgroup variabels rgroup <- c() n.rgroup <- c() output_data <- NULL for (varlabel in names(table_data)){ output_data <- rbind(output_data, table_data[[varlabel]]) rgroup <- c(rgroup, varlabel) n.rgroup <- c(n.rgroup, nrow(table_data[[varlabel]])) } # Add a column spanner for the death columns cgroup <- c("", "Death") n.cgroup <- c(2, 2) colnames(output_data) <- gsub("[ ]*death", "", colnames(output_data)) htmlTable(output_data, align="rrrr", rgroup=rgroup, n.rgroup=n.rgroup, rgroupCSSseparator="", cgroup = cgroup, n.cgroup = n.cgroup, rowlabel="", caption="Baseline characteristics", tfoot="<sup>†</sup> Age at the time of surgery. <br/><sup>‡</sup> Tumour thicknes, also known as Breslow thickness, measured in mm.", ctable=TRUE)
label(melanoma$sex) <- "Sex" label(melanoma$age) <- "Age" label(melanoma$ulcer) <- "Ulceration" label(melanoma$thickness) <- "Breslow thickness" # Setup needed for the rms coxph wrapper ddist <- datadist(melanoma) options(datadist = "ddist") # Do the cox regression model # for melanoma specific death msurv <- Surv(melanoma$time_years, melanoma$status=="Melanoma death") fit <- cph(msurv ~ sex + age + ulcer + thickness, data=melanoma) # Print the model printCrudeAndAdjustedModel(fit, desc_digits=0, caption="Adjusted and unadjusted estimates for melanoma specific death.", desc_column=TRUE, add_references=TRUE, ctable=TRUE) pvalues <- 1 - pchisq(coef(fit)^2/diag(vcov(fit)), df=1)
After adjusting for the three variables, age, sex, tumor thickness and ulceration, only the latter two remained significant (p-value r pvalueFormatter(pvalues["ulcer=Present"], sig.limit=10^-3)
and r pvalueFormatter(pvalues["thickness"], sig.limit=10^-3)
), see table r as.numeric(options("table_counter"))-1
and figure I.
# I've adjusted the coefficient for age to be by forestplotRegrObj(update(fit, .~.-age+I(age/10)), order.regexps=c("Female", "age", "ulc", "thi"), box.default.size=.25, xlog=TRUE, new_page=TRUE, clip=c(.5, 6), rowname.fn=function(x){ if (grepl("Female", x)) return("Female") if (grepl("Present", x)) return("Ulceration") if (grepl("age", x)) return("Age/10 years") return(capitalize(x)) })
There was no strong indication for non-linearity for any of the continuous variables although the impact of thickness did seem to lessen above 4 mm, see figure II.
plotHR_cap = paste0("The adjusted and unadjusted restricted cubic spline", " for tumor thickness. Solid line and confidence interval", " indicate the adjusted line while the dashed is", " the unadjusted line. The grey area at ", " the bottom indicates the density.") # Generate adjusted and anuadjusted regression models rcs_fit <- update(fit, .~.-thickness+rcs(thickness, 3)) rcs_fit_ua <- update(fit, .~+rcs(thickness, 3)) # Make sure the axes stay at the exact intended points par(xaxs="i", yaxs="i") plotHR(list(rcs_fit, rcs_fit_ua), col.dens="#00000033", lty.term=c(1, 2), col.term=c("blue", "#444444"), col.se = c("#0000FF44", "grey"), polygon_ci=c(TRUE, FALSE), term="thickness", xlab="Thickness (mm)", ylim=c(.1, 4), xlim=c(min(melanoma$thickness), 4), plot.bty="l", y.ticks=c(.1, .25, .5, 1, 2, 4)) legend(x=.1, y=1.1, legend=c("Adjusted", "Unadjusted"), fill=c("blue", "grey"), bty="n")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.