R/get_MS_predictions.R
In jhmigueles/MINISTOP: MINISTOP predictions
Defines functions
get_MS_predictions
Documented in
get_MS_predictions
#' get_MS_predictions
#'
#' @description Get predictions from SEM as presented in paper.
#' @return
#' @export
#'
#' @examples
get_MS_predictions = function() {
# packages
suppressMessages(
require(compositions)
)
# # basic data ----
# compo_names = c("VPA", "MPA", "LPA", "SB", "Sleep")
# mean composition and ilrs at baseline ----
m4 = c(10.44407, 55.27578, 347.65288, 511.51217, 515.11511)
sbp <- matrix(c( 1, -1, -1, -1, -1,
0, 1, -1, -1, -1,
0, 0, 1, -1, -1,
0, 0, 0, 1, -1),
ncol = 5, byrow = TRUE)
psi <- compositions::gsi.buildilrBase(t(sbp))
ilr.m4 = compositions::ilr(m4, V = psi)
# new compositions for predictions ----
# change in the ilr of interest (a = change in ilr of interest; b = change in the others to compensate a)
# Change in ILR 1 (from -10 to +30 min of VPA, min by min)
a = c(1, seq(m4[1] - 10, m4[1] + 30) / m4[1])
b = (1440 - (m4[1] * a)) / (m4[2] + m4[3] + m4[4] + m4[5])
ab_ilr1 = data.frame(ilr = c(0, rep(1, length(a) - 1)), a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*a, MPA4 = m4[2]*b, LPA4 = m4[3]*b,
SB4 = m4[4]*b, Sleep4 = m4[5]*b)
rm(a,b)
# Change in ILR 2 (from -30 to +30 min of MPA by 2 min)
# 1440 = 10.44407*(a*b*b*b)^(1/4) + 55.27578*a + 347.65288*b + 511.51217*b + 515.11511*b
# b values calculated in https://www.symbolab.com/solver/equation-calculator
a = seq(m4[2] - 30, m4[2] + 30, by = 2) / m4[2]
b = c(1.02307, 1.0215, 1.01994, 1.01838, 1.01683, 1.01528,
1.01374, 1.0122, 1.01067, 1.00913, 1.00761, 1.00608,
1.00608, 1.00304,1.00152, 1.000000000, 0.998486,0.996973,
0.995462, 0.993953, 0.992446, 0.99094, 0.989435, 0.987932,
0.98643, 0.98493, 0.983431, 0.981932, 0.980435, 0.978939,
0.977444)
ab_ilr2 = data.frame(ilr = 2, a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*(a*b^3)^(1/4), MPA4 = m4[2]*a, LPA4 = m4[3]*b,
SB4 = m4[4]*b, Sleep4 = m4[5]*b)
# Change in ILR 3 (from -90 to +90 min of LPA)
# 1440 = 10.44407(((a·b^2)^(1/3))^(1/4) · (a·b^2)^(1/4)) + 55.27578((a·b^2)^(1/3)) + 347.65288·a + 1026.627b
# b values calculated in https://www.symbolab.com/solver/equation-calculator
a = seq(m4[3] - 90, m4[3] + 90, by = 5) / m4[3]
b = c(1.09031, 1.08524, 1.08017, 1.07512, 1.07006, 1.06502, 1.05998,
1.05495, 1.04993, 1.04491, 1.03989, 1.03489, 1.02989, 1.02489,
1.0199, 1.01492, 1.00994, 1.00497, 1.00000, 0.995037,0.99008,
0.985127,0.980179,0.975236,0.970297,0.965363,0.960434,0.955509,
0.950588,0.945672,0.94076, 0.935853,0.930949,0.92605, 0.921154,
0.916263,0.911375)
ab_ilr3 = data.frame(ilr = 3, a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*((a*b*b)^(1/3))^(1/4) * (a*b*b)^(1/4), MPA4 = m4[2]*(a*b*b)^(1/3), LPA4 = m4[3]*a, SB4 = m4[4]*b, Sleep4 = m4[5]*b)
rm(a,b)
# Change in ILR 4 (from -300 to +300 min of SB)
# 1440=(10.44407(ab)^(1/4))(sqrt(ab)^(1/4))((sqrt(ab)(ab))^(1/3))^(1/4)+55.27578(sqrt(ab)·(ab))^(1/3)+347.65288·sqrt(ab)+511.51217·a+515.11511b
# b values calculated in https://www.symbolab.com/solver/equation-calculator
a = seq(m4[4] - 300, m4[4] + 300, by = 15) / m4[4]
b = c(1.71006670229299, 1.66640191455040, 1.62394854148532, 1.58261212859228, 1.54231162605595, 1.50297737901422, 1.46454833997031,
1.42697109098268, 1.39019849515969, 1.35418872980947, 1.31890450514219, 1.28431230994092, 1.25038235623840, 1.21708743618891,
1.18440303206502, 1.15230689069605, 1.12077876394383, 1.08980008292947, 1.05935418943891, 1.02942553932630, 1.00000000000000,
0.971064552474187,0.942607323564974,0.914617354422225,0.887084495238876,0.859999693825665,0.833354369296217,0.807140733529092,
0.781351647358753,0.755980577100591,0.731021471655306,0.706469067493431,0.682318347884886,0.658564858123805,0.635204596272442,
0.612233992576200,0.589649891480090,0.567449460973740,0.545630480697629,0.524190876360710,0.503129014405917)
ab_ilr4 = data.frame(ilr = 4, a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*(a*b)^(1/4)*(sqrt(a*b)^(1/4))*((sqrt(a*b)*(a*b))^(1/3))^(1/4),
MPA4 = m4[2]*(sqrt(a*b)*(a*b))^(1/3), LPA4 = m4[3]*sqrt(a*b),
SB4 = m4[4]*a, Sleep4 = m4[5]*b)
out = rbind(ab_ilr1, ab_ilr2, ab_ilr3, ab_ilr4)
# ilr for the new compositions
out[, 9:12] = compositions::ilr(out[,4:8], V = psi)
colnames(out)[9:12] = c("ilr1_bl", "ilr2_bl", "ilr3_bl", "ilr4_bl")
# clean environment
rm(a, b, ab_ilr1, ab_ilr2, ab_ilr3, ab_ilr4)
# -------------------------------------------------------------------------
# Predictions of outcomes at 9 years old ----------------------------------
# -------------------------------------------------------------------------
# Mean predictions at follow up (predicted from lavaan) ----
attach(dat)
for (i in 1:nrow(out)) {
ilr.m9_bc = predict_lavaan(fit = fCompo1_BC.sex.age,
newdata = cbind.data.frame(FMI4 = mean(FMI4),
FFMI4 = mean(FFMI4),
VPA4_MPA4.LPA4.SB4.SLEEP4 = out[i, "ilr1_bl"],
MPA4_LPA4.SB4.SLEEP4 = out[i, "ilr2_bl"],
LPA4_SB4.SLEEP4 = out[i, "ilr3_bl"],
SB4_SLEEP4 = out[i, "ilr4_bl"],
sex = mean(sex),
age4 = mean(age4),
age9 = mean(age9)))
ilr.m9_pf = predict_lavaan(fit = fCompo1_PF.sex.age,
newdata = cbind.data.frame(laps20m4 = mean(laps20m4),
Speed4 = mean(Speed4),
Str4 = mean(Str4),
VPA4_MPA4.LPA4.SB4.SLEEP4 = out[i, "ilr1_bl"],
MPA4_LPA4.SB4.SLEEP4 = out[i, "ilr2_bl"],
LPA4_SB4.SLEEP4 = out[i, "ilr3_bl"],
SB4_SLEEP4 = out[i, "ilr4_bl"],
sex = mean(sex),
age4 = mean(age4),
age9 = mean(age9)))
if (i == 1) {
preds9 = as.data.frame(matrix(NA, nrow = nrow(out), ncol = 14))
colnames(preds9) = c("ilr1_fu", "ilr2_fu", "ilr3_fu", "ilr4_fu",
"VPA9", "MPA9", "LPA9", "SB9", "Sleep9",
"FMI9", "FFMI9","laps20m9", "Speed9", "Str9")
}
# store output
preds9[i,] = c(as.numeric(ilr.m9_bc[1, c("VPA9_MPA9.LPA9.SB9.SLEEP9", "MPA9_LPA9.SB9.SLEEP9",
"LPA9_SB9.SLEEP9", "SB9_SLEEP9")]),
compositions::clo(compositions::ilrInv(ilr.m9_bc[1, c("VPA9_MPA9.LPA9.SB9.SLEEP9", "MPA9_LPA9.SB9.SLEEP9",
"LPA9_SB9.SLEEP9", "SB9_SLEEP9")], V = psi), total = 1440),
as.numeric(ilr.m9_bc[1, c("FMI9", "FFMI9")]),
as.numeric(ilr.m9_pf[1, c("laps20m9", "Speed9", "Str9")]))
}
detach(dat)
out = cbind(out, preds9)
out = round(out, 3)
for (i in 2:nrow(out)) out[i,4:8] = out[i,4:8] - out[1,4:8]
for (i in 2:nrow(out)) out[i,17:26] = out[i,17:26] - out[1,17:26]
# return
return(out)
}
jhmigueles/MINISTOP documentation built on June 16, 2022, 9:54 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions get_MS_predictions
Documented in get_MS_predictions
#' get_MS_predictions
#'
#' @description Get predictions from SEM as presented in paper.
#' @return
#' @export
#'
#' @examples
get_MS_predictions = function() {
# packages
suppressMessages(
require(compositions)
)
# # basic data ----
# compo_names = c("VPA", "MPA", "LPA", "SB", "Sleep")
# mean composition and ilrs at baseline ----
m4 = c(10.44407, 55.27578, 347.65288, 511.51217, 515.11511)
sbp <- matrix(c( 1, -1, -1, -1, -1,
0, 1, -1, -1, -1,
0, 0, 1, -1, -1,
0, 0, 0, 1, -1),
ncol = 5, byrow = TRUE)
psi <- compositions::gsi.buildilrBase(t(sbp))
ilr.m4 = compositions::ilr(m4, V = psi)
# new compositions for predictions ----
# change in the ilr of interest (a = change in ilr of interest; b = change in the others to compensate a)
# Change in ILR 1 (from -10 to +30 min of VPA, min by min)
a = c(1, seq(m4[1] - 10, m4[1] + 30) / m4[1])
b = (1440 - (m4[1] * a)) / (m4[2] + m4[3] + m4[4] + m4[5])
ab_ilr1 = data.frame(ilr = c(0, rep(1, length(a) - 1)), a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*a, MPA4 = m4[2]*b, LPA4 = m4[3]*b,
SB4 = m4[4]*b, Sleep4 = m4[5]*b)
rm(a,b)
# Change in ILR 2 (from -30 to +30 min of MPA by 2 min)
# 1440 = 10.44407*(a*b*b*b)^(1/4) + 55.27578*a + 347.65288*b + 511.51217*b + 515.11511*b
# b values calculated in https://www.symbolab.com/solver/equation-calculator
a = seq(m4[2] - 30, m4[2] + 30, by = 2) / m4[2]
b = c(1.02307, 1.0215, 1.01994, 1.01838, 1.01683, 1.01528,
1.01374, 1.0122, 1.01067, 1.00913, 1.00761, 1.00608,
1.00608, 1.00304,1.00152, 1.000000000, 0.998486,0.996973,
0.995462, 0.993953, 0.992446, 0.99094, 0.989435, 0.987932,
0.98643, 0.98493, 0.983431, 0.981932, 0.980435, 0.978939,
0.977444)
ab_ilr2 = data.frame(ilr = 2, a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*(a*b^3)^(1/4), MPA4 = m4[2]*a, LPA4 = m4[3]*b,
SB4 = m4[4]*b, Sleep4 = m4[5]*b)
# Change in ILR 3 (from -90 to +90 min of LPA)
# 1440 = 10.44407(((a·b^2)^(1/3))^(1/4) · (a·b^2)^(1/4)) + 55.27578((a·b^2)^(1/3)) + 347.65288·a + 1026.627b
# b values calculated in https://www.symbolab.com/solver/equation-calculator
a = seq(m4[3] - 90, m4[3] + 90, by = 5) / m4[3]
b = c(1.09031, 1.08524, 1.08017, 1.07512, 1.07006, 1.06502, 1.05998,
1.05495, 1.04993, 1.04491, 1.03989, 1.03489, 1.02989, 1.02489,
1.0199, 1.01492, 1.00994, 1.00497, 1.00000, 0.995037,0.99008,
0.985127,0.980179,0.975236,0.970297,0.965363,0.960434,0.955509,
0.950588,0.945672,0.94076, 0.935853,0.930949,0.92605, 0.921154,
0.916263,0.911375)
ab_ilr3 = data.frame(ilr = 3, a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*((a*b*b)^(1/3))^(1/4) * (a*b*b)^(1/4), MPA4 = m4[2]*(a*b*b)^(1/3), LPA4 = m4[3]*a, SB4 = m4[4]*b, Sleep4 = m4[5]*b)
rm(a,b)
# Change in ILR 4 (from -300 to +300 min of SB)
# 1440=(10.44407(ab)^(1/4))(sqrt(ab)^(1/4))((sqrt(ab)(ab))^(1/3))^(1/4)+55.27578(sqrt(ab)·(ab))^(1/3)+347.65288·sqrt(ab)+511.51217·a+515.11511b
# b values calculated in https://www.symbolab.com/solver/equation-calculator
a = seq(m4[4] - 300, m4[4] + 300, by = 15) / m4[4]
b = c(1.71006670229299, 1.66640191455040, 1.62394854148532, 1.58261212859228, 1.54231162605595, 1.50297737901422, 1.46454833997031,
1.42697109098268, 1.39019849515969, 1.35418872980947, 1.31890450514219, 1.28431230994092, 1.25038235623840, 1.21708743618891,
1.18440303206502, 1.15230689069605, 1.12077876394383, 1.08980008292947, 1.05935418943891, 1.02942553932630, 1.00000000000000,
0.971064552474187,0.942607323564974,0.914617354422225,0.887084495238876,0.859999693825665,0.833354369296217,0.807140733529092,
0.781351647358753,0.755980577100591,0.731021471655306,0.706469067493431,0.682318347884886,0.658564858123805,0.635204596272442,
0.612233992576200,0.589649891480090,0.567449460973740,0.545630480697629,0.524190876360710,0.503129014405917)
ab_ilr4 = data.frame(ilr = 4, a = a, b = b,
# time-use composition at 4 years old
VPA4 = m4[1]*(a*b)^(1/4)*(sqrt(a*b)^(1/4))*((sqrt(a*b)*(a*b))^(1/3))^(1/4),
MPA4 = m4[2]*(sqrt(a*b)*(a*b))^(1/3), LPA4 = m4[3]*sqrt(a*b),
SB4 = m4[4]*a, Sleep4 = m4[5]*b)
out = rbind(ab_ilr1, ab_ilr2, ab_ilr3, ab_ilr4)
# ilr for the new compositions
out[, 9:12] = compositions::ilr(out[,4:8], V = psi)
colnames(out)[9:12] = c("ilr1_bl", "ilr2_bl", "ilr3_bl", "ilr4_bl")
# clean environment
rm(a, b, ab_ilr1, ab_ilr2, ab_ilr3, ab_ilr4)
# -------------------------------------------------------------------------
# Predictions of outcomes at 9 years old ----------------------------------
# -------------------------------------------------------------------------
# Mean predictions at follow up (predicted from lavaan) ----
attach(dat)
for (i in 1:nrow(out)) {
ilr.m9_bc = predict_lavaan(fit = fCompo1_BC.sex.age,
newdata = cbind.data.frame(FMI4 = mean(FMI4),
FFMI4 = mean(FFMI4),
VPA4_MPA4.LPA4.SB4.SLEEP4 = out[i, "ilr1_bl"],
MPA4_LPA4.SB4.SLEEP4 = out[i, "ilr2_bl"],
LPA4_SB4.SLEEP4 = out[i, "ilr3_bl"],
SB4_SLEEP4 = out[i, "ilr4_bl"],
sex = mean(sex),
age4 = mean(age4),
age9 = mean(age9)))
ilr.m9_pf = predict_lavaan(fit = fCompo1_PF.sex.age,
newdata = cbind.data.frame(laps20m4 = mean(laps20m4),
Speed4 = mean(Speed4),
Str4 = mean(Str4),
VPA4_MPA4.LPA4.SB4.SLEEP4 = out[i, "ilr1_bl"],
MPA4_LPA4.SB4.SLEEP4 = out[i, "ilr2_bl"],
LPA4_SB4.SLEEP4 = out[i, "ilr3_bl"],
SB4_SLEEP4 = out[i, "ilr4_bl"],
sex = mean(sex),
age4 = mean(age4),
age9 = mean(age9)))
if (i == 1) {
preds9 = as.data.frame(matrix(NA, nrow = nrow(out), ncol = 14))
colnames(preds9) = c("ilr1_fu", "ilr2_fu", "ilr3_fu", "ilr4_fu",
"VPA9", "MPA9", "LPA9", "SB9", "Sleep9",
"FMI9", "FFMI9","laps20m9", "Speed9", "Str9")
}
# store output
preds9[i,] = c(as.numeric(ilr.m9_bc[1, c("VPA9_MPA9.LPA9.SB9.SLEEP9", "MPA9_LPA9.SB9.SLEEP9",
"LPA9_SB9.SLEEP9", "SB9_SLEEP9")]),
compositions::clo(compositions::ilrInv(ilr.m9_bc[1, c("VPA9_MPA9.LPA9.SB9.SLEEP9", "MPA9_LPA9.SB9.SLEEP9",
"LPA9_SB9.SLEEP9", "SB9_SLEEP9")], V = psi), total = 1440),
as.numeric(ilr.m9_bc[1, c("FMI9", "FFMI9")]),
as.numeric(ilr.m9_pf[1, c("laps20m9", "Speed9", "Str9")]))
}
detach(dat)
out = cbind(out, preds9)
out = round(out, 3)
for (i in 2:nrow(out)) out[i,4:8] = out[i,4:8] - out[1,4:8]
for (i in 2:nrow(out)) out[i,17:26] = out[i,17:26] - out[1,17:26]
# return
return(out)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.