-
Notifications
You must be signed in to change notification settings - Fork 0
Description
R code to test against:
Set up data as per https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_sect016.htm
and https://www.statsmodels.org/dev/examples/notebooks/generated/quasibinomial.html
raw <- c(0.05,0.00,1.25,2.50,5.50,1.00,5.00,5.00,17.50,
0.00,0.05,1.25,0.50,1.00,5.00,0.10,10.00,25.00,
0.00,0.05,2.50,0.01,6.00,5.00,5.00,5.00,42.50,
0.10,0.30,16.60,3.00,1.10,5.00,5.00,5.00,50.00,
0.25,0.75,2.50,2.50,2.50,5.00,50.00,25.00,37.50,
0.05,0.30,2.50,0.01,8.00,5.00,10.00,75.00,95.00,
0.50,3.00,0.00,25.00,16.50,10.00,50.00,50.00,62.50,
1.30,7.50,20.00,55.00,29.50,5.00,25.00,75.00,95.00,
1.50,1.00,37.50,5.00,20.00,50.00,50.00,75.00,95.00,
1.50,12.70,26.25,40.00,43.50,75.00,75.00,75.00,95.00)
df <- data.frame(blotch = raw / 100,
variety = as.factor(rep(1:9, times = length(raw) / length(1:9))),
site = as.factor(unlist(lapply(1:10, rep, times = length(raw) / length(1:10)))))
mod <- glm(blotch ~ variety + site, family = quasibinomial, data = df)
summary(mod)