Syntax for fit.contrast

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Syntax for fit.contrast

Sorkin, John
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 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)

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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast

David Winsemius

> On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

Sorkin, John
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 <[hidden email]>
Sent: Sunday, October 22, 2017 1:20 PM
To: Sorkin, John
Cc: [hidden email]
Subject: Re: [R] Syntax for fit.contrast


> On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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]]


______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

David Winsemius

> On Oct 22, 2017, at 3:56 PM, Sorkin, John <[hidden email]> wrote:
>
> 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):

I'm guessing that the second contrast you were hoping for was actually for "regional".

Contrasts, hence the name, are for differences between two levels (or more accurately between the means on the scale specified by the link parameter. In the absence of another level the only other reference point would be a value of zero or perhaps the value you specified by your offset term.

--
David


>
> 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 <[hidden email]>
> Sent: Sunday, October 22, 2017 1:20 PM
> To: Sorkin, John
> Cc: [hidden email]
> Subject: Re: [R] Syntax for fit.contrast
>  
>
> > On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
> >
> > ______________________________________________
> > [hidden email] 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

David Winsemius
Alameda, CA, USA

'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

Sorkin, John
David,

Again  you have my thanks!.

You are correct. What I want is not technically a contrast. What I want is the estimate for "regional" and its SE. I don't mind if I get these on the log scale; I can get the anti-log. Can you suggest  how I can get the point estimate and its SE for "regional"? The predict function will give the point estimate, but not (to my knowledge) the SE.

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 <[hidden email]>
Sent: Sunday, October 22, 2017 7:56 PM
To: Sorkin, John
Cc: [hidden email]
Subject: Re: [R] Syntax for fit.contrast (from package gmodels)


> On Oct 22, 2017, at 3:56 PM, Sorkin, John <[hidden email]> wrote:
>
> 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):
I'm guessing that the second contrast you were hoping for was actually for "regional".

Contrasts, hence the name, are for differences between two levels (or more accurately between the means on the scale specified by the link parameter. In the absence of another level the only other reference point would be a value of zero or perhaps the value you specified by your offset term.

--
David


>
> 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 <[hidden email]>
> Sent: Sunday, October 22, 2017 1:20 PM
> To: Sorkin, John
> Cc: [hidden email]
> Subject: Re: [R] Syntax for fit.contrast
>
>
> > On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
> >
> > ______________________________________________
> > [hidden email] 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
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]]


______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

David Winsemius

> On Oct 22, 2017, at 5:01 PM, Sorkin, John <[hidden email]> wrote:
>
> David,
> Again  you have my thanks!.
> You are correct. What I want is not technically a contrast. What I want is the estimate for "regional" and its SE.

There needs to be a reference value for the contrast. Contrasts are differences. I gave you the choice of two references (treatment constrast or the offset value you specified). Pick one or suggest an alternate value. Typical alternate values are the global mean or zero.

Read ?predict.glm

"se.fit    logical switch indicating if standard errors are required."


> I don't mind if I get these on the log scale; I can get the anti-log. Can you suggest  how I can get the point estimate and its SE for "regional"? The predict function will give the point estimate, but not (to my knowledge) the SE.


> 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 <[hidden email]>
> Sent: Sunday, October 22, 2017 7:56 PM
> To: Sorkin, John
> Cc: [hidden email]
> Subject: Re: [R] Syntax for fit.contrast (from package gmodels)
>  
>
> > On Oct 22, 2017, at 3:56 PM, Sorkin, John <[hidden email]> wrote:
> >
> > 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):
>
> I'm guessing that the second contrast you were hoping for was actually for "regional".
>
> Contrasts, hence the name, are for differences between two levels (or more accurately between the means on the scale specified by the link parameter. In the absence of another level the only other reference point would be a value of zero or perhaps the value you specified by your offset term.
>
> --
> David
>
>
> >
> > 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 <[hidden email]>
> > Sent: Sunday, October 22, 2017 1:20 PM
> > To: Sorkin, John
> > Cc: [hidden email]
> > Subject: Re: [R] Syntax for fit.contrast
> >  
> >
> > > On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
> > >
> > > ______________________________________________
> > > [hidden email] 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
>
> David Winsemius
> Alameda, CA, USA
>
> 'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law

David Winsemius
Alameda, CA, USA

'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

Sorkin, John
David,

predict.glm and se.fit were exactly what I was looking for.

Many  thanks!

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 <[hidden email]>
Sent: Sunday, October 22, 2017 8:15 PM
To: Sorkin, John
Cc: [hidden email]
Subject: Re: [R] Syntax for fit.contrast (from package gmodels)


> On Oct 22, 2017, at 5:01 PM, Sorkin, John <[hidden email]> wrote:
>
> David,
> Again  you have my thanks!.
> You are correct. What I want is not technically a contrast. What I want is the estimate for "regional" and its SE.

There needs to be a reference value for the contrast. Contrasts are differences. I gave you the choice of two references (treatment constrast or the offset value you specified). Pick one or suggest an alternate value. Typical alternate values are the global mean or zero.

Read ?predict.glm

"se.fit     logical switch indicating if standard errors are required."


> I don't mind if I get these on the log scale; I can get the anti-log. Can you suggest  how I can get the point estimate and its SE for "regional"? The predict function will give the point estimate, but not (to my knowledge) the SE.


> 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 <[hidden email]>
> Sent: Sunday, October 22, 2017 7:56 PM
> To: Sorkin, John
> Cc: [hidden email]
> Subject: Re: [R] Syntax for fit.contrast (from package gmodels)
>
>
> > On Oct 22, 2017, at 3:56 PM, Sorkin, John <[hidden email]> wrote:
> >
> > 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):
>
> I'm guessing that the second contrast you were hoping for was actually for "regional".
>
> Contrasts, hence the name, are for differences between two levels (or more accurately between the means on the scale specified by the link parameter. In the absence of another level the only other reference point would be a value of zero or perhaps the value you specified by your offset term.
>
> --
> David
>
>
> >
> > 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 <[hidden email]>
> > Sent: Sunday, October 22, 2017 1:20 PM
> > To: Sorkin, John
> > Cc: [hidden email]
> > Subject: Re: [R] Syntax for fit.contrast
> >
> >
> > > On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
> > >
> > > ______________________________________________
> > > [hidden email] 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
>
> David Winsemius
> Alameda, CA, USA
>
> 'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law
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]]


______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

David Winsemius

> On Oct 22, 2017, at 5:26 PM, Sorkin, John <[hidden email]> wrote:
>
> David,
> predict.glm and se.fit were exactly what I was looking for.

The default 'se' delivered for a listed contrast is not for that particular level per se, but rather for the difference between that level and the non-listed factor level (its mean) which is "embedded" as it were in the Intercept. The se for the Intercept is for the mean of all the reference levels jointly being zero.

Apologies for any mis-steps in this effort. I'm enjoying my fourth IPA. I usually do not drink and derive.

--
David.


> Many  thanks!
> 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 <[hidden email]>
> Sent: Sunday, October 22, 2017 8:15 PM
> To: Sorkin, John
> Cc: [hidden email]
> Subject: Re: [R] Syntax for fit.contrast (from package gmodels)
>  
>
> > On Oct 22, 2017, at 5:01 PM, Sorkin, John <[hidden email]> wrote:
> >
> > David,
> > Again  you have my thanks!.
> > You are correct. What I want is not technically a contrast. What I want is the estimate for "regional" and its SE.
>
> There needs to be a reference value for the contrast. Contrasts are differences. I gave you the choice of two references (treatment constrast or the offset value you specified). Pick one or suggest an alternate value. Typical alternate values are the global mean or zero.
>
> Read ?predict.glm
>
> "se.fit     logical switch indicating if standard errors are required."
>
>
> > I don't mind if I get these on the log scale; I can get the anti-log. Can you suggest  how I can get the point estimate and its SE for "regional"? The predict function will give the point estimate, but not (to my knowledge) the SE.
>
>
> > 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 <[hidden email]>
> > Sent: Sunday, October 22, 2017 7:56 PM
> > To: Sorkin, John
> > Cc: [hidden email]
> > Subject: Re: [R] Syntax for fit.contrast (from package gmodels)
> >  
> >
> > > On Oct 22, 2017, at 3:56 PM, Sorkin, John <[hidden email]> wrote:
> > >
> > > 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):
> >
> > I'm guessing that the second contrast you were hoping for was actually for "regional".
> >
> > Contrasts, hence the name, are for differences between two levels (or more accurately between the means on the scale specified by the link parameter. In the absence of another level the only other reference point would be a value of zero or perhaps the value you specified by your offset term.
> >
> > --
> > David
> >
> >
> > >
> > > 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 <[hidden email]>
> > > Sent: Sunday, October 22, 2017 1:20 PM
> > > To: Sorkin, John
> > > Cc: [hidden email]
> > > Subject: Re: [R] Syntax for fit.contrast
> > >  
> > >
> > > > On Oct 22, 2017, at 6:04 AM, Sorkin, John <[hidden email]> 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]]
> > > >
> > > > ______________________________________________
> > > > [hidden email] 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
> >
> > David Winsemius
> > Alameda, CA, USA
> >
> > 'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law
>
> David Winsemius
> Alameda, CA, USA
>
> 'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law

David Winsemius
Alameda, CA, USA

'Any technology distinguishable from magic is insufficiently advanced.'   -Gehm's Corollary to Clarke's Third Law

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Syntax for fit.contrast (from package gmodels)

Martin Maechler
In reply to this post by Sorkin, John
>>>>> Sorkin, John <[hidden email]>
>>>>>     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

______________________________________________
[hidden email] 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.