I’ve asked a couple of colleagues who are experienced in r, for help in converting SAS code to r. The challenge involves a multiway contingency table analysis, using a relatively simple data set. Although both colleague are veteran r users, they were independently, unable to generate output in r that even closely matched output from SAS — leaving them both perplexed. Below are the Data and SAS output and the SAS program. If anyone is up for a challenge in solving the r script for this…it would be much appreciated.
data new;
input prey $ stage $ habitat $ number;
cards ;
chr a l 69
chr j p 8
chr j l 43
chr a p 21
amp a l 29
amp j p 1
amp j l 19
amp a p 6
zoo a l 6
zoo j p 28
zoo j l 3
zoo a p 72
ost a l 10
ost j p 4
ost j l 6
ost a p 4
;
proc genmod data=new;
class prey stage habitat/param=effect;
model number=prey|stage|habitat / dist=poisson type3 wald;
output out=out p=p;
run;
proc sort data=out;
by prey stage habitat;
run;
proc print data=out;
run;
As mentioned previously, I’d asked two colleagues who are well-versed in r, to help solve a SAS-to-r conversion for a multiway contingency table analysis. Neither colleague was able to generate results using r that were even closely similar to SAS. I also conferred with a Senior Statistician at SAS, who corroborated my [SAS] results with his own analysis, where he used a number of different procedures. So the SAS output is correct.
Steve is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
8
Here’s an identical result to the SAS code, making use of the joint_tests() function from the R package emmeans (thanks as always to Russ Lenth for creating the package, it’s a lifesaver)! By the way, I was able to determine that SAS was doing joint tests by running the OP’s code in SAS which gives the following message:
Note: Under full-rank parameterizations, Type 3 effect tests are replaced by joint tests. The joint test for an effect is a test that all the parameters associated with that effect are zero. Such joint tests might not be equivalent to Type 3 effect tests under GLM parameterization.
library(emmeans)
dt <- data.frame(
prey = rep(c('chr','amp','zoo','ost'), each = 4),
stage = c('a','j','j','a'),
habitat = c('l', 'p'),
number = c(69, 8, 43, 21, 29, 1, 19, 6, 6, 28, 3, 72, 10, 4, 6, 4)
)
fit <- glm(number ~ prey * stage * habitat, family = poisson, data = dt)
joint_tests(fit)
We get this output:
model term df1 df2 F.ratio Chisq p.value
prey 3 Inf 16.844 50.532 <.0001
stage 1 Inf 12.103 12.103 0.0005
habitat 1 Inf 5.739 5.739 0.0166
prey:stage 3 Inf 0.556 1.668 0.6439
prey:habitat 3 Inf 28.608 85.824 <.0001
stage:habitat 1 Inf 0.923 0.923 0.3368
prey:stage:habitat 3 Inf 0.635 1.905 0.5925
SAS output screenshot for comparison:
3