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).
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.
Thank you,
John
John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric
Medicine
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
________________________________
From: David Winsemius <[email protected]>
Sent: Sunday, October 22, 2017 1:20 PM
To: Sorkin, John
Cc: [email protected]
Subject: Re: [R] Syntax for fit.contrast
> On Oct 22, 2017, at 6:04 AM, Sorkin, John <[email protected]> wrote:
>
> I have a model (run with glm) that has a factor, type. Type has two levels,
> "general" and "regional". I am trying to get estimates (and SEs) for the
> model with type="general" and type ="regional" using fit.contrast
?fit.contrast
No documentation for �fit.contrast� in specified packages and libraries:
you could try �??fit.contrast�
Perhaps the gmodels function of that name?
> but I can't get the syntax of the coefficients to use in fit.contrast
> correct. I hope someone can show me how to use fit.contrast, or some other
> method to get estimate with SEs. (I know I can use the variance co-variance
> matrix, but I would rather not have to code the linear contrast my self from
> the coefficients of the matrix)
>
I'm having trouble understanding what you are trying to extract. There are only
2 levels so there is really only one interesting contrast ("general" vs
"regional") , and it's magnitude would be reported by just typing `model`, and
it's SE would show up in output of `summary(model)`.
I'm thinking you should pick one of the examples in gmodels::fit.contrast that
most resembles your real problem, post it, and and then explain what
difficulties you are having with interpretation.
--
David.
> Thank you,
>
> John
>
>
> My model:
>
> model=glm(events~type,family=poisson(link=log),offset=log(SS),data=data)
>
>
> Model details:
>
>> summary(data$type)
>
> general regional
> 16 16
>
>> levels(data$type)
> [1] "general" "regional"
>
>> contrasts(data$type)
> regional
> general 0
> regional 1
>
>
> I have tried the following syntax for fit.contrast
>
> fit.contrast(model,type,c(1,0))
> and get an error:
> Error in `[[<-`(`*tmp*`, varname, value = cmat) :
> no such index at level 1
>
>
>> fit.contrast(model,type,c(0,1),showall=TRUE)
> and get an error:
> Error in `[[<-`(`*tmp*`, varname, value = cmat) :
> no such index at level 1
>
>
>
>> fit.contrast(model,type,c(1,-1),showall=TRUE)
> and get an error:
> Error in `[[<-`(`*tmp*`, varname, value = cmat) :
> no such index at level 1
>
>
>> fit.contrast(model,type,c(0))
> and get an error:
> Error in make.contrasts(coeff, ncol(coeff)) :
> Too many contrasts specified. Must be less than the number of factor levels
> (columns).
>
>> fit.contrast(model,type,c(1))
> Error in make.contrasts(coeff, ncol(coeff)) :
> and get an error
> Too many contrasts specified. Must be less than the number of factor levels
> (columns).
>
>
>
>
>
>
>
>
> John David Sorkin M.D., Ph.D.
> Professor of Medicine
> Chief, Biostatistics and Informatics
> University of Maryland School of Medicine Division of Gerontology and
> Geriatric Medicine
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [email protected] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius
Alameda, CA, USA
'Any technology distinguishable from magic is insufficiently advanced.'
-Gehm's Corollary to Clarke's Third Law
[[alternative HTML version deleted]]
______________________________________________
[email protected] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
