# tests/test1.R In dglm: Double Generalized Linear Models

```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.