Pilot Study 1
Analyze EQ5D scores
## source functions to analyze EQ5D scores
source("eq5d-mapping-functions.R")
## Loading required package: XLConnectJars
## XLConnect 0.2-11 by Mirai Solutions GmbH [aut],
## Martin Studer [cre],
## The Apache Software Foundation [ctb, cph] (Apache POI, Apache Commons
## Codec),
## Stephen Colebourne [ctb, cph] (Joda-Time Java library)
## http://www.mirai-solutions.com ,
## http://miraisolutions.wordpress.com
## define GBD disability weights
DW2010 <- c(0.003, 0.040, 0.127, 0.319, 0.756)[-4]
DW2013 <- c(0.005, 0.037, 0.128, NA, 0.778)[-4]
## scatterplot
par(mar = c(4, 4, 1, 1))
plot(DW2010, DW2013)
## load data
wb <- loadWorkbook("data-pilot-1.xlsx")
db <- readWorksheet(wb, sheet = 1)
## calculate utility for each observation
id <- seq(nrow(db))
utilities <- t(sapply(id, analyze_p1))
## remove health state 4 (epilepsy)
## and give this a new name to avoid confusion
utilities4 <- utilities[, -4]
## calculate mean utility per health state
## note: this sorts utilities by health state
utilities_mean <- apply(utilities4, 2, mean)
utilities_median <- apply(utilities4, 2, median)
utilities_sd <- apply(utilities4, 2, sd)
Plot utilities
## boxplot of utilities for each health state
df <-
data.frame(utility = c(utilities4),
health.state = factor(rep(1:4, each = nrow(utilities4))))
ggplot(df, aes(x = health.state, y = utility)) +
geom_boxplot() +
theme_bw()
## scatter plot of mean utility vs GBD DW
par(mfrow = c(1, 2), mar = c(4, 4, 1, 1) + .5)
plot_utilities1(utilities4, DW2010, ylab = "GBD 2010 DW")
plot_utilities1(utilities4, DW2013, ylab = "GBD 2013 DW")
Fit loess and predict DW from EQ5D utility - GBD 2010
## fit EQ5D vs GBD 2010 DW
fit_2010 <-
loess(logit(c(0.0001, DW2010, 0.9999)) ~
c(0.9999, utilities_mean, 0.0001))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## function to predict GBD 2010 DWs
predict_2010 <-
function(utility) {
expit(predict(fit_2010, utility))
}
## re-predict GBD 2010 DWs from utilities
DW2010_pred <- predict_2010(utilities_mean)
## compare utilities, predicted and true GBD 2010 DWs
round(
cbind(utility = utilities_mean,
"1-utility" = 1 - utilities_mean,
DW2010_pred,
DW2010,
pred.vs.true = DW2010_pred - DW2010),
3)
| 0.667 |
0.333 |
0.003 |
0.003 |
0 |
| 0.530 |
0.470 |
0.040 |
0.040 |
0 |
| 0.474 |
0.526 |
0.127 |
0.127 |
0 |
| 0.289 |
0.711 |
0.756 |
0.756 |
0 |
par(mar = c(4, 4, 1, 1))
plot(DW2010_pred, DW2010); abline(0, 1)
cor.test(DW2010_pred, DW2010, method = "spearman")
##
## Spearman's rank correlation rho
##
## data: DW2010_pred and DW2010
## S = 0, p-value = 0.08333
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 1
summary(lm(DW2010_pred ~ DW2010))
## Warning in summary.lm(lm(DW2010_pred ~ DW2010)): essentially perfect fit:
## summary may be unreliable
##
## Call:
## lm(formula = DW2010_pred ~ DW2010)
##
## Residuals:
## 1 2 3 4
## 8.946e-17 -7.701e-17 -1.944e-17 6.985e-18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.551e-17 5.312e-17 -1.045e+00 0.406
## DW2010 1.000e+00 1.384e-16 7.226e+15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.474e-17 on 2 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 5.221e+31 on 1 and 2 DF, p-value: < 2.2e-16
Fit loess and predict DW from EQ5D utility - GBD 2013
## fit EQ5D vs GBD 2013 DW
fit_2013 <-
loess(logit(c(0.0001, DW2013, 0.9999)) ~
c(0.9999, utilities_mean, 0.0001))
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : Chernobyl! trL>n 6
## function to predict GBD 2013 DWs
predict_2013 <-
function(utility) {
expit(predict(fit_2013, utility))
}
## re-predict GBD 2013 DWs from utilities
DW2013_pred <- predict_2013(utilities_mean)
## compare utilities, predicted and true GBD 2013 DWs
round(
cbind(utility = utilities_mean,
"1-utility" = 1 - utilities_mean,
DW2013_pred,
DW2013,
pred.vs.true = DW2013_pred - DW2013),
3)
| 0.667 |
0.333 |
0.005 |
0.005 |
0 |
| 0.530 |
0.470 |
0.037 |
0.037 |
0 |
| 0.474 |
0.526 |
0.128 |
0.128 |
0 |
| 0.289 |
0.711 |
0.778 |
0.778 |
0 |
par(mar = c(4, 4, 1, 1))
plot(DW2013_pred, DW2013); abline(0, 1)
cor.test(DW2013_pred, DW2013, method = "spearman")
##
## Spearman's rank correlation rho
##
## data: DW2013_pred and DW2013
## S = 0, p-value = 0.08333
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 1
summary(lm(DW2013_pred ~ DW2013))
## Warning in summary.lm(lm(DW2013_pred ~ DW2013)): essentially perfect fit:
## summary may be unreliable
##
## Call:
## lm(formula = DW2013_pred ~ DW2013)
##
## Residuals:
## 1 2 3 4
## -1.372e-17 -5.967e-18 2.312e-17 -3.432e-18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.551e-17 1.227e-17 4.523e+00 0.0456 *
## DW2013 1.000e+00 3.109e-17 3.216e+16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.963e-17 on 2 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.034e+33 on 1 and 2 DF, p-value: < 2.2e-16
R session info
sessionInfo()
## R version 3.2.3 (2015-12-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
##
## locale:
## [1] LC_COLLATE=Dutch_Belgium.1252 LC_CTYPE=Dutch_Belgium.1252
## [3] LC_MONETARY=Dutch_Belgium.1252 LC_NUMERIC=C
## [5] LC_TIME=Dutch_Belgium.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_2.1.0 XLConnect_0.2-11 XLConnectJars_0.2-9
## [4] printr_0.0.5 knitr_1.12.3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.3 digest_0.6.9 plyr_1.8.3 grid_3.2.3
## [5] gtable_0.2.0 formatR_1.3 magrittr_1.5 scales_0.4.0
## [9] evaluate_0.8.3 highr_0.5.1 stringi_1.0-1 rmarkdown_0.9.2
## [13] labeling_0.3 tools_3.2.3 stringr_1.0.0 munsell_0.4.3
## [17] yaml_2.1.13 colorspace_1.2-6 rJava_0.9-8 htmltools_0.3