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

[Ben Bolker]

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.
> 
> 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).

  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