[ADMB Users] Integrating R and ADMB to generalize TPL files

[dave fournier]

On 12-11-07 09:29 AM, Ben Bolker wrote:
> On 12-11-07 11:58 AM, dave fournier wrote:
>> Here is an example of why I dislike the R approach
>> so much.
>>
>> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q4/019349.html
>>
>> I was curious to find out why introducing to additional parameters for the
>> variance of the random effects would produce no improvement in the
>> model.  So first I wanted to run glmmadmb as seth had done.
>>
>> I put his data in a file and read it into R.
>>
>> Then I ran his model by grabbing his call to glmmadmb.
>> I was surprised to see that I got a different result from Seth.
>>
>> Examination of the design matrix and I realized that one of the columns of
>> the data needed to be changed to a factor.  If I had not got the
>> different result I would have thought that the analysis was correct.
>> That was the first gotcha.  With R you never really know what
>> you are doing.   This was for the simpler model with paramerter for
>> the variance of the RE's.
>    I did see that, and I agree there is a price to pay for convenience.

It is not convenient to effortless get to the wrong answer,
just a bit of a waste of time.
>> Then I used Seth's R script to run the more complicated model.
>> He was right that the likelihood did not change.  Strange!
>>
>> So I examined the design matrix for the RE's.
>> The problem was obvious.  rather than a different
>> parameter for each time period the design matrix had
>> an overall mean and differences from it.  While this
>> would be the same for fixed effects it is not equivalent
>> for random effects because it assumes a different covariance
>> structure.
>>
>> I hacked together the correct model by hand and ran it on his data.
>    I agree that the model Seth was specifying was that there is variation
> among trees in their time-1 response; variation in the difference
> between time 1 and time 2; and variation in the difference between time
> 1 and time 3 (unless I'm wrong, which has been known to happen).
All I know is that he said

The data suggest that variance may differ
sharply before and after treatment (factor variable 'time', values 1/2/3),

I don't see anything about differences just that the variance should 
depend on time.
>
>    It sounds like you fitted a model with variance around the time-1
> response, variance around the time-2 response, and variance around the
> time-3 response.  I *think* one could specify that via
> random=~time-1|subject, but I haven't tried it.
>
>    I agree that unless one specifies corStruct="full" these are not
> equivalent models.
>
>    Do we know which of these two models answers the question Seth wanted
> to ask?
>
>
>> Results are:
>>
>> # Number of parameters = 8  Objective function value = 255.508 Maximum
>> gradient component = 1.44753e-06
>> # pz:
>> 0.000100000000000
>> # beta:
>>   -10.7785676169 -2.55267985235 -4.72926437596 6.91095066201
>> # tmpL:
>>   0.352458239856 0.489012745656 0.496221176099
>> # tmpL1:
>>   0.000100000000000
>> # log_alpha:
>> 0.591420283490
>> the original fit is
>>
>>      # Number of parameters = 6  Objective function value = 255.758
>> Maximum gradient component = 2.23893e-05
>> # pz:
>> 0.000100000000000
>> # beta:
>>   -10.4061157821 -1.87212686310 -3.41511253652 8.12797421786
>> # tmpL:
>>   0.431551628285
>> # tmpL1:
>>   0.000100000000000
>> # log_alpha:
>> 0.593869268695
>>
>> So it appears that the variance is larger in the second and third time
>> period, but not
>> significantly so.
>>
>> I welcome your snide comments.
>>
>>
>>
>> _______________________________________________
>> Users mailing list
>> Users at admb-project.org
>> http://lists.admb-project.org/mailman/listinfo/users