Segmentation Analysis
Julie Wisch
September 12, 2018
I used segmentation analysis to try and capture that change that we saw in the ADRC data review.
I used maximum likelihood tests to optimize a mixed effect model. I included the following as covariates: sex, education, age at PET scan, race, APOE status, interaction of race and sex, interaction of race and age, and interaction of sex and PET scan, taking into account repeated measures for individuals.
## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ EDUC + PETage + Race + apoe4 + Race:PETage +
## Model: (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 8 1048.6 1087.0 -516.30 1032.6
## Model.null 11 1044.3 1097.2 -511.17 1022.3 10.261 3 0.01647 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + PETage + Race + apoe4 +
## Model: GENDER:Race + Race:PETage + GENDER:PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 10 1044.7 1092.7 -512.32 1024.7
## Model.null 11 1044.3 1097.2 -511.17 1022.3 2.3165 1 0.128## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ EDUC + GENDER + Race + apoe4 + GENDER:Race +
## Model: (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 8 1295.2 1333.6 -639.60 1279.2
## Model.null 11 1044.3 1097.2 -511.17 1022.3 256.88 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ EDUC + GENDER + Race + PETage +
## Model: GENDER:Race + Race:PETage + GENDER:PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 10 1125.4 1173.5 -552.72 1105.4
## Model.null 11 1044.3 1097.2 -511.17 1022.3 83.1 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ EDUC + GENDER + Race + PETage +
## Model: GENDER:Race + Race:PETage + apoe4 + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 10 1044.5 1092.5 -512.26 1024.5
## Model.null 11 1044.3 1097.2 -511.17 1022.3 2.1902 1 0.1389## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + apoe4 +
## Model: GENDER:PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 8 1041.3 1079.8 -512.67 1025.3
## Model.null 11 1044.3 1097.2 -511.17 1022.3 3.0059 3 0.3907## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model: apoe4 + GENDER:Race + GENDER:PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 10 1043.1 1091.1 -511.55 1023.1
## Model.null 11 1044.3 1097.2 -511.17 1022.3 0.769 1 0.3805## Data: as.data.frame(df[complete.cases(df$apoe4), ])
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model: apoe4 + Race:PETage + GENDER:PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + EDUC + PETage + Race +
## Model.null: apoe4 + GENDER:Race + Race:PETage + GENDER:PETage + (1 |
## Model.null: ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 10 1042.9 1090.9 -511.44 1022.9
## Model.null 11 1044.3 1097.2 -511.17 1022.3 0.5508 1 0.458The following parameters were significant: sex, age at PET scan and APOE status.
Model.final<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + PETage + apoe4 + (1|ID), data = as.data.frame(df[complete.cases(df$apoe4),]))
summary(Model.final)## Linear mixed model fit by REML ['lmerMod']
## Formula: PUP_fSUVR_rsf_TOT_CORTMEAN ~ GENDER + PETage + apoe4 + (1 | ID)
## Data: as.data.frame(df[complete.cases(df$apoe4), ])
##
## REML criterion at convergence: 1053
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1474 -0.2891 -0.0423 0.2745 4.2380
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.43279 0.6579
## Residual 0.02898 0.1702
## Number of obs: 900, groups: ID, 526
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.367049 0.150031 -9.112
## GENDERmale 0.152764 0.059866 2.552
## PETage 0.037132 0.002079 17.861
## apoe41 0.578385 0.061118 9.463
##
## Correlation of Fixed Effects:
## (Intr) GENDER PETage
## GENDERmale -0.090
## PETage -0.956 -0.072
## apoe41 -0.187 -0.032 0.045I used a broken stick mixed effect model to look for a breakpoint in amyloid accumulation (Muggeo, 2016). There’s a breakpoint at 63.6 years.
library(nlme)
df<-df[complete.cases(df$apoe4),]
df.male<-subset(df, df$GENDER == "male")
df.female<-subset(df, df$GENDER == "female")
o<-lme(PUP_fSUVR_rsf_TOT_CORTMEAN~PETage+apoe4, random=list(ID=pdDiag(~1+PETage)), data=df)
oseg<-segmented.lme(o, Z=PETage, random=list(ID=pdDiag(~1+PETage+U+G0)))
oseg## Segmented Linear mixed-effects model
## psi.link = identity
##
## Linear mixed-effects model fit by REML
## Data: mf
## Log-restricted-likelihood: -219.1503
## Fixed: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + U + G0
## (Intercept) PETage apoe41 U G0
## 1.233665919 -0.004689664 0.094203902 0.065520106 62.955182462
##
## Random effects:
## Formula: ~1 + PETage + U + G0 | id
## Structure: Diagonal
## (Intercept) PETage U G0 Residual
## StdDev: 0.001693854 0.002627061 0.05548306 5.214227 0.09434226
##
## Number of Observations: 900
## Number of Groups: 526confint(oseg)## (Intercept) PETage apoe41 U G0
## 2.5% 1.029205 -0.008138159 0.04469935 0.05735606 61.91863
## 97.5% 1.438127 -0.001241168 0.14370845 0.07368416 63.99173This breakpoint holds up for females, but not for males.
o.f<-lme(PUP_fSUVR_rsf_TOT_CORTMEAN~PETage+apoe4, random=list(ID=pdDiag(~1+PETage)), data=df.female)
oseg.f<-segmented.lme(o.f, Z=PETage, random=list(ID=pdDiag(~1+PETage+U+G0)))
o.f1<-lme(PUP_fSUVR_rsf_TOT_CORTMEAN~PETage+apoe4, random=list(ID=pdDiag(~1+PETage)),
data=as.data.frame(subset(df.female, df.female$PETage < 63.6)))
o.f1<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1|ID),
data=as.data.frame(subset(df.female, df.female$PETage < 63.6)))
#confint(o.f1)
o.f2<-lme(PUP_fSUVR_rsf_TOT_CORTMEAN~PETage+apoe4, random=list(ID=pdDiag(~1+PETage)),
data=as.data.frame(subset(df.female, df.female$PETage > 63.6)))
o.f2<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1|ID),
data=as.data.frame(subset(df.female, df.female$PETage > 63.6)))
#confint(o.f2)
o.m<-lme(PUP_fSUVR_rsf_TOT_CORTMEAN~PETage+apoe4, random=list(ID=pdDiag(~1+PETage+apoe4)), data=df.male)
oseg.m<-segmented.lme(o.m, Z=PETage, random=list(ID=pdDiag(~1+PETage+U+G0)), psi.link = "identity")
oseg.f## Segmented Linear mixed-effects model
## psi.link = identity
##
## Linear mixed-effects model fit by REML
## Data: mf
## Log-restricted-likelihood: -97.43765
## Fixed: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + U + G0
## (Intercept) PETage apoe41 U G0
## 0.859090781 0.002085661 0.086727805 0.052606673 63.609634605
##
## Random effects:
## Formula: ~1 + PETage + U + G0 | id
## Structure: Diagonal
## (Intercept) PETage U G0 Residual
## StdDev: 0.003936403 0.002458759 0.05754963 6.431482 0.08928796
##
## Number of Observations: 551
## Number of Groups: 309oseg.m## Segmented Linear mixed-effects model
## psi.link = identity
##
## Linear mixed-effects model fit by REML
## Data: mf
## Log-restricted-likelihood: -117.0129
## Fixed: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + U + G0
## (Intercept) PETage apoe41 U G0
## 1.683418e+06 7.690488e+02 5.483149e-01 -7.690075e+02 -2.450406e+03
##
## Random effects:
## Formula: ~1 + PETage + U + G0 | id
## Structure: Diagonal
## (Intercept) PETage U G0 Residual
## StdDev: 0.00266161 1.665838e-05 0.0002737582 3.172825e-06 0.1454577
##
## Number of Observations: 349
## Number of Groups: 217Model.null<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1|ID), data=as.data.frame(subset(df.female, df.female$PETage < 63.6)))
Model<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage+ (1|ID), data=as.data.frame(subset(df.female, df.female$PETage < 63.6)))
anova(Model.null, Model)## Data: as.data.frame(subset(df.female, df.female$PETage < 63.6))
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 4 13.3766 25.8962 -2.6883 5.3766
## Model.null 5 -6.8028 8.8467 8.4014 -16.8028 22.179 1 2.483e-06
##
## Model
## Model.null ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ apoe4+ (1|ID), data=as.data.frame(subset(df.female, df.female$PETage < 63.6)))
anova(Model.null, Model)## Data: as.data.frame(subset(df.female, df.female$PETage < 63.6))
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ apoe4 + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 4 20.6706 33.190 -6.3353 12.671
## Model.null 5 -6.8028 8.847 8.4014 -16.803 29.473 1 5.669e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model.null<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1|ID), data=as.data.frame(subset(df.female, df.female$PETage > 63.6)))
Model<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage+ (1|ID), data=as.data.frame(subset(df.female, df.female$PETage > 63.6)))
anova(Model.null, Model)## Data: as.data.frame(subset(df.female, df.female$PETage > 63.6))
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 4 549.19 564.97 -270.59 541.19
## Model.null 5 500.38 520.11 -245.19 490.38 50.808 1 1.018e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ apoe4+ (1|ID), data=as.data.frame(subset(df.female, df.female$PETage > 63.6)))
anova(Model.null, Model)## Data: as.data.frame(subset(df.female, df.female$PETage > 63.6))
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ apoe4 + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 4 589.44 605.22 -290.72 581.44
## Model.null 5 500.38 520.11 -245.19 490.38 91.06 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model.null<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1|ID), data=df.male)
Model<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage+ (1|ID), data=df.male)
anova(Model.null, Model)## Data: df.male
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 4 474.35 489.77 -233.18 466.35
## Model.null 5 448.36 467.64 -219.18 438.36 27.989 1 1.22e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model<-lmer(PUP_fSUVR_rsf_TOT_CORTMEAN ~ apoe4+ (1|ID), data=df.male)
anova(Model.null, Model)## Data: df.male
## Models:
## Model: PUP_fSUVR_rsf_TOT_CORTMEAN ~ apoe4 + (1 | ID)
## Model.null: PUP_fSUVR_rsf_TOT_CORTMEAN ~ PETage + apoe4 + (1 | ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## Model 4 568.13 583.55 -280.06 560.13
## Model.null 5 448.36 467.64 -219.18 438.36 121.76 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here’s a visualization of what that looks like. The average rate of amyloid accumulation triples for females after 63.6 (and is significantly different after 63.6 than before it). For males it is uniform, but closer to the higher rate of accumulation than the lower rate. Here’s a graph.
Confidence<-as.data.frame(cbind(as.data.frame(confint(o.f1)), as.data.frame(confint(o.f2)),as.data.frame(confint(Model.null))))
colnames(Confidence)<-c("o.f1_lower", "o.f1_upper", "o.f2_lower", "o.f2_upper", "o.m_lower", "o.m_upper")
x1<-seq(from = 40, to = 90, by = 1)
x2<-seq(from = 40, to = 63.6, length.out = 51)
x3<-seq(from = 63.6, to = 90, length.out = 51)
df.plot<-as.data.frame(cbind(seq(from = 40, to = 90, by = 1), 0.041193*x1-1.499763,
Confidence$o.m_lower[3]+Confidence$o.m_lower[4]*x1,
Confidence$o.m_upper[3]+Confidence$o.m_upper[4]*x1,
seq(from = 40, to = 63.6, length.out = 51),
0.02194053*x2-0.24517285,
Confidence$o.f1_lower[3]+Confidence$o.f1_lower[4]*x2,
Confidence$o.f1_upper[3]+Confidence$o.f1_upper[4]*x2,
seq(from = 63.6, to = 90, length.out = 51),
0.05075227*x3-2.39503348,
Confidence$o.f2_lower[3]+Confidence$o.f2_lower[4]*x3,
Confidence$o.f2_upper[3]+Confidence$o.f2_upper[4]*x3))
colnames(df.plot)<-c("x.m", "y.m", "ylow.m", "yup.m","x.f1", "y.f1", "ylow.f1", "yup.f1",
"x.f2", "y.f2", "ylow.f2", "yup.f2")
ggplot(data = df, aes(x = PETage, y = PUP_fSUVR_rsf_TOT_CORTMEAN, group = ID, colour = GENDER)) +
scale_color_manual(values=c("red", "blue"))+
scale_fill_manual(values=c("red", "blue"))+
geom_point(alpha = 0.3) +
#geom_line(data = df.plot.m,aes( x = x, y = y))
geom_abline(slope = 0.041193, intercept = -1.499763, colour = "blue") + #Males
# geom_ribbon(aes(ymin = Confidence$o.m_lower[3]+Confidence$o.m_lower[4]*x,
# ymax = Confidence$o.m_upper[3]+Confidence$o.m_upper[4]*x))
geom_segment(x = 40, y = 0.02194053*40-0.24517285, xend = 63.6, yend = 0.02194053*67.6-0.24517285, colour = "red")+
geom_segment(x = 63.6, y = 0.05075227*67.6-2.39503348, xend = 90, yend = 0.05075227*90-2.39503348, colour = "red")
#ggplot(data = df.plot, aes(x = x.m, y = y.m))+geom_line()+geom_ribbon(aes(ymin = ylow.m, ymax = yup.m), alpha = 0.2)+
# geom_line(data = df.plot, aes(x = x.f1, y = y.f1))+geom_ribbon(aes(ymin = ylow.f1, ymax = yup.f1), alpha = 0.2)+
#geom_line(data = df.plot, aes(x = x.f2, y = y.f2))+geom_ribbon(aes(ymin = ylow.f2, ymax = yup.f2), alpha = 0.2)