[R] Syntax for fit.contrast (from package gmodels)

[Martin Maechler]

>>>>> Sorkin, John <jsorkin at som.umaryland.edu>
>>>>>     on Sun, 22 Oct 2017 22:56:16 +0000 writes:

    > David,
    > Thank you for responding to my post.


    > Please consider the following output (typeregional is a factor having two levels, "regional" vs. "general"):

    > Call:
    > glm(formula = events ~ type, family = poisson(link = log), data = data,
    > offset = log(SS))

    > Deviance Residuals:
    > Min       1Q   Median       3Q      Max
    > -43.606  -17.295   -4.651    4.204   38.421

    > Coefficients:
    > Estimate Std. Error z value Pr(>|z|)
    > (Intercept)  -2.52830    0.01085 -233.13   <2e-16 ***
    > typeregional  0.33788    0.01641   20.59   <2e-16 ***


    > Let's forget for a moment that the model is a Poisson regression and pretend that the output is from a simple linear regression, e.g. from lm.


    > To get the estimate for "general" one simply needs to use the value of the intercept i.e. -2.5830. Similarly to get the 95% CI of general one simply needs to compute -2.52830-(1.96*0.01085) and -2.52830+(1.96*0.01085).

I'm pretty sure you can just use  (the base R) functions

  dummy.coef()

or
  model.tables()

possibly with SE=TRUE to get coefficients for all levels of a factor..
I'd like to have tried to show this here, but for that we'd have
wanted to see a "MRE" or "ReprEx" (minimal reproducible example) ..

    > To get the estimate for "regional" one needs to compute intercept + typeregional, i.e. -2.52830 + 0.33788. To get the 95% CI is somewhat more difficult as one needs to use results from the variance-covariance matix, specifically the variance of intercept, the variance of "regional", and the covariance of (intercept,"regional") which involves:

    > var =  var(intercept) + var(regional) +2*(covar(intercept,regional)),

    > and then get the SE of the variance

    > SE=sqrt(var)

    > 95% CI = intercept + regional - 1.95*SE and intercept + regional + 1.95*SE.


    > I was hoping that a contrast statement could be written that would give me the point estimate and SE for "general" and its SE and another contrast statement could be written that would give me the point estimate and SE for "general" and it SE without my having to work directly with the variance-covariance matrix. I tried doing this using the fit.contrast statements (from the gmodels package):


    > fit.contrast(model,type,c(1,0),showall=TRUE)

    > fit.contrast(model,type,c(0,1),showall=TRUE)


    > and received the error message,

    > Error in `[[<-`(`*tmp*`, varname, value = c(0, 1)) :
    > no such index at level 1


    > Perhaps fit.contrast is not the way to accomplish my goal. Perhaps my goal can be accomplished without a contrast statement, but I don't know how.

My guess is that "standard R" aka "base R" would be
sufficient to get what you'd want, notably if you'd consider
using  se.contrast() additionally.

Martin


    > Thank you,
    > John