size <- 10 # Number of trials; N in the notation above
nn <- 200
odata <- data.frame(pstr1 = logitlink( 0, inv = TRUE), # 0.50
mubin1 = logitlink(-1, inv = TRUE), # Binomial mean
svec = rep(size, length = nn),
x2 = runif(nn))
odata <- transform(odata,
mubin2 = logitlink(-1 + x2, inv = TRUE))
odata <- transform(odata,
y1 = roiposbinom(nn, svec, pr = mubin1, pstr1 = pstr1),
y2 = roiposbinom(nn, svec, pr = mubin2, pstr1 = pstr1))
with(odata, table(y1))
fit1 <- vglm(y1 / svec ~ 1, oiposbinomial, data = odata,
weights = svec, trace = TRUE, crit = "coef")
fit2 <- vglm(y2 / svec ~ x2, oiposbinomial, data = odata,
weights = svec, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1) # Useful for intercept-only models
head(fitted(fit1, type = "pobs1")) # Estimate of P(Y = 1)
head(fitted(fit1))
with(odata, mean(y1)) # Compare this with fitted(fit1)
summary(fit1)
Run the code above in your browser using DataLab