x <- seq(0,8,by=0.1) plot(x,dgamma(x,0.75), type="l", ylab="", xlab="", ylim=c(0,1.25), xaxs="i", yaxs="i") plot(x,dgamma(x,1.0), type="l", ylab="", xlab="", ylim=c(0,1.25), xaxs="i",yaxs="i") plot(x,dgamma(x,2.0), type="l", ylab="", xlab="", ylim=c(0,1.25), xaxs="i",yaxs="i") data(wafer, package="faraway") summary(wafer) llmdl <- lm(log(resist) ~ .^2, wafer) rlmdl <- step(llmdl) library(faraway) sumary(rlmdl) gmdl <- glm(resist ~ .^2, family=Gamma(link=log), wafer) rgmdl <- step(gmdl) sumary(rgmdl) sqrt(summary(rgmdl)$dispersion) library(MASS) gamma.dispersion(rgmdl) data(motorins, package="faraway") motori <- motorins[motorins$Zone == 1,] gl <- glm(Payment ~ offset(log(Insured)) + as.numeric(Kilometres) + Make+Bonus, family=Gamma(link=log), motori) sumary(gl) llg <- glm(log(Payment) ~ offset(log(Insured))+as.numeric(Kilometres)+Make+Bonus,family=gaussian , motori) sumary(llg) x <- seq(0,5,by=0.05) plot(x,dgamma(x,1/0.55597,scale=0.55597),type="l",ylab="", xlab="",yaxs="i",ylim=c(0,1)) plot(x,dlnorm(x,meanlog=-0.30551,sdlog=sqrt(0.61102)),type="l", ylab="",xlab="",yaxs="i",ylim=c(0,1)) x0 <- data.frame(Make="1",Kilometres=1,Bonus=1,Insured=100) predict(gl,new=x0,se=T,type="response") predict(llg,new=x0,se=T,type="response") c(exp(10.998),exp(10.998)*0.16145) library(SuppDists) x <- seq(0,8,by=0.1) plot(x,dinvGauss(x,1,0.5),type="l",ylab="",xlab="",ylim=c(0,1.5), xaxs="i",yaxs="i") plot(x,dinvGauss(x,1,1),type="l",ylab="",xlab="",ylim=c(0,1.5), xaxs="i",yaxs="i") plot(x,dinvGauss(x,1,5),type="l",ylab="",xlab="",ylim=c(0,1.5), xaxs="i",yaxs="i") data(cpd, package="faraway") lmod <- lm(actual ~ projected-1,cpd) summary(lmod) plot(actual ~ projected, cpd) abline(lmod) igmod <- glm(actual ~ projected-1, family=inverse.gaussian(link="identity"), cpd) sumary(igmod) abline(igmod,lty=2) plot(residuals(igmod) ~ log(fitted(igmod)),ylab="Deviance residuals", xlab=expression(log(hat(mu)))) abline(h=0) data(weldstrength, package="faraway") lmod <- lm(Strength ~ Drying + Material + Preheating, weldstrength) sumary(lmod) h <- influence(lmod)$hat d <- residuals(lmod)^2/(1-h) gmod <- glm(d ~ Material+Preheating,family=Gamma(link=log), weldstrength,weights=1-h) w <- 1/fitted(gmod) lmod <- lm(Strength ~ Drying + Material + Preheating, weldstrength, weights=w) sumary(lmod) sumary(gmod) gl <- glm(resist ~ x1 + x2 + x3 + x4, family=Gamma(link=log), wafer) glq <- glm(resist ~ x1 + x2 + x3 + x4, family=quasi(link=log,variance="mu^2"), wafer) c(deviance(gl),df.residual(gl),logLik(gl)) c(deviance(glq),df.residual(glq),logLik(glq)) linkpow <- seq(0,1,by=0.1) devpow <- numeric(length(linkpow)) for(i in seq(along=linkpow)){ glqq <- glm(resist ~ x1 + x2 + x3 + x4, family=quasi(link=power(linkpow[i]),variance="mu^2"), wafer) devpow[i] <- deviance(glqq) } plot(linkpow, devpow, type="l", xlab="Link power", ylab="Deviance") powfam <- quasi(link="log",variance="mu^2") varpow <- seq(1,3,by=0.1) devpow <- numeric(length(varpow)) for(i in seq(along=varpow)){ powfam[["variance"]] <- function(mu) mu^varpow[i] glqq <- glm(resist ~ x1 + x2 + x3 + x4, family=powfam, wafer) devpow[i] <- deviance(glqq) } plot(varpow, devpow, type="l", xlab="Variance power", ylab="Deviance") data(chredlin, package="faraway") sum(chredlin$involact == 0) library(mgcv) twmod <- gam(involact ~ race+fire+theft+age+log(income), family=tw(link="log"), data=chicago) summary(twmod) xgrid <- seq(1e-10, 1.25, len=100) p <- 1.108 phi <- 0.30847 (mu <- twmod$fit[1]) (poismean <- mu^(2-p)/(phi * (2- p))) (p0 <- exp(-poismean)) twden <- exp(ldTweedie(xgrid, mu, p=p, phi=phi)[,1]) plot(xgrid, twden*(1-p0), type="l",xlab="x",ylab="Density") dmax <- max(twden*(1-p0)) segments(0,0,0,dmax,lwd=2) text(0.05,dmax,paste0("p=",signif(p0,3)))