run with error when the nonlinear predictor with division and how to introduce random effect for parameter in nonlinear predictors

[ZSRNOG]

Hello!

I run into two tricky problems with inlabru

1, The nonlinear predictor includes division, e.g. 1/z 2, Introducing random for the parameter in nonlinear prector,

` library(inlabru)

Case 1, run correctly

z <- 3 z1 <- 2 input.df <- data.frame(x = cos(1:10)) input.df <- within(input.df, y <- z1+ exp(z)x + rnorm(10, mean = 0, sd = 0.1)) lik <- like( family = "gaussian", data = input.df, formula = y ~ z1 + exp(z) x + 0 ) fit <- bru(~ z(1) + z1(1) + 0 , lik) fit$summary.fixed

Case 2 include (1/z1) run with error "Error in rep(k, sum(nonzero)) : 'times'is incorrect"

z <- 3 z1 <- 2 input.df <- data.frame(x = cos(1:10)) input.df <- within(input.df, y <- 1/z1 + z1+ exp(z)x + rnorm(10, mean = 0, sd = 0.1)) lik <- like( family = "gaussian", data = input.df, formula = y ~ 1/z1 + z1 + exp(z) x + 0 ) fit <- bru(~ z(1) + z1(1) + 0 , lik) fit$summary.fixed

Case 3, how can I set the random effect for z

z_mean <- 3 z_sd <- 0.5 z1 <- 2 zs = rnorm(3, z_mean, z_sd) input.df <- data.frame(x = cos(1:10)) input.df1 <- within(input.df, y <- z1+ exp(zs[1])x + rnorm(10, mean = 0, sd = 0.1)) input.df2 <- within(input.df, y <- z1+ exp(zs[2])x + rnorm(10, mean = 0, sd = 0.1)) input.df3 <- within(input.df, y <- z1+ exp(zs[3])*x + rnorm(10, mean = 0, sd = 0.1))

input.df <- rbind(input.df1,input.df2,input.df3) input.df$unit <- rep(paste('unit',1:3),rep(10,3))

lik <- like( family = "gaussian", data = input.df, formula = y ~ z1 + exp(z) * x + f(unit,'iid') + 0 ) fit <- bru(~ z(1) + z1(1) + 0 , lik) fit$summary.fixed `

[finnlindgren]

3. To make an iid random effect for an indexed input z, use component specification z(z, model="iid").
4. For the predictor expression 1/z1 + z1 to do something reasonable, inlabru needs a better linearisation starting point than the default, 0. From the help text for bru_options:

   bru_initial
   An inla object returned from previous calls of INLA::inla, bru() or lgcp(), or a list
   of named vectors of starting values for the latent variables. This will be used as
   a starting point for further improvement of the approximate posterior.

   Something like bru(..., options = list(bru_initial = list(z1 = 1))) should help. Not however that since the posterior is multimodal and has singularities due to the potential of division by zero, this may not entirely solve the problem as stated, but I'm guessing this was just an example triggering the problem and not that actual model you wanted to run.

[ZSRNOG]

Thank you Finnlindgren. Following your suggestions, I solve the first problem by specifying the bru_initial. But, I'm still not sure how to add random effects.

case 1

` z_mean <- -2 z_sd <- 0.1 z1 <- 5 zs = rnorm(6, z_mean, z_sd) input.df <- data.frame(x = cos(1:20)) input.df1 <- within(input.df, y <- 1/z1+ exp(zs[1])x + rnorm(20, mean = 0, sd = 0.1)) input.df2 <- within(input.df, y <- 1/z1+ exp(zs[2])x + rnorm(20, mean = 0, sd = 0.1)) input.df3 <- within(input.df, y <- 1/z1+ exp(zs[3])x + rnorm(20, mean = 0, sd = 0.1)) input.df4 <- within(input.df, y <- 1/z1+ exp(zs[4])x + rnorm(20, mean = 0, sd = 0.1)) input.df5 <- within(input.df, y <- 1/z1+ exp(zs[5])x + rnorm(20, mean = 0, sd = 0.1)) input.df6 <- within(input.df, y <- 1/z1+ exp(zs[6])x + rnorm(20, mean = 0, sd = 0.1))

input.df <- rbind(input.df1,input.df2,input.df3,input.df4,input.df5,input.df6) input.df$unit <- rep(paste('unit',1:6),rep(20,6))

lik <- like( family = "gaussian", data = input.df, formula = y ~ 0+ 1/(z1+1e-5) + exp(z) * x, options = list(bru_initial = list(z1 = 5, z=-2)) ) fit <- bru(~ z1(1) + z(1) +0 , lik) fit$summary.fixed library(brinla) bri.hyperpar.summary(fit) `

I tried to introduce random effects through the following method, but the result seems implausible.

fit1 <- bru(~ z1(1) + z(1) +0 + f(unit, model='iid'), lik) fit1$summary.fixed library(brinla) bri.hyperpar.summary(fit1)

[finnlindgren]

library(inlabru)
#> Loading required package: ggplot2
#> Loading required package: sp
z_mean <- -2
z_sd <- 0.1
z1 <- 5
zs = rnorm(6, z_mean, z_sd)
input.df <- data.frame(x = rep(cos(1:20), times = 6),
                       unit = rep(1:6, each = 20))
input.df <- within(input.df, y <- 1/z1+ exp(zs[unit])*x + rnorm(20*6, mean = 0, sd = 0.1))

comp <- ~ 0 + z1(1) + z(unit, model = "iid")
lik <- like(
  family = "gaussian", data = input.df,
  formula = y ~ 1/z1 + exp(z) * x
)
fit <- bru(comp , lik, options = list(bru_initial = list(z1 = 5, z=-2)))
summary(fit)
#> inlabru version: 2.3.1.9000
#> INLA version: 21.04.16
#> Components:
#>   z1: Model types main='linear', group='exchangeable', replicate='iid'
#>   z: Model types main='iid', group='exchangeable', replicate='iid'
#> Likelihoods:
#>   Family: 'gaussian'
#>     Data class: 'data.frame'
#>     Predictor: y ~ 1/z1 + exp(z) * x
#> Time used:
#>     Pre = 0.181, Running = 0.187, Post = 0.0979, Total = 0.465 
#> Fixed effects:
#>     mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> z1 5.367 0.242       4.89    5.367      5.843 5.367   0
#> 
#> Random effects:
#>   Name     Model
#>     z IID model
#> 
#> Model hyperparameters:
#>                                            mean     sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 120.012 15.840     91.313   119.23
#> Precision for z                           0.298  0.154      0.088     0.27
#>                                         0.975quant    mode
#> Precision for the Gaussian observations    153.428 117.907
#> Precision for z                              0.676   0.209
#> 
#> Expected number of effective parameters(stdev): 6.88(0.07)
#> Number of equivalent replicates : 17.45 
#> 
#> Deviance Information Criterion (DIC) ...............: -222.80
#> Deviance Information Criterion (DIC, saturated) ....: -269.21
#> Effective number of parameters .....................: 8.06
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -222.56
#> Effective number of parameters .................: 7.68
#> 
#> Marginal log-Likelihood:  81.03 
#> Posterior marginals for the linear predictor and
#>  the fitted values are computed

^{Created on 2021-05-03 by the reprex package (v2.0.0)}

[finnlindgren]

Note that the iid model has mean zero, so your z_mean aspect to the random effects isn't captured properly by that. You may want to model the mean of zs explicitly instead for this example.