Testing get_draws_mixed_model() function

For this, I copied and pasted the code from the previous version, and changed ‘se_tau2’ to NA

# get_draws_mixed_model <- function(
#   n_draws, coefs, vcov, tau2, se_tau2) {
  
  dev <- TRUE
  
  if (dev) {
    n_draws <- 10000
    coefs <- c(intercept = -0.4238, age_mid = 0.0301) # must be named vector
    vcov <- rbind(c(0.0264985, -0.0003751), c(-0.0003751, 0.0000060))
    tau2 <- 0.0785
    se_tau2 <- NA
  }
  
  require(MASS)

## Loading required package: MASS

  # take draws from the variance covariance matrix to capture intercept and age effects
  draws_fe <- mvrnorm(n = n_draws, mu = coefs, Sigma = vcov)
  
  head(draws_fe) # look at the first six rows (there are 10000)

##       intercept    age_mid
## [1,] -0.2730491 0.02809247
## [2,] -0.5754017 0.03185617
## [3,] -0.2331068 0.02729583
## [4,] -0.3080840 0.02974443
## [5,] -0.1438975 0.02575644
## [6,] -0.6333532 0.03394269

# in this case, se_tau2 has the value of NA, so the 'if' statement is executed
  if (is.na(se_tau2)) {
    draws_tau2 <- rep(tau2, times = n_draws)
  } else {
    draws_tau2 <- rnorm(n = n_draws, mean = tau2, sd = se_tau2)
  }

  str(draws_tau2) # 'draws_tau2' is a n=10000 vector of tau-squared values

##  num [1:10000] 0.0785 0.0785 0.0785 0.0785 0.0785 0.0785 0.0785 0.0785 0.0785 0.0785 ...

  hist(draws_tau2) # all the values are the same

  draws_tau_tmp <- suppressWarnings({ sqrt(draws_tau2) })
  str(draws_tau_tmp) # 'draws_tau_tmp' is a n=10000 vector of tau values, with no NaN values

##  num [1:10000] 0.28 0.28 0.28 0.28 0.28 ...

  table(is.nan(draws_tau_tmp))

## 
## FALSE 
## 10000

  # convert NaN to zero (superfluous in this case)
  draws_tau_tmp2 <- ifelse(is.nan(draws_tau_tmp), 0, draws_tau_tmp)
  str(draws_tau_tmp2)

##  num [1:10000] 0.28 0.28 0.28 0.28 0.28 ...

  table(is.nan(draws_tau_tmp2))

## 
## FALSE 
## 10000

  # get final draws of tau
  draws_tau <- rnorm(n = n_draws, mean = 0, sd = draws_tau_tmp2)
  str(draws_tau)

##  num [1:10000] -0.0949 -0.6957 -0.1541 -0.0115 0.3084 ...

  # the first draw used the first element of 'draws_tau_tmp2', the second used the second, etc.
  # here's a check that the function works element-wise like I think it does
  rnorm(n = 5, mean = 0, sd = c(1, 999, 1, 999, 1))

## [1]    0.01343896 -313.21986232    0.82937636  -63.32913653   -0.51840231

  # the final output
  out <- as.data.frame(cbind(draws_fe, draws_tau))
  str(out)

## 'data.frame':    10000 obs. of  3 variables:
##  $ intercept: num  -0.273 -0.575 -0.233 -0.308 -0.144 ...
##  $ age_mid  : num  0.0281 0.0319 0.0273 0.0297 0.0258 ...
##  $ draws_tau: num  -0.0949 -0.6957 -0.1541 -0.0115 0.3084 ...

# }

# I add the levels of 'x_age' for which I want predictions
# I'll aggregate 5000 draws for each age in this example
out$x_age <- rep(c(25, 100), each = 5000)
str(out)

## 'data.frame':    10000 obs. of  4 variables:
##  $ intercept: num  -0.273 -0.575 -0.233 -0.308 -0.144 ...
##  $ age_mid  : num  0.0281 0.0319 0.0273 0.0297 0.0258 ...
##  $ draws_tau: num  -0.0949 -0.6957 -0.1541 -0.0115 0.3084 ...
##  $ x_age    : num  25 25 25 25 25 25 25 25 25 25 ...

require(dplyr)

## Loading required package: dplyr

## Warning: package 'dplyr' was built under R version 3.4.4

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# generate the predictions
out2 <- out %>%
  mutate(pred = intercept + x_age*age_mid + draws_tau)

## Warning: package 'bindrcpp' was built under R version 3.4.4

str(out2)

## 'data.frame':    10000 obs. of  5 variables:
##  $ intercept: num  -0.273 -0.575 -0.233 -0.308 -0.144 ...
##  $ age_mid  : num  0.0281 0.0319 0.0273 0.0297 0.0258 ...
##  $ draws_tau: num  -0.0949 -0.6957 -0.1541 -0.0115 0.3084 ...
##  $ x_age    : num  25 25 25 25 25 25 25 25 25 25 ...
##  $ pred     : num  0.334 -0.475 0.295 0.424 0.808 ...

# aggregate the draws to get point estimate and UI

out3 <- out2 %>%
  group_by(x_age) %>%
  summarize(
  pred_mean = mean(pred),
  pred_lo = quantile(pred, probs = 0.025),
  pred_hi = quantile(pred, probs = 0.975))

print(out3)

## # A tibble: 2 x 4
##   x_age pred_mean pred_lo pred_hi
##   <dbl>     <dbl>   <dbl>   <dbl>
## 1  25.0     0.326  -0.239   0.907
## 2 100       2.58    2.01    3.17