Closed hendersontrent closed 2 years ago
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)
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)