Loop code doesn't work properly

Javier9 · November 16, 2021, 11:05am

Hi, we developed a loop to obtain terciles of several variables. My ultimate goal is to create terciles according to one variable and then compare (using t.test of independent groups) the terciles created for the rest of the variables. The problem is I have too many variables to perform it manually, so I developed a loop to iterate through variables, creating terciles of all variables and comparing all terciles of the variables. I have done it manually with several variables, and the results dont`t match with those obtained by the loop code. I have gone through it many times, consult the statician of my research group, but we don't know what is happening

The small piece of the database

sample<-structure(list(paciente = structure(c(6363, 6052, 6519, 6371, 
6555, 6185, 6002, 6155, 6287, 6217), format.spss = "F5.0"), sexo_s1 = structure(c(1L, 
2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L), .Label = c("Hombre", "Mujer"
), label = "Sexo", class = "factor"), edad_s1 = structure(c(66, 
63, 72, 64, 72, 64, 70, 73, 63, 65), label = "Edad", format.spss = "F3.0"), 
    peso1_v00 = structure(c(84.4, 76.2, 88.3, 122, 72, 86.4, 
    78, 79, 65.5, 98.8), label = "Peso: 1a determinación", format.spss = "F5.1"), 
    cintura1_v00 = structure(c(104.5, 101, 107.5, 128.5, 107, 
    109.5, 105.5, 109, 97, 118.5), label = "Cintura: 1a determinación", format.spss = "F5.1"), 
    tasis2_e_v00 = structure(c(155, 131, 148, 147, 136, 154, 
    130, 154, 147, 139), label = "TA: tensión arterial 2: sistólica", format.spss = "F4.0"), 
    tadias2_e_v00 = structure(c(98, 76, 83, 84, 73, 80, 64, 80, 
    82, 78), label = "TA: tensión arterial 2: diastólica", format.spss = "F4.0"), 
    p17_total_v00 = structure(c(10, 8, 10, 10, 10, 5, 11, 9, 
    8, 5), label = "Cuestionario de 17 puntos: Suma de puntuación de P17", format.spss = "F3.0"), 
    geaf_tot_v00 = structure(c(1048.95, 4195.8, 4615.38, 3356.64, 
    839.16, 1608.39, 2958.04, 10209.79, 4335.66, 157.34), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2"), 
    glucosa_v00 = structure(c(97, 122, 109, 201, 143, 139, 95, 
    110, 101, 181), label = "Analítica: Glucosa en mg/dL", format.spss = "F4.0"), 
    albumi_v00 = structure(c(4.61, 4.52, 4.44, 4.48, 4.75, 4.87, 
    4.66, 4.86, 4.75, 4.99), label = "Analítica: Albúmina en g/dL", format.spss = "F6.2"), 
    coltot_v00 = structure(c(221, 218, 130, 261, 190, 221, 199, 
    185, 233, 232), label = "Analítica: Colesterol total en mg/dL", format.spss = "F4.0"), 
    hdl_v00 = structure(c(52, 54, 43, 42, 50, 67, 90, 50, 43, 
    51), label = "Analítica: Colesterol HDL en mg/dL", format.spss = "F4.0"), 
    ldl_calc_v00 = structure(c(148, 135, 75, 173, 128, 133, 84, 
    109, 134, 144), label = "Analítica: LDL calculado en mg/dL si trigli<=300", format.spss = "F4.0"), 
    trigli_v00 = structure(c(103, 144, 58, 232, 60, 106, 126, 
    131, 282, 186), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    hba1c_v00 = structure(c(5.87, NA, 5.76, 7.98, 6.38, 7.01, 
    NA, 6.51, 5.95, 9.01), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    peso1_v66 = structure(c(76.6, 74.2, 82.2, 115, 64, 73.5, 
    74.4, 75.5, 59.5, 100), label = "Peso: 1a determinación", format.spss = "F5.1"), 
    cintura1_v66 = structure(c(97.5, 99, 100, 122.5, 99, 101, 
    97, 104.5, 82.5, 119.5), label = "Cintura: 1a determinación", format.spss = "F5.1"), 
    tasis2_e_v66 = structure(c(133, 129, 144, 160, 122, 123, 
    130, 155, 139, 153), label = "TA: tensión arterial 2: sistólica", format.spss = "F4.0"), 
    tadias2_e_v66 = structure(c(87, 75, 77, 86, 68, 53, 64, 81, 
    73, 89), label = "TA: tensión arterial 2: diastólica", format.spss = "F4.0", display_width = 13L), 
    p17_total_v66 = structure(c(12, 12, 12, 13, 11, 15, 14, 9, 
    15, 12), label = "Cuestionario de 17 puntos: Suma total de P17", format.spss = "F3.0"), 
    geaf_tot_v66 = structure(c(4615.38, 0, 3461.54, 3356.64, 
    2097.9, 5132.87, 5118.88, 8111.89, 6772.03, 209.79), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2"), 
    glucosa_v66 = structure(c(96, 107, 111, 173, 107, 98, 88, 
    118, 97, 185), label = "Analítica: Glucosa en mg/dL", format.spss = "F4.0"), 
    albumi_v66 = structure(c(4.5, 4.93, 4.41, 4.42, 4.82, 4.54, 
    4.8, 4.46, 4.64, 4.84), label = "Analítica: Albúmina en g/dL", format.spss = "F6.2"), 
    coltot_v66 = structure(c(215, 226, 156, 235, 154, 210, 283, 
    182, 225, 171), label = "Analítica: Colesterol total en mg/dL", format.spss = "F4.0"), 
    hdl_v66 = structure(c(54, 65, 54, 40, 51, 73, 88, 58, 46, 
    37), label = "Analítica: Colesterol HDL en mg/dL", format.spss = "F4.0"), 
    ldl_calc_v66 = structure(c(147, 133, 89, 156, 94, 123, 175, 
    106, 137, 102), label = "Analítica: LDL calculado en mg/dL si trigli<=300", format.spss = "F4.0"), 
    trigli_v66 = structure(c(72, 138, 63, 197, 47, 72, 99, 89, 
    209, 160), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    hba1c_v66 = structure(c(5.67, NA, 5.54, 8.05, 5.5, 5.95, 
    NA, 6.17, 5.75, 8.98), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    peso1_v01 = structure(c(75.2, 72, 82.4, 116, 63, 73.5, 72.4, 
    79, 58.5, 100), label = "Peso: 1a determinación", format.spss = "F5.1"), 
    cintura1_v01 = structure(c(97.5, 98, 100, 122.5, 99, 101, 
    96, 106, 88.5, 119.5), label = "Cintura: 1a determinación", format.spss = "F5.1"), 
    tasis2_e_v01 = structure(c(130, 137, 137, 183, 122, 152, 
    138, 139, 147, 154), label = "TA: tensión arterial 2: sistólica", format.spss = "F4.0"), 
    tadias2_e_v01 = structure(c(84, 75, 72, 97, 71, 63, 72, 72, 
    60, 80), label = "TA: tensión arterial 2: diastólica", format.spss = "F4.0"), 
    p17_total_v01 = structure(c(14, 11, 12, 10, 14, 15, 14, 11, 
    13, 12), label = "Cuestionario de 17 puntos: Suma total de P17", format.spss = "F3.0"), 
    geaf_tot_v01 = structure(c(1678.32, 0, 4713.29, 559.44, 2769.23, 
    3212.12, 6853.15, 3776.22, 5841.49, 1048.95), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2"), 
    glucosa_v01 = structure(c(93, 116, 100, 200, 112, 109, 105, 
    118, 94, 242), label = "Analítica: Glucosa en mg/dL", format.spss = "F4.0"), 
    albumi_v01 = structure(c(4.57, 4.42, 4.68, 4.36, 4.98, 4.74, 
    4.87, 4.8, 4.81, 4.88), label = "Analítica: Albúmina en g/dL", format.spss = "F6.2"), 
    coltot_v01 = structure(c(198, 236, 158, 270, 181, 187, 213, 
    204, 226, 192), label = "Analítica: Colesterol total en mg/dL", format.spss = "F4.0"), 
    hdl_v01 = structure(c(58, 59, 48, 74, 60, 60, 87, 52, 49, 
    37), label = "Analítica: Colesterol HDL en mg/dL", format.spss = "F4.0"), 
    ldl_calc_v01 = structure(c(128, 160, 99, 168, 107, 109, 105, 
    130, 147, 125), label = "Analítica: LDL calculado en mg/dL si trigli<=300", format.spss = "F4.0"), 
    trigli_v01 = structure(c(62, 83, 53, 139, 71, 90, 105, 110, 
    148, 151), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    hba1c_v01 = structure(c(5.61, 6.48, 5.53, 8.26, 5.86, 6.4, 
    5.26, 6.5, 5.7, 9.62), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    i_hucpeptide_v00 = structure(c(704.96, NA, 675.43, 913.16, 
    1325.4, 840.98, NA, 1932.23, 459.83, 422.15), label = "Hu C-peptide (72) IMIM S'han substituit en les següents var els codis de inf i sup a limit de detecció per el limit inf i sup de detecció", format.spss = "F9.2", display_width = 13L), 
    i_hucpeptide_v66 = structure(c(510.07, NA, 824.39, 926.34, 
    1199.51, 488.01, NA, 1461.11, 346.92, 679.09), label = "Hu C-peptide (72) IMIM", format.spss = "F9.2", display_width = 13L), 
    i_hucpeptide_v01 = structure(c(432.38, NA, 737.66, 707.55, 
    1057.83, 699.08, NA, 1512.69, 345.69, 356.67), label = "Hu C-peptide (72) IMIM", format.spss = "F9.2", display_width = 12L), 
    i_hughrelin_v00 = structure(c(1823.83, NA, 1050.11, 424.25, 
    198.06, 534.27, NA, 709.73, 117.69, 420.37), label = "Hu Ghrelin (26) IMIM", format.spss = "F7.2", display_width = 10L), 
    i_hughrelin_v66 = structure(c(1407.04, NA, 746.93, 423.13, 
    207.63, 464.17, NA, 728.57, 113.23, 463.65), label = "Hu Ghrelin (26) IMIM", format.spss = "F7.2", display_width = 13L), 
    i_hughrelin_v01 = structure(c(1133.96, NA, 670.23, 405.11, 
    260.24, 418.79, NA, 533.78, 122.07, 220.22), label = "Hu Ghrelin (26) IMIM", format.spss = "F7.2", display_width = 11L), 
    i_hugip_v00 = structure(c(2.67, NA, 2.67, 2.67, 2.67, 2.67, 
    NA, 2.67, 2.67, 2.67), label = "Hu GIP (14) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_hugip_v66 = structure(c(2.67, NA, 2.67, 2.67, 256.21, 2.67, 
    NA, 2.67, 2.67, 2.67), label = "Hu GIP (14) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_hugip_v01 = structure(c(2.67, NA, 2.67, 24.74, 1165.16, 
    2.67, NA, 2.67, 2.67, 2.67), label = "Hu GIP (14) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huglp1_v00 = structure(c(14.14, NA, 14.14, 216.2, 116.16, 
    228.14, NA, 359.48, 14.14, 219.02), label = "Hu GLP-1 (27) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huglp1_v66 = structure(c(14.14, NA, 14.14, 202.5, 92.49, 
    278.55, NA, 400.16, 14.14, 274.62), label = "Hu GLP-1 (27) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huglp1_v01 = structure(c(14.14, NA, 14.14, 202.5, 56.41, 
    303.77, NA, 451.91, 14.14, 58.55), label = "Hu GLP-1 (27) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huglucagon_v00 = structure(c(273.5, NA, 231.41, 491.4, 
    542.39, 489.35, NA, 525.39, 241.55, 362.72), label = "Hu Glucagon (15) IMIM", format.spss = "F9.2", display_width = 11L), 
    i_huglucagon_v66 = structure(c(94.97, NA, 276.26, 456.39, 
    306.42, 288.62, NA, 523.49, 143.14, 440.22), label = "Hu Glucagon (15) IMIM", format.spss = "F9.2", display_width = 11L), 
    i_huglucagon_v01 = structure(c(182.57, NA, 233.95, 522.65, 
    10.48, 409.99, NA, 511.95, 216.47, 326.51), label = "Hu Glucagon (15) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huinsulin_v00 = structure(c(97.94, NA, 171.33, 286.66, 
    390.29, 221.06, NA, 668.14, 125.36, 349.99), label = "Hu Insulin (12) IMIM", format.spss = "F7.2", display_width = 10L), 
    i_huinsulin_v66 = structure(c(64.73, NA, 236.01, 255.05, 
    284.71, 73.24, NA, 392.19, 64.75, 381.42), label = "Hu Insulin (12) IMIM", format.spss = "F7.2", display_width = 9L), 
    i_huinsulin_v01 = structure(c(52.52, NA, 213.4, 703.51, 247.48, 
    147.22, NA, 521.07, 98.66, 298.72), label = "Hu Insulin (12) IMIM", format.spss = "F7.2"), 
    i_huleptin_v00 = structure(c(3493.7, NA, 1965.55, 4767.39, 
    5122.91, 12320.55, NA, 5367.22, 5217.35, 4682.33), label = "Hu Leptin (78) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huleptin_v66 = structure(c(1779.33, NA, 1410.07, 3977.58, 
    3645.73, 3608.76, NA, 4489.67, 3499.88, 5136.43), label = "Hu Leptin (78) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_huleptin_v01 = structure(c(1865.4, NA, 1312.5, 8371.63, 
    2128.98, 6921.89, NA, 3754.42, 3092.53, 3921.64), label = "Hu Leptin (78) IMIM", format.spss = "F9.2", display_width = 9L), 
    i_hupai1_v00 = structure(c(1581.8, NA, 1442.88, 3209.36, 
    2349.08, 3202.1, NA, 3177.94, 1463.04, 1701.4), label = "Hu PAI-1 (61) IMIM", format.spss = "F7.2"), 
    i_hupai1_v66 = structure(c(1625.45, NA, 1093.24, 2152.24, 
    2083.31, 982.35, NA, 2245.81, 1611.27, 2645.45), label = "Hu PAI-1 (61) IMIM", format.spss = "F7.2"), 
    i_hupai1_v01 = structure(c(1726.35, NA, 1166.53, 2511.45, 
    2268.52, 1592.08, NA, 2560.71, 1936.33, 2500.51), label = "Hu PAI-1 (61) IMIM", format.spss = "F7.2"), 
    i_huresistin_v00 = structure(c(4292.62, NA, 3951.76, 4101.48, 
    6430.17, 5599.94, NA, 4855.32, 2144.19, 2421.1), label = "Hu Resistin (65) IMIM", format.spss = "F8.2", display_width = 9L), 
    i_huresistin_v66 = structure(c(3201.72, NA, 4774.83, 4500.78, 
    7574.37, 4403.32, NA, 3224.09, 2102.65, 2003.5), label = "Hu Resistin (65) IMIM", format.spss = "F8.2", display_width = 7L), 
    i_huresistin_v01 = structure(c(3872.84, NA, 4595.27, 3581.62, 
    9521.7, 4225.4, NA, 3150.62, 2093.1, 2048.76), label = "Hu Resistin (65) IMIM", format.spss = "F8.2", display_width = 6L), 
    i_huvisfatin_v00 = structure(c(8.64, NA, 2.06, 560.32, 1498.58, 
    1356.01, NA, 632.07, 315.62, 461.86), label = "Hu Visfatin (22) IMIM", format.spss = "F9.2", display_width = 6L), 
    i_huvisfatin_v66 = structure(c(8.64, NA, 8.64, 472.64, 683.91, 
    8.64, NA, 486.94, 8.64, 477.56), label = "Hu Visfatin (22) IMIM", format.spss = "F9.2", display_width = 6L), 
    i_huvisfatin_v01 = structure(c(8.64, NA, 8.64, 2113.05, 1415.55, 
    800.08, NA, 155.08, 8.64, 108), label = "Hu Visfatin (22) IMIM", format.spss = "F9.2", display_width = 10L), 
    col_rema_v00 = structure(c(21, 29, 12, 46, 12, 21, 25, 26, 
    56, 37), format.spss = "F8.2", display_width = 14L), col_rema_v66 = structure(c(14, 
    28, 13, 39, 9, 14, 20, 18, 42, 32), format.spss = "F8.2", display_width = 14L), 
    col_rema_v01 = structure(c(12, 17, 11, 28, 14, 18, 21, 22, 
    30, 30), format.spss = "F8.2", display_width = 14L), homa_v00 = structure(c(422.230222222222, 
    NA, 829.998666666667, 2560.82933333333, 2480.50977777778, 
    1365.65955555556, NA, 3266.46222222222, 562.727111111111, 
    2815.47511111111), format.spss = "F8.2", display_width = 10L), 
    homa_v66 = structure(c(276.181333333333, NA, 1164.316, 1961.05111111111, 
    1353.95422222222, 319.000888888889, NA, 2056.81866666667, 
    279.144444444444, 3136.12), format.spss = "F8.2", display_width = 10L), 
    homa_v01 = structure(c(217.082666666667, NA, 948.444444444444, 
    6253.42222222222, 1231.90044444444, 713.199111111111, NA, 
    2732.72266666667, 412.179555555556, 3212.89955555556), format.spss = "F8.2", display_width = 10L), 
    d_homa_v66 = structure(c(-146.048888888889, NA, 334.317333333333, 
    -599.778222222222, -1126.55555555556, -1046.65866666667, 
    NA, -1209.64355555556, -283.582666666667, 320.644888888889
    ), format.spss = "F8.2", display_width = 12L), d_homa_v01 = structure(c(-205.147555555556, 
    NA, 118.445777777778, 3692.59288888889, -1248.60933333333, 
    -652.460444444444, NA, -533.739555555555, -150.547555555556, 
    397.424444444444), format.spss = "F8.2", display_width = 12L), 
    d_hughrelin_v66 = structure(c(-416.79, NA, -303.18, -1.12, 
    9.56999999999999, -70.1, NA, 18.84, -4.45999999999999, 43.28
    ), format.spss = "F8.2", display_width = 18L), d_hughrelin_v01 = structure(c(-689.87, 
    NA, -379.88, -19.14, 62.18, -115.48, NA, -175.95, 4.38, -200.15
    ), format.spss = "F8.2", display_width = 18L), d_huinsulin_v66 = structure(c(-33.21, 
    NA, 64.68, -31.61, -105.58, -147.82, NA, -275.95, -60.61, 
    31.43), format.spss = "F8.2", display_width = 17L), d_huinsulin_v01 = structure(c(-45.42, 
    NA, 42.07, 416.85, -142.81, -73.84, NA, -147.07, -26.7, -51.27
    ), format.spss = "F8.2", display_width = 17L), d_hucpeptide_v66 = structure(c(-194.89, 
    NA, 148.96, 13.1800000000001, -125.89, -352.97, NA, -471.12, 
    -112.91, 256.94), format.spss = "F8.2", display_width = 18L), 
    d_hucpeptide_v01 = structure(c(-272.58, NA, 62.23, -205.61, 
    -267.57, -141.9, NA, -419.54, -114.14, -65.48), format.spss = "F8.2", display_width = 18L), 
    d_huglucagon_v66 = structure(c(-178.53, NA, 44.85, -35.01, 
    -235.97, -200.73, NA, -1.89999999999998, -98.41, 77.5), format.spss = "F8.2", display_width = 18L), 
    d_huglucagon_v01 = structure(c(-90.93, NA, 2.53999999999999, 
    31.25, -531.91, -79.36, NA, -13.44, -25.08, -36.21), format.spss = "F8.2", display_width = 18L), 
    d_huleptin_v66 = structure(c(-1714.37, NA, -555.48, -789.81, 
    -1477.18, -8711.79, NA, -877.55, -1717.47, 454.1), format.spss = "F8.2", display_width = 16L), 
    d_huleptin_v01 = structure(c(-1628.3, NA, -653.05, 3604.24, 
    -2993.93, -5398.66, NA, -1612.8, -2124.82, -760.69), format.spss = "F8.2", display_width = 16L), 
    d_huresistin_v66 = structure(c(-1090.9, NA, 823.07, 399.3, 
    1144.2, -1196.62, NA, -1631.23, -41.54, -417.6), format.spss = "F8.2", display_width = 18L), 
    d_huresistin_v01 = structure(c(-419.78, NA, 643.51, -519.86, 
    3091.53, -1374.54, NA, -1704.7, -51.0900000000001, -372.34
    ), format.spss = "F8.2", display_width = 18L), d_huvisfatin_v66 = structure(c(0, 
    NA, 6.58, -87.6800000000001, -814.67, -1347.37, NA, -145.13, 
    -306.98, 15.7), format.spss = "F8.2", display_width = 18L), 
    d_huvisfatin_v01 = structure(c(0, NA, 6.58, 1552.73, -83.03, 
    -555.93, NA, -476.99, -306.98, -353.86), format.spss = "F8.2", display_width = 18L), 
    d_glucosa_v66 = structure(c(-1, -15, 2, -28, -36, -41, -7, 
    8, -4, 4), format.spss = "F8.2", display_width = 15L), d_glucosa_v01 = structure(c(-4, 
    -6, -9, -1, -31, -30, 10, 8, -7, 61), format.spss = "F8.2", display_width = 15L), 
    d_coltot_v66 = structure(c(-6, 8, 26, -26, -36, -11, 84, 
    -3, -8, -61), format.spss = "F8.2", display_width = 14L), 
    d_coltot_v01 = structure(c(-23, 18, 28, 9, -9, -34, 14, 19, 
    -7, -40), format.spss = "F8.2", display_width = 14L), d_hdl_v66 = structure(c(2, 
    11, 11, -2, 1, 6, -2, 8, 3, -14), format.spss = "F8.2", display_width = 11L), 
    d_hdl_v01 = structure(c(6, 5, 5, 32, 10, -7, -3, 2, 6, -14
    ), format.spss = "F8.2", display_width = 11L), d_ldl_calc_v66 = structure(c(-1, 
    -2, 14, -17, -34, -10, 91, -3, 3, -42), format.spss = "F8.2", display_width = 16L), 
    d_ldl_calc_v01 = structure(c(-20, 25, 24, -5, -21, -24, 21, 
    21, 13, -19), format.spss = "F8.2", display_width = 16L), 
    d_col_rema_v66 = structure(c(-7, -1, 1, -7, -3, -7, -5, -8, 
    -14, -5), format.spss = "F8.2", display_width = 16L), d_col_rema_v01 = structure(c(-9, 
    -12, -1, -18, 2, -3, -4, -4, -26, -7), format.spss = "F8.2", display_width = 16L), 
    d_trigli_v66 = structure(c(-31, -6, 5, -35, -13, -34, -27, 
    -42, -73, -26), format.spss = "F8.2", display_width = 14L), 
    d_trigli_v01 = structure(c(-41, -61, -5, -93, 11, -16, -21, 
    -21, -134, -35), format.spss = "F8.2", display_width = 14L), 
    d_hba1c_v66 = structure(c(-0.2, NA, -0.22, 0.0700000000000003, 
    -0.88, -1.06, NA, -0.34, -0.2, -0.0299999999999994), format.spss = "F8.2", display_width = 13L), 
    d_hba1c_v01 = structure(c(-0.26, NA, -0.23, 0.279999999999999, 
    -0.52, -0.609999999999999, NA, -0.00999999999999979, -0.25, 
    0.609999999999999), format.spss = "F8.2", display_width = 13L), 
    d_tasis2_e_v66 = structure(c(-22, -2, -4, 13, -14, -31, 0, 
    1, -8, 14), format.spss = "F8.2", display_width = 16L), d_tasis2_e_v01 = structure(c(-25, 
    6, -11, 36, -14, -2, 8, -15, 0, 15), format.spss = "F8.2", display_width = 16L), 
    d_tadias2_e_v66 = structure(c(-11, -1, -6, 2, -5, -27, 0, 
    1, -9, 11), format.spss = "F8.2", display_width = 17L), d_tadias2_e_v01 = structure(c(-14, 
    -1, -11, 13, -2, -17, 8, -8, -22, 2), format.spss = "F8.2", display_width = 17L), 
    d_peso1_v66 = structure(c(-7.80000000000001, -2, -6.09999999999999, 
    -7, -8, -12.9, -3.59999999999999, -3.5, -6, 1.2), format.spss = "F8.2", display_width = 13L), 
    d_peso1_v01 = structure(c(-9.2, -4.2, -5.89999999999999, 
    -6, -9, -12.9, -5.59999999999999, 0, -7, 1.2), format.spss = "F8.2", display_width = 13L), 
    d_cintura1_v66 = structure(c(-7, -2, -7.5, -6, -8, -8.5, 
    -8.5, -4.5, -14.5, 1), format.spss = "F8.2", display_width = 16L), 
    d_cintura1_v01 = structure(c(-7, -3, -7.5, -6, -8, -8.5, 
    -9.5, -3, -8.5, 1), format.spss = "F8.2", display_width = 16L), 
    d_geaf_tot_v66 = structure(c(3566.43, -4195.8, -1153.84, 
    0, 1258.74, 3524.48, 2160.84, -2097.9, 2436.37, 52.45), format.spss = "F8.2", display_width = 16L), 
    d_geaf_tot_v01 = structure(c(629.37, -4195.8, 97.9099999999999, 
    -2797.2, 1930.07, 1603.73, 3895.11, -6433.57, 1505.83, 891.61
    ), format.spss = "F8.2", display_width = 16L), d_p17_total_v66 = structure(c(2, 
    4, 2, 3, 1, 10, 3, 0, 7, 7), format.spss = "F8.2", display_width = 11L), 
    d_p17_total_v01 = structure(c(4, 3, 2, 0, 4, 10, 3, 2, 5, 
    7), format.spss = "F8.2"), d_hupai1_v66 = structure(c(43.6500000000001, 
    NA, -349.64, -1057.12, -265.77, -2219.75, NA, -932.13, 148.23, 
    944.05), format.spss = "F8.2", display_width = 13L), d_hupai1_v01 = structure(c(144.55, 
    NA, -276.35, -697.91, -80.5599999999999, -1610.02, NA, -617.23, 
    473.29, 799.11), format.spss = "F8.2", display_width = 13L), 
    d_hugip_v66 = structure(c(0, NA, 0, 0, 253.54, 0, NA, 0, 
    0, 0), format.spss = "F8.2", display_width = 13L), d_hugip_v01 = structure(c(0, 
    NA, 0, 22.07, 1162.49, 0, NA, 0, 0, 0), format.spss = "F8.2", display_width = 13L), 
    d_huglp1_v66 = structure(c(0, NA, 0, -13.7, -23.67, 50.41, 
    NA, 40.68, 0, 55.6), format.spss = "F8.2", display_width = 13L), 
    d_huglp1_v01 = structure(c(0, NA, 0, -13.7, -59.75, 75.63, 
    NA, 92.43, 0, -160.47), format.spss = "F8.2", display_width = 13L), 
    ln_trigli_v00 = structure(c(4.63472898822964, 4.969813299576, 
    4.06044301054642, 5.44673737166631, 4.0943445622221, 4.66343909411207, 
    4.83628190695148, 4.87519732320115, 5.64190707093811, 5.2257466737132
    ), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    ln_trigli_v66 = structure(c(4.27666611901606, 4.92725368515721, 
    4.14313472639153, 5.28320372873799, 3.85014760171006, 4.27666611901606, 
    4.59511985013459, 4.48863636973214, 5.34233425196481, 5.07517381523383
    ), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    ln_trigli_v01 = structure(c(4.12713438504509, 4.4188406077966, 
    3.97029191355212, 4.93447393313069, 4.26267987704132, 4.49980967033027, 
    4.65396035015752, 4.70048036579242, 4.99721227376411, 5.01727983681492
    ), label = "Analítica: Triglicéridos en mg/dL", format.spss = "F5.0"), 
    ln_homa_v00 = structure(c(6.04555071556718, NA, 6.72142409436365, 
    7.84808644334377, 7.81621937359094, 7.21939278180313, NA, 
    8.09146278899802, 6.33279480567865, 7.94288630470217), format.spss = "F8.2", display_width = 10L), 
    ln_homa_v66 = structure(c(5.62105765481488, NA, 7.05988906911122, 
    7.58123588965634, 7.21078464361132, 5.76519388926651, NA, 
    7.62891573149717, 5.63172936987266, 8.05074164604432), format.spss = "F8.2", display_width = 10L), 
    ln_homa_v01 = structure(c(5.38027823337748, NA, 6.85482321564537, 
    8.74088414843216, 7.11631333274751, 6.56976063964933, NA, 
    7.91305370498932, 6.02145906886522, 8.07492909673972), format.spss = "F8.2", display_width = 10L), 
    ln_hba1c_v00 = structure(c(1.76985463384001, NA, 1.7509374747078, 
    2.07693841146172, 1.8531680973567, 1.9473377010465, NA, 1.87333945622048, 
    1.78339121955754, 2.19833507162025), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    ln_hba1c_v66 = structure(c(1.73518911773966, NA, 1.71199450075919, 
    2.08567209143047, 1.70474809223843, 1.78339121955754, NA, 
    1.8196988379173, 1.74919985480926, 2.19499988231411), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    ln_hba1c_v01 = structure(c(1.72455071953461, 1.86872051036418, 
    1.71018781553424, 2.11142458753289, 1.76814960358892, 1.85629799036563, 
    1.66013102674962, 1.87180217690159, 1.7404661748405, 2.26384426467762
    ), label = "Analítica: Hemoglobina glicosilada (HbA1c %)", format.spss = "F5.2"), 
    ln_geaf_tot_v00 = structure(c(6.95554494281799, 8.34183930393788, 
    8.4371494837422, 8.11869575262367, 6.73240139150378, 7.38298895764493, 
    7.99228216582975, 9.23110234287667, 8.37462912676087, 5.0584090689011
    ), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2"), 
    ln_geaf_tot_v66 = structure(c(8.4371494837422, -Inf, 8.14946885573527, 
    8.11869575262367, 7.64869212337793, 8.5434202359197, 8.54069094410428, 
    9.00108616558123, 8.82055617324912, 5.34610703038388), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2"), 
    ln_geaf_tot_v01 = structure(c(7.42554857206372, -Inf, 8.45814145696367, 
    6.32693628339561, 7.92632458219889, 8.07468643426916, 8.83246367957041, 
    8.23647878828005, 8.67274118026667, 6.95554494281799), label = "AF: Gasto energético en actividad física total (MET•min/sem)", format.spss = "F8.2")), row.names = c(NA, 
-10L), class = c("tbl_df", "tbl", "data.frame"))

The manual method is to create terciles of one variable d_cintura1_v66 = difference of waist circumference between baseline and 6 months. Afterwards, we need to filter and do the analysis (imagine doing it for all variables all the combinations: 1st tercile vs 2nd, 1st vs 3rd, 2nd vs 3rd)

sample %>% 
dplyr::select(paciente, matches(c("_v00", "_v01", "_v66")))  %>%
  mutate(terciles_d_cintura1_v66 = ntile(.$d_cintura1_v66,3)) %>% 
    dplyr::select(paciente, matches("_v66")) %>% 
  pivot_longer(!c(paciente, terciles_d_cintura1_v66)) %>% 
  dplyr::filter(terciles_d_cintura1_v66 == "1" | terciles_d_cintura1_v66 == "2") %>% 
  dplyr::filter(str_detect(name, "d_")) %>% 
   mutate(terciles_d_cintura1_v66 = factor(terciles_d_cintura1_v66)) %>%
  pivot_wider(names_from= terciles_d_cintura1_v66, values_from=value) %>% 
  select(-paciente) %>% 
  group_by(name) %>% 
  do(tidy(t.test(.$"1", .$"2", na.rm = TRUE))) %>% 
  dplyr::ungroup()

My loop has the appearance below, I have noticed that changing the order of colvars and rowvars has signficance value (I wast told that this did not affect because iterate searching string)

Clarification: I ommited all variables with _v01 in the loop which are related the same way but 12 months intervention but I need to do it as well
Clarification 2: the temp object I don't really understand its function. It is something thought up by my colleague

colvars <- c("d_peso1_v66","d_homa_v66", "d_hughrelin_v66", "d_huinsulin_v66", "d_hucpeptide_v66", "d_huglucagon_v66",  "d_huleptin_v66", "d_huresistin_v66", "d_huvisfatin_v66", "d_hupai1_v66","d_huglp1_v66", "d_hugip_v66",  "d_glucosa_v66",  "d_coltot_v66","d_hdl_v66", "d_ldl_calc_v66","d_col_rema_v66","d_trigli_v66","d_hba1c_v66","d_tasis2_e_v66","d_tadias2_e_v66","d_cintura1_v66","d_geaf_tot_v66","d_p17_total_v66")

# Variables which should be added to the table. 
# The "for" loop works with the visit determined by colvars[i].
rowvars <- c("peso1", "homa", "i_hughrelin", "i_huinsulin", "i_hucpeptide", "i_huglucagon", "i_huleptin", "i_huresistin", "i_huvisfatin", "i_hupai1", "i_huglp1", "i_hugip", "glucosa", "coltot", "hdl", "ldl_calc", "col_rema", "trigli", "hba1c", "tasis2_e", "tadias2_e", "cintura1", "geaf_tot", "p17_total")

#loop
for(i in 1:length(colvars)){
  
  # Settle variables to work with.
  var <- colvars[i]
  visita <- substr(var, nchar(var)-3, nchar(var))
  rvars <- paste0("d_", rowvars, visita)
  rvars <- sub("d_i_", "d_", rvars)


P1 <- sample%>% 
    dplyr::select(paciente, all_of(var), rvars)  %>%
    mutate(temp = ntile(eval(parse(text=var)) ,3)) %>% 
    pivot_longer(!c(paciente, temp)) %>% 
    dplyr::filter(temp == "1" | temp == "2") %>% 
    mutate(temp = factor(temp)) %>%
    pivot_wider(names_from= temp, values_from=value) %>% 
    select(-paciente) %>% 
    group_by(name) %>% 
    do(tidy(t.test(.$"1", .$"2", na.rm = TRUE))) %>% 
    dplyr::ungroup() %>%
    pull(p.value)


  P1 <- ifelse(P1 < 2.2e-16, "< 2.2e-16",
        ifelse(P1 < 0.001, "< 0.001",
        sprintf("%.3f", round(P1, digits = 3))))
  

  # Groups to be compared 1st tercile vs 3rd tercile
  P2 <- sample%>% 
    dplyr::select(paciente, all_of(var), rvars)  %>%
    mutate(temp = ntile(eval(parse(text=var)) ,3)) %>% 
    pivot_longer(!c(paciente, temp)) %>% 
    dplyr::filter(temp == "1" | temp == "3") %>% 
    mutate(temp = factor(temp)) %>%
    pivot_wider(names_from= temp, values_from=value) %>% 
    select(-paciente) %>% 
    group_by(name) %>% 
    do(tidy(t.test(.$"1", .$"3", na.rm = TRUE))) %>% 
    dplyr::ungroup() %>%
    pull(p.value)
   

  
  P2 <- ifelse(P2 < 2.2e-16, "< 2.2e-16",
        ifelse(P2 < 0.001, "< 0.001",
        sprintf("%.3f", round(P2, digits = 3))))
  
  # Groups to be compared 2nd tercile vs 3rd tercile
  
  P3 <- sample%>% 
    dplyr::select(paciente, all_of(var), rvars)  %>%
    mutate(temp = ntile(eval(parse(text=var)) ,3)) %>% 
    pivot_longer(!c(paciente, temp)) %>% 
    dplyr::filter(temp == "2" | temp == "3") %>% 
    mutate(temp = factor(temp)) %>%
    pivot_wider(names_from= temp, values_from=value) %>% 
    select(-paciente) %>% 
    group_by(name) %>% 
    do(tidy(t.test(.$"2", .$"3", na.rm = TRUE))) %>% 
    dplyr::ungroup() %>%
    pull(p.value)
  
  P3 <- ifelse(P3 < 2.2e-16, "< 2.2e-16",
        ifelse(P3 < 0.001, "< 0.001",
        sprintf("%.3f", round(P3, digits = 3))))
  
  P<- cbind(P1, P2, P3)
  
  
  # Makes the table
  tab <- cbind(P)
  rownames(tab) <- rvars
  colnames(tab) <- c("1st vs 2nd","1st vs 3rd","2nd vs 3rd")
  
  
  # Exports the table to excel
  write.csv2(tab, paste0("terciles of sample ",var,".csv"))
}

So, this is it. I have made comparisons and don't match. Don't guess where the loop falls apart
Thanks in advance

nirgrahamuk · November 16, 2021, 3:18pm

pretty sure you arent accounting for the order of things so when you do

 rownames(tab) <- rvars

its somewhat arbitrary what you would end up with. i.e. your calculations are correct but your loop mislabels

consider adding both the following modifications:

sorting the colvars and rowvars that drive your loop

colvars <-sort( c("d_peso1_v66","d_homa_v66", "d_hughrelin_v66", "d_huinsulin_v66", "d_hucpeptide_v66", "d_huglucagon_v66",  "d_huleptin_v66", "d_huresistin_v66", "d_huvisfatin_v66", "d_hupai1_v66","d_huglp1_v66", "d_hugip_v66",  "d_glucosa_v66",  "d_coltot_v66","d_hdl_v66", "d_ldl_calc_v66","d_col_rema_v66","d_trigli_v66","d_hba1c_v66","d_tasis2_e_v66","d_tadias2_e_v66","d_cintura1_v66","d_geaf_tot_v66","d_p17_total_v66"))
rowvars <- sort(c("peso1", "homa", "i_hughrelin", "i_huinsulin", "i_hucpeptide", "i_huglucagon", "i_huleptin", "i_huresistin", "i_huvisfatin", "i_hupai1", "i_huglp1", "i_hugip", "glucosa", "coltot", "hdl", "ldl_calc", "col_rema", "trigli", "hba1c", "tasis2_e", "tadias2_e", "cintura1", "geaf_tot", "p17_total"))

sort your t.test results

 P1 <- sample%>% 
   dplyr::select(paciente, all_of(var), rvars)  %>%
   mutate(temp = ntile(eval(parse(text=var)) ,3)) %>% 
   pivot_longer(!c(paciente, temp)) %>% 
   dplyr::filter(temp == "1" | temp == "2") %>% 
   mutate(temp = factor(temp)) %>%
   pivot_wider(names_from= temp, values_from=value) %>% 
   select(-paciente) %>% 
   group_by(name) %>% 
   arrange(name) %>%   # <------------*** sorting
   do(tidy(t.test(.$"1", .$"2", na.rm = TRUE))) %>% 
   dplyr::ungroup() %>%
   pull(p.value)

I think those would be enough to correct the mislabeling.
Other advice might be to make a function to do what you need to do to make P1,P2,P3 because they seem similar but for the 1/2,2/3,1/3 combo's which could be passed as params.

Javier9 · November 16, 2021, 8:02pm

Totally right again, thank you.
The select function sort variables in alphabetical order. I was assuming that the order remained based on the definition of the object "rowvars"!

I am not an expert in functions, in fact I am kind of scared of syntax involving functions. Do you have a little suggestion?

system · November 23, 2021, 8:02pm

This topic was automatically closed 7 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.