Pilot Study 2
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.004, 0.005, 0.008, 0.009, 0.016,
0.021, 0.029, 0.033, 0.040, 0.051, 0.065,
0.087, 0.114, 0.127, 0.225, 0.263, 0.293,
0.352, 0.390, 0.433, 0.492, 0.519, 0.547,
0.625, 0.707, 0.756)
DW2013 <-
c(0.005, 0.003, 0.004, 0.006, 0.015, 0.014,
0.019, 0.027, 0.026, 0.037, 0.088, 0.100,
0.050, 0.117, 0.128, 0.231, 0.267, 0.295,
0.369, 0.279, 0.441, 0.501, 0.569, 0.582,
0.637, 0.719, 0.778)
## scatterplot
par(mar = c(4, 4, 1, 1))
plot(DW2010, DW2013)
## load data
wb <- loadWorkbook("data-pilot-2.xlsx")
db <- readWorksheet(wb, sheet = 1)
## calculate utility for each observation
id <- seq(nrow(db))
utilities <- sapply(id, analyze_p2)
## calculate mean utility per health state
## note: this sorts utilities by health state
utilities_mean <- tapply(utilities, db$HealthstatID, mean)
utilities_median <- tapply(utilities, db$HealthstatID, median)
utilities_sd <- tapply(utilities, db$HealthstatID, sd)
Plot utilities
## boxplot of utilities for each health state
df <-
data.frame(utility = utilities,
health.state = factor(db$HealthstatID))
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_utilities2(utilities, db$HealthstatID, DW2010, ylab = "GBD 2010 DW")
plot_utilities2(utilities, db$HealthstatID, 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))
## 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)
| 1 |
0.694 |
0.306 |
0.010 |
0.003 |
0.007 |
| 2 |
0.826 |
0.174 |
0.002 |
0.004 |
-0.002 |
| 5 |
0.813 |
0.187 |
0.003 |
0.005 |
-0.002 |
| 8 |
0.650 |
0.350 |
0.017 |
0.008 |
0.009 |
| 11 |
0.809 |
0.191 |
0.003 |
0.009 |
-0.006 |
| 22 |
0.644 |
0.356 |
0.018 |
0.016 |
0.002 |
| 29 |
0.550 |
0.450 |
0.091 |
0.021 |
0.070 |
| 37 |
0.639 |
0.361 |
0.019 |
0.029 |
-0.010 |
| 47 |
0.615 |
0.385 |
0.027 |
0.033 |
-0.006 |
| 55 |
0.350 |
0.650 |
0.278 |
0.040 |
0.238 |
| 60 |
0.304 |
0.696 |
0.316 |
0.051 |
0.265 |
| 77 |
0.553 |
0.447 |
0.088 |
0.065 |
0.023 |
| 96 |
0.389 |
0.611 |
0.252 |
0.087 |
0.165 |
| 107 |
0.474 |
0.526 |
0.191 |
0.114 |
0.077 |
| 111 |
0.486 |
0.514 |
0.177 |
0.127 |
0.050 |
| 149 |
0.455 |
0.545 |
0.215 |
0.225 |
-0.010 |
| 155 |
0.391 |
0.609 |
0.250 |
0.263 |
-0.013 |
| 161 |
0.193 |
0.807 |
0.538 |
0.293 |
0.245 |
| 174 |
0.331 |
0.669 |
0.294 |
0.352 |
-0.058 |
| 186 |
0.106 |
0.894 |
0.934 |
0.390 |
0.544 |
| 193 |
0.548 |
0.452 |
0.094 |
0.433 |
-0.339 |
| 200 |
0.449 |
0.551 |
0.221 |
0.492 |
-0.271 |
| 203 |
0.317 |
0.683 |
0.306 |
0.519 |
-0.213 |
| 206 |
0.360 |
0.640 |
0.272 |
0.547 |
-0.275 |
| 214 |
0.200 |
0.800 |
0.504 |
0.625 |
-0.121 |
| 219 |
0.144 |
0.856 |
0.803 |
0.707 |
0.096 |
| 220 |
0.341 |
0.659 |
0.285 |
0.756 |
-0.471 |
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 = 630, p-value = 2.041e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8076923
summary(lm(DW2010_pred ~ DW2010))
##
## Call:
## lm(formula = DW2010_pred ~ DW2010)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27163 -0.08724 -0.03970 0.03099 0.60469
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.08701 0.04939 1.762 0.090375 .
## DW2010 0.62125 0.14843 4.185 0.000307 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1854 on 25 degrees of freedom
## Multiple R-squared: 0.412, Adjusted R-squared: 0.3885
## F-statistic: 17.52 on 1 and 25 DF, p-value: 0.0003074
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))
## 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)
| 1 |
0.694 |
0.306 |
0.011 |
0.005 |
0.006 |
| 2 |
0.826 |
0.174 |
0.002 |
0.003 |
-0.001 |
| 5 |
0.813 |
0.187 |
0.003 |
0.004 |
-0.001 |
| 8 |
0.650 |
0.350 |
0.017 |
0.006 |
0.011 |
| 11 |
0.809 |
0.191 |
0.003 |
0.015 |
-0.012 |
| 22 |
0.644 |
0.356 |
0.018 |
0.014 |
0.004 |
| 29 |
0.550 |
0.450 |
0.092 |
0.019 |
0.073 |
| 37 |
0.639 |
0.361 |
0.019 |
0.027 |
-0.008 |
| 47 |
0.615 |
0.385 |
0.028 |
0.026 |
0.002 |
| 55 |
0.350 |
0.650 |
0.290 |
0.037 |
0.253 |
| 60 |
0.304 |
0.696 |
0.331 |
0.088 |
0.243 |
| 77 |
0.553 |
0.447 |
0.088 |
0.100 |
-0.012 |
| 96 |
0.389 |
0.611 |
0.254 |
0.050 |
0.204 |
| 107 |
0.474 |
0.526 |
0.189 |
0.117 |
0.072 |
| 111 |
0.486 |
0.514 |
0.175 |
0.128 |
0.047 |
| 149 |
0.455 |
0.545 |
0.213 |
0.231 |
-0.018 |
| 155 |
0.391 |
0.609 |
0.251 |
0.267 |
-0.016 |
| 161 |
0.193 |
0.807 |
0.545 |
0.295 |
0.250 |
| 174 |
0.331 |
0.669 |
0.309 |
0.369 |
-0.060 |
| 186 |
0.106 |
0.894 |
0.931 |
0.279 |
0.652 |
| 193 |
0.548 |
0.452 |
0.095 |
0.441 |
-0.346 |
| 200 |
0.449 |
0.551 |
0.219 |
0.501 |
-0.282 |
| 203 |
0.317 |
0.683 |
0.322 |
0.569 |
-0.247 |
| 206 |
0.360 |
0.640 |
0.281 |
0.582 |
-0.301 |
| 214 |
0.200 |
0.800 |
0.512 |
0.637 |
-0.125 |
| 219 |
0.144 |
0.856 |
0.801 |
0.719 |
0.082 |
| 220 |
0.341 |
0.659 |
0.299 |
0.778 |
-0.479 |
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 = 648, p-value = 2.222e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8021978
summary(lm(DW2013_pred ~ DW2013))
##
## Call:
## lm(formula = DW2013_pred ~ DW2013)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25506 -0.10049 -0.06975 0.03662 0.67219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.10210 0.05184 1.970 0.06007 .
## DW2013 0.56119 0.15261 3.677 0.00113 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1956 on 25 degrees of freedom
## Multiple R-squared: 0.351, Adjusted R-squared: 0.3251
## F-statistic: 13.52 on 1 and 25 DF, p-value: 0.001129
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