[R-sig-ME] Likelihood drops on adding random effect

[Murray Jorgensen]

I guess we may have a case of a sort of spike in the loglikelihood that 
indicates that the MLE (or at least the lmer estimate) for the model 
without the 'tree' effect is spurious. (Probably the sort of thing that 
would be smoothed away by any reasonable prior.) I really was looking 
more closely at the other fixed effect parameters, regarding the 
intercept as a bit of a nuisance. I hope to say more about the data set 
later. Maybe I will do a sort of profile around the intercept value and 
just fit the other parameters.

I would like to ask the list for some more help with the 'start =' 
parameter though because I don't seem to be quite using it correctly.

 > ip2 = c(0, 8.29931, 2.56368e-06,  0.00000)
 > f0128bit = lmer(fincr ~ icfac + (1|gap) + (1|box)+ (1|gap:box) + 
(1|tree), family = binomial,
+                  verbose = TRUE, start = ip2)
   0:     775.67797:  0.00000  8.29931 2.56368e-06  0.00000 -3.22238 
1.44743 -2.20257  1.02515
   1:     747.52247:  0.00000  5.04306 6.75518e-05  0.00000 -4.34098 
2.12417 -3.74289 -0.121576
   2:     736.84946:  0.00000  4.55878 0.000215665  0.00000 -4.39655 
1.52182 -3.63273 -0.115459
...
  45:     671.01076: 4.57910e-09 0.287720 0.350164 0.711054 -3.43481 
1.45346 -2.57398 0.293897
  46:     671.01076:  0.00000 0.287721 0.350162 0.711064 -3.43483 
1.45346 -2.57397 0.293955

That works well for initialising the variance parameters only.

 > f0128bit = lmer(fincr ~ icfac + (1|gap) + (1|box)+ (1|gap:box) + 
(1|tree), family = binomial,
+               verbose = TRUE, start = list( fixef = list( -9.28405, 
2.81300, -4.75935,  2.91080),
+ ST = list(0, 8.29931, 2.56368e-06,  0.00000))
+ )
Error: class(STnew[[i]]) == class(ST[[i]]) is not TRUE
In addition: Warning message:
In sort(names(start)) == sort(names(FL)) :
   longer object length is not a multiple of shorter object length

A failed attempt to initialise both fixed and random parameters.

 > f0128bit = lmer(fincr ~ icfac + (1|gap) + (1|box)+ (1|gap:box) + 
(1|tree), family = binomial,
+               verbose = TRUE, start = list( fixef = c( -9.28405, 
2.81300, -4.75935,  2.91080),
+ ST = c(0, 8.29931, 2.56368e-06,  0.00000))
+ )
Error: is.list(STnew) is not TRUE
In addition: Warning message:
In sort(names(start)) == sort(names(FL)) :
   longer object length is not a multiple of shorter object length

Another failure.

 > initpar = c(0, 8.29931, 2.56368e-06,  0.00000, -9.28405,  2.81300, 
-4.75935,  2.91080)
 > f0128bit = lmer(fincr ~ icfac + (1|gap) + (1|box)+ (1|gap:box) + 
(1|tree), family = binomial,
+                  verbose = TRUE, start = initpar)
   0:     696.35176:  1.15470 0.222222 0.182574 0.157135 -3.22238 
1.44743 -2.20257  1.02515
   1:     687.88422: 0.646498 0.349241 0.199295 0.375671 -4.02981 
1.29061 -2.23330  1.01588
  ...
  26:     671.01076:  0.00000 0.287718 0.350162 0.711060 -3.43482 
1.45346 -2.57397 0.293937
  27:     671.01076:  0.00000 0.287721 0.350163 0.711061 -3.43482 
1.45346 -2.57398 0.293936

Runs, but does not seem to use the supplied starting parameters.

Thanks, for your comments.

Murray

On 18/05/2012 5:48 a.m., Douglas Bates wrote:
> Look at the values of the coefficients and standard deviations that
> you are "converging" to.  Your intercept is -9.28, which, with a
> binomial family, corresponds to probabilities below 1e-4.  With icfac
> = fem the linear predictor is -9.28 - 4.76 = -14.04 corresponding to a
> probability of 8e-07.
>
> You are going to need to look at the data and the proportions of
> positives for different levels of icfac to see what would make sense.
> This problem will create a very ill-defined likelihood surface because
> the fitted values will lose sensitivity to the parameters when the
> probabilities are so extreme.  If you start extreme values you will
> never be able to converge.
>
> On Thu, May 17, 2012 at 12:37 AM, Murray Jorgensen<maj at waikato.ac.nz>  wrote:
>> PS  I also tried
>>
>> start = list( fixef = c( -9.28405,  2.81300, -4.75935,  2.91080),
>>         ST = c(0, 8.29931, 2.56368e-06,  0.00000))
>>
>> and
>>
>> start = list( fixef = list( -9.28405,  2.81300, -4.75935,  2.91080),
>>        ST = list(0, 8.29931, 2.56368e-06,  0.00000))
>>
>> to no avail.
>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

-- 
Dr Murray Jorgensen      http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: maj at waikato.ac.nz      majmurr at gmail.com         Fax 7 838 4155
Phone  +64 7 838 4773 wk    Home +64 7 825 0441   Mobile 021 0200 8350