ndata <- data.frame(x2 = runif(nn <- 2000))
# Note that coeff1 + coeff2 + coeff5 == 1. So try "multilogitlink".
myoffset <- 10
ndata <- transform(ndata,
coeff1 = 0.25, # "multilogitlink"
coeff2 = 0.25, # "multilogitlink"
coeff3 = exp(-0.5), # "loglink"
# "logofflink" link:
coeff4 = logofflink(+0.5, offset = myoffset, inverse = TRUE),
coeff5 = 0.50, # "multilogitlink"
coeff6 = 1.00, # "identitylink"
v2 = runif(nn),
v3 = runif(nn),
v4 = runif(nn),
v5 = rnorm(nn),
v6 = rnorm(nn))
ndata <- transform(ndata,
Coeff1 = 0.25 - 0 * x2,
Coeff2 = 0.25 - 0 * x2,
Coeff3 = logitlink(-0.5 - 1 * x2, inverse = TRUE),
Coeff4 = logloglink( 0.5 - 1 * x2, inverse = TRUE),
Coeff5 = 0.50 - 0 * x2,
Coeff6 = 1.00 + 1 * x2)
ndata <- transform(ndata,
y1 = coeff1 * 1 +
coeff2 * v2 +
coeff3 * v3 +
coeff4 * v4 +
coeff5 * v5 +
coeff6 * v6 + rnorm(nn, sd = exp(0)),
y2 = Coeff1 * 1 +
Coeff2 * v2 +
Coeff3 * v3 +
Coeff4 * v4 +
Coeff5 * v5 +
Coeff6 * v6 + rnorm(nn, sd = exp(0)))
# An intercept-only model
fit1 <- vglm(y1 ~ 1,
form2 = ~ 1 + v2 + v3 + v4 + v5 + v6,
normal.vcm(link.list = list("(Intercept)" = "multilogitlink",
"v2" = "multilogitlink",
"v3" = "loglink",
"v4" = "logofflink",
"(Default)" = "identitylink",
"v5" = "multilogitlink"),
earg.list = list("(Intercept)" = list(),
"v2" = list(),
"v4" = list(offset = myoffset),
"v3" = list(),
"(Default)" = list(),
"v5" = list()),
zero = c(1:2, 6)),
data = ndata, trace = TRUE)
coef(fit1, matrix = TRUE)
summary(fit1)
# This works only for intercept-only models:
multilogitlink(rbind(coef(fit1, matrix = TRUE)[1, c(1, 2)]), inverse = TRUE)
# A model with covariate x2 for the regression coefficients
fit2 <- vglm(y2 ~ 1 + x2,
form2 = ~ 1 + v2 + v3 + v4 + v5 + v6,
normal.vcm(link.list = list("(Intercept)" = "multilogitlink",
"v2" = "multilogitlink",
"v3" = "logitlink",
"v4" = "logloglink",
"(Default)" = "identitylink",
"v5" = "multilogitlink"),
earg.list = list("(Intercept)" = list(),
"v2" = list(),
"v3" = list(),
"v4" = list(),
"(Default)" = list(),
"v5" = list()),
zero = c(1:2, 6)),
data = ndata, trace = TRUE)
coef(fit2, matrix = TRUE)
summary(fit2)
Run the code above in your browser using DataLab