Looping/sapply through nlme function

General
,
1

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

2

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))
}
3

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"))
4

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]])

5

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

6

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 ?

7

With my full data I get this instead

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

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.

9

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

10

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.