thr3ads.net - R help - [R] Poisson Regression [Oct 2010]

If this information is useful, please help other people find it:
Share via:

Antonio Paredes

2010-Oct-13 20:50 UTC

[R] Poisson Regression

Hello everyone,

I wanted to ask if there is an R-package to fit the following Poisson
regression model

log(\lambda_{ijk}) = \phi_{i} + \alpha_{j} + \beta_{k}
i=1,\cdots,N (subjects)
j=0,1 (two levels)
k=0,1 (two levels)

treating the \phi_{i} as nuinsance parameters.

Thank you very much
-- 
-Tony

	[[alternative HTML version deleted]]

David Winsemius

2010-Oct-14 00:21 UTC

head link

[R] Poisson Regression

On Oct 13, 2010, at 4:50 PM, Antonio Paredes wrote:
> Hello everyone,
>
> I wanted to ask if there is an R-package to fit the following Poisson
> regression model
>
> log(\lambda_{ijk}) = \phi_{i} + \alpha_{j} + \beta_{k}
> i=1,\cdots,N (subjects)
> j=0,1 (two levels)
> k=0,1 (two levels)
>
> treating the \phi_{i} as nuinsance parameters.
If I am reading this piece correctly there should be no difference  
between a conditional treatment of phi_i in  that model and results  
from the unconditional model one would get from fitting with

glm(lambda ~ phi + alpha + beta  ,family="poisson").

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.6.9679&rep=rep1&type=pdf

(But I am always looking for corrections to my errors.)

-- 
David Winsemius, MD
West Hartford, CT

Bill.Venables at csiro.au

2010-Oct-14 00:44 UTC

head link

[R] Poisson Regression

One possible way to treat parameters as "nuisance parameters" is to
model them as random.  This gives allows them to have a reduced parametric load.

There are many packages with funcitons to fit glmms.  One you may wish to look
at is lme4, which has the lmer fitting function

library(lme4)
fm <- glmer(Y ~ A + B + (1|Subject), family = poisson, data = pData)

for example, may be a useful alternative to a fully fixed effects approach.

W. 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of David Winsemius
Sent: Thursday, 14 October 2010 10:22 AM
To: Antonio Paredes
Cc: r-help at r-project.org
Subject: Re: [R] Poisson Regression

On Oct 13, 2010, at 4:50 PM, Antonio Paredes wrote:
> Hello everyone,
>
> I wanted to ask if there is an R-package to fit the following Poisson
> regression model
>
> log(\lambda_{ijk}) = \phi_{i} + \alpha_{j} + \beta_{k}
> i=1,\cdots,N (subjects)
> j=0,1 (two levels)
> k=0,1 (two levels)
>
> treating the \phi_{i} as nuinsance parameters.
If I am reading this piece correctly there should be no difference  
between a conditional treatment of phi_i in  that model and results  
from the unconditional model one would get from fitting with

glm(lambda ~ phi + alpha + beta  ,family="poisson").

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.6.9679&rep=rep1&type=pdf

(But I am always looking for corrections to my errors.)

-- 
David Winsemius, MD
West Hartford, CT

______________________________________________
R-help at r-project.org mailing list
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.

Charles C. Berry

2010-Oct-14 00:51 UTC

head link

[R] Poisson Regression

On Wed, 13 Oct 2010, David Winsemius wrote:
>
> On Oct 13, 2010, at 4:50 PM, Antonio Paredes wrote:
>
>> Hello everyone,
>> 
>> I wanted to ask if there is an R-package to fit the following Poisson
>> regression model
>> 
>> log(\lambda_{ijk}) = \phi_{i} + \alpha_{j} + \beta_{k}
>> i=1,\cdots,N (subjects)
>> j=0,1 (two levels)
>> k=0,1 (two levels)
>> 
>> treating the \phi_{i} as nuinsance parameters.
>
> If I am reading this piece correctly there should be no difference between
a
> conditional treatment of phi_i in  that model and results from the 
> unconditional model one would get from fitting with
>
> glm(lambda ~ phi + alpha + beta  ,family="poisson").
Right.

But if N is large, the model.matrix will be huge and there may be problems 
with memory and elapsed time.

loglin() and loglm() will fit the same model without need for a 
model.matrix (modulo having enough data to actually fit that model), and 
large values of N are no big deal.

HTH,

Chuck>
>
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.6.9679&rep=rep1&type=pdf
>
> (But I am always looking for corrections to my errors.)
>
> -- 
> David Winsemius, MD
> West Hartford, CT
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
>
Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

Michael Friendly

2010-Oct-14 12:55 UTC

head link

[R] Poisson Regression

On 10/13/2010 4:50 PM, Antonio Paredes wrote:> Hello everyone,
>
> I wanted to ask if there is an R-package to fit the following Poisson
> regression model
>
> log(\lambda_{ijk}) = \phi_{i} + \alpha_{j} + \beta_{k}
> i=1,\cdots,N (subjects)
> j=0,1 (two levels)
> k=0,1 (two levels)
>
> treating the \phi_{i} as nuinsance parameters.
>
> Thank you very much
You can use the gnm() function in the gnm package with eliminate= to get 
rid of parameters for each subject.
Something like

gnm( lambda ~ Row + Col, eliminate = Subject, family=poisson, data=myData)

where Row, Col, Subject are suitable *factors*.  This is equivalent to

gnm( lambda ~ -1 + Subject Row + Col, family=poisson)
except that Subject parameters aren't estimated explicitly.  See 
vignette("gnmOverview")


-- 
Michael Friendly     Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    Web:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA

avsha38

2010-Nov-10 07:43 UTC

head link

[R] poisson regression

Hello All,

I have a question regarding using Poission Regression, I would like to Model
the number of hospitalizations by a set of covariates. 

The issue I ran into is "lack of fit" even after I tried to solve the
"overdispersion" problem with negetive binomial. 

What would you suggest?
1. Is the Poission Regression the right Model?
2. is there another Model that is better for my needs?

Any ideas will be more the welcome.
Thank you,
Avi

The File is attached and the code is below:
Please note that the "offset" is for dealingwith different followup
times of
the subjects.

pr2<-read.csv("c:/rfiles/pr/pr_na.csv",header=TRUE)
head(pr2)
# distribution of the Response variable "zb44" , number of
hospitalization.
tab <- table(pr2$zb44)
tab
x <- as.numeric(names(tab))
plot(x, tab, type='h', xlab='Number of hospitalizations',
ylab='Frequency',
main="Distribution of hospitalization")
points(x, tab, pch=16)
mean(pr2$zb44)
# find the "best model" according to AIC
fit.pr0<- glm(zb44 ~1+ offset(log(timem)),data=pr2, family=poisson) 
fit.full<- glm(zb44 ~
ESHKOL+factor(Phys_activity)+factor(PTCAwin45d)+factor(White_collar)+factor
(Asia_Africa)+factor(Sex)+Age+ factor(SteadyPartner)
         + factor(RelIncome)+ factor(Employment)+
Education+factor(CA_blockers)+Health_Prob+factor(PerceivHealth)+factor(Diuretics)+
factor(ACE)
         + factor(DiabTreat)+ factor(Insulin)+ factor(LLD)+ factor(Aspirin)+
factor(Smoking)+factor(CHD_duration)+ factor(HTN)+ factor(Diabetes)
         + factor(Hyperchol)+ factor(Intens_care)+ factor(Thromb)+
factor(AP)+ factor(QWave)+factor(AntMI)+hosp_days+factor(COPD)+
factor(B_blockers)
         + factor(Cancer)+factor(Ulcer)+factor(kilip2)+ factor(HF_index)+
factor(CVA)+Charlson_cat
         + factor(Cardiac_death)+
factor(CABGwin45d)+offset(log(timem)),data=pr2,family=poisson) 
Forw <- step(fit.pr0, scope = list(lower =fit.pr0, upperfit.full),direction =
"both")
summary(Forw)
anova.glm(Forw ,test = "Chisq") # display significate factors
# checking model fit
1-pchisq(Forw$dev, Forw$df.residual)
# p value is nearly zero, There is a lack of fit!!!

# Residuals plots
dev.new()
par(mfrow=c(2,2))
plot(Forw , pch=23, bg='red', cex=2)

# Examine (over) dispersion
# Pearson?s residuals
(dp <- sum(resid(Forw , type = "pearson")^2)/Forw$df.resid)
summary(Forw,dispersion=dp)
# Deviance
(dp1<-Forw$deviance/Forw$df.resid)
summary(Forw,dispersion=dp1)
# Model is over dispersed

# Consider a negative binomial  to deal with overdispersion

NB.model<- step(glm.nb(hosp_total~
ESHKOL+factor(Phys_activity)+factor(PTCAwin45d)+factor(White_collar)+factor(Asia_Africa)+factor(Sex)+Age+
factor(SteadyPartner)
         + factor(RelIncome)+ factor(Employment)+
Education+factor(CA_blockers)+Health_Prob+factor(PerceivHealth)+factor(Diuretics)+
factor(ACE)
         + factor(DiabTreat)+ factor(Insulin)+ factor(LLD)+ factor(Aspirin)+
factor(Smoking)+factor(CHD_duration)+ factor(HTN)+ factor(Diabetes)
         + factor(Hyperchol)+ factor(Intens_care)+ factor(Thromb)+
factor(AP)+ factor(QWave)+factor(AntMI)+hosp_days+factor(COPD)+
factor(B_blockers)
         + factor(Cancer)+factor(Ulcer)+factor(kilip2)+ factor(HF_index)+
factor(CVA)+Charlson_cat
         + factor(Cardiac_death)+
factor(CABGwin45d)+offset(log(timem)),data=pr2)) 
summary(NB.model)
1 - pchisq(NB.model$dev, NB.model$df.residual)
par(mfrow=c(2,2))
plot(NB.model, pch=23, bg='red', cex=2)
# p value is nearly zero, There is still a lack of fit!!!
http://r.789695.n4.nabble.com/file/n3035613/pr_na.csv pr_na.csv 
-- 
View this message in context:
http://r.789695.n4.nabble.com/poisson-regression-tp3035613p3035613.html
Sent from the R help mailing list archive at Nabble.com.

Reasonably Related Threads

Search for more maybe matching threads

R help - Oct 2010 - Poisson Regression

[R] Poisson Regression

[R] Poisson Regression

[R] Poisson Regression

[R] Poisson Regression

[R] Poisson Regression

[R] poisson regression

Reasonably Related Threads