Random effects standard deviations

[John Sibert]

Is it reasonable to ask about the estimated standard deviations of estimated random effects? If so, how does one compute the estimated standard deviations of estimated random effects?

Thanks,

John

[Gavin Fay]

Hi John,

You can extract them very easily from the sdreport().

e.g. if the result of MakADFun() is "obj", then:

sdr <- sdreport(obj)
summary(sdr,select="random")

will give you the best estimates for the random effects plus their std. errors.

Cheers,

Gavin

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/8ca81aa8-c588-4561-9bc3-54365696abcao%40googlegroups.com.

--

Gavin Fay, Ph.D.
Fisheries & Ecosystem management, Ecological modeling, Stock assessment
gfa...@gmail.com
@gavin_fay
www.gavinfay.com

[John Sibert]

Hi Gavin,

Thanks much for the tip. It got me a long way down the road.

I'm still curious how these standard errors are calculated and how they should be interpreted.

This is all part of my corona virus social isolation activities.

Cheers,

John

On Wednesday, August 5, 2020 at 5:16:52 PM UTC-7, Gavin Fay wrote:

Hi John,

You can extract them very easily from the sdreport().

e.g. if the result of MakADFun() is "obj", then:

sdr <- sdreport(obj)
summary(sdr,select="random")

will give you the best estimates for the random effects plus their std. errors.

Cheers,

Gavin

On Wed, Aug 5, 2020 at 7:24 PM John Sibert <johnr...@gmail.com> wrote:

Is it reasonable to ask about the estimated standard deviations of estimated random effects? If so, how does one compute the estimated standard deviations of estimated random effects?

Thanks,

John

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+unsubscribe@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/8ca81aa8-c588-4561-9bc3-54365696abcao%40googlegroups.com.

[Kasper Kristensen]

Hi John,

The standard errors are calculated as described under 'sdreport' here:

https://cran.r-project.org/web/packages/TMB/TMB.pdf

A joint covariance matrix of estimation error is constructed and all marginal standard errors are derived from this. The approach is the same as that of ADMB.
Note that standard errors of random effects are based on V(hat(u) - u), *not* V(hat(u)) as wrongly documented in TMB versions <= 1.7.16.

A standard frequentist interpretation is valid. However, because several approximations are assumed to derive the covariance formula, you may want to perform a simulation study to check the coverage degree:

1. Fix the true parameter vector theta
2. repeat {
  * Simulate random effects and data
  * Estimate parameters and random effects
  * Construct 95% confidence intervals of selected random effects based on sdreport output
}
Coverage should be close to 95%.

Enjoy your time in isolation :)

Cheers,
Kasper

[Mark Maunder]

How do you simulate the random effects?

I was trying to find an easy way to sample from the random effects in TMB. The SIMULATE option allows you to simulate the data. However, we are trying to implement a Frequentist version of posterior predictive checks which requires also sampling the parameters from their uncertainty distribution. In this case we are using a normal approximation using the MLEs and the variance-covariance matrix. We have it implemented for a fixed effects model, but now want to apply it to a random effects (state-space) model. We can call the $simulate(par=sample_pars) function with a random set of parameters. However, we need to generate this set of random parameters, which includes the random effects (actually all that is needed is the random effect representing biomass over time in a theta-logistic type population dynamics model). Do you have an example of sampling from the random effects or is there a function in TMB that automates this.

Kasper replied (to my e-mail, which is copied above)

Please have a look at:

http://kaskr.github.io/adcomp/_book/Simulation.html#overview-of-simulation-methods-in-tmb

It sounds to me like section 13.2.3 has the information you ask for?

Simulating the data is easy (where PARAMETER_VECTOR(logB) and B=exp(logB);)

SIMULATE {

  Type mu=0;

    for(int year=0;year<Nyears;year++){

        Y[year]=B[year]*q*exp(rnorm(mu, sd)-0.5*sd*sd);

    }

  REPORT(C);

  REPORT(Y);

  }

The question is how to simulate from the random effect logB.

You could supply the parameter values to the simulate call obj$simulate(par=par_sample)

which would require sampling the parameters from a multivariate normal using the MLEs and the variance covariance matrix

pars_sample <- MASS::mvrnorm(1, mu = pl, Sigma = cov)

sdr<-sdreport(obj)

pl <- as.list(sdr,"Est")

****** But how do we get cov? *******

A possible alternative is in the SIMULATE functionaity of TMB

Under 13.2.3 Advanced examples in http://kaskr.github.io/adcomp/_book/Simulation.html#overview-of-simulation-methods-in-tmb

You might be able to use 

MVNORM_t<Type> N_0_Sigma(Sigma); 

N_0_Sigma.simulate()

in SIMULATE {} 

***** but how do you access Sigma with a TMB program? ****** 

I would appreciate any help on working out how to do these simulations.

Automating the simulation of the parameters and random effects within the SIMULATE function using a multivariate normal would be a great way to produce Frequentist versions of posterior predictive checks. 

Thanks,

Mark

[Mark Maunder]

The two questions I am looking for answers from in the previous post are

1. How do you get the covariance matrix for the random effects (or the random effects and the fixed effects combined, or the derived parameters) from a converged model?

 

2. How do you access and use the covariance matrix for the random effects (or ....) within the SIMULATE{} function in a TMB model?

Thanks,

Mark

On Wednesday, August 5, 2020 at 4:24:08 PM UTC-7, John Sibert wrote:

[James Thorson]

Mark (and anyone interested),

We have examples of how to do both of those in FishStatsUtils::simulate_data(.), code here: https://github.com/James-Thorson-NOAA/FishStatsUtils/blob/master/R/simulate_data.R

Jim

--

To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/9410f3fb-a8a9-45e5-9bb0-90bdab172927o%40googlegroups.com.

[Mark Maunder]

Thanks Jim,

I think this is what I need

newpar = rmvnorm_prec( mu=Obj$env$last.par.best, prec=fit$parameter_estimates$SD$jointPrecision, n.sims=1, random_seed=random_seed )[,1]

    # Simulate

    Obj$env$data$Options_list$Options['simulate_random_effects'] = FALSE

    Return = Obj$simulate( par=newpar, complete=TRUE )

rmvnorm_prec <- function(mu, prec, n.sims, random_seed ) {

    set.seed( random_seed )

    z = matrix(rnorm(length(mu) * n.sims), ncol=n.sims)

    L = Matrix::Cholesky(prec, super=TRUE)

    z = Matrix::solve(L, z, system = "Lt") ## z = Lt^-1 %*% z

    z = Matrix::solve(L, z, system = "Pt") ## z = Pt    %*% z

    z = as.matrix(z)

    return(mu + z)

  }

On Wednesday, August 5, 2020 at 4:24:08 PM UTC-7, John Sibert wrote:

[Mark Maunder]

I think I have to also do this

sdr<-sdreport(obj,getJointPrecision = TRUE)

newpar = rmvnorm_prec( mu=obj$env$last.par.best, prec=sdr$jointPrecision, n.sims=1, random_seed=12345 )[,1]

On Wednesday, August 5, 2020 at 4:24:08 PM UTC-7, John Sibert wrote:

[Craig Marsh]

Another resource is the R file found at  https://github.com/kaskr/adcomp/blob/master/tmb_examples/validation/randomwalkvalidation.R which follows the goodness of fit paper by Thygesen et. al (2017). Which shows how to simulate from the prior [Y,X] for the single sample from the posterior residuals.

Can I ask (although maybe not appropriate for this forum), why you would use the joint covariance of (fixed effects & random effects) for predictions of Y? Instead of using the fixed effect covariance for simulating parameter uncertainty for  just fixed effects and simulating the random effects from the assumed superpopulation distribution? Looking at Jim's R code, it seems the former is type=3. But the latter is not available, I was thinking the R-code would like something like,

# simulate Fixed effect parameters to account for fixed effect parameter error
Hess = Obj$env$spHess(par=Obj$env$last.par.best, random=F)
newfixed = rmvnorm_prec( mu=rep(0,nrow(Hess)), prec=Hess, n.sims=1, random_seed=random_seed )[,1]
# Set internal model to simulate Random effects, condition on simulated fixed effects
Obj$env$data$Options_list$Options['simulate_random_effects'] = TRUE
Return = Obj$simulate( par = newfixed, complete=TRUE )

Chairs,

C

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/CAMDbYva1RBeutPzdVUS%3Dx08JSPYu_NnpOZHFp23p%2BgYz_vgj%2Bg%40mail.gmail.com.

[Craig Marsh]

Sorry there would be a slight error in that implementation in that code, as obj$simulate() requires both fixed and random effects, not just fixed effects as I wrote above.

pars = Obj$env$last.par.best 

# simulate Fixed effect parameters to account for fixed effect parameter error
Hess = Obj$env$spHess(par=Obj$env$last.par.best, random=F)
newfixed = rmvnorm_prec( mu=rep(0,nrow(Hess)), prec=Hess, n.sims=1, random_seed=random_seed )[,1]

pars[!Obj$env$random] = newfixed # you would want to double check the correct mapping here!

# Set internal model to simulate Random effects, which will override the values in pars (depending on set up of SIMULATE{} call
Obj$env$data$Options_list$Options['simulate_random_effects'] = TRUE

Return = Obj$simulate( par = pars, complete=TRUE )  

[Mark Maunder]

Thanks Craig,

If I follow what you are describing. 

Rather than sampling from the joint fixed and random effects distribution, you are sampling from the fixed effects distribution first, then conditional on that random sample from the fixed effects you are sampling from the random effects distribution. I am not sure if these will produce the same outcome.

The quantity of interest is a function of both the random and fixed effects. I am simulating CPUE from the product of the (log) biomass which is a random effect and the catchability which is a fixed effect. Thats why I am using Jim's method type=3 that uses the precision matrix for both the fixed and random effects in the multivariate normal. 

On Wednesday, August 5, 2020 at 4:24:08 PM UTC-7, John Sibert wrote:

[Craig Marsh]

Yeah that was what I was describing, It was more a general thought/query (always dangerous!!), rather than expecting them to yield the same outcome.  

Just as a thought experiment. The standard errors and confidence intervals for fixed effects parameters from a frequentist perspective relate to estimates falling within an interval from hypothetical replicate data. So could we say that when we simulate values to capture this variability, we are in essence mimicking a replicate data experiment? 

If yes to the above, if you follow that thought, regarding a simulation as representing a replicate data experiment onto the random effects. I would not expect a realisation of the random variable from a replicated data experiment to be associated with previously estimated values (i.e. drawn from the joint fixed and random effect distribution), but a random draw from the distribution assumed in the model (i.e. B_t ~ lognormal(B_t-1, sd)

That was my thinking around the nuances of simulating data from mixed effect models. Would love to know peoples thoughts? Or if I am embarrassingly misinformed, it wouldn't be the first time and won't be the last time :)

--

To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/aa128ab3-7469-478a-91d7-0ac543c34720o%40googlegroups.com.

[Mark Maunder]

Craig,

What we are trying to do is apply a Frequentist version of posterior predictive checks.

See 

Besbeas, P., Morgan, B.J.T., 2014. Goodness-of-fit of integrated population models using calibrated simulation. Methods Ecol. Evol. 5, 1373-1382. https://doi.org/10.1111/2041-210X.12279

Bayesian posterior predictive checks are reviewed in 

Conn, P.B., Johnson, D.S., Williams, P.J., Melin, S.R., Hooten, M.B., 2018. A guide to Bayesian model checking for ecologists. Ecol. Monogr. 88, 526-542. https://doi.org/10.1002/ecm.1314

My understanding is that given the model is correct and the MLE estimates, we are evaluating the probability that the model produced the observed data. If we simulate the data from the MLE estimates many times, then the observed data should fall in the 95% distribution of the simulated results. We can actually do this for each data component separately so it is good for integrated models that use different data types. We just need a measure of the discrepancy between the model and each data type.

Even if the model is correct, the MLE may not be the "true" parameters, therefore we need to simulate the data from all possibly parameter values. So, we sample from the posterior distribution of the parameters in a Bayesian context and from the multivariate normal approximation or bootstrap distribution in a Frequentist context for each data set created. Some calibration may also be needed.  

Just like the uncertainty in the fixed effects, different realizations of the random effects may be the truth. Sampling from the Bayesian posterior automatically does this. We probably can't do la place approximation in the simulations to implement the Frequentist approach, so I think sampling them in the joint multivariate normal makes sense.   

Obviously, we are still working through this, so some of the logic might be a bit off.

 

On Wednesday, August 5, 2020 at 4:24:08 PM UTC-7, John Sibert wrote:

[James Thorson]

Mark,

I believe that this is what the oneStepPredict(.) feature in TMB is designed to test.  I recommend exploring that work-flow, which is designed to account for a variety of issues (nonlinearity, leverage defined with a mixed-modelling framework, etc.) that arise when doing goodness-of-fit tests in a hierarchical model.  FWIW, I have also started using package DHARMa to automate the visualization of output from diagnostics such as this.  

--

To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/69ac8202-f00c-414f-a57a-842a3011dabbo%40googlegroups.com.