[R-sig-ME] Unidentifiable model in lmer

Takahiro Fushimi taka.6765 at gmail.com
Fri Sep 18 17:08:28 CEST 2015


Sorry about it. Does this work?

dataset <- data.frame(y = c(60625.19, 51088.38, 67283.03, 56550.54, 
NA,   NA, 58720.44,63939.01, 66298.32, 62662.15, 53604.41, 88338.55, 
62818.78, 68696.71,57433.98, 59105.99, 44801.97, 68655.48, 62645.16, 
69409.15, 43568.70,55693.21, 57322.63, 51095.55, 57660.41, 64128.63, 
55216.42, 59945.17,67284.31, 55401.13, 60471.59, NA, 63633.07, 
65939.14,NA, 69488.38, 43950.50, 52782.99, 55674.94, 70130.25, 72560.25, 
69297.95, 53188.98, 75687.53, 68242.23, 60696.15, 65087.04, 59377.71, 
55055.45, 66673.40, 65933.37, 62636.56, 49606.80, 58668.13, 56694.13, 
NA, NA, 61928.95,63208.84, 83786.71, 72727.53, 78873.45, 47482.32, 
65967.84, 53228.25, 65699.59, 61792.81, 44816.72, 71048.06, 59671.84, 
89513.43, 56972.30, 45439.25, 51357.83, NA, 60146.86, 56845.88, 
63181.13, 63704.80, 49140.22, 54538.14, 49411.32, 66059.73, 73147.55, 
55665.96, 66821.27, 66935.12, 54350.63, 48116.25, 67664.66, 64278.54, 
64555.03, 62463.60, 55831.53, 57392.61, 75727.00, 64839.50, 51009.78, 
65274.59, 63339.91, 62276.49, 66159.42, 58333.76, 60275.09) , A = 
as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)) , B = as.factor(c(1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1 ,1 ,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0)), ID = as.factor(c(72,  73,  74,  75,  76, 78,  79,  86,  88,  
89,  90, 91,  92,  93,  94,  95,  96,  121, 122, 123, 124, 77,  80,  81, 
82,  83,  84,  98,  116, 117, 118, 119, 120, 125, 128, 100, 101, 102, 
103, 104, 105, 106, 107, 108, 109, 110, 112, 113, 114, 115, 126, 127, 
72, 73, 74, 75, 76, 78, 79, 86, 88, 89, 90, 91, 92, 93, 94, 95, 96, 121, 
122, 123, 124, 77, 80, 81, 82, 83, 84, 98, 116, 117, 118, 119, 120, 125, 
128, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 112, 113, 
114, 115, 126, 127)))

Best regards,
Takahiro

On 9/18/15 10:41 AM, Thierry Onkelinx wrote:
> Dear Takahiro,
>
> Please send the output of dput(dataset). That's a lot easier to
> copy-paste into an R session.
>
> Best regards,
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be no
> more than asking him to perform a post-mortem examination: he may be
> able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does
> not ensure that a reasonable answer can be extracted from a given body
> of data. ~ John Tukey
>
>
> 2015-09-18 16:40 GMT+02:00 Takahiro Fushimi <taka.6765 at gmail.com>:
>> The data set is as follows:
>>> dataset
>>             y A B  ID
>> 1   60625.19 0 1  72
>> 2   51088.38 0 1  73
>> 3   67283.03 0 1  74
>> 4   56550.54 0 1  75
>> 5         NA 0 1  76
>> 6         NA 0 1  78
>> 7   58720.44 0 1  79
>> 8   63939.01 0 1  86
>> 9   66298.32 0 1  88
>> 10  62662.15 0 1  89
>> 11  53604.41 0 1  90
>> 12  88338.55 0 1  91
>> 13  62818.78 0 1  92
>> 14  68696.71 0 1  93
>> 15  57433.98 0 1  94
>> 16  59105.99 0 1  95
>> 17  44801.97 0 1  96
>> 18  68655.48 0 1 121
>> 19  62645.16 0 1 122
>> 20  69409.15 0 1 123
>> 21  43568.70 0 1 124
>> 22  55693.21 1 1  77
>> 23  57322.63 1 1  80
>> 24  51095.55 1 1  81
>> 25  57660.41 1 1  82
>> 26  64128.63 1 1  83
>> 27  55216.42 1 1  84
>> 28  59945.17 1 1  98
>> 29  67284.31 1 1 116
>> 30  55401.13 1 1 117
>> 31  60471.59 1 1 118
>> 32        NA 1 1 119
>> 33  63633.07 1 1 120
>> 34  65939.14 1 1 125
>> 35        NA 1 1 128
>> 36  69488.38 2 1 100
>> 37  43950.50 2 1 101
>> 38  52782.99 2 1 102
>> 39  55674.94 2 1 103
>> 40  70130.25 2 1 104
>> 41  72560.25 2 1 105
>> 42  69297.95 2 1 106
>> 43  53188.98 2 1 107
>> 44  75687.53 2 1 108
>> 45  68242.23 2 1 109
>> 46  60696.15 2 1 110
>> 47  65087.04 2 1 112
>> 48  59377.71 2 1 113
>> 49  55055.45 2 1 114
>> 50  66673.40 2 1 115
>> 51  65933.37 2 1 126
>> 52  62636.56 2 1 127
>> 53  49606.80 0 0  72
>> 54  58668.13 0 0  73
>> 55  56694.13 0 0  74
>> 56        NA 0 0  75
>> 57        NA 0 0  76
>> 58  61928.95 0 0  78
>> 59  63208.84 0 0  79
>> 60  83786.71 0 0  86
>> 61  72727.53 0 0  88
>> 62  78873.45 0 0  89
>> 63  47482.32 0 0  90
>> 64  65967.84 0 0  91
>> 65  53228.25 0 0  92
>> 66  65699.59 0 0  93
>> 67  61792.81 0 0  94
>> 68  44816.72 0 0  95
>> 69  71048.06 0 0  96
>> 70  59671.84 0 0 121
>> 71  89513.43 0 0 122
>> 72  56972.30 0 0 123
>> 73  45439.25 0 0 124
>> 74  51357.83 1 0  77
>> 75        NA 1 0  80
>> 76  60146.86 1 0  81
>> 77  56845.88 1 0  82
>> 78  63181.13 1 0  83
>> 79  63704.80 1 0  84
>> 80  49140.22 1 0  98
>> 81  54538.14 1 0 116
>> 82  49411.32 1 0 117
>> 83  66059.73 1 0 118
>> 84  73147.55 1 0 119
>> 85  55665.96 1 0 120
>> 86  66821.27 1 0 125
>> 87  66935.12 1 0 128
>> 88  54350.63 2 0 100
>> 89  48116.25 2 0 101
>> 90  67664.66 2 0 102
>> 91  64278.54 2 0 103
>> 92  64555.03 2 0 104
>> 93  62463.60 2 0 105
>> 94  55831.53 2 0 106
>> 95  57392.61 2 0 107
>> 96  75727.00 2 0 108
>> 97  64839.50 2 0 109
>> 98  51009.78 2 0 110
>> 99  65274.59 2 0 112
>> 100 63339.91 2 0 113
>> 101 62276.49 2 0 114
>> 102 66159.42 2 0 115
>> 103 58333.76 2 0 126
>> 104 60275.09 2 0 127
>>
>> Best regards,
>> Takahiro
>>
>>
>> On 9/18/15 3:23 AM, Thierry Onkelinx wrote:
>>> The formula shouldn't be the problem. model.matrix() is clever enough
>>> to handle this. Have a look at the examples below.
>>>
>>> model.matrix(~ A + B + A * B, data= data.frame(A = 1, B = 1))
>>> model.matrix(~ A + B + A * B + B + A:B + B:A, data= data.frame(A = 1, B =
>>> 1))
>>>
>>> I suggest to give use a reproducible example of the problem so that we
>>> can invesigate what when wrong.
>>> ir. Thierry Onkelinx
>>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>>> and Forest
>>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>>> Kliniekstraat 25
>>> 1070 Anderlecht
>>> Belgium
>>>
>>> To call in the statistician after the experiment is done may be no
>>> more than asking him to perform a post-mortem examination: he may be
>>> able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
>>> The plural of anecdote is not data. ~ Roger Brinner
>>> The combination of some data and an aching desire for an answer does
>>> not ensure that a reasonable answer can be extracted from a given body
>>> of data. ~ John Tukey
>>>
>>>
>>> 2015-09-18 7:23 GMT+02:00 Ché Lucero <chelucero at uchicago.edu>:
>>>> I don't know if this is causing the problem you're seeing, but A*B
>>>> expands
>>>> to A + B + A:B, so your model right now is R ~ A + B + A + B + A:B +
>>>> (1|S).
>>>>
>>>> On Thu, Sep 17, 2015 at 5:53 PM Takahiro Fushimi <taka.6765 at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi everyone,
>>>>>
>>>>> I have been working on a linear mixed effect model by using lmer()
>>>>> function, but I got an error saying that the fitted model is not
>>>>> identifiable.
>>>>>
>>>>> The data set includes the following variables:
>>>>> y = a numeric variable
>>>>> factorA = a 3-level categorical variable
>>>>> factorB = a 2-level categorical variable
>>>>> subjectID = subject id number. 2 measurements of y for each subject
>>>>>
>>>>> R code and output are as follows:
>>>>>    > (result <- lmer(y ~ factorA + factorB + factorA*factorB +
>>>>> (1|subjectID)))
>>>>> Linear mixed model fit by REML ['merModLmerTest']
>>>>> Formula: y ~ factorA + factorB + factorA * factorB + (1 | subjectID)
>>>>> REML criterion at convergence: 1928.966
>>>>> Random effects:
>>>>>     Groups    Name        Std.Dev.
>>>>>     subjectID (Intercept) 4711
>>>>>     Residual              7688
>>>>> Number of obs: 97, groups:  subjectID, 51
>>>>> Fixed Effects:
>>>>>          (Intercept)           factorA1           factorA2 factorB1
>>>>> factorA1:factorB1  factorA2:factorB1
>>>>>                62411              -2700              -1124
>>>>> -1037               1279               2482
>>>>>    > summary(result)
>>>>> Model is not identifiable...
>>>>> summary from lme4 is returned
>>>>> some computational error has occurred in lmerTest
>>>>> Linear mixed model fit by REML ['lmerMod']
>>>>> Formula: y ~ factorA + factorB + factorA * factorB + (1 | subjectID)
>>>>>
>>>>> REML criterion at convergence: 1929
>>>>>
>>>>> Scaled residuals:
>>>>>         Min       1Q   Median       3Q      Max
>>>>> -1.93423 -0.62611  0.01837  0.48887  2.73380
>>>>>
>>>>> Random effects:
>>>>>     Groups    Name        Variance Std.Dev.
>>>>>     subjectID (Intercept) 22194074 4711
>>>>>     Residual              59108495 7688
>>>>> Number of obs: 97, groups:  subjectID, 51
>>>>>
>>>>> Fixed effects:
>>>>>                      Estimate Std. Error t value
>>>>> (Intercept)          62411       2065  30.227
>>>>> factorA1             -2700       3238  -0.834
>>>>> factorA2             -1124       3008  -0.374
>>>>> factorB1             -1037       2512  -0.413
>>>>> factorA1:factorB1     1279       4013   0.319
>>>>> factorA2:factorB1     2482       3642   0.681
>>>>>
>>>>> Correlation of Fixed Effects:
>>>>>                (Intr) fctrA1 fctrA2 fctrB1 fA1:B1
>>>>> factorA1    -0.638
>>>>> factorA2    -0.687  0.438
>>>>> factorB1    -0.608  0.388  0.418
>>>>> fctrA1:fcB1  0.381 -0.601 -0.262 -0.626
>>>>> fctrA2:fcB1  0.420 -0.268 -0.606 -0.690  0.432
>>>>>    >
>>>>>
>>>>> Could anyone give me some idea of why this unidentifiability problem
>>>>> happens and how to fix it?
>>>>> Any help would be appreciated.
>>>>>
>>>>> Best regards,
>>>>> Takahiro
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>
>>>>           [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> --
>> Takahiro Fushimi
>> Columbia University in the City of New York
>> Master of Arts in Statistics
>>

-- 
Takahiro Fushimi
Columbia University in the City of New York
Master of Arts in Statistics



More information about the R-sig-mixed-models mailing list