Hierarchical Generalized Linear Models: The R Package HGLMMM

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Hierarchical Generalized Linear Models: The R Package HGLMMM

The R package HGLMMM has been developed to fit generalized linear models with random effects using the h-likelihood approach. The response variable is allowed to follow a binomial, Poisson, Gaussian or gamma distribution. The distribution of random effects can be specified as Gaussian, gamma, inverse-gamma or beta. Complex structures as multi-membership design or multilevel designs can be handled. Further, dispersion parameters of random components and the residual dispersion (overdispersion) can be modeled as a function of covariates. Overdispersion parameter can be fixed or estimated. Fixed effects in the mean structure can be estimated using extended likelihood or a first order Laplace approximation to the marginal likelihood. Dispersion parameters are estimated using first order adjusted profile likelihood
4.pdf (562.16 KB, 需要: 1 个论坛币)

2013-12-10 12:40:17 上传

Hierarchical Generalized Linear Models: The R Package HGLMMM
需要: 1 个论坛币  [购买]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏2 回帖

关键词：Hierarchical Generalized Generalize General package developed specified Complex package design

library("HGLMMM")

###########################
##### Salamander Data #####
###########################

RSal<-data.frame(int=rep(1,60))

modBin <- HGLMfit (DistResp = "Binomial", DistRand = c("Normal", "Normal"),
                Link = "Logit", LapFix = TRUE, ODEst = FALSE, ODEstVal = c(0),
                formulaMain = Mate ~ TypeF + TypeM + TypeF * TypeM +
                    (1|Female) + (1|Male),
                formulaOD = ~ 1, formulaRand = list(one = ~ 1, two = ~ 1),
                DataMain = salamander,DataRand = list(RSal, RSal),
                BinomialDen = rep(1, 360), INFO = TRUE, DEBUG = FALSE)

modBin
summary(modBin, V=TRUE)

################
# Normal Model #
################

cake$repbatch<-100*cake$Replicate+cake$Batch
R1Cake<-data.frame(int=rep(1,15))
R2Cake<-data.frame(int=rep(1,45))
modCake1<-HGLMfit(DistResp="Normal",DistRand=c("Normal","Normal"),Link="Identity",LapFix=FALSE,ODEst=TRUE,ODEstVal=c(0),
            formulaMain=Angle~as.factor(Recipe)+as.factor(Temperature)+as.factor(Recipe)*as.factor(Temperature)+(1|Replicate)+(1|repbatch),
            ,formulaOD=~1,formulaRand=list(one=~1,two=~1),
            DataMain=cake,DataRand=list(R1Cake,R2Cake),Offset=NULL,BinomialDen=rep(1,360),
            ,INFO=TRUE,DEBUG=FALSE)     

# Creating residual plots for the normal-normal-normal model #

# postscript("V:\\HomeDir\\576790\\EUR\\ZIPHGLM_Package\\Paper_JSS\\graphs\\cakeG1.eps",horiz=TRUE)
op<-par(mfrow=c(2,2),
    oma = c(1,1,2,1),
    mar = c(3,3,4,1) +1.2      # Margins
    )

res<-modCake1$Details$StdDevianceResidualY
# Scaled fitted values vs residuals #
xmat<-model.matrix(~as.factor(Recipe)+as.factor(Temperature)+as.factor(Recipe)*as.factor(Temperature),data=cake)
mu<-xmat%*%modCake1$Results$Beta
plot(mu,res,pch=18,col="black",xlab="Scaled Fitted Values",ylab="Deviance Residuals")
los<-loess.smooth(mu, res, span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)
   
# Scaled fitted values vs absolute residuals #
plot(mu,abs(res),pch=18,col="black",xlab="Scaled Fitted Values",ylab="Absolute Deviance Residuals")
los<-loess.smooth(mu, abs(res), span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)

# qqplot #
qqnorm(res,pch=18)
abline(h=0,v=0)
abline(a=0,b=1,lty=3)

# histogram of residuals #
breaks<-seq(-4,4,by=0.8)
hist(res,breaks=breaks,xlab="Deviance Residuals",main="Histogram")
par(op)
mtext("Diagnostics for Cake Model Normal",
      side=3, line=1.5, font=2, cex=2, col='black')
par(op)
# dev.off()

############################
##### Cake model Gamma #####
############################
cake$repbatch<-100*cake$Replicate+cake$Batch
R1Cake<-data.frame(int=rep(1,15))
R2Cake<-data.frame(int=rep(1,45))
modCake2<-HGLMfit(DistResp="Gamma",DistRand=c("Normal","Normal"),Link="Log",LapFix=FALSE,ODEst=TRUE,ODEstVal=c(0),
            formulaMain=Angle~as.factor(Recipe)+as.factor(Temperature)+as.factor(Recipe)*as.factor(Temperature)+(1|Replicate)+(1|repbatch),
            ,formulaOD=~1,formulaRand=list(one=~1,two=~1),
            DataMain=cake,DataRand=list(R1Cake,R2Cake),Offset=NULL,BinomialDen=rep(1,360),
            ,INFO=TRUE,DEBUG=FALSE)     

# Creating residual plots for the gamma-normal-normal model #

#postscript("V:\\HomeDir\\576790\\EUR\\ZIPHGLM_Package\\Paper_JSS\\graphs\\cakeG2.eps",horiz=TRUE)
op<-par(mfrow=c(2,2),
    oma = c(1,1,2,1),
    mar = c(3,3,4,1) +1.2      # Margins
    )

res<-modCake2$Details$StdDevianceResidualY
# Scaled fitted values vs residuals #
xmat<-model.matrix(~as.factor(Recipe)+as.factor(Temperature)+as.factor(Recipe)*as.factor(Temperature),data=cake)
eta<-xmat%*%modCake2$Results$Beta
mu<-exp(eta)
mu<-log(mu)
plot(mu,res,pch=18,col="black",xlab="Scaled Fitted Values",ylab="Deviance Residuals")
los<-loess.smooth(mu, res, span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)
   
# Scaled fitted values vs absolute residuals #
plot(mu,abs(res),pch=18,col="black",xlab="Scaled Fitted Values",ylab="Absolute Deviance Residuals")
los<-loess.smooth(mu, abs(res), span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)

# qqplot #
qqnorm(res,pch=18)
abline(h=0,v=0)
abline(a=0,b=1,lty=3)

# histogram of residuals #
breaks<-seq(-4,4,by=0.8)
hist(res,breaks=breaks,xlab="Deviance Residuals",main="Histogram")
par(op)
mtext("Diagnostics for Cake Model Gamma",
      side=3, line=1.5, font=2, cex=2, col='black')
par(op)
#dev.off()

##### Cake 3 model - gamma without interaction #####
modCake3<-HGLMfit(DistResp="Gamma",DistRand=c("Normal","Normal"),Link="Log",LapFix=FALSE,ODEst=TRUE,ODEstVal=c(0),
            formulaMain=Angle~as.factor(Recipe)+as.factor(Temperature)+(1|Replicate)+(1|repbatch),
            ,formulaOD=~1,formulaRand=list(one=~1,two=~1),
            DataMain=cake,DataRand=list(R1Cake,R2Cake),Offset=NULL,BinomialDen=rep(1,360),
            ,INFO=TRUE,DEBUG=FALSE)     
HGLMLRTest(modCake3,modCake2)

#postscript("V:\\HomeDir\\576790\\EUR\\ZIPHGLM_Package\\Paper_JSS\\graphs\\cakeG3.eps",horiz=TRUE)
op<-par(mfrow=c(2,2),
    oma = c(1,1,2,1),
    mar = c(3,3,4,1) +1.2      # Margins
    )

res<-modCake3$Details$StdDevianceResidualY
# Scaled fitted values vs residuals #
xmat<-model.matrix(~as.factor(Recipe)+as.factor(Temperature),data=cake)
eta<-xmat%*%modCake3$Results$Beta
mu<-exp(eta)
mu<-log(mu)
plot(mu,res,pch=18,col="black",xlab="Scaled Fitted Values",ylab="Deviance Residuals")
los<-loess.smooth(mu, res, span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)
   
# Scaled fitted values vs absolute residuals #
plot(mu,abs(res),pch=18,col="black",xlab="Scaled Fitted Values",ylab="Absolute Deviance Residuals")
los<-loess.smooth(mu, abs(res), span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)

# qqplot #
qqnorm(res,pch=18)
abline(h=0,v=0)
abline(a=0,b=1,lty=3)

# histogram of residuals #
breaks<-seq(-4,4,by=0.8)
hist(res,breaks=breaks,xlab="Deviance Residuals",main="Histogram")
par(op)
mtext("Diagnostics for Cake Model Gamma(Final)",
      side=3, line=1.5, font=2, cex=2, col='black')
par(op)
#dev.off()

##### Gamma models with no receipe effect #####
modCake4<-HGLMfit(DistResp="Gamma",DistRand=c("Normal","Normal"),Link="Log",LapFix=FALSE,ODEst=TRUE,ODEstVal=c(0),
            formulaMain=Angle~as.factor(Temperature)+(1|Replicate)+(1|repbatch),
            ,formulaOD=~1,formulaRand=list(one=~1,two=~1),
            DataMain=cake,DataRand=list(R1Cake,R2Cake),Offset=NULL,BinomialDen=rep(1,360),
            ,INFO=TRUE,DEBUG=FALSE)     
HGLMLRTest(modCake4,modCake3)

postscript("V:\\HomeDir\\576790\\EUR\\ZIPHGLM_Package\\Paper_JSS\\graphs\\cakeG4.eps",horiz=TRUE)
op<-par(mfrow=c(2,2),
    oma = c(1,1,2,1),
    mar = c(3,3,4,1) +1.2      # Margins
    )

res<-modCake4$Details$StdDevianceResidualY
# Scaled fitted values vs residuals #
xmat<-model.matrix(~as.factor(Temperature),data=cake)
eta<-xmat%*%modCake4$Results$Beta
mu<-exp(eta)
mu<-log(mu)
plot(mu,res,pch=18,col="black",xlab="Scaled Fitted Values",ylab="Deviance Residuals")
los<-loess.smooth(mu, res, span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)
   
# Scaled fitted values vs absolute residuals #
plot(mu,abs(res),pch=18,col="black",xlab="Scaled Fitted Values",ylab="Absolute Deviance Residuals")
los<-loess.smooth(mu, abs(res), span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)

# qqplot #
qqnorm(res,pch=18)
abline(h=0,v=0)
abline(a=0,b=1,lty=3)

# histogram of residuals #
breaks<-seq(-4,4,by=0.8)
hist(res,breaks=breaks,xlab="Deviance Residuals",main="Histogram")
par(op)
mtext("Diagnostics for Cake Model Gamma(Final)",
      side=3, line=1.5, font=2, cex=2, col='black')
par(op)
dev.off()

###########################
##### Rat data models #####
###########################

Rrat<-data.frame(WBC=tapply(rat$WhiteBloodCells,rat$Subject,mean),RBC=tapply(rat$RedBloodCells,rat$Subject,mean))
modRat2<-HGLMfit(DistResp="Poisson",DistRand=c("Normal"),Link="Log",LapFix=FALSE,ODEst=TRUE,ODEstVal=c(0),formulaMain=Y~WhiteBloodCells+RedBloodCells+as.factor(Drug)+(1|Subject),
            ,formulaOD=~1,formulaRand=list(one=~WBC+I(WBC^2)),
            DataMain=rat,DataRand=list(Rrat),
            ,INFO=TRUE,DEBUG=FALSE)  

##################################
# Residual component diagnostics #
##################################

#postscript("V:\\HomeDir\\576790\\EUR\\ZIPHGLM_Package\\Paper_JSS\\graphs\\ratQP1.eps",horiz=TRUE)
op<-par(mfrow=c(2,2),
    oma = c(1,1,2,1),
    mar = c(3,3,4,1) +1.4      # Margins
    )

res<-modRat2$Details$StdDevianceResidualY
# Scaled fitted values vs residuals #
xmat<-model.matrix(~WhiteBloodCells+RedBloodCells+as.factor(Drug),data=rat)
eta<-xmat%*%modRat2$Results$Beta
mu<-exp(eta)
mu<-sqrt(mu)
plot(mu,res,pch=18,col="black",xlab="Scaled Fitted Values",ylab="Deviance Residuals")
los<-loess.smooth(mu, res, span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)
   
# Scaled fitted values vs absolute residuals #
plot(mu,abs(res),pch=18,col="black",xlab="Scaled Fitted Values",ylab="Absolute Deviance Residuals")
los<-loess.smooth(mu, abs(res), span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)

# qqplot #
qqnorm(res,pch=18)
abline(h=0,v=0)
abline(a=0,b=1,lty=3)

# histogram of residuals #
breaks<-seq(-4,4,by=0.8)
hist(res,breaks=breaks,xlab="Deviance Residuals",main="Histogram")
par(op)
mtext("Diagnostics for Rat Model Quasi-Poisson (y|v)",
      side=3, line=1.5, font=2, cex=2, col='black')
par(op)
#dev.off()

################################
# Random component diagnostics #
################################

#postscript("V:\\HomeDir\\576790\\EUR\\ZIPHGLM_Package\\Paper_JSS\\graphs\\ratQP2.eps",horiz=TRUE)
op<-par(mfrow=c(2,2),
    oma = c(1,1,2,1),
    mar = c(3,3,4,1) +1.4      # Margins
    )

res<-modRat2$Details$StdDevianceResidualDisp[[1]]
# Scaled fitted values vs residuals #
xmat<-model.matrix(~~WBC+I(WBC^2),data=Rrat)
eta<-xmat%*%modRat2$Results$Dispersion
mu<-exp(eta)
mu<-log(mu)
mu<-as.vector(mu)
plot(mu,res,pch=18,col="black",xlab="Scaled Fitted Values",ylab="Deviance Residuals")
los<-loess.smooth(mu, res, span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)
   
# Scaled fitted values vs absolute residuals #
plot(mu,abs(res),pch=18,col="black",xlab="Scaled Fitted Values",ylab="Absolute Deviance Residuals")
los<-loess.smooth(mu, abs(res), span = 1/2, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50)
lines(los$x,los$y,col="black",lwd=2)

# qqplot #
qqnorm(res,pch=18)
abline(h=0,v=0)
abline(a=0,b=1,lty=3)

# histogram of residuals #
breaks<-seq(-3,2,by=1)
hist(res,breaks=breaks,xlab="Deviance Residuals",main="Histogram")
par(op)
mtext("Diagnostics for Rat Model Quasi-Poisson (v)",
      side=3, line=1.5, font=2, cex=2, col='black')
par(op)
#dev.off()

modRat1<-HGLMfit(DistResp="Poisson",DistRand=c("Normal"),Link="Log",LapFix=FALSE,ODEst=FALSE,ODEstVal=c(0),formulaMain=Y~WhiteBloodCells+RedBloodCells+as.factor(Drug)+(1|Subject),
            ,formulaOD=~1,formulaRand=list(one=~1),
            DataMain=rat,DataRand=list(Rrat),
            ,INFO=TRUE,DEBUG=FALSE)  

BootstrapEnvelopeHGLM(modRat1,19,9999)

收藏一下。