Alan Mishler
2013-Aug-07 16:01 UTC
[R] lmer( ) parameter estimates changing depending on dummy-coded reference level
Hi everyone, I'm writing with a question about lmer( ) in the lme4 package. I've searched around for answers and done quite a bit of experimentation with toy data sets to figure out my issue, and I haven't been able to resolve it. I'm running linear mixed effects models on a large, sparse dataset in which I'm regressing reaction time (a continuous variable) on several categorical factors: Block (Block1/Block2/Block3), Group (monolingual/bilingual), and Type (target/nontarget). As a way of examining simple effects, I am dummy-coding specific factors, setting each level of a given factor as the reference level in turn. For example, I generate three models with each of the three levels of Block coded as the reference level, without changing the codings of the other factors: ## Model with Block 1 as reference level contrasts(nback.low$Group) <- c(1, 0) # monoling ref contrasts(nback.low$Type) <- c(1, -1) contrasts(nback.low$Block) <- matrix(c(0, 1, 0, 0, 0, 1), ncol=2) #B1 ref, B2: 1, B3: 2 glmerNL.RS.SI_RTgxb1 <- glmer(WinRTs~(Group*Block*Type) + (1+Block+Type|Subject) + (1+Group+Block|Item),data=nback.low) ## Model with Block 2 as reference level contrasts(nback.low$Group) <- c(1, 0) # monoling ref contrasts(nback.low$Type) <- c(1, -1) contrasts(nback.low$Block) <- matrix(c(1, 0, 0, 0, 0, 1), ncol=2) #B2 ref, B1: 1, B3: 2 glmerNL.RS.SI_RTgxb2 <- glmer(WinRTs~(Group*Block*Type) + (1+Block+Type|Subject) + (1+Group+Block|Item),data=nback.low) ## Model with Block 3 as reference level contrasts(nback.low$Group) <- c(1, -1) # monoling 'ref' contrasts(nback.low$Type) <- c(1, 0) # target ref contrasts(nback.low$Block) <- matrix(c(1, 0, 0, 0, 1, 0), ncol=2) #B3 ref, B1: 1, B2: 2 glmerNL.RS.SI_RTbxt3 <- glmer(WinRTs~(Group*Block*Type) + (1+Block+Type|Subject) + (1+Group+Block|Item),data=nback.low) summary(glmerNL.RS.SI_RTbxt3) The issue I'm having is that contrasts that I believe should be identical are not. Below are summaries of the three models. You can see that the estimate of the fixed effect of Block1 (the contrast between Block 1 and Block 2) is -117.98 in the first model and 118.478 in the second model. To my understanding, they should be identical except for the sign. Similar discrepancies can be seen in the other Block contrasts. There are two subjects who have no data at Block 1, so I removed them and re-ran the models, but the same issue occurred. Separately, I removed the random effects for Item, without removing those two subjects, and when I did that the discrepancies disappeared. I have a feeling this means that my models are too complex for my data, but I'm not sure what I should look at to (dis)confirm this hunch or how exactly to proceed if that is the case. (As an example of the sparseness of the data, items are repeated across subjects, but each subject has only one data point per item, or zero data points per item for trials where they didn't respond correctly. However, I didn't get any warnings about model convergence, or any warnings at all.) Any clues as to why I'm getting this results would be very much appreciated. Thanks in advance, Alan Mishler Research Assistant University of Maryland -- ## Model 1 output: Block 1 as reference level ##> glmerNL.RS.SI_RTgxb1Linear mixed model fit by REML Formula: WinRTs ~ (Group * Block * Type) + (1 + Block + Type | Subject) + (1 + Group + Block | Item) Data: nback.low AIC BIC logLik deviance REMLdev 181965 182210 -90949 181990 181899 Random effects: Groups Name Variance Std.Dev. Corr Item (Intercept) 6134.64 78.324 Group1 2336.31 48.335 -1.000 Block1 1759.56 41.947 1.000 -1.000 Block2 846.06 29.087 0.771 -0.771 0.771 Subject (Intercept) 132462.22 363.954 Block1 6011.83 77.536 -0.034 Block2 10883.42 104.324 -0.127 0.915 Type1 13653.97 116.850 -0.194 0.260 0.109 Residual 97048.54 311.526 Number of obs: 12640, groups: Item, 288; Subject, 52 Fixed effects: Estimate Std. Error t value (Intercept) 1016.73 73.39 13.855 Group1 12.39 101.67 0.122 Block1 -117.98 18.97 -6.220 Block2 -175.96 23.60 -7.455 Type1 -136.36 24.96 -5.463 Group1:Block1 46.43 26.05 1.782 Group1:Block2 76.70 32.53 2.358 Group1:Type1 20.55 34.15 0.602 Block1:Type1 62.16 10.47 5.934 Block2:Type1 96.66 10.46 9.243 Group1:Block1:Type1 -18.39 14.05 -1.309 Group1:Block2:Type1 -39.94 14.13 -2.826 Correlation of Fixed Effects: (Intr) Group1 Block1 Block2 Type1 Gr1:B1 Gr1:B2 Gr1:T1 Bl1:T1 Bl2:T1 G1:B1: Group1 -0.721 Block1 -0.065 0.049 Block2 -0.146 0.106 0.815 Type1 -0.186 0.134 0.219 0.106 Grop1:Blck1 0.053 -0.074 -0.715 -0.587 -0.160 Grop1:Blck2 0.108 -0.150 -0.586 -0.721 -0.077 0.818 Group1:Typ1 0.136 -0.189 -0.160 -0.078 -0.721 0.226 0.110 Block1:Typ1 0.012 -0.009 -0.086 -0.038 -0.163 0.063 0.028 0.131 Block2:Typ1 0.013 -0.009 -0.051 -0.070 -0.189 0.037 0.051 0.144 0.533 Grp1:Bl1:T1 -0.009 0.013 0.064 0.028 0.153 -0.090 -0.043 -0.216 -0.702 -0.371 Grp1:Bl2:T1 -0.010 0.014 0.038 0.052 0.156 -0.055 -0.072 -0.218 -0.371 -0.717 0.527 ## Model 2 output: Block 2 as reference level ##> glmerNL.RS.SI_RTgxb2 [Block 2 as reference level]Linear mixed model fit by REML Formula: WinRTs ~ (Group * Block * Type) + (1 + Block + Type | Subject) + (1 + Group + Block | Item) Data: nback.low AIC BIC logLik deviance REMLdev 181931 182177 -90933 181957 181865 Random effects: Groups Name Variance Std.Dev. Corr Item (Intercept) 14663.50 121.093 Group1 2400.09 48.991 -1.000 Block1 6561.73 81.004 -0.582 0.582 Block2 590.58 24.302 -0.677 0.677 -0.205 Subject (Intercept) 136638.85 369.647 Block1 5868.05 76.603 -0.172 Block2 2128.70 46.138 -0.146 -0.388 Type1 13788.61 117.425 -0.140 -0.259 -0.188 Residual 95743.04 309.424 Number of obs: 12640, groups: Item, 288; Subject, 52 Fixed effects: Estimate Std. Error t value (Intercept) 898.925 74.610 12.048 Group1 58.779 103.098 0.570 Block1 118.478 19.267 6.149 Block2 -58.059 13.680 -4.244 Type1 -74.323 25.553 -2.909 Group1:Block1 -46.515 25.799 -1.803 Group1:Block2 30.232 18.734 1.614 Group1:Type1 2.166 34.141 0.063 Block1:Type1 -64.291 11.246 -5.717 Block2:Type1 33.892 10.045 3.374 Group1:Block1:Type1 19.392 13.995 1.386 Group1:Block2:Type1 -21.160 13.634 -1.552 Correlation of Fixed Effects: (Intr) Group1 Block1 Block2 Type1 Gr1:B1 Gr1:B2 Gr1:T1 Bl1:T1 Bl2:T1 G1:B1: Group1 -0.720 Block1 -0.183 0.126 Block2 -0.152 0.108 -0.033 Type1 -0.132 0.096 -0.174 -0.095 Grop1:Blck1 0.126 -0.176 -0.699 0.019 0.130 Grop1:Blck2 0.106 -0.147 0.021 -0.722 0.070 -0.030 Group1:Typ1 0.099 -0.137 0.130 0.071 -0.714 -0.188 -0.102 Block1:Typ1 0.009 -0.007 -0.079 -0.049 -0.239 0.059 0.036 0.148 Block2:Typ1 0.010 -0.008 -0.036 -0.107 -0.219 0.027 0.079 0.152 0.416 Grp1:Bl1:T1 -0.008 0.010 0.063 0.040 0.137 -0.092 -0.050 -0.194 -0.655 -0.345 Grp1:Bl2:T1 -0.008 0.010 0.026 0.079 0.141 -0.036 -0.101 -0.200 -0.315 -0.722 0.482 ## Model 3 output: Block 3 as reference level ##> glmerNL.RS.SI_RTgxb3 [Block 3 as reference level]Linear mixed model fit by REML Formula: WinRTs ~ (Group * Block * Type) + (1 + Block + Type | Subject) + (1 + Group + Block | Item) Data: nback.low AIC BIC logLik deviance REMLdev 181932 182177 -90933 181958 181866 Random effects: Groups Name Variance Std.Dev. Corr Item (Intercept) 11296.69 106.286 Group1 2419.00 49.183 -1.000 Block1 7418.92 86.133 -0.522 0.522 Block2 494.56 22.239 0.596 -0.596 0.373 Subject (Intercept) 133759.20 365.731 Block1 10713.51 103.506 -0.153 Block2 2123.58 46.082 0.021 0.733 Type1 13764.21 117.321 -0.165 -0.107 0.189 Residual 95762.39 309.455 Number of obs: 12640, groups: Item, 288; Subject, 52 Fixed effects: Estimate Std. Error t value (Intercept) 840.86 73.76 11.399 Group1 89.07 102.02 0.873 Block1 176.58 23.92 7.381 Block2 58.07 13.66 4.251 Type1 -40.40 25.30 -1.597 Group1:Block1 -76.71 32.29 -2.376 Group1:Block2 -30.30 18.72 -1.618 Group1:Type1 -18.88 34.08 -0.554 Block1:Type1 -98.00 11.46 -8.555 Block2:Type1 -33.93 10.03 -3.382 Group1:Block1:Type1 40.12 14.06 2.853 Group1:Block2:Type1 21.06 13.63 1.545 Correlation of Fixed Effects: (Intr) Group1 Block1 Block2 Type1 Gr1:B1 Gr1:B2 Gr1:T1 Bl1:T1 Bl2:T1 G1:B1: Group1 -0.720 Block1 -0.170 0.119 Block2 -0.031 0.024 0.595 Type1 -0.156 0.113 -0.073 0.138 Grop1:Blck1 0.119 -0.164 -0.706 -0.434 0.054 Grop1:Blck2 0.026 -0.034 -0.429 -0.723 -0.101 0.603 Group1:Typ1 0.116 -0.161 0.054 -0.102 -0.717 -0.079 0.142 Block1:Typ1 0.009 -0.006 -0.063 -0.045 -0.231 0.047 0.034 0.148 Block2:Typ1 0.009 -0.007 -0.032 -0.107 -0.177 0.024 0.079 0.138 0.461 Grp1:Bl1:T1 -0.007 0.009 0.052 0.037 0.141 -0.073 -0.048 -0.195 -0.651 -0.357 Grp1:Bl2:T1 -0.007 0.009 0.025 0.079 0.144 -0.030 -0.101 -0.200 -0.324 -0.723 0.490
Bert Gunter
2013-Aug-07 17:14 UTC
[R] lmer( ) parameter estimates changing depending on dummy-coded reference level
I think it's fair to say that this should really be posted on the mixed-model specific (especially using lme4) list r-sig-mixed-models and not here. Cheers, Bert On Wed, Aug 7, 2013 at 9:01 AM, Alan Mishler <amishler at umd.edu> wrote:> Hi everyone, > > I'm writing with a question about lmer( ) in the lme4 package. I've > searched around for > answers and done quite a bit of experimentation with toy data sets to > figure out my issue, and I haven't been able to resolve it. > > I'm running linear mixed effects models on a large, sparse dataset in > which I'm regressing reaction time (a continuous variable) on several > categorical factors: Block (Block1/Block2/Block3), Group > (monolingual/bilingual), and Type (target/nontarget). > > As a way of examining simple effects, I am dummy-coding specific > factors, setting each level of a given factor as the reference level > in turn. For example, I generate three models with each of the three > levels of Block coded as the reference level, without changing the > codings of the other factors: > > ## Model with Block 1 as reference level > contrasts(nback.low$Group) <- c(1, 0) # monoling ref > contrasts(nback.low$Type) <- c(1, -1) > contrasts(nback.low$Block) <- matrix(c(0, 1, 0, 0, 0, 1), ncol=2) #B1 > ref, B2: 1, B3: 2 > glmerNL.RS.SI_RTgxb1 <- glmer(WinRTs~(Group*Block*Type) + > (1+Block+Type|Subject) + (1+Group+Block|Item),data=nback.low) > > ## Model with Block 2 as reference level > contrasts(nback.low$Group) <- c(1, 0) # monoling ref > contrasts(nback.low$Type) <- c(1, -1) > contrasts(nback.low$Block) <- matrix(c(1, 0, 0, 0, 0, 1), ncol=2) #B2 > ref, B1: 1, B3: 2 > glmerNL.RS.SI_RTgxb2 <- glmer(WinRTs~(Group*Block*Type) + > (1+Block+Type|Subject) + (1+Group+Block|Item),data=nback.low) > > ## Model with Block 3 as reference level > contrasts(nback.low$Group) <- c(1, -1) # monoling 'ref' > contrasts(nback.low$Type) <- c(1, 0) # target ref > contrasts(nback.low$Block) <- matrix(c(1, 0, 0, 0, 1, 0), ncol=2) #B3 > ref, B1: 1, B2: 2 > glmerNL.RS.SI_RTbxt3 <- glmer(WinRTs~(Group*Block*Type) + > (1+Block+Type|Subject) + (1+Group+Block|Item),data=nback.low) > summary(glmerNL.RS.SI_RTbxt3) > > The issue I'm having is that contrasts that I believe should be > identical are not. Below are summaries of the three > models. You can see that the estimate of the fixed effect of Block1 > (the contrast between Block 1 and Block 2) is -117.98 in the first > model and 118.478 in the second model. To my understanding, they > should be identical except for the sign. Similar discrepancies can be > seen in the other Block contrasts. > > There are two subjects who have no data at Block 1, so I removed them > and re-ran the models, but the same issue occurred. Separately, I > removed the random effects for Item, without removing those two > subjects, and when I did that the discrepancies disappeared. I have a > feeling this means that my models are too complex for my data, but I'm > not sure what I should look at to (dis)confirm this hunch or how > exactly to proceed if that is the case. (As an example of the > sparseness of the data, items are repeated across subjects, but each > subject has only one data point per item, or zero data points per item > for trials where they didn't respond correctly. However, I didn't get > any warnings about model convergence, or any warnings at > all.) > > Any clues as to why I'm getting this results would be > very much appreciated. > > Thanks in advance, > > Alan Mishler > Research Assistant > University of Maryland > -- > ## Model 1 output: Block 1 as reference level ## >> glmerNL.RS.SI_RTgxb1 > Linear mixed model fit by REML > Formula: WinRTs ~ (Group * Block * Type) + (1 + Block + Type | > Subject) + (1 + Group + Block | Item) > Data: nback.low > AIC BIC logLik deviance REMLdev > 181965 182210 -90949 181990 181899 > > Random effects: > Groups Name Variance Std.Dev. Corr > Item (Intercept) 6134.64 78.324 > Group1 2336.31 48.335 -1.000 > Block1 1759.56 41.947 1.000 -1.000 > Block2 846.06 29.087 0.771 -0.771 0.771 > Subject (Intercept) 132462.22 363.954 > Block1 6011.83 77.536 -0.034 > Block2 10883.42 104.324 -0.127 0.915 > Type1 13653.97 116.850 -0.194 0.260 0.109 > Residual 97048.54 311.526 > Number of obs: 12640, groups: Item, 288; Subject, 52 > > Fixed effects: > Estimate Std. Error t value > (Intercept) 1016.73 73.39 13.855 > Group1 12.39 101.67 0.122 > Block1 -117.98 18.97 -6.220 > Block2 -175.96 23.60 -7.455 > Type1 -136.36 24.96 -5.463 > Group1:Block1 46.43 26.05 1.782 > Group1:Block2 76.70 32.53 2.358 > Group1:Type1 20.55 34.15 0.602 > Block1:Type1 62.16 10.47 5.934 > Block2:Type1 96.66 10.46 9.243 > Group1:Block1:Type1 -18.39 14.05 -1.309 > Group1:Block2:Type1 -39.94 14.13 -2.826 > > Correlation of Fixed Effects: > (Intr) Group1 Block1 Block2 Type1 Gr1:B1 Gr1:B2 Gr1:T1 > Bl1:T1 Bl2:T1 G1:B1: > Group1 -0.721 > Block1 -0.065 0.049 > Block2 -0.146 0.106 0.815 > Type1 -0.186 0.134 0.219 0.106 > Grop1:Blck1 0.053 -0.074 -0.715 -0.587 -0.160 > Grop1:Blck2 0.108 -0.150 -0.586 -0.721 -0.077 0.818 > Group1:Typ1 0.136 -0.189 -0.160 -0.078 -0.721 0.226 0.110 > Block1:Typ1 0.012 -0.009 -0.086 -0.038 -0.163 0.063 0.028 0.131 > Block2:Typ1 0.013 -0.009 -0.051 -0.070 -0.189 0.037 0.051 0.144 > 0.533 > Grp1:Bl1:T1 -0.009 0.013 0.064 0.028 0.153 -0.090 -0.043 -0.216 > -0.702 -0.371 > Grp1:Bl2:T1 -0.010 0.014 0.038 0.052 0.156 -0.055 -0.072 -0.218 > -0.371 -0.717 0.527 > > ## Model 2 output: Block 2 as reference level ## >> glmerNL.RS.SI_RTgxb2 [Block 2 as reference level] > Linear mixed model fit by REML > Formula: WinRTs ~ (Group * Block * Type) + (1 + Block + Type | > Subject) + (1 + Group + Block | Item) > Data: nback.low > AIC BIC logLik deviance REMLdev > 181931 182177 -90933 181957 181865 > > Random effects: > Groups Name Variance Std.Dev. Corr > Item (Intercept) 14663.50 121.093 > Group1 2400.09 48.991 -1.000 > Block1 6561.73 81.004 -0.582 0.582 > Block2 590.58 24.302 -0.677 0.677 -0.205 > Subject (Intercept) 136638.85 369.647 > Block1 5868.05 76.603 -0.172 > Block2 2128.70 46.138 -0.146 -0.388 > Type1 13788.61 117.425 -0.140 -0.259 -0.188 > Residual 95743.04 309.424 > Number of obs: 12640, groups: Item, 288; Subject, 52 > > Fixed effects: > Estimate Std. Error t value > (Intercept) 898.925 74.610 12.048 > Group1 58.779 103.098 0.570 > Block1 118.478 19.267 6.149 > Block2 -58.059 13.680 -4.244 > Type1 -74.323 25.553 -2.909 > Group1:Block1 -46.515 25.799 -1.803 > Group1:Block2 30.232 18.734 1.614 > Group1:Type1 2.166 34.141 0.063 > Block1:Type1 -64.291 11.246 -5.717 > Block2:Type1 33.892 10.045 3.374 > Group1:Block1:Type1 19.392 13.995 1.386 > Group1:Block2:Type1 -21.160 13.634 -1.552 > > Correlation of Fixed Effects: > (Intr) Group1 Block1 Block2 Type1 Gr1:B1 Gr1:B2 Gr1:T1 > Bl1:T1 Bl2:T1 G1:B1: > Group1 -0.720 > Block1 -0.183 0.126 > Block2 -0.152 0.108 -0.033 > Type1 -0.132 0.096 -0.174 -0.095 > Grop1:Blck1 0.126 -0.176 -0.699 0.019 0.130 > Grop1:Blck2 0.106 -0.147 0.021 -0.722 0.070 -0.030 > Group1:Typ1 0.099 -0.137 0.130 0.071 -0.714 -0.188 -0.102 > Block1:Typ1 0.009 -0.007 -0.079 -0.049 -0.239 0.059 0.036 0.148 > Block2:Typ1 0.010 -0.008 -0.036 -0.107 -0.219 0.027 0.079 0.152 > 0.416 > Grp1:Bl1:T1 -0.008 0.010 0.063 0.040 0.137 -0.092 -0.050 -0.194 > -0.655 -0.345 > Grp1:Bl2:T1 -0.008 0.010 0.026 0.079 0.141 -0.036 -0.101 -0.200 > -0.315 -0.722 0.482 > > ## Model 3 output: Block 3 as reference level ## >> glmerNL.RS.SI_RTgxb3 [Block 3 as reference level] > Linear mixed model fit by REML > Formula: WinRTs ~ (Group * Block * Type) + (1 + Block + Type | > Subject) + (1 + Group + Block | Item) > Data: nback.low > AIC BIC logLik deviance REMLdev > 181932 182177 -90933 181958 181866 > > Random effects: > Groups Name Variance Std.Dev. Corr > Item (Intercept) 11296.69 106.286 > Group1 2419.00 49.183 -1.000 > Block1 7418.92 86.133 -0.522 0.522 > Block2 494.56 22.239 0.596 -0.596 0.373 > Subject (Intercept) 133759.20 365.731 > Block1 10713.51 103.506 -0.153 > Block2 2123.58 46.082 0.021 0.733 > Type1 13764.21 117.321 -0.165 -0.107 0.189 > Residual 95762.39 309.455 > Number of obs: 12640, groups: Item, 288; Subject, 52 > > Fixed effects: > Estimate Std. Error t value > (Intercept) 840.86 73.76 11.399 > Group1 89.07 102.02 0.873 > Block1 176.58 23.92 7.381 > Block2 58.07 13.66 4.251 > Type1 -40.40 25.30 -1.597 > Group1:Block1 -76.71 32.29 -2.376 > Group1:Block2 -30.30 18.72 -1.618 > Group1:Type1 -18.88 34.08 -0.554 > Block1:Type1 -98.00 11.46 -8.555 > Block2:Type1 -33.93 10.03 -3.382 > Group1:Block1:Type1 40.12 14.06 2.853 > Group1:Block2:Type1 21.06 13.63 1.545 > > Correlation of Fixed Effects: > (Intr) Group1 Block1 Block2 Type1 Gr1:B1 Gr1:B2 Gr1:T1 > Bl1:T1 Bl2:T1 G1:B1: > Group1 -0.720 > Block1 -0.170 0.119 > Block2 -0.031 0.024 0.595 > Type1 -0.156 0.113 -0.073 0.138 > Grop1:Blck1 0.119 -0.164 -0.706 -0.434 0.054 > Grop1:Blck2 0.026 -0.034 -0.429 -0.723 -0.101 0.603 > Group1:Typ1 0.116 -0.161 0.054 -0.102 -0.717 -0.079 0.142 > Block1:Typ1 0.009 -0.006 -0.063 -0.045 -0.231 0.047 0.034 0.148 > Block2:Typ1 0.009 -0.007 -0.032 -0.107 -0.177 0.024 0.079 0.138 > 0.461 > Grp1:Bl1:T1 -0.007 0.009 0.052 0.037 0.141 -0.073 -0.048 -0.195 > -0.651 -0.357 > Grp1:Bl2:T1 -0.007 0.009 0.025 0.079 0.144 -0.030 -0.101 -0.200 > -0.324 -0.723 0.490 > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.-- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm