tests/test1.R

library(dglm)

n <- 1000

TestFunc <- function(mean.formula, var.formula, model.df) {
  
  my.glm <- glm(formula = mean.formula,
                data = model.df)
  
  x <- predict(my.glm, se.fit = TRUE)
  
  my.dglm <- dglm(formula = mean.formula,
                  dformula = var.formula,
                  data = model.df)
  
  y <- predict(my.dglm, se.fit = TRUE)
  
  z <- predict(my.dglm$dispersion.fit, se.fit = TRUE)
  
  return(list(glm = x, dglm.mean = y, dglm.var = z))
}


a <- runif(n)
b <- runif(n)
c <- runif(n)
d <- runif(n)

betas <- c(runif(n = 3, min = -10, max = 10), runif(n = 3, min = -1, max = 1))

mu <- betas[1] + a*betas[2] + b*betas[3]
var <- exp(betas[4] + c*betas[5] + d*betas[6])
y <- mu + rnorm(n = n, sd = sd(var))

my.df <- data.frame(y, a, b, c, d)

mean.formula <- as.formula('y ~ a + b')
var.formula <- as.formula('~ c + d')

l <- TestFunc(mean.formula = mean.formula, var.formula = var.formula, model.df = my.df)

par(mfrow = c(2, 2))
plot(mu, l$dglm.mean$fit); abline(0, 1)
legend(x = 'topleft', legend = paste(c('mu:', 'a', 'b'), round(betas[1:3], 1)), bty = 'n')
plot(var, exp(l$dglm.var$fit/l$dglm.var$residual.scale)); abline(0, 1)
legend(x = 'topleft', legend = paste(c('mu:', 'c:', 'd:'), round(betas[4:6], 2)))

my.dglm <- dglm(formula = mean.formula, dformula = var.formula, data = my.df, method = 'reml')

plot(mu, predict(my.dglm))
plot(sqrt(var), predict(my.dglm$dispersion.fit))


anova(my.dglm)

Try the dglm package in your browser

Any scripts or data that you put into this service are public.

dglm documentation built on May 29, 2017, 10:42 p.m.