Hippo therapy evaluation /assessment tool (HEAT)_Dr.M.Shotwell
D. Weissman-Miller, Sc.D.
November 6, 2014, modified January 6, 2015
Object of this analysis
HEAT is hippo therapy evaluation /assessment tool designed to measure ongoing progress and potentially outcomes for participants in hippo therapy in 4 domains:
- Static posture, 2) Dynamic motor behavior, 3) Sensory processing and 4) Psycho-social performance. HEAT supplements other measures in this field to provide support for assessing factors that may contribute to a client’s health.
The motivation for this statistical analysis is to develop criterion validity for the means of the HEATTOT variable with respect to Disability or Non-Disability in the multi-year assessment database. A logit (logistical) regression was chosen and the outcome variable re-coded as 0 to 1. The logistic regression coefficients give the change in the log odds of the outcome for a one unit increase in the predictor variable.
In addition, the HEATTOT was renamed HEATRANK (in means) to reflect the rank order of this variable. The mean of 74 for the HEAT variables equals the disabled participants. The mean of 80 for the HEAT variables equals the total mean for all participants. The mean of 88 for the HEAT variables equals the typical paticipants. The modified spreadsheet was used with only the variables of interest (Static Posture, Dynamic Motor Behavior, Sensory Processing and Psychsocial behavior.)
Prepare the R environment
library(knitr)
opts_chunk$set(echo=TRUE, results= 'hold')Load required libraries
library(data.table)
library(ggplot2)Load the data
HEAT <- read.table('HeatDataDwm_final.csv', header=TRUE, sep= ",")HEAT dataset (first 6 rows)
head(HEAT)## DISvsNON HEATSTAT HEATDYN HEATSENS HEATPSYCH HEATRANK
## 1 0 9 17 6 5 37
## 2 0 17 34 10 3 64
## 3 0 20 43 17 10 90
## 4 0 18 38 12 6 74
## 5 0 12 20 11 5 48
## 6 0 18 33 14 14 79Summary of logit regression
summary(HEAT)## DISvsNON HEATSTAT HEATDYN HEATSENS
## Min. :0.0000 Min. : 0.00 Min. : 0.00 Min. : 6.00
## 1st Qu.:0.0000 1st Qu.:16.25 1st Qu.:32.25 1st Qu.:15.25
## Median :0.0000 Median :17.00 Median :36.00 Median :19.00
## Mean :0.4429 Mean :16.49 Mean :34.39 Mean :17.26
## 3rd Qu.:1.0000 3rd Qu.:18.75 3rd Qu.:39.00 3rd Qu.:20.00
## Max. :1.0000 Max. :20.00 Max. :45.00 Max. :20.00
## HEATPSYCH HEATRANK
## Min. : 0.00 Min. : 6.00
## 1st Qu.:10.00 1st Qu.: 79.00
## Median :14.00 Median : 86.50
## Mean :11.89 Mean : 80.17
## 3rd Qu.:15.00 3rd Qu.: 88.75
## Max. :15.00 Max. :100.00This is a plot of ‘pairs’ in the reduced HEAT dataset
pairs(HEAT)
Here is a logit regression model of disability Vs. non-disability from HEAT predictors
fit <- glm(HEAT$DISvsNON ~ HEAT$HEATSTAT + HEAT$HEATDYN + HEAT$HEATSENS + HEAT$HEATPSYCH +
HEAT$HEATRANK, data=HEAT, family = "binomial")## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(fit)##
## Call:
## glm(formula = HEAT$DISvsNON ~ HEAT$HEATSTAT + HEAT$HEATDYN +
## HEAT$HEATSENS + HEAT$HEATPSYCH + HEAT$HEATRANK, family = "binomial",
## data = HEAT)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.07958 -0.01873 0.00000 0.24964 2.52740
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -46.2821 20.5117 -2.256 0.024 *
## HEAT$HEATSTAT -1.7430 1.2244 -1.424 0.155
## HEAT$HEATDYN 0.3072 1.0140 0.303 0.762
## HEAT$HEATSENS 1.7491 1.2721 1.375 0.169
## HEAT$HEATPSYCH 0.7895 1.0765 0.733 0.463
## HEAT$HEATRANK 0.2441 1.0092 0.242 0.809
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 96.124 on 69 degrees of freedom
## Residual deviance: 24.678 on 64 degrees of freedom
## AIC: 36.678
##
## Number of Fisher Scoring iterations: 9Now the confidence intervals are predicted for this model
confint(fit)## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) -99.4871840 -21.106140
## HEAT$HEATSTAT -6.3395573 -0.295648
## HEAT$HEATDYN -3.0721772 1.595216
## HEAT$HEATSENS -1.4928049 4.190024
## HEAT$HEATPSYCH -2.6569220 2.280009
## HEAT$HEATRANK -0.8002501 3.942860Confidence intervals using SE
confint.default(fit)## 2.5 % 97.5 %
## (Intercept) -86.4842970 -6.0799094
## HEAT$HEATSTAT -4.1428550 0.6568333
## HEAT$HEATDYN -1.6801836 2.2946671
## HEAT$HEATSENS -0.7441617 4.2423827
## HEAT$HEATPSYCH -1.3204521 2.8994877
## HEAT$HEATRANK -1.7338646 2.2220601The coefficients are exponentiated and interpreted as odds ratios
## odds ratios only
exp(coef(fit))## (Intercept) HEAT$HEATSTAT HEAT$HEATDYN HEAT$HEATSENS HEAT$HEATPSYCH
## 7.942148e-21 1.749927e-01 1.359670e+00 5.749486e+00 2.202334e+00
## HEAT$HEATRANK
## 1.276469e+00## odds ratios and 95% CI
exp(cbind(OR = coef(fit), confint(fit)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 7.942148e-21 6.212480e-44 6.818985e-10
## HEAT$HEATSTAT 1.749927e-01 1.765083e-03 7.440493e-01
## HEAT$HEATDYN 1.359670e+00 4.632020e-02 4.929393e+00
## HEAT$HEATSENS 5.749486e+00 2.247414e-01 6.602435e+01
## HEAT$HEATPSYCH 2.202334e+00 7.016385e-02 9.776773e+00
## HEAT$HEATRANK 1.276469e+00 4.492166e-01 5.156585e+01Calculate the mean for each HEAT variable when disability == 0
newdata0 <- subset(HEAT, DISvsNON == 0)
newdata1 <- with(newdata0, data.frame(HEATSTAT=mean(HEATSTAT), HEATDYN=mean(HEATDYN),
HEATRANK=as.numeric(74:88)))
newdata1## HEATSTAT HEATDYN HEATRANK
## 1 15.79487 32.41026 74
## 2 15.79487 32.41026 75
## 3 15.79487 32.41026 76
## 4 15.79487 32.41026 77
## 5 15.79487 32.41026 78
## 6 15.79487 32.41026 79
## 7 15.79487 32.41026 80
## 8 15.79487 32.41026 81
## 9 15.79487 32.41026 82
## 10 15.79487 32.41026 83
## 11 15.79487 32.41026 84
## 12 15.79487 32.41026 85
## 13 15.79487 32.41026 86
## 14 15.79487 32.41026 87
## 15 15.79487 32.41026 88Calculate the mean for each variable when disability ==1
newdata0 <- subset(HEAT, DISvsNON == 1)
newdata1 <- with(newdata0, data.frame(HEATSTAT=mean(HEATSTAT), HEATDYN=mean(HEATDYN),
HEATRANK=as.numeric(74:88)))
newdata1## HEATSTAT HEATDYN HEATRANK
## 1 17.35484 36.87097 74
## 2 17.35484 36.87097 75
## 3 17.35484 36.87097 76
## 4 17.35484 36.87097 77
## 5 17.35484 36.87097 78
## 6 17.35484 36.87097 79
## 7 17.35484 36.87097 80
## 8 17.35484 36.87097 81
## 9 17.35484 36.87097 82
## 10 17.35484 36.87097 83
## 11 17.35484 36.87097 84
## 12 17.35484 36.87097 85
## 13 17.35484 36.87097 86
## 14 17.35484 36.87097 87
## 15 17.35484 36.87097 88Code to generate predicted probabilities
newdata2 <- with(newdata1, data.frame(HEATRANK = seq(21,90, by= 1), 70)) Combine code
newdata3 <- cbind(newdata2, predict(fit, newdata = newdata2, type = "link",
se = TRUE))
newdata3 <- within(newdata3, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})View final dataset including predicted probabilities for HEATRANK
library(knitr)
kable(newdata3[,1:8], format= "markdown")##
##
## | HEATRANK| X70| fit| se.fit| residual.scale| UL| LL| PredictedProb|
## |--------:|---:|-----------:|----------:|--------------:|---------:|---------:|-------------:|
## | 21| 70| -33.2722222| 13.4327035| 1| 0.0009634| 0.0000000| 0.0000000|
## | 22| 70| -29.9851535| 9.5156029| 1| 0.0000120| 0.0000000| 0.0000000|
## | 23| 70| -8.3320706| 2.7287559| 1| 0.0481673| 0.0000011| 0.0002406|
## | 24| 70| -22.1914455| 7.1449815| 1| 0.0002781| 0.0000000| 0.0000000|
## | 25| 70| -26.1489016| 9.4178975| 1| 0.0004573| 0.0000000| 0.0000000|
## | 26| 70| -12.6928020| 4.9024568| 1| 0.0437726| 0.0000000| 0.0000031|
## | 27| 70| -2.1353367| 1.5696166| 1| 0.7193242| 0.0054224| 0.1057094|
## | 28| 70| -18.8664486| 7.0956404| 1| 0.0069711| 0.0000000| 0.0000000|
## | 29| 70| -2.6977630| 2.4915714| 1| 0.8989595| 0.0005097| 0.0631055|
## | 30| 70| 0.9775104| 0.8569072| 1| 0.9344441| 0.3313628| 0.7266140|
## | 31| 70| -3.3852113| 2.4737710| 1| 0.8120491| 0.0002655| 0.0327609|
## | 32| 70| -6.2358605| 2.1181900| 1| 0.1106412| 0.0000308| 0.0019541|
## | 33| 70| -8.4871946| 3.3557351| 1| 0.1289860| 0.0000003| 0.0002060|
## | 34| 70| -10.4063801| 4.4304825| 1| 0.1515322| 0.0000000| 0.0000302|
## | 35| 70| -20.7616144| 7.8175968| 1| 0.0043245| 0.0000000| 0.0000000|
## | 36| 70| -2.8248217| 1.9195858| 1| 0.7186078| 0.0013760| 0.0559975|
## | 37| 70| -14.4717602| 4.6274100| 1| 0.0044872| 0.0000000| 0.0000005|
## | 38| 70| -9.2707323| 9.7413177| 1| 0.9999458| 0.0000000| 0.0000941|
## | 39| 70| -9.2296597| 4.3210722| 1| 0.3185685| 0.0000000| 0.0000981|
## | 40| 70| 2.0401050| 1.5958563| 1| 0.9943362| 0.2520378| 0.8849440|
## | 41| 70| -2.6056122| 1.5678657| 1| 0.6147693| 0.0034066| 0.0687781|
## | 42| 70| 0.1850422| 0.8346177| 1| 0.8606732| 0.1898772| 0.5461290|
## | 43| 70| -0.0610645| 1.5583003| 1| 0.9522662| 0.0424792| 0.4847386|
## | 44| 70| -12.4785519| 4.0146730| 1| 0.0098558| 0.0000000| 0.0000038|
## | 45| 70| -0.2802553| 1.1577855| 1| 0.8796378| 0.0724591| 0.4303912|
## | 46| 70| -3.2950786| 1.4527081| 1| 0.3898909| 0.0021452| 0.0357404|
## | 47| 70| -31.6231817| 13.8222423| 1| 0.0106490| 0.0000000| 0.0000000|
## | 48| 70| -14.4548003| 4.9005169| 1| 0.0077692| 0.0000000| 0.0000005|
## | 49| 70| -34.3228535| 16.1666873| 1| 0.0668480| 0.0000000| 0.0000000|
## | 50| 70| -12.4045291| 4.1378816| 1| 0.0134627| 0.0000000| 0.0000041|
## | 51| 70| -6.8852425| 2.3565855| 1| 0.0939497| 0.0000101| 0.0010217|
## | 52| 70| -5.2432431| 2.7450962| 1| 0.5342327| 0.0000243| 0.0052553|
## | 53| 70| -2.3595055| 1.2989860| 1| 0.5464921| 0.0073512| 0.0863132|
## | 54| 70| -10.2313272| 3.9021736| 1| 0.0702362| 0.0000000| 0.0000360|
## | 55| 70| -10.2844982| 6.2826587| 1| 0.8838611| 0.0000000| 0.0000342|
## | 56| 70| -18.1260467| 6.4747248| 1| 0.0043398| 0.0000000| 0.0000000|
## | 57| 70| -33.4893127| 11.6151491| 1| 0.0000220| 0.0000000| 0.0000000|
## | 58| 70| -15.3783543| 5.2309060| 1| 0.0059065| 0.0000000| 0.0000002|
## | 59| 70| -2.4045493| 1.6917267| 1| 0.7132528| 0.0032679| 0.0828264|
## | 60| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 61| 70| 2.6605079| 1.3166862| 1| 0.9947340| 0.5199402| 0.9346557|
## | 62| 70| 3.7632055| 1.3219807| 1| 0.9982637| 0.7635286| 0.9773172|
## | 63| 70| 3.7632055| 1.3219807| 1| 0.9982637| 0.7635286| 0.9773172|
## | 64| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 65| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 66| 70| 2.4193792| 0.7350687| 1| 0.9793690| 0.7268392| 0.9182932|
## | 67| 70| 3.9183109| 1.7159654| 1| 0.9993124| 0.6352992| 0.9805127|
## | 68| 70| 1.9540817| 0.8601439| 1| 0.9744183| 0.5666508| 0.8758910|
## | 69| 70| 1.9540817| 0.8601439| 1| 0.9744183| 0.5666508| 0.8758910|
## | 70| 70| 0.8344242| 0.7431857| 1| 0.9081346| 0.3492768| 0.6972896|
## | 71| 70| 2.1091871| 1.0229721| 1| 0.9839239| 0.5260169| 0.8917929|
## | 72| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 73| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 74| 70| 2.4193792| 0.7350687| 1| 0.9793690| 0.7268392| 0.9182932|
## | 75| 70| 3.2118474| 1.3449504| 1| 0.9971225| 0.6400877| 0.9612777|
## | 76| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 77| 70| 1.5408691| 0.8117080| 1| 0.9581862| 0.4874830| 0.8235910|
## | 78| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 79| 70| 3.6428183| 5.9671742| 1| 0.9999998| 0.0003181| 0.9744894|
## | 80| 70| 1.9911971| 1.5357760| 1| 0.9933159| 0.2652371| 0.8798697|
## | 81| 70| 3.7632055| 1.3219807| 1| 0.9982637| 0.7635286| 0.9773172|
## | 82| 70| -0.1131494| 0.9154815| 1| 0.8430626| 0.1292627| 0.4717428|
## | 83| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 84| 70| 3.6081002| 1.0351886| 1| 0.9964494| 0.8290814| 0.9736119|
## | 85| 70| -3.1519924| 1.5458330| 1| 0.4694980| 0.0020624| 0.0410128|
## | 86| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 87| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 88| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|
## | 89| 70| 0.4261709| 0.8011435| 1| 0.8804199| 0.2415738| 0.6049590|
## | 90| 70| 3.4529948| 0.9573319| 1| 0.9951762| 0.8287259| 0.9693203|Plot the results
It can also be helpful to use graphs of predicted probabilities to understand the model. The ggplot2 package has been used for graphing. The variable HEATRANK is plotted against the Predicted Probabilities with 95% confidence intervals (CI). In this plot, the minimum y value is the lower CI while the maximum y values are given by the upper CI values.
heatplot <- ggplot(newdata3, aes(x = HEATRANK, y = PredictedProb)) + geom_smooth(aes(ymin = LL,
ymax = UL, line = HEATRANK), alpha = 0.3, size = 1.5) + theme(panel.background=element_rect(fill='lightblue')) +
ggtitle("HEATRANK vs. Probability including Confidence intervals") +
geom_line(aes(colour = HEATRANK), size = 1.5, )
print(heatplot)