inst/doc/agridat_mixed_model_example.R

In agridat: Agricultural Datasets

## ----setup------------------------------------------------------------------------------
library(agridat)
data(crowder.seeds)
dat <- crowder.seeds

## ----plot-------------------------------------------------------------------------------
libs(lattice)
dotplot(germ/n~gen|extract, dat, main="crowder.seeds", xlab="gen")

## ----brms, eval=FALSE-------------------------------------------------------------------
#  if(require(brms)){
#    m1.brms <- brms::brm( germ|trials(n)~ gen*extract,
#       data = dat,
#       family = binomial,
#       chains=3, iter=3000, warmup=1000)
#    summary(m1.brms)
#    # round(  summary(m1.brms)$fixed[,1:4] , 2)
#    #                        Estimate Est.Error l-95% CI u-95% CI
#    # Intercept                 -0.42      0.18    -0.77    -0.06
#    # genO75                    -0.14      0.22    -0.56     0.29
#    # extractcucumber            0.55      0.25     0.07     1.05
#    # genO75:extractcucumber     0.77      0.30     0.18     1.36
#  }

## ----glm,eval=FALSE---------------------------------------------------------------------
#  # ----- GLM.
#  # family=binomial() fixes dispersion at 1
#  # family=quasibinomial() estimates dispersion, had larger std errors
#  m1.glm <- glm(cbind(germ,n-germ) ~ gen*extract,
#                data=dat,
#                #family="binomial",
#                family=quasibinomial()
#                )
#  summary(m1.glm)
#  ## round(summary(m1.glm)$coef,2)
#  ##                        Estimate Std. Error t value Pr(>|t|)
#  ## (Intercept)               -0.41       0.25   -1.64     0.12
#  ## genO75                    -0.15       0.30   -0.48     0.64
#  ## extractcucumber            0.54       0.34    1.58     0.13
#  ## genO75:extractcucumber     0.78       0.42    1.86     0.08

## ----eval=FALSE-------------------------------------------------------------------------
#    # ----- Stan using pre-built models from rstanarm
#    libs(tidyverse, rstan, rstanarm,bayesplot)
#    set.seed(42)
#    m1.stan <- stan_glm( cbind(germ,n-germ) ~ gen*extract,
#                        data=dat,
#                        family = binomial(link="logit")  )
#    summary(m1.stan)
#    ## round(posterior_interval(m1.stan, prob=.90),3)
#    #                            5%    95%
#    # (Intercept)            -0.728 -0.115
#    # genO75                 -0.506  0.243
#    # extractcucumber         0.133  0.947
#    # genO75:extractcucumber  0.255  1.267
#  
#    libs(bayesplot)
#    mcmc_areas(m1.stan, prob = 0.9) +
#      ggtitle("Posterior distributions",
#              "with medians and 95 pct intervals")
#  

## ----asreml,eval=FALSE------------------------------------------------------------------
#  if(require(asreml)){
#    m1.asreml <- asreml(germ ~ gen*extract,
#                        data=dat,
#                        random= ~ plate,
#                        family=asr_binomial(dispersion=1, total=n))
#    summary(m1.asreml)
#    ##
#    ##                          effect
#    ## (Intercept)               -0.47
#    ## gen_O73                    0.00
#    ## gen_O75                   -0.08
#    ## extract_bean               0.00
#    ## extract_cucumber           0.51
#    ## gen_O73:extract_bean       0.00
#    ## gen_O73:extract_cucumber   0.00
#    ## gen_O75:extract_bean       0.00
#    ## gen_O75:extract_cucumber   0.83
#  
#  }

## ----glmmpql,eval=FALSE-----------------------------------------------------------------
#  # --- GLMM.  Assumes Gaussian random effects
#  libs(MASS)
#  m1.glmm <- glmmPQL(cbind(germ, n-germ) ~ gen*extract,
#                     random= ~1|plate,
#                     family=binomial(), data=dat)
#  summary(m1.glmm)
#  ## round(summary(m1.glmm)$tTable,2)
#  ##                        Value Std.Error DF t-value p-value
#  ## (Intercept)            -0.44      0.25 17   -1.80    0.09
#  ## genO75                 -0.10      0.31 17   -0.34    0.74
#  ## extractcucumber         0.52      0.34 17    1.56    0.14
#  ## genO75:extractcucumber  0.80      0.42 17    1.88    0.08
#  

## ----glmmtmb,eval=FALSE-----------------------------------------------------------------
#  libs(glmmTMB)
#  m1.glmmtmb <- glmmTMB(cbind(germ, n-germ) ~ gen*extract + (1|plate),
#                        data=dat,
#                        family=binomial)
#  round(summary(m1.glmmtmb)$coefficients$cond , 2)
#  ##                        Estimate Std. Error z value Pr(>|z|)
#  ## (Intercept)               -0.45       0.22   -2.03     0.04
#  ## genO75                    -0.10       0.28   -0.35     0.73
#  ## extractcucumber            0.53       0.30    1.74     0.08
#  ## genO75:extractcucumber     0.81       0.38    2.11     0.04

## ----hglm,eval=FALSE--------------------------------------------------------------------
#  # ----- HGML package. Beta-binomial with beta-distributed random effects
#  if(require(hglm)){
#    m1.hglm <- hglm(fixed= germ/n ~ I(gen=="O75")*extract, weights=n, data=dat,
#                    random=~1|plate, family=binomial(), rand.family=Beta(),
#                    fix.disp=1)
#    summary(m1.hglm)
#    # round(summary(m1.hglm)$FixCoefMat,2)
#    ##                                     Estimate Std. Error t-value Pr(>|t|)
#    ## (Intercept)                            -0.47       0.24   -1.92     0.08
#    ## I(gen == "O75")TRUE                    -0.08       0.31   -0.25     0.81
#    ## extractcucumber                         0.51       0.33    1.53     0.16
#    ## I(gen == "O75")TRUE:extractcucumber     0.83       0.43    1.92     0.08
#  }

## ----inla,eval=FALSE--------------------------------------------------------------------
#  if(require(INLA)){
#    #gen,extract are fixed. plate is a random effect
#    #Priors for hyper parameters. See: inla.doc("pc.prec")
#    hyper1 = list(theta = list(prior="pc.prec", param=c(1,0.01)))
#    m1.inla = inla(germ ~ gen*extract + f(plate, model="iid", hyper=hyper1),
#                   data=crowder.seeds,
#                   family="binomial", Ntrials=n,
#                   control.family=list(control.link=list(model="logit")))
#    round( summary(m1.inla)$fixed, 2)
#    ##                         mean   sd 0.025quant 0.5quant 0.975quant  mode kld
#    ## (Intercept)            -0.47 0.24      -0.96    -0.46       0.00 -0.46   0
#    ## genO75                 -0.08 0.31      -0.68    -0.09       0.54 -0.09   0
#    ## extractcucumber         0.53 0.33      -0.13     0.53       1.18  0.53   0
#    ## genO75:extractcucumber  0.82 0.43      -0.01     0.82       1.69  0.82   0
#  }

## ----rjags,eval=FALSE-------------------------------------------------------------------
#  
#  # JAGS/BUGS.  See https://mathstat.helsinki.fi/openbugs/Examples/Seeds.html
#  # Germination rate depends on p, which is a logit of a linear predictor
#  # based on genotype and extract, plus random deviation to intercept
#  
#  # To match the output on the BUGS web page, use: dat$gen=="O73".
#  # We use dat$gen=="O75" to compare with the parameterization above.
#  
#  jdat =list(germ = dat$germ, n = dat$n,
#  root = as.numeric(dat$extract=="cucumber"),
#  gen = as.numeric(dat$gen=="O75"),
#  nobs = nrow(dat))
#  
#  jinit = list(int = 0, genO75 = 0, extcuke = 0, g75ecuke = 0, tau = 10)
#  
#  # Use logical names (unlike BUGS documentation)
#  mod.bug =
#  "model {
#    for(i in 1:nobs) {
#      germ[i] ~ dbin(p[i], n[i])
#      b[i] ~ dnorm(0.0, tau)
#      logit(p[i]) <- int + genO75 * gen[i] + extcuke * root[i] +
#                     g75ecuke * gen[i] * root[i] + b[i]
#    }
#    int ~ dnorm(0.0, 1.0E-6)
#    genO75 ~ dnorm(0.0, 1.0E-6)
#    extcuke ~ dnorm(0.0, 1.0E-6)
#    g75ecuke ~ dnorm(0.0, 1.0E-6)
#    tau ~ dgamma(0.001, 0.001)
#    sigma <- 1 / sqrt(tau)
#  }"
#  
#  libs(rjags)
#  oo <- textConnection(mod.bug)
#  j1 <- jags.model(oo, data=jdat, inits=jinit, n.chains=1)
#  close(oo)
#  
#  c1 <- coda.samples(j1, c("int","genO75","g75ecuke","extcuke","sigma"),
#                     n.iter=20000)
#  summary(c1) # Medians are very similar to estimates from hglm
#  # libs(lucid)
#  # print(vc(c1),3)
#  ##             Mean    SD    2.5%  Median   97.5%
#  ## extcuke   0.543  0.331 -0.118   0.542   1.2
#  ## g75ecuke  0.807  0.436 -0.0586  0.802   1.7
#  ## genO75   -0.0715 0.309 -0.665  -0.0806  0.581
#  ## int      -0.479  0.241 -0.984  -0.473  -0.0299
#  ## sigma     0.289  0.142  0.0505  0.279   0.596
#  
#  # Plot observed data with HPD intervals for germination probability
#  c2 <- coda.samples(j1, c("p"), n.iter=20000)
#  hpd <- HPDinterval(c2)[[1]]
#  med <- summary(c2, quantiles=.5)$quantiles
#  fit <- data.frame(med, hpd)
#  
#  libs(latticeExtra)
#  obs <- dotplot(1:21 ~ germ/n, dat,
#                 main="crowder.seeds", ylab="plate",
#                 col=as.numeric(dat$gen), pch=substring(dat$extract,1))
#  obs + segplot(1:21 ~ lower + upper, data=fit, centers=med)
#  

## ----R2jags,eval=FALSE------------------------------------------------------------------
#  libs("agridat")
#  libs("R2jags")
#  dat <- crowder.seeds
#  
#  # To match the output on the BUGS web page, use: dat$gen=="O73".
#  # We use dat$gen=="O75" to compare with the parameterization above.
#  jdat =list(germ = dat$germ, n = dat$n,
#             root = as.numeric(dat$extract=="cucumber"),
#             gen = as.numeric(dat$gen=="O75"),
#             nobs = nrow(dat))
#  
#  jinit = list(list(int = 0, genO75 = 0, extcuke = 0, g75ecuke = 0, tau = 10))
#  
#  mod.bug = function() {
#    for(i in 1:nobs) {
#      germ[i] ~ dbin(p[i], n[i])
#      b[i] ~ dnorm(0.0, tau)
#      logit(p[i]) <- int + genO75 * gen[i] + extcuke * root[i] +
#        g75ecuke * gen[i] * root[i] + b[i]
#    }
#    int ~ dnorm(0.0, 1.0E-6)
#    genO75 ~ dnorm(0.0, 1.0E-6)
#    extcuke ~ dnorm(0.0, 1.0E-6)
#    g75ecuke ~ dnorm(0.0, 1.0E-6)
#    tau ~ dgamma(0.001, 0.001)
#    sigma <- 1 / sqrt(tau)
#  }
#  
#  parms <- c("int","genO75","g75ecuke","extcuke","sigma")
#  
#  j1 <- jags(data=jdat, inits=jinit, parms, model.file=mod.bug,
#             n.iter=20000, n.chains=1)
#  print(j1)
#  ##          mu.vect sd.vect   2.5%    25%     50%     75%   97.5%
#  ## extcuke    0.519   0.325 -0.140  0.325   0.531   0.728   1.158
#  ## g75ecuke   0.834   0.429 -0.019  0.552   0.821   1.101   1.710
#  ## genO75    -0.096   0.305 -0.670 -0.295  -0.115   0.089   0.552
#  ## int       -0.461   0.236 -0.965 -0.603  -0.455  -0.312   0.016
#  ## sigma      0.255   0.148  0.033  0.140   0.240   0.352   0.572
#  ## deviance 103.319   7.489 90.019 98.010 102.770 108.689 117.288
#  
#  traceplot(as.mcmc(j1))
#  densityplot(as.mcmc(j1))
#  HPDinterval(as.mcmc(j1))
#  
#  }