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