Looping/sapply through nlme function

I am trying to execute a loop with mixed-model effects with response variable changing. I came from here and here. I should say that I have tried sthg creating a function and then sapply or lapply (wihtout success)

I provide a small dataset (really small) just to represent my original database (much larger and similar to those of longitudinal studies)

data<- structure(list(paciente = structure(c(6527, 6302, 6454, 6302, 
6238), label = "Paciente", format.spss = "F6.0"), edad_s1 = structure(c(63, 
60, 64, 60, 66), label = "Edad", format.spss = "F3.0"), sexo_s1 = structure(c(1L, 
1L, 1L, 1L, 2L), .Label = c("Hombre", "Mujer"), label = "Sexo", class = "factor"), 
    peso1 = c(90.5, 110.2, 92, 107.7, 85), time = structure(c(1L, 
    3L, 3L, 1L, 3L), .Label = c("0", "1", "2"), class = "factor"), 
    cintura1 = c(114, 121, 113, 121, 118), tasis2_e = c(121, 
    142, 123, 138, 146), tadias2_e = c(56, 79, 76, 85, 81), p17_total = c(4, 
    10, 11, 7, 10), geaf_tot = c(559.44, 923.08, 1272.73, 1384.62, 
    1258.74), glucosa = c(90, 149, 90, 126, 185), albumi = c(4.48, 
    4.65, 4.59, 4.65, 4.72), coltot = c(250, 241, 280, 211, 207
    ), hdl = c(49, 55, 53, 45, 56), ldl_calc = c(178, 145, NA, 
    137, 114), trigli = c(117, 203, 414, 143, 183), hba1c = c(5.9, 
    6.38, 5.24, 5.86, 8.02), i_hucpeptide = c(988.17, 1046.24, 
    1250.17, 1021.56, 1548.18), i_hughrelin = c(1292.73, 375.44, 
    727.79, 406.6, 859.89), i_hugip = c(2.67, 2.67, 2.67, 2.67, 
    2.67), i_huglp1 = c(60.29, 374.77, 213.82, 258.84, 192.61
    ), i_huglucagon = c(379.54, 781.85, 503.48, 642.2, 554.8), 
    i_huinsulin = c(214.73, 532.7, 391.98, 518.59, 650.66), i_huleptin = c(6049.44, 
    6423.4, 13185.8, 7678.36, 8906.6), i_hupai1 = c(1925.44, 
    3374.85, 2507.28, 2723.42, 2844.02), i_huresistin = c(4502.68, 
    4481.43, 5895.57, 4570.61, 2417.61), i_huvisfatin = c(816.18, 
    1266.7, 2321.13, 1276.66, 1196.37), col_rema = c(23, 41, 
    NA, 29, 37), homa = c(7.953, 32.663, 14.518, 26.89, 49.536
    ), i_pcr = c(NA, 0.14, 0.21, 0.67, 0.17)), row.names = c(NA, 
-5L), class = c("tbl_df", "tbl", "data.frame"))

Afterwards I am defining my iteration and my variables database

ex<- subset(data[, 6:30])

for (i in 1:length(ex)) {
  var_1 <- ex[,i]
  var_1 <- unlist(var_1)
  lme_1 <- lme(var_1 ~ sexo_s1*peso1 + edad_s1 + p17_total + poly(time, 2)*grupo_int_v00,
             random = ~ poly(time, 2)|paciente, control=lmeControl(opt="optim"),
           data = dat_longer, subset = !is.na(var_1))

Error in model.frame.default(formula = ~time + var_1 + sexo_s1 + peso1 +  : 
  invalid type (list) for variable 'var_1'

I have tried unlisting/as.data.frame in before running the loop

for (i in 1:length(data)) {
  var_1 <- data[,i]
  var_1 <- unlist(var_1) #or as.data.frame(var_1)
    lme_1 <- lme(var_1 ~ sexo_s1*peso1 + edad_s1 + p17_total + poly(time, 2)*grupo_int_v00,
             random = ~ poly(time, 2)|paciente, control=lmeControl(opt="optim"),
           data = dat_longer, subset = !is.na(var_1))
}
Error in model.frame.default(formula = ~time + var_1 + sexo_s1 + peso1 +  : 
  variable lengths differ (found for 'var_1')

I have also tried to develop a new function to iterate over

lme_z <- function(z){
out <-  lme(z ~ sexo_s1*peso1 + edad_s1 + p17_total + poly(time, 2)*grupo_int_v00,
             random = ~ poly(time, 2)|paciente, control=lmeControl(opt="optim"),
           data = dat_longer, subset = !is.na(z))
  
}
Error

If there is some contribution to iterate in the response variable (I know Ben Bolker is an expert) Thanks in advance

I tried to attempt to help you , but your code expects a dat_longer that hasnt been provided.
I think you are confusing some differing concepts... the lme() call you are seemingly trying to combine a raw vector of values with symbols that relate to columns of a dataframe. you should probably be consistent and just have a formula which is symbolic and relating to the columns of the data provided in.
I think that for you, constructing the formula as if were a string would be your path of least resistance.
something like



ex<- subset(data[, 6:30])

for (i in 1:length(ex)) {
  var_1 <- names(ex[,i])
  print(paste0("now ", var_1))
  
  lme_1 <- lme(as.formula(paste0(var_1," ~ sexo_s1*peso1 + edad_s1 + p17_total + poly(time, 2)*grupo_int_v00")),
               random = ~ poly(time, 2)|paciente, control=lmeControl(opt="optim"),
               data = data, subset = !is.na(data[var_1]))
  print(summary(lm_1))
}
  

I have updated the database to include grupo_int_v00. It is true that I did not provide dat_longer.

Thank you for the contribution
However, I did encounter the next error:

Error in na.fail.default(list(time = c(0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1,  : 
  missing values in object

Running the model alone (for just one variable) it worked properly
I attach however the database with a extra observations


structure(list(paciente = structure(c(6195, 6149, 6133, 6457, 
6248, 6368, 6001, 6533, 6470, 6506, 6060, 6438, 6394, 6259, 6531, 
6378, 6551, 6227, 6138, 6400, 6040, 6279, 6316, 6137, 6274, 6432, 
6206, 6228, 6550, 6326), label = "Paciente", format.spss = "F6.0"), 
    edad_s1 = structure(c(63, 60, 65, 61, 68, 71, 72, 58, 67, 
    69, 66, 61, 68, 65, 64, 64, 72, 69, 62, 71, 69, 72, 68, 65, 
    70, 71, 63, 71, 63, 70), label = "Edad", format.spss = "F3.0"), 
    sexo_s1 = structure(c(1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 
    2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 
    2L, 1L, 1L, 1L, 1L, 2L), .Label = c("Hombre", "Mujer"), label = "Sexo", class = "factor"), 
    grupo_int_v00 = structure(c(1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 
    2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 
    2L, 2L, 2L, 2L, 1L, 2L, 2L), .Label = c("A", "B"), label = "Grupo de intervención", class = "factor"), 
    time = c(0, 2, 0, 2, 0, 0, 2, 1, 0, 2, 0, 1, 1, 0, 2, 0, 
    2, 2, 1, 2, 2, 2, 1, 2, 0, 1, 1, 2, 1, 1), peso1 = c(127.6, 
    79.6, 70.1, 95, 85.5, 84.1, 78.3, 131.3, 92.6, 72, 91.9, 
    79.6, 81, 81.5, 71, 98.7, 88.3, 84.8, 82, 85, 73.7, 86.9, 
    66, 65.2, 90.3, 72.8, 87.5, 103.3, 82.2, 59.1), cintura1 = c(128, 
    103.5, 103, 103.5, 110, 97.5, 110, 128, 115, 106.5, 117, 
    109, 100, 107, 98, 113.5, 106, 103.5, 108.5, 103.5, 104, 
    104.5, 105, 98, 115, 96, 102.5, 120.5, 95.5, 95), tasis2_e = c(146, 
    162, 139, 131, 156, 139, 141, 146, 167, 148, 139, 135, 128, 
    164, 134, 145, 116, 138, 124, 114, 132, 146, 94, 130, 161, 
    142, 131, 131, 115, 138), tadias2_e = c(73, 96, 73, 80, 82, 
    78, 63, 73, 76, 77, 68, 77, 78, 81, 68, 71, 65, 72, 74, 53, 
    65, 69, 60, 80, 78, 72, 77, 65, 57, 75), p17_total = c(4, 
    9, 7, 5, 11, 9, 13, 13, 9, 7, 3, 11, 13, 8, 9, 11, 10, 11, 
    11, 12, 10, 10, 14, 14, 7, 15, 7, 8, 15, 12), geaf_tot = c(3356.64, 
    3202.8, 2531.47, 1230.77, 1706.29, 223.78, 6965.03, 2587.41, 
    839.16, 3356.64, 1678.32, 1006.99, 2909.09, 4461.54, 2895.1, 
    1678.32, 4979.02, 3020.98, 2657.34, 4090.91, 4615.38, 9328.67, 
    1958.04, 4727.27, 5687.65, 7776.22, 979.02, 4876.46, 5361.31, 
    1258.74), glucosa = c(99, 104, 135, 95, 103, 156, 132, 145, 
    119, 98, 93, 130, 133, 106, 119, 119, 98, 103, 126, 95, 101, 
    106, 111, 102, 137, 102, 135, 101, 109, 116), albumi = c(4.83, 
    5.25, 4.84, 4.78, 4.74, 4.74, 4.58, 4.72, 4.59, 4.32, 4.59, 
    4.85, 4.38, 4.25, 5.1, 4.52, 4.57, 4.91, 4.78, 4.17, 4.59, 
    4.54, 4.73, 4.65, 4.76, 4.62, 4.42, 4.64, 4.33, 4.65), coltot = c(260, 
    217, 263, 276, 241, 275, 195, 132, 248, 278, 178, 188, 201, 
    232, 254, 188, 175, 248, 241, 260, 206, 195, 321, 275, 325, 
    157, 172, 211, 187, 237), hdl = c(50, 39, 57, 57, 50, 72, 
    86, 43, 49, 75, 43, 63, 70, 43, 70, 40, 42, 50, 52, 57, 63, 
    59, 56, 65, 46, 56, 46, 66, 60, 48), ldl_calc = c(167, 134, 
    182, 204, 153, 185, 93, 67, 155, 184, 101, 103, 118, 164, 
    160, 122, 121, 157, 154, 171, 112, 122, NA, 183, 246, 84, 
    97, 127, 105, 147), trigli = c(216, 220, 122, 74, 189, 88, 
    80, 109, 219, 94, 169, 112, 64, 124, 120, 129, 60, 203, 173, 
    161, 155, 68, 311, 137, 165, 86, 143, 89, 112, 208), hba1c = c(5.81, 
    5.97, NA, 5.68, 6.49, 6.27, 6.28, 6.44, 6.11, 5.57, NA, 6.66, 
    5.8, 6.03, 5.97, 5.87, 5.71, 6.09, 5.83, 5.78, 5.87, 5.53, 
    5.9, 5.98, 6.98, 5.62, 6.91, 5.69, 5.48, 6.22), i_hucpeptide = c(764.89, 
    1528.01, NA, 466.64, 1131.82, 1928.51, NA, 1648.74, 847.89, 
    314.53, NA, 975.24, 667.79, 830.84, 990.65, 1351.61, 545.1, 
    1307.36, NA, 138.36, NA, 578.8, 286.46, NA, 813.83, 832.27, 
    1638.12, 649.19, 819.65, 861.62), i_hughrelin = c(276.53, 
    770.63, NA, 410.97, 763.9, 1476.06, NA, 1338.17, 453, 1101.59, 
    NA, 841.87, 1361.9, 629.31, 583.35, 424.99, 442.24, 478.57, 
    NA, 397.89, NA, 313.16, 1603.63, NA, 1064.98, 1097.91, 230.18, 
    456.28, 527.01, 1639.29), i_hugip = c(2.67, 2.67, NA, 2.67, 
    2.67, 2.67, NA, 2.67, 2.67, 2.67, NA, 2.67, 2.67, 2.67, 2.67, 
    2.67, 2.67, 2.67, NA, 2.67, NA, 2.67, 2.67, NA, 2.67, 2.67, 
    2.67, 2.67, 2.67, 2.67), i_huglp1 = c(162.96, 350.07, NA, 
    127.66, 118.74, 193.61, NA, 186.67, 200.13, 65.09, NA, 183.17, 
    64.39, 114.85, 14.14, 45.47, 113.69, 14.14, NA, 43.04, NA, 
    14.14, 14.14, NA, 14.14, 14.14, 414.37, 141.23, 234.84, 95.8
    ), i_huglucagon = c(543.93, 483.91, NA, 333.79, 470, 563.09, 
    NA, 376.14, 726.99, 469.21, NA, 532.28, 374.23, 471.29, 410.69, 
    549.54, 239.55, 473.43, NA, 434.59, NA, 241.55, 151.01, NA, 
    277.81, 375.12, 616.31, 186.88, 610.5, 403.96), i_huinsulin = c(191.5, 
    338.31, NA, 129.24, 324.97, 591.08, NA, 583.01, 299.75, 58.99, 
    NA, 295.51, 150.46, 301.65, 203.19, 295.91, 64.12, 326.06, 
    NA, 121.68, NA, 126.33, 73.11, NA, 154.8, 169.85, 810.69, 
    232.17, 267.54, 176.84), i_huleptin = c(14162.44, 4831.07, 
    NA, 3992.49, 9441.46, 10447.88, NA, 6989.59, 8409.76, 7571.15, 
    NA, 6537.04, 9322.29, 9754.34, 7940.63, 3.9, 4053.58, 10703.55, 
    NA, 9771.57, NA, 4261.31, 4194.34, NA, 5613.27, 1944.27, 
    2521.25, 2160.27, 2965.92, 6397.05), i_hupai1 = c(3478.99, 
    3593.14, NA, 997.29, 2506.14, 1564.74, NA, 2523, 3085.25, 
    2061.57, NA, 4635.7, 1952.29, 2415.79, 2093.76, 2956.08, 
    994.61, 4735.39, NA, 3329.24, NA, 1215.78, 944.76, NA, 1327.49, 
    1638.69, 1603.74, 1625.27, 1730.55, 2228.95), i_huresistin = c(6391.27, 
    9663.25, NA, 3044.48, 2529.7, 6077.57, NA, 4043.41, 3221.72, 
    2866.02, NA, 7641.27, 2871.47, 3652.36, 2930.08, 4435.12, 
    3659.3, 2746.74, NA, 4855.17, NA, 2187.33, 2536.88, NA, 1156.05, 
    3998.54, 5520.42, 2994.15, 5170.95, 2932.3), i_huvisfatin = c(2277.03, 
    452.48, NA, 302.3, 1684.98, 1497.28, NA, 627.63, 2989.72, 
    1791.18, NA, 1628.96, 1144.52, 1341.28, 1487.32, 2155, 390.72, 
    1664.39, NA, 1367.63, NA, 334.56, 284.77, NA, 384.15, 8.64, 
    8.64, 8.64, 8.64, 1205.34), col_rema = c(43, 44, 24, 15, 
    38, 18, 16, 22, 44, 19, 34, 22, 13, 25, 24, 26, 12, 41, 35, 
    32, 31, 14, NA, 27, 33, 17, 29, 18, 22, 42), homa = c(7.802, 
    14.479, NA, 5.053, 13.774, 37.946, NA, 34.789, 14.679, 2.379, 
    NA, 15.809, 8.235, 13.158, 9.95, 14.491, 2.586, 13.821, NA, 
    4.757, NA, 5.511, 3.34, NA, 8.727, 7.13, 45.038, 9.65, 12.001, 
    8.442)), row.names = c(NA, -30L), class = c("tbl_df", "tbl", 
"data.frame"))

do you encounter this error when you run not your real script, but specifically your reprex as you provided it? I ask because I don't get that error, perhaps the reprex doesnt capture the relevant aspects of what you are doing...
that said, there may be a subtletly to fix. Given the differences between tibbles and dataframes I would probably change subset = !is.na(data[var_1]) to subset = !is.na(data[[var_1]])

The reprex, running the reprex I got that error.

I have made the modification. Right now what I encounter is:

Error in lme.formula(as.formula(paste0(var_1, " ~ sexo_s1*peso1 + edad_s1 + p17_total + poly(time, 2)*grupo_int_v00")), : fewer observations than random effects in all level 1 groups

Not pretty sure I'm missing sthg

That message is describing that your model is too complex for the data.
Does it go away when you run it on your full data rather than this smaller sample ?

With my full data I get this instead

[1] "now paciente"
Error in MEestimate(lmeSt, grps) : 
  Singularity in backsolve at level 0, block 1

I believe we are away from a looping issue and into more purely an lme issue.
I'm not a frequent lme user, so I'll be unlikely to help you further in this matter.

Thanks for your help, this one and the previous. I am looking for help now

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.