12/21/2020

Data Introduction

  • SOF is a tool for screen Frailty Status with total 3 items, the score range from 0-3, higher score means more frailty
  • KCL is another tool for screen Risk of Frailty with total 25 items, he score range from 0-25, higher score means more frailty
  • The dataset also included demographic data such as Age, Gender, Education, Living Status…etc
  • The measurement between pre and post test is Functional Fittness, contains 6 components

Data Structure

dta<- read_excel("project.xlsx", sheet=1)
dim(dta)
## [1] 740  56
names(dta)
##  [1] "Year"             "Area"             "Class"            "ID"              
##  [5] "Age"              "Gender"           "Education"        "Living"          
##  [9] "everfall"         "activity"         "Rate"             "ever_participate"
## [13] "Time"             "SOF_1"            "SOF_2"            "SOF_3"           
## [17] "SOF_S"            "SOF_def"          "K_Q1"             "K_Q2"            
## [21] "K_Q3"             "K_Q4"             "K_Q5"             "K_Q6"            
## [25] "K_Q7"             "K_Q8"             "K_Q9"             "K_Q10"           
## [29] "K_Q11"            "K_Q12"            "K_Q13"            "K_Q14"           
## [33] "K_Q15"            "K_Q16"            "K_Q17"            "K_Q18"           
## [37] "K_Q19"            "K_Q20"            "K_Q21"            "K_Q22"           
## [41] "K_Q23"            "K_Q24"            "K_Q25"            "KCLs"            
## [45] "KCL_def"          "Height"           "Weight"           "BMI"             
## [49] "WC"               "HIPC"             "STAND_30sec"      "Grip"            
## [53] "UpReach"          "DownReach"        "STEP_2min"        "WALK"
# data wrangling
dta <- dta %>% mutate(
        KCL_def=factor(KCL_def, levels=0:1, labels=c("No_risk", "At_Risk")), 
        activity=factor(activity, levels=0:2, labels=c("Low", "Medium", "High")),
        SOF_def=factor(SOF_def, levels = c("Robust", "Pre-frailty", "Frailty")), 
        Gender=factor(ifelse(Gender=="男", "Male", "Female")))
dta$Rate<-as.numeric(dta$Rate)      
dta$WC<-as.numeric(dta$WC)
dta$HIPC<-as.numeric(dta$HIPC)
dta$STAND_30sec<-as.numeric(dta$STAND_30sec)
dta$Grip=as.numeric(dta$Grip)
dta$UpReach=as.numeric(dta$UpReach)
dta$DownReach=as.numeric(dta$DownReach)
dta$STEP_2min=as.numeric(dta$STEP_2min)
dta$WALK=as.numeric(dta$WALK)

Descriptive statistics

Class by Time

addmargins(with(dta, table(Class, Time)))
##      Time
## Class Baseline Intervention Sum
##   B         17           13  30
##   C         27           24  51
##   D         49           46  95
##   E         28           28  56
##   F         36           36  72
##   G         37           37  74
##   H         48           48  96
##   I         36           36  72
##   J         40           34  74
##   K         44           38  82
##   L         38            0  38
##   Sum      400          340 740

Class L only got Baseline Data

dta<-subset(dta, Class!="L")

Delete Class L data

Gender by Time

addmargins((with(dta, table(Gender, Time))))
##         Time
## Gender   Baseline Intervention Sum
##   Female      268          254 522
##   Male         94           86 180
##   Sum         362          340 702
addmargins(with(dta, table(Education, Class)))
##          Class
## Education   B   C   D   E   F   G   H   I   J   K Sum
##    大學     0   0  16   0   0  14  28  18  18  17 111
##    不識字   9  13  11  18  25  14   0  10  12  14 126
##    研究所   0   0   2   0   0   0   4   4   1   0  11
##    缺       0   0   0   4   0   0   0   0   0   0   4
##    高中職   5   2  18   2   9   6  22  12  27  15 118
##    國小    15  32  28  32  37  22   8  12  12  26 224
##    國中     1   4  20   0   1  18  34  16   4  10 108
##    Sum     30  51  95  56  72  74  96  72  74  82 702

Frailty by Time

addmargins(with(dta, table(SOF_def, Time)))
##              Time
## SOF_def       Baseline Intervention Sum
##   Robust           272          238 510
##   Pre-frailty       70           58 128
##   Frailty           18            9  27
##   Sum              360          305 665
a<-as.data.frame(with(dta, ftable(SOF_def,Time)))  
summary(xtabs(Freq~., a))  #Chisq 
## Call: xtabs(formula = Freq ~ ., data = a)
## Number of cases in table: 665 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 1.8555, df = 2, p-value = 0.3954

More Informations

summarySE(dta, measurevar="Age", groupvars=c("Gender", "SOF_def","Time"), na.rm=TRUE)
##    Gender     SOF_def         Time   N      Age        sd        se        ci
## 1  Female      Robust     Baseline 198 72.85859  7.726295 0.5490839  1.082837
## 2  Female      Robust Intervention 173 73.92486  7.328672 0.5571886  1.099808
## 3  Female Pre-frailty     Baseline  54 75.07407  8.377535 1.1400381  2.286627
## 4  Female Pre-frailty Intervention  44 74.79545  8.679244 1.3084453  2.638731
## 5  Female     Frailty     Baseline  14 75.35714 10.322545 2.7588163  5.960060
## 6  Female     Frailty Intervention   8 73.37500 11.070392 3.9139745  9.255079
## 7  Female        <NA>     Baseline   2 73.00000  1.414214 1.0000000 12.706205
## 8  Female        <NA> Intervention  29 67.20690  6.061467 1.1255861  2.305659
## 9    Male      Robust     Baseline  74 75.01351  6.829616 0.7939267  1.582294
## 10   Male      Robust Intervention  65 75.78462  6.936075 0.8603142  1.718675
## 11   Male Pre-frailty     Baseline  16 78.50000 10.250203 2.5625508  5.461948
## 12   Male Pre-frailty Intervention  14 75.35714  9.739745 2.6030564  5.623561
## 13   Male     Frailty     Baseline   4 74.75000  6.238322 3.1191612  9.926563
## 14   Male     Frailty Intervention   1 76.00000        NA        NA       NaN
## 15   Male        <NA> Intervention   6 71.50000  8.871302 3.6216939  9.309861

summarySE(dta, measurevar="Grip", groupvars=c("SOF_def","Time"), na.rm=TRUE)
##       SOF_def         Time   N     Grip       sd        se         ci
## 1      Robust     Baseline 266 25.36429 8.106292 0.4970286  0.9786276
## 2      Robust Intervention 194 25.19691 7.594728 0.5452697  1.0754526
## 3 Pre-frailty     Baseline  70 23.52000 6.387411 0.7634416  1.5230243
## 4 Pre-frailty Intervention  47 22.30851 6.292433 0.9178456  1.8475274
## 5     Frailty     Baseline  18 23.64444 7.582505 1.7872135  3.7706909
## 6     Frailty Intervention   7 22.25714 8.421175 3.1829050  7.7882880
## 7        <NA>     Baseline   2 16.25000 1.767767 1.2500000 15.8827559
## 8        <NA> Intervention   3 23.30000 7.483983 4.3208795 18.5912441

Data Visulization

Take Grip for Example

Overall difference between Before and After Intervention

diffPlot(Grip~Time, data=dta, ylim = c(23, 26), 
         xlab="Intervention Effect", ylab="Mean Grip(kg)")

Relationships between Age, Frailty and Grip strength

Age by Grip

By Gender

Different Baseline of Grip gains by Different Frailty Status

Gain and Baseline Grip

Class gain grip means and CIs

Different Class of Grip gains

sdta <- dtagW %>% group_by(Class) %>% na.omit() %>%
  dplyr::summarise(gain_ave=mean(gain),
                   gain_sd2=var(gain), 
                   n=n())
m_g<-lme4::lmer(gain~(1|Class), data=dtagW)
sdta<-data.frame(Class=row.names(coef(m_g)$Class),
                 yhat=coef(m_g)$Class[,1]) %>%
  inner_join(., sdta, by="Class") %>% select(-2) %>%
  data.frame()
arrange(sdta, -gain_ave)
##   Class   gain_ave  gain_sd2  n
## 1     B  3.0700000 16.957889 10
## 2     E  2.7043478 11.456798 23
## 3     F  1.3750000  5.164667 16
## 4     D  1.3178571 16.305966 28
## 5     C  0.2666667 13.129697 12
## 6     G -0.1562500 13.470625 16
## 7     K -0.5275862  7.430640 29
## 8     H -1.1950000  6.716289 20
## 9     J -2.7095238 10.646905 21
  • Improvement on Grip: B>E>F>D>C>G>K>H>J

By Gender and Class

By Time and Frailty Status

Individuals differences

Analysis

Paired t test for Pre and Post Grip data

with(dtagL,(t.test(Baseline, Intervention, paired = TRUE)))
## 
##  Paired t-test
## 
## data:  Baseline and Intervention
## t = -1.6491, df = 349, p-value = 0.1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.70666040  0.06208898
## sample estimates:
## mean of the differences 
##              -0.3222857

NO significant difference between Pre and Post data

Random Class Intercept

m0<-lme4::lmer(gain~Class+(1|Class), data=dtagL)
m0_a<-lme4::lmer(gain~Class+Time+(1|Class), data=dtagL)
anova(m0, m0_a)
## Data: dtagL
## Models:
## m0: gain ~ Class + (1 | Class)
## m0_a: gain ~ Class + Time + (1 | Class)
##      npar    AIC    BIC logLik deviance Chisq Df Pr(>Chisq)
## m0     11 1834.4 1876.8 -906.2   1812.4                    
## m0_a   12 1836.4 1882.7 -906.2   1812.4     0  1          1

Models Comparison

m1<-lmer(gain~Baseline+Class+(1|Class), data=dtagL)
m2<-lmer(gain~Baseline+Gender+Class+(1|Class), data=dtagL)
m3<-lmer(gain~Baseline+Gender+Age+Class+(1|Class), data=dtagL)
m4<-lmer(gain~Baseline+Gender*Age+Class+(1|Class), data=dtagL) #Interaction
anova(m0, m1, m2, m3, m4)
## Data: dtagL
## Models:
## m0: gain ~ Class + (1 | Class)
## m1: gain ~ Baseline + Class + (1 | Class)
## m2: gain ~ Baseline + Gender + Class + (1 | Class)
## m3: gain ~ Baseline + Gender + Age + Class + (1 | Class)
## m4: gain ~ Baseline + Gender * Age + Class + (1 | Class)
##    npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## m0   11 1834.4 1876.8 -906.20   1812.4                          
## m1   12 1814.4 1860.7 -895.21   1790.4 21.9941  1  2.735e-06 ***
## m2   13 1795.1 1845.2 -884.54   1769.1 21.3328  1  3.861e-06 ***
## m3   14 1788.2 1842.2 -880.12   1760.2  8.8471  1   0.002935 ** 
## m4   15 1789.5 1847.4 -879.77   1759.5  0.6841  1   0.408177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(m2, m3)
  gain gain
Predictors Estimates CI p Estimates CI p
(Intercept) 7.45 3.73 – 11.17 <0.001 13.84 8.23 – 19.46 <0.001
Baseline -0.22 -0.28 – -0.15 <0.001 -0.25 -0.32 – -0.19 <0.001
Gender [Male] 2.53 1.46 – 3.61 <0.001 3.10 1.97 – 4.23 <0.001
Class [C] -2.55 -7.45 – 2.35 0.308 -2.27 -7.08 – 2.55 0.356
Class [D] -1.47 -6.28 – 3.34 0.549 -1.46 -6.18 – 3.26 0.545
Class [E] -1.08 -5.91 – 3.74 0.660 -0.93 -5.68 – 3.81 0.699
Class [F] -1.66 -6.52 – 3.19 0.502 -1.51 -6.28 – 3.26 0.535
Class [G] -2.52 -7.39 – 2.35 0.311 -2.72 -7.51 – 2.06 0.265
Class [H] -3.28 -8.12 – 1.56 0.184 -3.08 -7.84 – 1.68 0.204
Class [J] -5.13 -9.96 – -0.30 0.037 -5.13 -9.88 – -0.38 0.034
Class [K] -3.13 -7.94 – 1.68 0.203 -2.96 -7.68 – 1.77 0.220
Age -0.08 -0.13 – -0.03 0.003
Random Effects
σ2 9.47 9.27
τ00 2.69 Class 2.59 Class
ICC 0.22 0.22
N 9 Class 9 Class
Observations 350 350
Marginal R2 / Conditional R2 0.255 / 0.420 0.271 / 0.430

Age effect with Random Slope

m5<-update(m3,.~(Age|Class))
m5a<-update(m3,.~(Time|Class))
m5b<-update(m3, .~(Baseline+Age|Class))
m5c<-update(m3, .~(Baseline+Time|Class))
anova(m5, m5a, m5b, m5c)
## Data: dtagL
## Models:
## m5: gain ~ (Age | Class)
## m5a: gain ~ (Time | Class)
## m5b: gain ~ (Baseline + Age | Class)
## m5c: gain ~ (Baseline + Time | Class)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m5     5 1850.5 1869.8 -920.23   1840.5                         
## m5a    5 1852.7 1872.0 -921.36   1842.7  0.000  0               
## m5b    8 1828.5 1859.4 -906.25   1812.5 30.224  3  1.238e-06 ***
## m5c    8 1832.6 1863.4 -908.28   1816.6  0.000  0               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compare to m3 model

anova(m3, m5b)
## Data: dtagL
## Models:
## m5b: gain ~ (Baseline + Age | Class)
## m3: gain ~ Baseline + Gender + Age + Class + (1 | Class)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m5b    8 1828.5 1859.4 -906.25   1812.5                         
## m3    14 1788.2 1842.2 -880.12   1760.2 52.273  6  1.643e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

After model comparing, it seems NO need to consider random slope

tab_model(m3)
  gain
Predictors Estimates CI p
(Intercept) 13.84 8.23 – 19.46 <0.001
Baseline -0.25 -0.32 – -0.19 <0.001
Gender [Male] 3.10 1.97 – 4.23 <0.001
Age -0.08 -0.13 – -0.03 0.003
Class [C] -2.27 -7.08 – 2.55 0.356
Class [D] -1.46 -6.18 – 3.26 0.545
Class [E] -0.93 -5.68 – 3.81 0.699
Class [F] -1.51 -6.28 – 3.26 0.535
Class [G] -2.72 -7.51 – 2.06 0.265
Class [H] -3.08 -7.84 – 1.68 0.204
Class [J] -5.13 -9.88 – -0.38 0.034
Class [K] -2.96 -7.68 – 1.77 0.220
Random Effects
σ2 9.27
τ00 Class 2.59
ICC 0.22
N Class 9
Observations 350
Marginal R2 / Conditional R2 0.271 / 0.430

Interpretation

  • The estimate mean of gain Grip for all Class is 13.84 with as variance of 2.59
  • The estimate mean of gain Grip for all Case is 13.84 with as variance of 2.59+9.27=11.86
  • On average, there is a drop of 0.25 kg on gain Grip from an increase of one kg in Baseline average Grip performance
  • Male Grip performance averagely higher than Female, about 3.1kg
  • And Age, there is a drop of 0.08kg on one year growth
  • After adjusting for school difference in Baseline Grip, correlation between gain scores among case attending the same class is 0.22
  • No significant difference between Pre and Post data

Residuals

plot(m3)

Normality

qqnorm(resid(m3))
qqline(resid(m3))

Generalized Linear Models

At Frailty Risk?

Descriptive statistics

dtak<-dta %>% select(Class, ID,KCL_def, Age, Gender, Time, Education, Living)
dtaf <- with(dtak, ftable(Class, Time, KCL_def)) %>% as.data.frame
addmargins(with(dtak, table(Time, KCL_def)))
##               KCL_def
## Time           No_risk At_Risk Sum
##   Baseline         321      32 353
##   Intervention     269      11 280
##   Sum              590      43 633
addmargins(with(dtak, table(Gender, KCL_def)))
##         KCL_def
## Gender   No_risk At_Risk Sum
##   Female     432      33 465
##   Male       158      10 168
##   Sum        590      43 633

Add Variables

dtakW <-dtaf %>% spread(KCL_def, Freq) %>% 
  inner_join(., dtak, by=c("Time", "Class")) %>%
  mutate(Total=No_risk+At_Risk)
g0<-glm(cbind(At_Risk, No_risk)~Time+Gender*Age+Education+Living, data=dtakW, 
        family=binomial(link =logit ))
g1<-glm(cbind(At_Risk, No_risk)~Time+Gender*Age+Education, data=dtakW, 
        family=binomial(link =logit ))
anova(g0, g1, test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: cbind(At_Risk, No_risk) ~ Time + Gender * Age + Education + Living
## Model 2: cbind(At_Risk, No_risk) ~ Time + Gender * Age + Education
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       690     608.25                     
## 2       691     608.57 -1 -0.32146   0.5707

tab_model(g1)
  cbind(At_Risk, No_risk)
Predictors Odds Ratios CI p
(Intercept) 0.03 0.02 – 0.06 <0.001
Time [Intervention] 0.36 0.32 – 0.41 <0.001
Gender [Male] 2.54 0.74 – 8.65 0.137
Age 1.01 1.00 – 1.02 0.007
Education [攼㸹㠼㹡㤼㹡㤼㹡謢戼㸶謖戼㸳 1.45 1.18 – 1.77 <0.001
Education [攼㹥㤼㹤憐<89>挼㸲€] 0.84 0.50 – 1.33 0.479
Education [攼㸸璁㤼㹣 2.55 1.30 – 4.59 0.003
Education [攼㸶㤼㸳朴攼㸸㤼㸱㠼㸹蹓 1.22 1.01 – 1.48 0.037
Education [<9c>攼㹦㤼㸵㤼㹤攼㹥㠼㹦㠼㹢 1.49 1.26 – 1.77 <0.001
Education [<9c>攼㹦㤼㸵㤼㸱 0.97 0.79 – 1.19 0.770
Gender [Male] * Age 0.99 0.97 – 1.01 0.212
Observations 702
R2 Tjur 0.000

Goodness of fit

(g1$null.deviance - g1$deviance)/g1$null.deviance
## [1] 0.373979
1 - pchisq(g1$deviance, g1$df.residual)
## [1] 0.9891091

Residual

qqnorm(resid(g1))
qqline(resid(g1))

Interpretation

  • The Odds of after intervention being at frailty risk is between 0.32-0.41
  • The Odds of Higher Education(大學、研究所) at frailty risk is respectively between 0.56-0.84 and 0.34-0.93

Limitation

##              Time
## SOF_def       Baseline Intervention Sum
##   Robust           272          238 510
##   Pre-frailty       70           58 128
##   Frailty           18            9  27
##   Sum              360          305 665
  • The dropout rate roughly calculated is (360-305)/360=15.2%
  • The SOF_def group sample size is significantly different, hard to infer the effects on different Frailty Status, which is one of major research question.