6
$\begingroup$

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 
  link=logit
  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, 
                  family=quasibinomial(link="logit"), data=testdata)

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Eric Jul 24 '15 at 13:42
  • 1
    $\begingroup$ @Eric Calculation of estimates in glm is iterative, but none of the usual algorithms is random. $\endgroup$ – Glen_b -Reinstate Monica Jul 25 '15 at 2:44
13
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Kabau Jul 24 '15 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.