Hi all,
I was trying to fit a Gamma hierarchical model using "glmer", but got
weird error message that I could not understand. On the other hand, a similar
call to the glmmPQL leads to results that are close to what I expect. I also
tried to change tha "nAGQ" argument in "glmer", but it did
not solve the problem. The model I was fitting has a simple structure - one
hierarchy and one offset, and I did not expect "glmer" to fail on this
simple data set. Did I specify something wrong here? I appreciate any help.
Thanks in advance.
#############################################################
## Code to fit the model
#############################################################
library(MASS)
library(lme4)
address <-
"http://instruction.bus.wisc.edu/jfrees/jfreesbooks/Longitudinal%20and%20Panel%20Data/Book/Data/TXTData/WorkerC.txt"
wc <- read.table(address,header=T)
# clean data
cl <- unique(wc$CL[wc$LOSS==0])
wc2 <- subset(wc, !(CL %in% cl) )
# fit model
fit1 <- glmmPQL(LOSS~ YR+offset(log(PR)),
random=~1|CL,family=Gamma(link="log"),
data=wc2)
fit2 <- glmer(LOSS~ (1|CL) + YR+offset(log(PR)),
family=Gamma(link="log"),
data=wc2)
#############################################################
## Results of fit1 and fit2
#############################################################
> fit1
Linear mixed-effects model fit by maximum likelihood
Data: wc2
Log-likelihood: NA
Fixed: LOSS ~ YR + offset(log(PR))
(Intercept) YR
-4.2830507147 0.0005085944
Random effects:
Formula: ~1 | CL
(Intercept) Residual
StdDev: 0.8313193 0.5346455
Variance function:
Structure: fixed weights
Formula: ~invwt
Number of Observations: 700
Number of Groups: 100> fit2
Error in asMethod(object) : matrix is not symmetric [1,2]
Wayne (Yanwei) Zhang
Statistical Research
CNA
Email: Yanwei.Zhang@cna.com<mailto:Yanwei.Zhang@cna.com>
NOTICE: This e-mail message, including any attachments and appended messages,
is for the sole use of the intended recipients and may contain confidential and
legally privileged information.
If you are not the intended recipient, any review, dissemination, distribution,
copying, storage or other use of all or any portion of this message is strictly
prohibited.
If you received this message in error, please immediately notify the sender by
reply e-mail and delete this message in its entirety.
[[alternative HTML version deleted]]
Hi: I believe that the problem with the glmer() call is that you are invoking the Gamma family, which appears not to be implemented at present. See http://r-project.markmail.org/search/?q=r-sig-mixed-models#query:r-sig-mixed-models%20date%3A201105%20+page:8+mid:q2b2uoupejdzmvtl+state:results Dennis On Mon, May 16, 2011 at 3:31 PM, Zhang,Yanwei <Yanwei.Zhang at cna.com> wrote:> Hi all, > > I was trying to fit a Gamma hierarchical model using "glmer", but got weird error message that I could not understand. On the other hand, a similar call to the glmmPQL leads to results that are close to what I expect. I also tried to change tha "nAGQ" argument in "glmer", but it did not solve the problem. The model I was fitting has a simple structure - one hierarchy and one offset, and I did not expect "glmer" to fail on this simple data set. Did I specify something wrong here? I appreciate any help. Thanks in advance. > > ############################################################# > ## Code to fit the model > ############################################################# > library(MASS) > library(lme4) > address <- "http://instruction.bus.wisc.edu/jfrees/jfreesbooks/Longitudinal%20and%20Panel%20Data/Book/Data/TXTData/WorkerC.txt" > wc <- read.table(address,header=T) > # clean data > cl <- unique(wc$CL[wc$LOSS==0]) > wc2 <- subset(wc, !(CL %in% cl) ) > # fit model > fit1 <- glmmPQL(LOSS~ YR+offset(log(PR)), random=~1|CL,family=Gamma(link="log"), > ? ? ? ? ?data=wc2) > fit2 <- glmer(LOSS~ (1|CL) + YR+offset(log(PR)), family=Gamma(link="log"), > ? ? ? ? ?data=wc2) > > > ############################################################# > ## Results of fit1 and fit2 > ############################################################# > >> fit1 > Linear mixed-effects model fit by maximum likelihood > ?Data: wc2 > ?Log-likelihood: NA > ?Fixed: LOSS ~ YR + offset(log(PR)) > ?(Intercept) ? ? ? ? ? ?YR > -4.2830507147 ?0.0005085944 > > Random effects: > ?Formula: ~1 | CL > ? ? ? ?(Intercept) ?Residual > StdDev: ? 0.8313193 0.5346455 > > Variance function: > ?Structure: fixed weights > ?Formula: ~invwt > Number of Observations: 700 > Number of Groups: 100 >> fit2 > Error in asMethod(object) : matrix is not symmetric [1,2] > > > > > Wayne (Yanwei) Zhang > Statistical Research > CNA > Email: Yanwei.Zhang at cna.com<mailto:Yanwei.Zhang at cna.com> > > > > > NOTICE: ?This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. > If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. > If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety. > > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > 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. >
Hi Dennis, Thanks for your input. I am quite shocked that the Gamma model is not implemented in glmer yet, but do you happen to know if there is other packages in R that could fit Gamma hierarchical models besides MASS? The glmmPQL results are a bit off compared to what have been published. Thanks. Wayne (Yanwei) Zhang Statistical Research>CNA-----Original Message----- From: Dennis Murphy [mailto:djmuser at gmail.com] Sent: Tuesday, May 17, 2011 3:05 AM To: Zhang,Yanwei Cc: r-help at r-project.org Subject: Re: [R] Problem on glmer Hi: I believe that the problem with the glmer() call is that you are invoking the Gamma family, which appears not to be implemented at present. See http://r-project.markmail.org/search/?q=r-sig-mixed-models#query:r-sig-mixed-models%20date%3A201105%20+page:8+mid:q2b2uoupejdzmvtl+state:results Dennis On Mon, May 16, 2011 at 3:31 PM, Zhang,Yanwei <Yanwei.Zhang at cna.com> wrote:> Hi all, > > I was trying to fit a Gamma hierarchical model using "glmer", but got weird error message that I could not understand. On the other hand, a similar call to the glmmPQL leads to results that are close to what I expect. I also tried to change tha "nAGQ" argument in "glmer", but it did not solve the problem. The model I was fitting has a simple structure - one hierarchy and one offset, and I did not expect "glmer" to fail on this simple data set. Did I specify something wrong here? I appreciate any help. Thanks in advance. > > ############################################################# > ## Code to fit the model > ############################################################# > library(MASS) > library(lme4) > address <- "http://instruction.bus.wisc.edu/jfrees/jfreesbooks/Longitudinal%20and%20Panel%20Data/Book/Data/TXTData/WorkerC.txt" > wc <- read.table(address,header=T) > # clean data > cl <- unique(wc$CL[wc$LOSS==0]) > wc2 <- subset(wc, !(CL %in% cl) ) > # fit model > fit1 <- glmmPQL(LOSS~ YR+offset(log(PR)), random=~1|CL,family=Gamma(link="log"), > ? ? ? ? ?data=wc2) > fit2 <- glmer(LOSS~ (1|CL) + YR+offset(log(PR)), family=Gamma(link="log"), > ? ? ? ? ?data=wc2) > > > ############################################################# > ## Results of fit1 and fit2 > ############################################################# > >> fit1 > Linear mixed-effects model fit by maximum likelihood > ?Data: wc2 > ?Log-likelihood: NA > ?Fixed: LOSS ~ YR + offset(log(PR)) > ?(Intercept) ? ? ? ? ? ?YR > -4.2830507147 ?0.0005085944 > > Random effects: > ?Formula: ~1 | CL > ? ? ? ?(Intercept) ?Residual > StdDev: ? 0.8313193 0.5346455 > > Variance function: > ?Structure: fixed weights > ?Formula: ~invwt > Number of Observations: 700 > Number of Groups: 100 >> fit2 > Error in asMethod(object) : matrix is not symmetric [1,2] > > > > > Wayne (Yanwei) Zhang > Statistical Research > CNA > Email: Yanwei.Zhang at cna.com<mailto:Yanwei.Zhang at cna.com> > > > > > NOTICE: ?This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. > If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. > If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety. > > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > 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. >NOTICE: This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety.
Hi: This might be helpful: http://glmm.wikidot.com/faq It has some commentary and potential workarounds re gamma GLMMs toward the bottom of the page. I might also suggest that r-sig-mixed-models might be a better forum for these types of questions; you can subscribe by going to http://www.r-project.org/mail.html and scrolling down a bit to where the special interest groups are listed. Same process as subscribing to R-help. Dennis On Tue, May 17, 2011 at 7:17 AM, Zhang,Yanwei <Yanwei.Zhang at cna.com> wrote:> Hi Dennis, > > Thanks for your input. I am quite shocked that the Gamma model is not implemented in glmer yet, but do you happen to know if there is other packages in R that could fit Gamma hierarchical models besides MASS? The glmmPQL results are a bit off compared to what have been published. ?Thanks. > > > > Wayne (Yanwei) Zhang > Statistical Research >>CNA > > -----Original Message----- > From: Dennis Murphy [mailto:djmuser at gmail.com] > Sent: Tuesday, May 17, 2011 3:05 AM > To: Zhang,Yanwei > Cc: r-help at r-project.org > Subject: Re: [R] Problem on glmer > > Hi: > > I believe that the problem with the glmer() call is that you are > invoking the Gamma family, which appears not to be implemented at > present. See > > http://r-project.markmail.org/search/?q=r-sig-mixed-models#query:r-sig-mixed-models%20date%3A201105%20+page:8+mid:q2b2uoupejdzmvtl+state:results > > Dennis > > On Mon, May 16, 2011 at 3:31 PM, Zhang,Yanwei <Yanwei.Zhang at cna.com> wrote: >> Hi all, >> >> I was trying to fit a Gamma hierarchical model using "glmer", but got weird error message that I could not understand. On the other hand, a similar call to the glmmPQL leads to results that are close to what I expect. I also tried to change tha "nAGQ" argument in "glmer", but it did not solve the problem. The model I was fitting has a simple structure - one hierarchy and one offset, and I did not expect "glmer" to fail on this simple data set. Did I specify something wrong here? I appreciate any help. Thanks in advance. >> >> ############################################################# >> ## Code to fit the model >> ############################################################# >> library(MASS) >> library(lme4) >> address <- "http://instruction.bus.wisc.edu/jfrees/jfreesbooks/Longitudinal%20and%20Panel%20Data/Book/Data/TXTData/WorkerC.txt" >> wc <- read.table(address,header=T) >> # clean data >> cl <- unique(wc$CL[wc$LOSS==0]) >> wc2 <- subset(wc, !(CL %in% cl) ) >> # fit model >> fit1 <- glmmPQL(LOSS~ YR+offset(log(PR)), random=~1|CL,family=Gamma(link="log"), >> ? ? ? ? ?data=wc2) >> fit2 <- glmer(LOSS~ (1|CL) + YR+offset(log(PR)), family=Gamma(link="log"), >> ? ? ? ? ?data=wc2) >> >> >> ############################################################# >> ## Results of fit1 and fit2 >> ############################################################# >> >>> fit1 >> Linear mixed-effects model fit by maximum likelihood >> ?Data: wc2 >> ?Log-likelihood: NA >> ?Fixed: LOSS ~ YR + offset(log(PR)) >> ?(Intercept) ? ? ? ? ? ?YR >> -4.2830507147 ?0.0005085944 >> >> Random effects: >> ?Formula: ~1 | CL >> ? ? ? ?(Intercept) ?Residual >> StdDev: ? 0.8313193 0.5346455 >> >> Variance function: >> ?Structure: fixed weights >> ?Formula: ~invwt >> Number of Observations: 700 >> Number of Groups: 100 >>> fit2 >> Error in asMethod(object) : matrix is not symmetric [1,2] >> >> >> >> >> Wayne (Yanwei) Zhang >> Statistical Research >> CNA >> Email: Yanwei.Zhang at cna.com<mailto:Yanwei.Zhang at cna.com> >> >> >> >> >> NOTICE: ?This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. >> If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. >> If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety. >> >> ? ? ? ?[[alternative HTML version deleted]] >> >> ______________________________________________ >> 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. >> > > NOTICE: ?This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. > If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. > If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety. >
Hglm does the work! Thanks! Also, I find that the developing version of lme4, called lme4a, has the capability to fit Gamma models. And both lme4a and hglm produce results consistent with the published ones. Problems solved! Wayne (Yanwei) Zhang Statistical Research>CNA-----Original Message----- From: Dennis Murphy [mailto:djmuser at gmail.com] Sent: Tuesday, May 17, 2011 12:46 PM To: Zhang,Yanwei Cc: r-help at r-project.org Subject: Re: [R] Problem on glmer Hi: This might be helpful: http://glmm.wikidot.com/faq It has some commentary and potential workarounds re gamma GLMMs toward the bottom of the page. I might also suggest that r-sig-mixed-models might be a better forum for these types of questions; you can subscribe by going to http://www.r-project.org/mail.html and scrolling down a bit to where the special interest groups are listed. Same process as subscribing to R-help. Dennis On Tue, May 17, 2011 at 7:17 AM, Zhang,Yanwei <Yanwei.Zhang at cna.com> wrote:> Hi Dennis, > > Thanks for your input. I am quite shocked that the Gamma model is not implemented in glmer yet, but do you happen to know if there is other packages in R that could fit Gamma hierarchical models besides MASS? The glmmPQL results are a bit off compared to what have been published. ?Thanks. > > > > Wayne (Yanwei) Zhang > Statistical Research >>CNA > > -----Original Message----- > From: Dennis Murphy [mailto:djmuser at gmail.com] > Sent: Tuesday, May 17, 2011 3:05 AM > To: Zhang,Yanwei > Cc: r-help at r-project.org > Subject: Re: [R] Problem on glmer > > Hi: > > I believe that the problem with the glmer() call is that you are > invoking the Gamma family, which appears not to be implemented at > present. See > > http://r-project.markmail.org/search/?q=r-sig-mixed-models#query:r-sig-mixed-models%20date%3A201105%20+page:8+mid:q2b2uoupejdzmvtl+state:results > > Dennis > > On Mon, May 16, 2011 at 3:31 PM, Zhang,Yanwei <Yanwei.Zhang at cna.com> wrote: >> Hi all, >> >> I was trying to fit a Gamma hierarchical model using "glmer", but got weird error message that I could not understand. On the other hand, a similar call to the glmmPQL leads to results that are close to what I expect. I also tried to change tha "nAGQ" argument in "glmer", but it did not solve the problem. The model I was fitting has a simple structure - one hierarchy and one offset, and I did not expect "glmer" to fail on this simple data set. Did I specify something wrong here? I appreciate any help. Thanks in advance. >> >> ############################################################# >> ## Code to fit the model >> ############################################################# >> library(MASS) >> library(lme4) >> address <- "http://instruction.bus.wisc.edu/jfrees/jfreesbooks/Longitudinal%20and%20Panel%20Data/Book/Data/TXTData/WorkerC.txt" >> wc <- read.table(address,header=T) >> # clean data >> cl <- unique(wc$CL[wc$LOSS==0]) >> wc2 <- subset(wc, !(CL %in% cl) ) >> # fit model >> fit1 <- glmmPQL(LOSS~ YR+offset(log(PR)), random=~1|CL,family=Gamma(link="log"), >> ? ? ? ? ?data=wc2) >> fit2 <- glmer(LOSS~ (1|CL) + YR+offset(log(PR)), family=Gamma(link="log"), >> ? ? ? ? ?data=wc2) >> >> >> ############################################################# >> ## Results of fit1 and fit2 >> ############################################################# >> >>> fit1 >> Linear mixed-effects model fit by maximum likelihood >> ?Data: wc2 >> ?Log-likelihood: NA >> ?Fixed: LOSS ~ YR + offset(log(PR)) >> ?(Intercept) ? ? ? ? ? ?YR >> -4.2830507147 ?0.0005085944 >> >> Random effects: >> ?Formula: ~1 | CL >> ? ? ? ?(Intercept) ?Residual >> StdDev: ? 0.8313193 0.5346455 >> >> Variance function: >> ?Structure: fixed weights >> ?Formula: ~invwt >> Number of Observations: 700 >> Number of Groups: 100 >>> fit2 >> Error in asMethod(object) : matrix is not symmetric [1,2] >> >> >> >> >> Wayne (Yanwei) Zhang >> Statistical Research >> CNA >> Email: Yanwei.Zhang at cna.com<mailto:Yanwei.Zhang at cna.com> >> >> >> >> >> NOTICE: ?This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. >> If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. >> If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety. >> >> ? ? ? ?[[alternative HTML version deleted]] >> >> ______________________________________________ >> 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. >> > > NOTICE: ?This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. > If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. > If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety. >NOTICE: This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information. If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited. If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety.