Difference in output between SAS's proc genmod and R's glm

I'm trying to replicate a model fitted in SAS in R but the fit I'm getting gives me slightly different coefficients and standard errors.

Data:

testdata <- data.frame(matrix(c("f","Test", 1.75,   16, 0,  16, 0,  1,  1,
"m",    "Test", 1.75,   15, 1,  16, 6.25,   1,  0,
"f",    "Test", 2.75,   4,  12, 16, 75, 1,  1,
"m",    "Test", 2.75,   9,  6,  15, 40, 1,  0,
"f",    "WHO",  1.75,   15, 1,  16, 6.25,   0,  1,
"m",    "WHO",  1.75,   14, 2,  16, 12.5,   0,  0,
"f",    "WHO",  2.75,   2,  13, 15, 86.6667,    0,  1,
"m",    "WHO",  2.75,   3,  13, 16, 81.25,  0,  0
), ncol=9, byrow=TRUE))
names(testdata) <- c("sex", "vaccine", "dose", "not_p", "para", "n", "pct",
"vacnum", "sexno")


SAS:

proc genmod data=model_data;
class sex;
model para/n  = dose sex vacnum

/dist=bin
type3;
run;

Analysis Of Maximum Likelihood Parameter Estimates
Parameter   DF Estimate Std Error Wald 95% Conf Lim Wald Chi-Square Pr > ChiSq
Intercept   1 -9.4020   1.6220    -12.5810  -6.2230  33.60           <.0001
dose        1  3.9208   0.6460      2.6546   5.1870  36.83           <.0001
sex f       1  0.5574   0.5184     -0.4587   1.5735   1.16           0.2823
sex m       0  0.0000   0.0000      0.0000   0.0000    .              .
vacnum      1 -1.3221   0.5483     -2.3967  -0.2475   5.81           0.0159
Scale       0  1.0000   0.0000      1.0000   1.0000


R:

testdata$sexno <- as.factor(testdata$sexno)
a <- contr.treatment(2, base = 1, contrasts = TRUE)

contrasts(testdata\$sexno) <- a

fitreduced <- glm(para/n ~ dose + as.factor(sex) + vacnum,

coef(summary(fitreduced))

Estimate Std. Error   t value    Pr(>|t|)
(Intercept)     -9.4013750  1.7613982 -5.337450 0.005935450
dose             3.9173794  0.7001133  5.595351 0.005007179
as.factor(sex)1  0.5704671  0.5568436  1.024466 0.363525300
vacnum          -1.3336100  0.5887552 -2.265135 0.086189704


I believe I have the right contrasts to give me a type III SS but there is a small discrepency in values, have a missed something here?

• The GLM process is iterative and dependent on random numbers. If you could set the same seed in each, I believe, you would get the exact same results. – Eric Jul 24 '15 at 13:42
• @Eric Calculation of estimates in glm is iterative, but none of the usual algorithms is random. – Glen_b Jul 25 '15 at 2:44

I notice several things here.

First, when you enter your data via matrix, all the data have to be the same type. Thus, they are coerced to be the most inclusive type, strings, which in turn are coerced to be factors by default. Note:

testdata <- data.frame(matrix(c("f","Test", 1.75,   16, 0,  16, 0,  1,  1,
...
sapply(testdata, class)
#      sex  vaccine     dose    not_p     para        n      pct   vacnum    sexno
# "factor" "factor" "factor" "factor" "factor" "factor" "factor" "factor" "factor"


Try using read.table(text='...', sep=",") instead:

testdata <- read.table(text='"f", "Test", 1.75,   16,   0,  16,  0,      1,  1
"m", "Test", 1.75,   15,   1,  16,  6.25,   1,  0
"f", "Test", 2.75,    4,  12,  16, 75,      1,  1
"m", "Test", 2.75,    9,   6,  15, 40,      1,  0
"f", "WHO",  1.75,   15,   1,  16,  6.25,   0,  1
"m", "WHO",  1.75,   14,   2,  16, 12.5,    0,  0
"f", "WHO",  2.75,    2,  13,  15, 86.6667, 0,  1
"m", "WHO",  2.75,    3,  13,  16, 81.25,   0,  0', sep=",")
names(testdata) <- c("sex", "vaccine", "dose", "not_p", "para", "n", "pct",
"vacnum", "sexno")
sapply(testdata, class)
#      sex   vaccine      dose     not_p      para         n       pct    vacnum
# "factor"  "factor" "numeric" "integer" "integer" "integer" "numeric" "integer"
#     sexno
# "integer"


(That was small potatoes.) The next trap to worry about is that SAS and R code logistic regression for binomial data differently. SAS uses "events over trials", but R uses the odds, successes/failures. Thus, your model formula should be:

form <- as.formula("cbind(para, n-para) ~ dose + sex + vacnum")


Finally, you specified family=quasibinomial (i.e., the quasibinomial) in your R code, but \DIST=BIN (i.e., the binomial) in your SAS code. To match the SAS output, use the binomial instead. Thus, your final model is:

fitreduced <- glm(form, family=binomial(link="logit"), data=testdata)
coef(summary(fitreduced))
#               Estimate Std. Error   z value     Pr(>|z|)
# (Intercept) -9.4020028  1.6219570 -5.796703 6.763131e-09
# dose         3.9207805  0.6460193  6.069138 1.285986e-09
# sexf         0.5574087  0.5184112  1.075225 2.822741e-01
# vacnum      -1.3221011  0.5482645 -2.411430 1.589012e-02


This seems to match the SAS estimates and standard errors.

• That's great. Apologies for the data input not being correct, i actually inputted it through .csv myself, i just wrote that to make a reproducible example but forget to check the classes so thanks. The SAS events/trials and R using success/failures is very informative. I never suspect the notation would be difference. Thank-you very much. – Kabau Jul 24 '15 at 14:16