CODEWITHSUNDEEP
rcurve-fitting
SunWuKung
SunWuKung

Reputation: 557

comparing segmented models in R

I have the hypothesis that the automatically found model below is not significantly better than one where the middle segment is horizontal. How could I test that?

df<-structure(list(ageThen=c(9,10,11,12,13,14,15,16,17,
18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,
34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,
50,51,52,53,54,55,56,57,58,59,60,61,62),mTh=c(-0.057,
-0.253,-0.345,-0.185,-0.155,-0.013,0.285,0.16,0.197,0.199,
0.215,0.288,0.401,0.363,0.387,0.37,0.387,0.28,0.571,
0.383,0.297,0.366,0.36,0.25,0.269,0.235,0.273,0.336,
0.354,0.286,0.331,0.21,0.32,0.278,0.195,0.257,0.259,
0.251,0.222,0.206,0.214,-0.072,-0.123,-0.043,-0.003,0.116,
-0.193,-0.218,-0.278,-0.265,-0.218,-0.541,-0.76,-0.401
),n=c(64L,524L,20595L,2504L,795L,704L,1700L,1239L,
1273L,1149L,1011L,1122L,1031L,814L,717L,667L,462L,414L,
405L,313L,256L,305L,187L,255L,240L,221L,262L,227L,230L,
239L,199L,290L,201L,246L,217L,215L,273L,229L,213L,193L,
199L,204L,159L,207L,148L,121L,115L,89L,87L,78L,68L,
85L,55L,80L)),class=c("tbl_df","tbl","data.frame"),row.names=c(NA,
-54L))

library(segmented)
m1<-lm(mTh~ageThen,data=df,weights=n)##initialfit
s2<-segmented(m1,psi=c(20,50))##twobreakpoints,estimatedstartingvalues
plot(mTh~ageThen,data=df)
lines(df$ageThen,predict(s2),col=2,lwd=2)

Upvotes: 1

Views: 82

Answers (1)

SunWuKung
SunWuKung

Reputation: 557

This workflow seems to solve my problem. The literature that I saw suggests that the there is a significant increase in the predictive value of the model if the elpd_diff/se_diff (in absolute terms) is larger than 2, in other places 5. In my case it isn't, so unconstrained model is not much better predictor than the constrained (theoretical). A citable source would be much appreciated to corroborate this.

https://lindeloev.github.io/mcp/articles/comparison.html#what-is-loo-cv.

https://discourse.mc-stan.org/t/interpreting-elpd-diff-loo-package/1628

library(mcp)

df<-structure(list(ageThen=c(9,10,11,12,13,14,15,16,17,
18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,
34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,
50,51,52,53,54,55,56,57,58,59,60,61,62),mTh=c(-0.057,
-0.253,-0.345,-0.185,-0.155,-0.013,0.285,0.16,0.197,0.199,
0.215,0.288,0.401,0.363,0.387,0.37,0.387,0.28,0.571,
0.383,0.297,0.366,0.36,0.25,0.269,0.235,0.273,0.336,
0.354,0.286,0.331,0.21,0.32,0.278,0.195,0.257,0.259,
0.251,0.222,0.206,0.214,-0.072,-0.123,-0.043,-0.003,0.116,
-0.193,-0.218,-0.278,-0.265,-0.218,-0.541,-0.76,-0.401
),n=c(64L,524L,20595L,2504L,795L,704L,1700L,1239L,
1273L,1149L,1011L,1122L,1031L,814L,717L,667L,462L,414L,
405L,313L,256L,305L,187L,255L,240L,221L,262L,227L,230L,
239L,199L,290L,201L,246L,217L,215L,273L,229L,213L,193L,
199L,204L,159L,207L,148L,121L,115L,89L,87L,78L,68L,
85L,55L,80L)),class=c("tbl_df","tbl","data.frame"),row.names=c(NA,
-54L))

modelC = list(
  mTh | weights(n) ~ 1 + ageThen,  # slope
  ~ 0 ,                            # joined plateau
  ~ 0 + ageThen                    # joined slope
)

fitC<-mcp(modelC, data = df, prior = list(cp_1 = "dunif(12, 25)", cp_2 = "dunif(35, 60)"))
plot(fitC)
summary(fitC)

modelUc = list(
  mTh | weights(n)  ~ 1 + ageThen,  # slope
  ~ 0 + ageThen,                    # joined slope
  ~ 0 + ageThen                     # joined slope
)
fitUc<-mcp(modelUc, data = df, prior = list(cp_1 = "dunif(12, 25)", cp_2 = "dunif(35, 60)"))
plot(fitUc)
summary(fitUc)

fitUc$loo = loo(fitUc)
fitC$loo = loo(fitC)

library(tidyverse)
as.data.frame(loo::loo_compare(fitUc$loo,fitC$loo))%>%
  mutate(diffRatio=elpd_diff/se_diff, .keep='used')

Upvotes: 1

Related Questions