Hierarchical Generalized Linear Models: The R Package HGLMMM

Lisrelchen

4303

收藏 2013-12-10

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 个论坛币马上下载

扫码加我拉你入群

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

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

全部回复

Lisrelchen

2013-12-10 12:41:55

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)

扫码加我拉你入群

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

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

Lisrelchen

2013-12-10 12:42:35

################
# 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()

扫码加我拉你入群

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

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

Lisrelchen

2013-12-10 12:43:18

############################
##### 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()

扫码加我拉你入群

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

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

Lisrelchen

2013-12-10 12:44:27

##### 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()

扫码加我拉你入群

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

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

Lisrelchen

2013-12-10 12:44:53

##### 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()

扫码加我拉你入群

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

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

点击查看更多内容…

Lisrelchen

2013-12-10 12:45:22

###########################
##### 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)

扫码加我拉你入群

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

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

Lisrelchen

2013-12-10 12:45:45

##################################
# 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()

扫码加我拉你入群

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

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

Lisrelchen

2013-12-10 12:46:07

################################
# 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)

扫码加我拉你入群

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

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

henryyhl

2014-4-2 14:12:20

收藏一下。

扫码加我拉你入群

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

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

witj

2015-2-12 11:23:18

已收藏，谢谢分享。

扫码加我拉你入群

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

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

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群

扫码加我拉你入群