# R code for CSPC tutorial # install.packages(c("support.CEs","censReg","lmtest","sandwich"), dependencies=TRUE) # Remove the leading hash symbol to install the packages used in this code # Generate a simple fractional factorial experimental design in R # U0=initial utility; LE=initial life expectancy; U1=final utility; nPats=patients treated; LYg=life years gained per patient library(support.CEs) CSPC.design <- Lma.design(attribute.names = list(Age=c("10", "40", "70"), U0=c("0.1", "0.5", "0.9"), LE=c("1m", "5yrs", "10yrs"), U1=c("0.1", "0.5", "0.9"), nPats=c("500", "2500", "5000"), LYg=c("1yr", "5yrs", "10yrs")), nalternatives=2, nblocks=3, row.renames=TRUE, seed = 123) CSPC.design # Design details questionnaire(CSPC.design, quote=FALSE) # Final questionnaire # Read sample dataset CSPC.Data <- read.csv("http://files.figshare.com/1411469/CSPC_data.csv", header=TRUE) # Read sample dataset str(CSPC.Data) # Summarize dataset # Histogram of budget allocations hist(CSPC.Data$BudgetB, breaks=101, freq=FALSE, col="grey", ylim=c(0,0.10), xlab="Budget to Program B", main="Histogram of Budget Allocations") # Calculate budget differences as allocation to Program B less allocation to Program A CSPC.Data$dBudget <- CSPC.Data$BudgetB - (100-CSPC.Data$BudgetB) # Histogram of budget differences hist(CSPC.Data$dBudget, breaks=201, freq=FALSE, col="grey", ylim=c(0,0.10), xlab="Budget to Program B", main="Histogram of Budget Differences") # Pooled tobit model library(censReg); library(lmtest); library(sandwich) pooled.tobit <- censReg(dBudget ~ dAge + dU0 + dLE + dU1 + I(dnPats/1000) + dLYg, left=-100, right=100, data=CSPC.Data) coeftest(pooled.tobit, vcov=sandwich) # Summary of results with robust SE # Pooled linear model library(lmtest); library(sandwich) pooled.linear <- glm(dBudget ~ dAge + dU0 + dLE + dU1 + I(dnPats/1000) + dLYg, data=CSPC.Data) coeftest(pooled.linear, vcov=sandwich) # Summary of results with robust SE # Calculate MRS, using LYg as numeraire, for discrete changes in attribute levels # Excludes LYg; by definition MRS(LYg:LYg)=1.0 MRS <- (coef(pooled.tobit)[2:6]/coef(pooled.tobit)["dLYg"])*c(-30,-0.4,-5,0.4,1000) MRS