5 votos

Examen de las tendencias con interacciones y con estratificación: obtención de resultados discordantes

Estoy examinando el efecto de los ingresos (clasificados en quintiles) sobre una variable de respuesta durante diferentes años (de 2003 a 2014). Ajusto por algunas otras covariables y tengo mediciones repetidas sobre el mismo individuo.

El problema: Para obtener el efecto de la renta sobre la variable de respuesta de cada año, incluyo un término de interacción entre el año y el quintil de renta. Esto me da resultados inverosímiles, ya que los que son ricos tienden a tener valores ajustados más altos, lo que es poco probable. Así que rehice el análisis, ajustando un modelo de regresión distinto para cada año, y eso me dio resultados plausibles; es decir, los ricos tenían valores ajustados más bajos.

¿No deberían los dos métodos dar resultados algo similares?

Mis cálculos:

    > # START
    > 
    > # Converting the year variable into a factor
    > data$year <- as.factor(data$year)
    > 
    > # Fitting a model with random effects for person
> require(lme4)
> fit  <- lmer(response ~ income*year + sex + age + education + biomarker + (1 | id), data=data)
> # Using lsmeans to predict means (95% confidence intervals)
> require(lsmeans)
> results1 <- lsmeans(fit, ~ income*year)
> summary(results1)
 income     year   lsmean        SE       df lower.CL upper.CL
 Quintile 1 2003 64.95472 0.2826723 190216.2 64.40069 65.50875
 Quintile 2 2003 65.49504 0.2716893 189962.8 64.96254 66.02755
 Quintile 3 2003 65.39961 0.2713204 189648.3 64.86782 65.93139
 Quintile 4 2003 65.51872 0.2715941 189734.3 64.98640 66.05103
 Quintile 5 2003 65.47502 0.2744592 190425.5 64.93709 66.01296
 Quintile 1 2004 64.03440 0.2541402 191982.7 63.53630 64.53251
 Quintile 2 2004 65.04364 0.2472738 191720.0 64.55899 65.52829
 Quintile 3 2004 64.98069 0.2425866 191644.4 64.50523 65.45616
 Quintile 4 2004 65.14219 0.2459067 191477.1 64.66022 65.62416
 Quintile 5 2004 65.25515 0.2496510 191947.6 64.76584 65.74446
 Quintile 1 2005 63.62453 0.2345767 192988.5 63.16476 64.08429
 Quintile 2 2005 64.18179 0.2294655 192853.2 63.73205 64.63154
 Quintile 3 2005 64.00250 0.2308264 192458.4 63.55008 64.45491
 Quintile 4 2005 63.79620 0.2276547 192775.7 63.35000 64.24240
 Quintile 5 2005 64.70116 0.2344171 193212.2 64.24171 65.16062
 Quintile 1 2006 64.17463 0.2228096 193482.8 63.73793 64.61134
 Quintile 2 2006 64.37025 0.2158588 193425.3 63.94717 64.79333
 Quintile 3 2006 64.30762 0.2141169 193366.8 63.88795 64.72728
 Quintile 4 2006 64.38315 0.2168376 193329.6 63.95815 64.80814
 Quintile 5 2006 64.22019 0.2198772 193526.5 63.78923 64.65114
 Quintile 1 2007 63.84188 0.2185725 193602.5 63.41349 64.27028
 Quintile 2 2007 63.92346 0.2104634 193527.6 63.51096 64.33597
 Quintile 3 2007 63.67954 0.2094212 193511.7 63.26908 64.09000
 Quintile 4 2007 64.08718 0.2103912 193478.5 63.67482 64.49955
 Quintile 5 2007 64.03627 0.2113704 193649.9 63.62199 64.45055
 Quintile 1 2008 64.65587 0.2043770 193372.1 64.25529 65.05644
 Quintile 2 2008 64.30141 0.1957655 193499.2 63.91772 64.68511
 Quintile 3 2008 65.04454 0.1955247 193589.4 64.66131 65.42776
 Quintile 4 2008 64.94073 0.1947715 193488.8 64.55898 65.32247
 Quintile 5 2008 65.00096 0.1979129 193119.2 64.61305 65.38886
 Quintile 1 2009 64.74611 0.1941592 191979.5 64.36556 65.12666
 Quintile 2 2009 64.68663 0.1872505 192801.9 64.31962 65.05363
 Quintile 3 2009 64.89048 0.1858611 192919.0 64.52620 65.25476
 Quintile 4 2009 65.19469 0.1848586 192734.9 64.83238 65.55701
 Quintile 5 2009 65.18344 0.1896058 192025.8 64.81182 65.55506
 Quintile 1 2010 64.99407 0.1863001 188874.1 64.62893 65.35922
 Quintile 2 2010 65.14159 0.1758624 190272.3 64.79691 65.48628
 Quintile 3 2010 65.21003 0.1740553 190655.7 64.86889 65.55118
 Quintile 4 2010 65.65492 0.1731845 190157.8 65.31548 65.99435
 Quintile 5 2010 65.42802 0.1774762 189156.3 65.08017 65.77587
 Quintile 1 2011 65.77711 0.1807193 179035.8 65.42290 66.13131
 Quintile 2 2011 65.81373 0.1719076 186787.1 65.47679 66.15066
 Quintile 3 2011 66.25649 0.1692650 187174.9 65.92473 66.58824
 Quintile 4 2011 66.24046 0.1672252 185968.9 65.91271 66.56822
 Quintile 5 2011 66.31895 0.1717513 184631.5 65.98232 66.65558
 Quintile 1 2012 66.45412 0.1815487 178604.9 66.09829 66.80995
 Quintile 2 2012 66.39743 0.1691640 185466.1 66.06587 66.72899
 Quintile 3 2012 66.61760 0.1674829 185934.4 66.28934 66.94586
 Quintile 4 2012 66.69909 0.1672601 185774.9 66.37126 67.02692
 Quintile 5 2012 66.74211 0.1697808 183081.3 66.40934 67.07487
 Quintile 1 2013 66.20088 0.1804738 177511.5 65.84716 66.55461
 Quintile 2 2013 66.05285 0.1710873 185354.9 65.71752 66.38817
 Quintile 3 2013 65.57061 0.1667456 185103.9 65.24379 65.89742
 Quintile 4 2013 65.96563 0.1669031 184335.3 65.63851 66.29276
 Quintile 5 2013 66.19121 0.1723801 183746.1 65.85335 66.52907
 Quintile 1 2014 65.37137 0.3060358 191891.8 64.77155 65.97120
 Quintile 2 2014 66.37503 0.2882805 188873.8 65.81001 66.94006
 Quintile 3 2014 65.49851 0.2876551 188265.5 64.93471 66.06231
 Quintile 4 2014 66.08503 0.2867954 188729.4 65.52291 66.64714
 Quintile 5 2014 65.98435 0.2901786 189169.5 65.41561 66.55309

Results are averaged over the levels of: sex, education 
Confidence level used: 0.95 
> # As you can see above, the richest (quintile 5) has higher values of a harmful biomarker, this is not likelybiomarker, this is not likely
    > 
    > # Redo the analysis by stratification, i.e one regression each year (I'll give 4 examples, years 2005, 2008, 2010 and 2014)
    > fit2005  <- lmer(response ~ income + sex + age + education + biomarker + (1 | id), data=data[data$year=="2005",])
        > fit2008  <- lmer(response ~ income + sex + age + education + biomarker + (1 | id), data=data[data$year=="2008",])
    > fit2010  <- lmer(response ~ income + sex + age + education + biomarker + (1 | id), data=data[data$year=="2010",])
        > fit2014  <- lmer(response ~ income + sex + age + education + biomarker + (1 | id), data=data[data$year=="2014",])
    > 
    > lsmeans(fit2005, "income")
     income       lsmean        SE      df lower.CL upper.CL
     Quintile 1 64.57677 0.3644627 7104.87 63.86232 65.29123
     Quintile 2 63.65426 0.3715929 7114.04 62.92582 64.38269
     Quintile 3 63.97948 0.3773103 7119.28 63.23984 64.71912
     Quintile 4 63.46368 0.3727073 7117.62 62.73306 64.19429
     Quintile 5 63.57556 0.3853513 7117.67 62.82016 64.33096

    Results are averaged over the levels of: sex, education 
    Confidence level used: 0.95 
    > lsmeans(fit2008, "income")
     income       lsmean        SE      df lower.CL upper.CL
     Quintile 1 64.27986 0.3268165 9880.81 63.63924 64.92049
     Quintile 2 65.07827 0.3288624 9866.53 64.43363 65.72291
     Quintile 3 64.98781 0.3265577 9859.02 64.34769 65.62793
     Quintile 4 64.89630 0.3305190 9868.49 64.24842 65.54419
     Quintile 5 63.66509 0.3428045 9898.35 62.99313 64.33706

    Results are averaged over the levels of: sex, education 
    Confidence level used: 0.95 
    > lsmeans(fit2010, "income")
     income       lsmean        SE       df lower.CL upper.CL
     Quintile 1 65.39530 0.2961521 11897.22 64.81480 65.97581
     Quintile 2 65.61791 0.2911321 11892.26 65.04724 66.18858
     Quintile 3 65.48423 0.2947050 11892.60 64.90656 66.06190
     Quintile 4 65.14303 0.2914349 11925.62 64.57177 65.71429
     Quintile 5 63.89145 0.3030998 11935.15 63.29733 64.48558

    Results are averaged over the levels of: sex, education 
    Confidence level used: 0.95 
    > lsmeans(fit2014, "income")
     income       lsmean        SE      df lower.CL upper.CL
     Quintile 1 67.05015 0.4888236 4523.42 66.09182 68.00849
     Quintile 2 65.41100 0.4801968 4523.32 64.46958 66.35242
     Quintile 3 65.35658 0.4740396 4525.63 64.42723 66.28592
     Quintile 4 65.03556 0.4793119 4526.56 64.09587 65.97524
     Quintile 5 64.94690 0.5001898 4529.88 63.96628 65.92751

    Results are averaged over the levels of: sex, education 
    Confidence level used: 0.95 

Debo haber entendido mal algo (importante) aquí...

¿Algún consejo?

0 votos

Su MODELO no incluye la interacción entre los ingresos y el año - sólo la llamada lsmeans. Reajuste el modelo y haga la lsmeans de nuevo y creo que verá un patrón diferente.

0 votos

Gracias por tomarse el tiempo de comentar. Ahora veo que el modelo descrito anteriormente no incluía la interacción, pero te aseguro que la llamada original en mi ordenador sí la incluye, así que he actualizado el post anterior para incluir la interacción y los nuevos resultados - el problema persiste. He subido los datos de muestra ahora. Creo que el problema es mi llamada a lsmeans en el modelo que tiene la interacción, porque sólo uso la interacción en la llamada a lsmeans...

0 votos

Creo que esto es una especie de paradoja de Simpson. Cuando se ajustan modelos separados por año, se están permitiendo diferentes efectos de todos los predictores para cada año. En esencia, esto incluye interacciones entre el año y los otros predictores, en comparación con su modelo combinado que asume que los efectos de esos predictores son los mismos cada año. Si incluye las interacciones, obtendrá estimaciones comparables (pero no errores estándar comparables porque se utiliza un término de error combinado)

3voto

anand Puntos 199

No tenía la id variable en el momento de preparar esta respuesta por lo que no reproduzco sus resultados. Pero la explicación se puede hacer utilizando lm se ajusta a estos datos.

En primer lugar, su modelo sin el id agrupación. Mostraré el modelo, las medias de LS y las comparaciones con el quintil de ingresos más bajo sólo para 2005:

> data$year = factor(data$year)
> fit.lm  <- lm(response ~ income*year + sex + age + education + biomarker, 
    data=data)
> lsmeans(fit.lm, trt.vs.ctrl1 ~ income, at = list(year = "2005"))

$lsmeans
 income       lsmean        SE    df lower.CL upper.CL
 Quintile 1 64.20122 1.0008221 19262 62.23952 66.16292
 Quintile 2 64.45945 1.0060940 19262 62.48742 66.43148
 Quintile 3 64.37015 0.9685111 19262 62.47178 66.26851
 Quintile 4 63.88540 0.9971111 19262 61.93098 65.83983
 Quintile 5 63.36031 1.0319215 19262 61.33765 65.38297

$contrasts
 contrast                  estimate       SE    df t.ratio p.value
 Quintile 2 - Quintile 1  0.2582246 1.417633 19262   0.182  0.9937
 Quintile 3 - Quintile 1  0.1689240 1.390806 19262   0.121  0.9976
 Quintile 4 - Quintile 1 -0.3158223 1.411699 19262  -0.224  0.9898
 Quintile 5 - Quintile 1 -0.8409129 1.437786 19262  -0.585  0.9013

(Nota: he eliminado algunas anotaciones en la salida.) Ahora, aquí están los mismos resultados, utilizando un modelo ajustado sólo a los datos de 2005:

> fit.lm2005  <- lm(response ~ income + sex + age + education + biomarker, 
    data=data, subset = (year=="2005"))  
> lsmeans(fit.lm2005, trt.vs.ctrl1 ~ income)

$lsmeans
 income       lsmean        SE  df lower.CL upper.CL
 Quintile 1 64.37016 0.9884950 998 62.43039 66.30993
 Quintile 2 64.31077 1.0036776 998 62.34121 66.28033
 Quintile 3 64.35983 0.9747897 998 62.44696 66.27270
 Quintile 4 64.05203 1.0154760 998 62.05932 66.04474
 Quintile 5 63.43306 1.0753143 998 61.32292 65.54320

$contrasts
 contrast                   estimate       SE  df t.ratio p.value
 Quintile 2 - Quintile 1 -0.05939015 1.385639 998  -0.043  0.9998
 Quintile 3 - Quintile 1 -0.01033274 1.357894 998  -0.008  1.0000
 Quintile 4 - Quintile 1 -0.31813180 1.400778 998  -0.227  0.9894
 Quintile 5 - Quintile 1 -0.93710083 1.465919 998  -0.639  0.8790

Ahora, considere la posibilidad de cambiar el modelo general para que los ingresos interactúen con cada otro predictor:

> fit.lmint  <- lm(response ~ year*(income + sex + age + education + biomarker), 
    data=data)
> lsmeans(fit.lmint, trt.vs.ctrl1 ~ income, at = list(year = "2005"))

$lsmeans
 income       lsmean       SE    df lower.CL upper.CL
 Quintile 1 64.33859 1.027770 19172 62.32408 66.35311
 Quintile 2 64.27920 1.034344 19172 62.25180 66.30661
 Quintile 3 64.32826 1.007009 19172 62.35444 66.30209
 Quintile 4 64.02046 1.047664 19172 61.96695 66.07398
 Quintile 5 63.40149 1.106308 19172 61.23303 65.56996

$contrasts
 contrast                   estimate       SE    df t.ratio p.value
 Quintile 2 - Quintile 1 -0.05939015 1.432060 19172  -0.041  0.9998
 Quintile 3 - Quintile 1 -0.01033274 1.403386 19172  -0.007  1.0000
 Quintile 4 - Quintile 1 -0.31813180 1.447706 19172  -0.220  0.9902
 Quintile 5 - Quintile 1 -0.93710083 1.515029 19172  -0.619  0.8878

Las medias de LS son ligeramente diferentes de las de fit.lm2005 pero las comparaciones son idénticas. La única razón por la que las medias LS están desplazadas es porque se utilizan diferentes valores de referencia para age y biomarker . En el siguiente código, muestro que obtenemos los mismos resultados (excepto por los errores de redondeo) de fit.lmint si utilizamos los mismos niveles de referencia:

> ref.grid(fit.lm2005)
'ref.grid' object with variables:
    income = Quintile 1, Quintile 2, Quintile 3, Quintile 4, Quintile 5
    sex = Female, Male
    age = 41.511
    education = <10 years education, 10 to 12 years education, College/University
    biomarker = 26.133

> lsmeans(fit.lmint, ~ income, at = list(year = "2005", 
    age = 41.511, biomarker = 26.133))
NOTE: Results may be misleading due to involvement in interactions
 income       lsmean       SE    df lower.CL upper.CL
 Quintile 1 64.37009 1.021611 19172 62.36765 66.37254
 Quintile 2 64.31070 1.037303 19172 62.27750 66.34391
 Quintile 3 64.35976 1.007447 19172 62.38508 66.33445
 Quintile 4 64.05196 1.049493 19172 61.99486 66.10906
 Quintile 5 63.43299 1.111337 19172 61.25468 65.61131

Mensaje para llevar a casa: Cuando se ajustan modelos por separado por algún factor, los ajustes diferirán de los de un modelo combinado, a menos que el modelo combinado incluya todas las interacciones con el factor que se utilizó para los ajustes por separado.

Además, tenga en cuenta que lsmeans sólo resume los resultados del modelo. Si utiliza una especificación como ~A*B en el lsmeans llamada, que hace no añadir un A:B al modelo, si esa interacción no estaba presente.

Nota adicional añadida posteriormente:

Con un modelo mixto, si el id es la variable cruzado con year entonces es probable que no obtener resultados idénticos con sólo incluir las interacciones de efectos fijos. esto se debe a que se estarían estimando diferentes id efectos en cada año también. Habría que hacer la parte aleatoria (1|id/year) .

0 votos

Es una hermosa manera de explicar esto. Aprecio sus esfuerzos.

0 votos

Sólo una pequeña nota: lsmeans(fit.lmint, trt.vs.ctrl1 ~ income, at = list(year = "1996")) y lsmeans(fit.lmint, trt.vs.ctrl1 ~ income, at = list(year = "2005")) Produce las mismas estimaciones, al igual que cualquier otro año...

0 votos

He probado esto y sí no obtener resultados idénticos con estos comandos para cada año. De hecho, no necesité probarlo porque la cuestión es que se pueden reproducir los resultados de un año determinado ajustados por separado utilizando fit.lsmint . Debe haber cometido un error. Por cierto, sin relación con esto, añadí un comentario sobre los modelos mixtos al final de la respuesta.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X