【独家发布】Extending the Linear Model with R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

Extending the Linear Model with RJulian Faraway

• The book is published by CRC press.
• Preface and Table of Contents
• The book refers to many datasets that can be found in the “faraway” package that needs to be added to R. Windows and Macintosh users will find it most convenient to select the “Install packages from CRAN” option under the “Package” menu while running R, and then choose “faraway” from the list of packages. A source version of the package can be obtained from the R website.
• The data from the faraway package in ASCII format. These are provided for the convenience of users of statistical software other than R. Users of R should install the faraway package.
• The R commands used in text.
• The Errata.
• The lme4 package, used to fit models with random effects and used extensively in Chapters 8 and 9, has undergone substantial changes since the publication of the book. Please read about these changes along with some suggested solutions. (last updated October 2014).
  

扫码加我 拉你入群

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

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

分享0 收藏2 回帖

关键词：Extending Linear Extend ending model

1. library(faraway)
   
2. data(gavote)
   
3. help(gavote)
   
4. help(quantile)
   
5. help.search("quantiles")
   
6. help.start()
   
7. gavote
   
8. head(gavote)
   
9. summary(gavote)
   
10. gavote$undercount <- (gavote$ballots-gavote$votes)/gavote$ballots
    
11. summary(gavote$undercount)
    
12. sum(gavote$ballots-gavote$votes)/sum(gavote$ballots)
    
13. hist(gavote$undercount,main="Undercount",xlab="Percent Undercount")
    
14. plot(density(gavote$undercount),main="Undercount")
    
15. rug(gavote$undercount)
    
16. pie(table(gavote$equip),col=gray(0:4/4))
    
17. barplot(sort(table(gavote$equip),decreasing=TRUE),las=2)
    
18. gavote$pergore <- gavote$gore/gavote$votes
    
19. plot(pergore ~ perAA, gavote, xlab="Proportion African American",ylab="Proportion for Gore")
    
20. plot(undercount ~ equip, gavote, xlab="", las=3)
    
21. xtabs(~ atlanta + rural, gavote)
    
22. nix <- c(3,10,11,12)
    
23. cor(gavote[,nix])
    
24. lmod <- lm(undercount ~ pergore + perAA, gavote)
    
25. coef(lmod)
    
26. predict(lmod)
    
27. residuals(lmod)
    
28. deviance(lmod)
    
29. df.residual(lmod)
    
30. nrow(gavote)-length(coef(lmod))
    
31. sqrt(deviance(lmod)/df.residual(lmod))
    
32. lmodsum <- summary(lmod)
    
33. lmodsum$sigma
    
34. lmodsum$r.squared
    
35. cor(predict(lmod),gavote$undercount)^2
    
36. lmodsum$adj.r.squared
    
37. summary(lmod)
    
38. contr.treatment(5)
    
39. gavote$cpergore <- gavote$pergore - mean(gavote$pergore)
    
40. gavote$cperAA <- gavote$perAA - mean(gavote$perAA)
    
41. lmodi <- lm(undercount ~ cperAA+cpergore*rural+equip, gavote)
    
42. summary(lmodi)
    
43. anova(lmod,lmodi)
    
44. drop1(lmodi,test="F")
    
45. confint(lmodi)
    
46. plot(lmodi)
    
47. gavote[cooks.distance(lmodi) > 0.1,]
    
48. halfnorm(influence(lmodi)$hat)
    
49. gavote[influence(lmodi)$hat>0.3,]
    
50. termplot(lmodi,partial=TRUE,terms=1)
    
51. library(MASS)
    
52. rlmodi <- rlm(undercount ~ cperAA+cpergore*rural+equip, gavote)
    
53. summary(rlmodi)
    
54. wlmodi <- lm(undercount ~ cperAA+cpergore*rural+equip, gavote, weights=ballots)
    
55. sqrt(0.035*(1-0.035)/881)
    
56. plmodi <- lm(undercount ~ poly(cperAA,4)+cpergore*rural+equip, gavote)
    
57. summary(plmodi)
    
58. termplot(plmodi,partial=TRUE,terms=1)
    
59. library(splines)
    
60. blmodi <- lm(undercount ~ cperAA+bs(cpergore,4)+rural+equip, gavote)
    
61. termplot(blmodi,partial=TRUE,terms=2)
    
62. biglm <- lm(undercount ~ (equip+econ+rural+atlanta)^2+(equip+econ+rural+atlanta)*(perAA+pergore), gavote)
    
63. smallm <- step(biglm,trace=F)
    
64. drop1(smallm,test="F")
    
65. finalm <- lm(undercount~equip + econ  + perAA + equip:econ + equip:perAA, gavote)
    
66. summary(finalm)
    
67. pdf <- data.frame(econ=rep(levels(gavote$econ),5),equip=rep(levels(gavote$equip),rep(3,5)),perAA=0.233)
    
68. pp <- predict(finalm,new=pdf)
    
69. xtabs(round(pp,3) ~ econ + equip, pdf)
    
70. pdf <- data.frame(econ=rep("middle",15),equip=rep(levels(gavote$equip),rep(3,5)),perAA=rep(c(.11,0.23,0.35),5))
    
71. pp <- predict(finalm,new=pdf)
    
72. propAA <- gl(3,1,15,labels=c("low","medium","high"))
    
73. xtabs(round(pp,3) ~ propAA + equip,pdf)

复制代码

1. library(faraway)
   
2. data(orings)
   
3. plot(damage/6 ~ temp, orings, xlim=c(25,85), ylim = c(0,1), xlab="Temperature", ylab="Prob of damage")
   
4. lmod <- lm(damage/6 ~ temp, orings)
   
5. abline(lmod)
   
6. logitmod <- glm(cbind(damage,6-damage) ~ temp, family=binomial, orings)
   
7. summary(logitmod)
   
8. plot(damage/6 ~ temp, orings, xlim=c(25,85), ylim = c(0,1), xlab="Temperature", ylab="Prob of damage")
   
9. x <- seq(25,85,1)
   
10. lines(x,ilogit(11.6630-0.2162*x))
    
11. probitmod <- glm(cbind(damage,6-damage) ~ temp, family=binomial(link=probit), orings)
    
12. summary(probitmod)
    
13. lines(x,pnorm(5.5915-0.1058*x),lty=2)
    
14. ilogit(11.6630-0.2162*31)
    
15. pnorm(5.5915-0.1058*31)
    
16. pchisq(deviance(logitmod),df.residual(logitmod),lower=FALSE)
    
17. pchisq(38.9,22,lower=FALSE)
    
18. pchisq(38.9-16.9,1,lower=FALSE)
    
19. c(-0.2162-1.96*0.0532,-0.2162+1.96*0.0532)
    
20. library(MASS)
    
21. confint(logitmod)
    
22. data(babyfood)
    
23. xtabs(disease/(disease+nondisease)~sex+food,babyfood)
    
24. mdl <- glm(cbind(disease,nondisease) ~ sex + food, family=binomial,babyfood)
    
25. summary(mdl)
    
26. drop1(mdl,test="Chi")
    
27. exp(-0.669)
    
28. exp(c(-0.669-1.96*0.153,-0.669+1.96*0.153))
    
29. library(MASS)
    
30. exp(confint(mdl))
    
31. 1-exp(-0.173)
    
32. babyfood[c(1,3),]
    
33. data(bliss)
    
34. bliss
    
35. modl <- glm(cbind(dead,alive) ~ conc, family=binomial, data=bliss)
    
36. modp <- glm(cbind(dead,alive) ~ conc, family=binomial(link=probit), data=bliss)
    
37. modc <- glm(cbind(dead,alive) ~ conc, family=binomial(link=cloglog), data=bliss)
    
38. fitted(modl)
    
39. coef(modl)[1]+coef(modl)[2]*bliss$conc
    
40. ilogit(modl$lin)
    
41. cbind(fitted(modl),fitted(modp),fitted(modc))
    
42. x <- seq(-2,8,0.2)
    
43. pl <- ilogit(modl$coef[1]+modl$coef[2]*x)
    
44. pp <- pnorm(modp$coef[1]+modp$coef[2]*x)
    
45. pc <- 1-exp(-exp((modc$coef[1]+modc$coef[2]*x)))
    
46. plot(x,pl,type="l",ylab="Probability",xlab="Dose")
    
47. lines(x,pp,lty=2)
    
48. lines(x,pc,lty=5)
    
49. matplot(x,cbind(pp/pl,(1-pp)/(1-pl)),type="l",xlab="Dose",ylab="Ratio")
    
50. matplot(x,cbind(pc/pl,(1-pc)/(1-pl)),type="l",xlab="Dose",ylab="Ratio")
    
51. data(hormone)
    
52. plot(estrogen ~ androgen,data=hormone,pch=as.character(orientation))
    
53. modl <- glm(orientation ~ estrogen + androgen, hormone, family=binomial)
    
54. summary(modl)
    
55. abline(-84.5/90.2,100.9/90.2)
    
56. library(brlr)
    
57. modb <- brlr(orientation ~ estrogen + androgen, hormone, family=binomial)
    
58. summary(modb)
    
59. abline(-3.65/3.586,4.074/3.586,lty=2)
    
60. modl <- glm(cbind(dead,alive) ~ conc, family=binomial, data=bliss)
    
61. sum(residuals(modl,type="pearson")^2)
    
62. deviance(modl)
    
63. (1-exp((modl$dev-modl$null)/150))/(1-exp(-modl$null/150))
    
64. data(bliss)
    
65. modl <- glm(cbind(dead,alive) ~ conc, family=binomial,data=bliss)
    
66. lmodsum <- summary(modl)
    
67. x0 <- c(1,2.5)
    
68. eta0 <- sum(x0*coef(modl))
    
69. ilogit(eta0)
    
70. (cm <- lmodsum$cov.unscaled)
    
71. se <- sqrt( t(x0) %*% cm %*% x0)
    
72. ilogit(c(eta0-1.96*se,eta0+1.96*se))
    
73. predict(modl,newdata=data.frame(conc=2.5),se=T)
    
74. ilogit(c(0.58095-1.96*0.2263,0.58095+1.96*0.2263))
    
75. x0 <- c(1,-5)
    
76. se <- sqrt( t(x0) %*% cm %*% x0)
    
77. eta0 <- sum(x0*lmod$coef)
    
78. ilogit(c(eta0-1.96*se,eta0+1.96*se))
    
79. (ld50 <- -lmod$coef[1]/lmod$coef[2])
    
80. dr <- c(-1/lmod$coef[2],lmod$coef[1]/lmod$coef[2]^2)
    
81. sqrt(dr %*% lmodsum$cov.un %*% dr)[,]
    
82. c(2-1.96*0.178,2+1.96*0.178)
    
83. ed90 <- (logit(0.9)-lmod$coef[1])/lmod$coef[2]
    
84. ed90
    
85. library(MASS)
    
86. dose.p(lmod,p=c(0.5,0.9))
    
87. data(troutegg)
    
88. ftable(xtabs(cbind(survive,total) ~ location+period, troutegg))
    
89. bmod <- glm(cbind(survive,total-survive) ~ location+period, family=binomial,troutegg)
    
90. bmod
    
91. halfnorm(residuals(bmod))
    
92. elogits <- log((troutegg$survive+0.5)/(troutegg$total-troutegg$survive+0.5))
    
93. with(troutegg,interaction.plot(period,location,elogits))
    
94. (sigma2 <- sum(residuals(bmod,type="pearson")^2)/12)
    
95. drop1(bmod,scale=sigma2,test="F")
    
96. summary(bmod,dispersion=sigma2)
    
97. data(amlxray)
    
98. head(amlxray)
    
99. amlxray[amlxray$downs=="yes",1:4]
    
100. (ii <- which(amlxray$downs=="yes"))
     
101. ramlxray <- amlxray[-c(ii,ii+1),]
     
102. library(survival)
     
103. cmod <- clogit(disease ~ Sex+Mray+Fray+CnRay+strata(ID),ramlxray)
     
104. summary(cmod)
     
105. cmodr <- clogit(disease ~ Fray+unclass(CnRay)+strata(ID),ramlxray)
     
106. summary(cmodr)
     
107. gmod <- glm(disease ~ Fray+unclass(CnRay),family=binomial,ramlxray)
     
108. summary(gmod)
     
109. data(esoph)
     
110. help(esoph)

复制代码

1. library(faraway)
   
2. data(gala)
   
3. gala <- gala[,-2]
   
4. modl <- lm(Species ~ . , gala)
   
5. plot(predict(modl),residuals(modl),xlab="Fitted",ylab="Residuals")
   
6. modt <- lm(sqrt(Species) ~ . , gala)
   
7. plot(predict(modt),residuals(modt),xlab="Fitted",ylab="Residuals")
   
8. summary(modt)
   
9. modp <- glm(Species ~ .,family=poisson,gala)
   
10. summary(modp)
    
11. halfnorm(residuals(modp))
    
12. plot(log(fitted(modp)),log((gala$Species-fitted(modp))^2),xlab=expression(hat(mu)),ylab=expression((y-hat(mu))^2))
    
13. abline(0,1)
    
14. (dp <- sum(residuals(modp,type="pearson")^2)/modp$df.res)
    
15. summary(modp,dispersion=dp)
    
16. drop1(modp,test="F")
    
17. data(dicentric)
    
18. round(xtabs(ca/cells ~ doseamt+doserate, dicentric),2)
    
19. with(dicentric,interaction.plot(doseamt,doserate,ca/cells))
    
20. lmod <- lm(ca/cells ~ log(doserate)*factor(doseamt), dicentric)
    
21. summary(lmod)$adj
    
22. plot(residuals(lmod) ~ fitted(lmod),xlab="Fitted",ylab="Residuals")
    
23. abline(h=0)
    
24. dicentric$dosef <- factor(dicentric$doseamt)
    
25. pmod <- glm(ca ~ log(cells)+log(doserate)*dosef, family=poisson,dicentric)
    
26. summary(pmod)
    
27. rmod <- glm(ca ~ offset(log(cells))+log(doserate)*dosef, family=poisson,dicentric)
    
28. summary(rmod)
    
29. data(solder)
    
30. modp <- glm(skips ~ . , family=poisson, data=solder)
    
31. deviance(modp)
    
32. df.residual(modp)
    
33. modp2  <- glm(skips ~ (Opening +Solder + Mask + PadType + Panel)^2 , family=poisson, data=solder)
    
34. deviance(modp2)
    
35. pchisq(deviance(modp2),df.residual(modp2),lower=FALSE)
    
36. library(MASS)
    
37. modn <- glm(skips ~ .,negative.binomial(1),solder)
    
38. modn
    
39. modn <- glm.nb(skips ~ .,solder)
    
40. summary(modn)

复制代码

1. library(faraway)
   
2. y <- c(320,14,80,36)
   
3. particle <- gl(2,1,4,labels=c("no","yes"))
   
4. quality <- gl(2,2,labels=c("good","bad"))
   
5. wafer <- data.frame(y,particle,quality)
   
6. wafer
   
7. (ov <- xtabs(y ~ quality+particle))
   
8. modl <- glm(y ~ particle+quality, poisson)
   
9. summary(modl)
   
10. drop1(modl,test="Chi")
    
11. (t(model.matrix(modl)) %*% y)[,]
    
12. (pp <- prop.table( xtabs(y ~ particle)))
    
13. (qp <- prop.table( xtabs(y ~ quality)))
    
14. (fv <- outer(qp,pp)*450)
    
15. 2*sum(ov*log(ov/fv))
    
16. sum( (ov-fv)^2/fv)
    
17. prop.test(ov)
    
18. (m <- matrix(y,nrow=2))
    
19. modb <- glm(m ~ 1, family=binomial)
    
20. deviance(modb)
    
21. fisher.test(ov)
    
22. (320*36)/(14*80)
    
23. data(haireye)
    
24. haireye
    
25. (ct <- xtabs(y ~ hair + eye, haireye))
    
26. summary(ct)
    
27. dotchart(ct)
    
28. mosaicplot(ct,color=T,main="")
    
29. modc <- glm(y ~ hair+eye,family=poisson,haireye)
    
30. summary(modc)
    
31. z <- xtabs(residuals(modc,type="pearson")~hair+eye,haireye)
    
32. svdz <- svd(z,2,2)
    
33. leftsv <- svdz$u %*% diag(sqrt(svdz$d[1:2]))
    
34. rightsv <- svdz$v %*% diag(sqrt(svdz$d[1:2]))
    
35. ll <- 1.1*max(abs(rightsv),abs(leftsv))
    
36. plot(rbind(leftsv,rightsv),asp=1,xlim=c(-ll,ll),ylim=c(-ll,ll),xlab="SV1",ylab="SV2",type="n")
    
37. abline(h=0,v=0)
    
38. text(leftsv,dimnames(z)[[1]])
    
39. text(rightsv,dimnames(z)[[2]])
    
40. data(eyegrade)
    
41. (ct <- xtabs(y ~ right+left, eyegrade))
    
42. summary(ct)
    
43. (symfac <- factor(apply(eyegrade[,2:3],1,function(x) paste(sort(x),collapse="-"))))
    
44. mods <- glm(y ~ symfac, eyegrade, family=poisson)
    
45. c(deviance(mods),df.residual(mods))
    
46. pchisq(deviance(mods),df.residual(mods),lower=F)
    
47. round(xtabs(residuals(mods) ~ right+left, eyegrade),3)
    
48. margin.table(ct,1)
    
49. margin.table(ct,2)
    
50. modq <- glm(y ~ right+left+symfac, eyegrade, family=poisson)
    
51. pchisq(deviance(modq),df.residual(modq),lower=F)
    
52. anova(mods,modq,test="Chi")
    
53. modqi <- glm(y ~ right+left, eyegrade, family=poisson, subset=-c(1,6,11,16))
    
54. pchisq(deviance(modqi),df.residual(modqi),lower=F)
    
55. data(femsmoke)
    
56. femsmoke
    
57. (ct <- xtabs(y ~ smoker+dead,femsmoke))
    
58. prop.table(ct,1)
    
59. summary(ct)
    
60. (cta <- xtabs(y ~ smoker+dead,femsmoke, subset=(age=="55-64")))
    
61. prop.table(cta,1)
    
62. prop.table(xtabs(y ~ smoker+age, femsmoke),2)
    
63. fisher.test(cta)
    
64. ct3 <- xtabs(y ~ smoker+dead+age,femsmoke)
    
65. apply(ct3, 3, function(x) (x[1,1]*x[2,2])/(x[1,2]*x[2,1]))
    
66. mantelhaen.test(ct3,exact=TRUE)
    
67. summary(ct3)
    
68. modi <- glm(y ~ smoker + dead + age, femsmoke, family=poisson)
    
69. c(deviance(modi),df.residual(modi))
    
70. (coefsmoke <- exp(c(0,coef(modi)[2])))
    
71. coefsmoke/sum(coefsmoke)
    
72. prop.table(xtabs(y ~ smoker, femsmoke))
    
73. modj <- glm(y ~ smoker*dead + age, femsmoke, family=poisson)
    
74. c(deviance(modj),df.residual(modj))
    
75. modc <- glm(y ~ smoker*age + age*dead, femsmoke, family=poisson)
    
76. c(deviance(modc),df.residual(modc))
    
77. modu <- glm(y ~ (smoker+age+dead)^2, femsmoke, family=poisson)
    
78. ctf <- xtabs(fitted(modu) ~ smoker+dead+age,femsmoke)
    
79. apply(ctf, 3, function(x) (x[1,1]*x[2,2])/(x[1,2]*x[2,1]))
    
80. exp(coef(modu)['smokerno:deadno'])
    
81. modsat <- glm(y ~ smoker*age*dead, femsmoke, family=poisson)
    
82. drop1(modsat,test="Chi")
    
83. drop1(modu,test="Chi")
    
84. ybin <- matrix(femsmoke$y,ncol=2)
    
85. modbin <- glm(ybin ~ smoker*age, femsmoke[1:14,], family=binomial)
    
86. drop1(modbin,test="Chi")
    
87. modbinr <- glm(ybin ~ smoker+age, femsmoke[1:14,], family=binomial)
    
88. drop1(modbinr,test="Chi")
    
89. deviance(modu)
    
90. deviance(modbinr)
    
91. exp(-coef(modbinr)[2])
    
92. modbinull <- glm(ybin ~ 1, femsmoke[1:14,], family=binomial)
    
93. deviance(modbinull)
    
94. modj <- glm(y ~ smoker*age + dead, femsmoke, family=poisson)
    
95. deviance(modj)
    
96. data(nes96)
    
97. xtabs( ~ PID + educ, nes96)
    
98. (partyed <- as.data.frame.table(xtabs( ~ PID + educ, nes96)))
    
99. nomod <- glm(Freq ~ PID + educ, partyed, family= poisson)
    
100. pchisq(deviance(nomod),df.residual(nomod),lower=F)
     
101. partyed$oPID <- unclass(partyed$PID)
     
102. partyed$oeduc <- unclass(partyed$educ)
     
103. ormod <- glm(Freq ~ PID + educ + I(oPID*oeduc), partyed, family= poisson)
     
104. anova(nomod,ormod,test="Chi")
     
105. summary(ormod)$coef['I(oPID * oeduc)',]
     
106. apid <- c(1,2,5,6,7,10,11)
     
107. aedu <- c(1,1,1,2,2,3,3)
     
108. ormoda <- glm(Freq ~ PID + educ + I(apid[oPID]*aedu[oeduc]), partyed, family= poisson)
     
109. anova(nomod,ormoda,test="Chi")
     
110. round(xtabs(predict(ormod,type="response") ~ PID + educ, partyed),2)
     
111. log(39.28*28.85/(47.49*23.19))
     
112. round(xtabs(residuals(ormod,type="response") ~ PID + educ, partyed),2)
     
113. cmod <- glm(Freq ~ PID + educ + educ:oPID, partyed, family= poisson)
     
114. anova(nomod,cmod,test="Chi")
     
115. summary(cmod)$coef[14:19,]
     
116. anova(ormod,cmod,test="Chi")
     
117. aedu <- c(1,1,2,2,2,2,2)
     
118. ormodb <- glm(Freq ~ PID + educ + I(oPID*aedu[oeduc]), partyed, family= poisson)
     
119. deviance(ormodb)
     
120. deviance(ormod)

复制代码

1. library(faraway)
   
2. data(nes96)
   
3. sPID <-  nes96$PID
   
4. levels(sPID) <- c("Democrat","Democrat","Independent","Independent","Independent","Republican","Republican")
   
5. summary(sPID)
   
6. inca <- c(1.5,4,6,8,9.5,10.5,11.5,12.5,13.5,14.5,16,18.5,21,23.5,27.5,32.5,37.5,42.5,47.5,55,67.5,82.5,97.5,115)
   
7. nincome <- inca[unclass(nes96$income)]
   
8. summary(nincome)
   
9. table(nes96$educ)
   
10. matplot(prop.table(table(nes96$educ,sPID),1),type="l",xlab="Education",ylab="Proportion",lty=c(1,2,5))
    
11. cutinc <- cut(nincome,7)
    
12. il <- c(8,26,42,58,74,90,107)
    
13. matplot(il,prop.table(table(cutinc,sPID),1),lty=c(1,2,5),type="l",ylab="Proportion",xlab="Income")
    
14. cutage <- cut(nes96$age,7)
    
15. al <- c(24,34,44,54,65,75,85)
    
16. matplot(al,prop.table(table(cutage,sPID),1),lty=c(1,2,5),type="l",ylab="Proportion",xlab="Age")
    
17. library(nnet)
    
18. mmod <- multinom(sPID ~ age + educ + nincome, nes96)
    
19. mmodi <- step(mmod)
    
20. mmode <- multinom(sPID ~ age +  nincome, nes96)
    
21. deviance(mmode) - deviance(mmod)
    
22. pchisq(16.206,mmod$edf-mmode$edf,lower=F)
    
23. predict(mmodi,data.frame(nincome=il),type="probs")
    
24. predict(mmodi,data.frame(nincome=il))
    
25. summary(mmodi)
    
26. cc <- c(0,-1.17493,-0.95036)
    
27. exp(cc)/sum(exp(cc))
    
28. predict(mmodi,data.frame(nincome=0),type="probs")
    
29. (pp <- predict(mmodi,data.frame(nincome=c(0,1)),type="probs"))
    
30. log(pp[1,1]*pp[2,2]/(pp[1,2]*pp[2,1]))
    
31. log(pp[1,1]*pp[2,3]/(pp[1,3]*pp[2,1]))
    
32. sPID[1:4]
    
33. cm <- diag(3)[unclass(sPID),]
    
34. cm[1:4,]
    
35. y <- as.numeric(t(cm))
    
36. resp.factor <- gl(944,3)
    
37. cat.factor <- gl(3,1,3*944,labels=c("D","I","R"))
    
38. rnincome <- rep(nincome,each=3)
    
39. head(data.frame(y,resp.factor,cat.factor,rnincome))
    
40. nullmod <- glm(y ~ resp.factor + cat.factor, family=poisson)
    
41. glmod <- glm(y ~ resp.factor + cat.factor + cat.factor:rnincome, family=poisson)
    
42. deviance(glmod)
    
43. deviance(mmodi)
    
44. coef(glmod)[c(1,945:949)]
    
45. coef(mmodi)
    
46. 0.016087-0.017665
    
47. data(cns)
    
48. cns
    
49. cns$CNS <- cns$An+cns$Sp+cns$Other
    
50. plot(log(CNS/NoCNS) ~ Water, cns, pch=as.character(Work))
    
51. binmodw <- glm(cbind(CNS,NoCNS) ~ Water + Work, cns, family=binomial)
    
52. binmoda <- glm(cbind(CNS,NoCNS) ~ Area + Work, cns, family=binomial)
    
53. anova(binmodw,binmoda,test="Chi")
    
54. halfnorm(residuals(binmodw))
    
55. summary(binmodw)
    
56. exp(-0.339058)
    
57. cmmod <- multinom(cbind(An,Sp,Other) ~ Water + Work, cns)
    
58. nmod <- step(cmmod)
    
59. nmod
    
60. cc <- c(0,0.28963,-0.98083)
    
61. names(cc) <- c("An","Sp","Other")
    
62. exp(cc)/sum(exp(cc))
    
63. multinom(cbind(NoCNS,An,Sp,Other) ~ Water + Work, cns)
    
64. library(MASS)
    
65. pomod <- polr(sPID ~ age + educ + nincome, nes96)
    
66. c(deviance(pomod),pomod$edf)
    
67. c(deviance(mmod),mmod$edf)
    
68. pomodi <- step(pomod)
    
69. deviance(pomodi)-deviance(pomod)
    
70. pchisq(11.151,pomod$edf-pomodi$edf,lower=F)
    
71. pim <- prop.table(table(nincome,sPID),1)
    
72. logit(pim[,1])-logit(pim[,1]+pim[,2])
    
73. summary(pomodi)
    
74. ilogit(0.209)
    
75. ilogit(1.292)-ilogit(0.209)
    
76. predict(pomodi,data.frame(nincome=il,row.names=il),type="probs")
    
77. x <- seq(-4,4,by=0.05)
    
78. plot(x,dlogis(x),type="l")
    
79. abline(v=c(0.209,1.292))
    
80. abline(v=c(0.209,1.292)-50*0.013120,lty=2)
    
81. abline(v=c(0.209,1.292)-100*0.013120,lty=5)
    
82. opmod <- polr(sPID ~ nincome, method="probit")
    
83. summary(opmod)
    
84. dems <- pnorm(0.128-il*0.008182)
    
85. demind <- pnorm(0.798-il*0.008182)
    
86. cbind(dems,demind-dems,1-demind)
    
87. polr(sPID ~ nincome, method="cloglog")

复制代码

1. library(faraway)
   
2. data(bliss)
   
3. modl <- glm(cbind(dead,alive) ~ conc, family=binomial, bliss)
   
4. summary(modl)$coef
   
5. y <- bliss$dead/30;  mu <- y
   
6. eta <- logit(mu)
   
7. z <- eta + (y-mu)/(mu*(1-mu))
   
8. w <- 30*mu*(1-mu)
   
9. lmod <- lm(z  ~ conc, weights=w, bliss)
   
10. coef(lmod)
    
11. for(i in 1:5){
    
12. eta <- lmod$fit
    
13. mu <- ilogit(eta)
    
14. z <- eta + (y-mu)/(mu*(1-mu))
    
15. w <- 30*mu*(1-mu)
    
16. lmod <- lm(z  ~ bliss$conc, weights=w)
    
17. cat(i,coef(lmod),"\n")
    
18. }
    
19. summary(lmod)
    
20. xm <- model.matrix(lmod)
    
21. wm <- diag(w)
    
22. sqrt(diag(solve(t(xm) %*% wm %*% xm)))
    
23. summary(lmod)$coef[,2]/summary(lmod)$sigma
    
24. summary(modl)
    
25. 1-pchisq(deviance(modl),df.residual(modl))
    
26. anova(modl,test="Chi")
    
27. modl2 <- glm(cbind(dead,alive) ~ conc+I(conc^2), family=binomial,bliss)
    
28. anova(modl,modl2,test="Chi")
    
29. anova(modl2,test="Chi")
    
30. residuals(modl)
    
31. residuals(modl,"pearson")
    
32. residuals(modl,"response")
    
33. bliss$dead/30 - fitted(modl)
    
34. residuals(modl,"working")
    
35. modl$residuals
    
36. residuals(lmod)
    
37. influence(modl)$hat
    
38. rstudent(modl)
    
39. influence(modl)$coef
    
40. cooks.distance(modl)
    
41. data(gala)
    
42. gala <- gala[,-2]
    
43. modp <- glm(Species ~ .,family=poisson,gala)
    
44. plot(residuals(modp) ~ predict(modp,type="response"),xlab=expression(hat(mu)),ylab="Deviance residuals")
    
45. plot(residuals(modp) ~ predict(modp,type="link"),xlab=expression(hat(eta)),ylab="Deviance residuals")
    
46. plot(residuals(modp,type="response") ~ predict(modp,type="link"),xlab=expression(hat(eta)),ylab="Response residuals")
    
47. plot(Species ~ Area, gala)
    
48. plot(Species ~ log(Area), gala)
    
49. mu <- predict(modp,type="response")
    
50. z <- predict(modp)+(gala$Species-mu)/mu
    
51. plot(z ~ log(Area), gala,ylab="Linearized Response")
    
52. modpl <- glm(Species ~ log(Area) + log(Elevation) + log(Nearest) + log(Scruz+0.1) + log(Adjacent), family=poisson, gala)
    
53. c(deviance(modp),deviance(modpl))
    
54. mu <- predict(modpl,type="response")
    
55. u <- (gala$Species-mu)/mu + coef(modpl)[2]*log(gala$Area)
    
56. plot(u ~ log(Area), gala,ylab="Partial Residual")
    
57. abline(0,coef(modpl)[2])
    
58. z <- predict(modpl)+(gala$Species-mu)/mu
    
59. plot(z ~ predict(modpl), xlab="Linear predictor", ylab="Linearized Response")
    
60. halfnorm(rstudent(modpl))
    
61. gali <- influence(modpl)
    
62. halfnorm(gali$hat)
    
63. halfnorm(cooks.distance(modpl))
    
64. plot(gali$coef[,5],ylab="Change in Scruz coef",xlab="Case no.")
    
65. modplr <- glm(Species ~ log(Area) + log(Elevation) + log(Nearest) + log(Scruz+0.1) + log(Adjacent), family=poisson, gala, subset=-25)
    
66. cbind(coef(modpl),coef(modplr))
    
67. modpla <- glm(Species ~ log(Area)+log(Adjacent), family=poisson, gala)
    
68. dp <- sum(residuals(modpla,type="pearson")^2)/modpla$df.res
    
69. summary(modpla,dispersion=dp)

复制代码

1. library(faraway)
   
2. x <- seq(0,8,by=0.1)
   
3. plot(x,dgamma(x,0.75),type="l",ylab="",xlab="",ylim=c(0,1.25),xaxs="i",yaxs="i")
   
4. plot(x,dgamma(x,1.0),type="l",ylab="",xlab="",ylim=c(0,1.25),xaxs="i",yaxs="i")
   
5. plot(x,dgamma(x,2.0),type="l",ylab="",xlab="",ylim=c(0,1.25),xaxs="i",yaxs="i")
   
6. data(wafer)
   
7. summary(wafer)
   
8. llmdl <- lm(log(resist) ~ .^2, wafer)
   
9. rlmdl <- step(llmdl)
   
10. summary(rlmdl)
    
11. gmdl <- glm(resist ~ .^2, family=Gamma(link=log), wafer)
    
12. rgmdl <- step(gmdl)
    
13. summary(rgmdl)
    
14. sqrt(0.0045249)
    
15. library(MASS)
    
16. gamma.dispersion(rgmdl)
    
17. data(motorins)
    
18. motori <- motorins[motorins$Zone == 1,]
    
19. gl <- glm(Payment ~ offset(log(Insured))+as.numeric(Kilometres)+Make+Bonus , family=Gamma(link=log), motori)
    
20. summary(gl)
    
21. llg <- glm(log(Payment) ~ offset(log(Insured))+as.numeric(Kilometres)+Make+Bonus,family=gaussian ,  motori)
    
22. summary(llg)
    
23. x <- seq(0,5,by=0.05)
    
24. plot(x,dgamma(x,1/0.55597,scale=0.55597),type="l",ylab="",xlab="",yaxs="i",ylim=c(0,1))
    
25. plot(x,dlnorm(x,meanlog=-0.30551,sdlog=sqrt(0.61102)),type="l",ylab="",xlab="",yaxs="i",ylim=c(0,1))
    
26. x0 <- data.frame(Make="1",Kilometres=1,Bonus=1,Insured=100)
    
27. predict(gl,new=x0,se=T,type="response")
    
28. predict(llg,new=x0,se=T,type="response")
    
29. c(exp(10.998),exp(10.998)*0.16145)
    
30. library(SuppDists)
    
31. x <- seq(0,8,by=0.1)
    
32. plot(x,dinvGauss(x,1,0.5),type="l",ylab="",xlab="",ylim=c(0,1.5),xaxs="i",yaxs="i")
    
33. plot(x,dinvGauss(x,1,1),type="l",ylab="",xlab="",ylim=c(0,1.5),xaxs="i",yaxs="i")
    
34. plot(x,dinvGauss(x,1,5),type="l",ylab="",xlab="",ylim=c(0,1.5),xaxs="i",yaxs="i")
    
35. data(cpd)
    
36. lmod <- lm(actual ~ projected-1,cpd)
    
37. summary(lmod)
    
38. plot(actual ~ projected, cpd)
    
39. abline(lmod)
    
40. igmod <- glm(actual ~ projected-1, family=inverse.gaussian(link="identity"), cpd)
    
41. summary(igmod)
    
42. abline(igmod,lty=2)
    
43. plot(residuals(igmod) ~ log(fitted(igmod)),ylab="Deviance residuals",xlab=expression(log(hat(mu))))
    
44. abline(h=0)
    
45. data(weldstrength)
    
46. lmod <- lm(Strength ~ Drying + Material + Preheating, weldstrength)
    
47. summary(lmod)
    
48. h <- influence(lmod)$hat
    
49. d <- residuals(lmod)^2/(1-h)
    
50. gmod <- glm(d ~ Material+Preheating,family=Gamma(link=log),weldstrength,weights=1-h)
    
51. w <- 1/fitted(gmod)
    
52. lmod <- lm(Strength ~ Drying + Material + Preheating, weldstrength, weights=w)
    
53. summary(lmod)
    
54. summary(gmod)
    
55. data(mammalsleep)
    
56. mammalsleep$pdr <- with(mammalsleep, dream/sleep)
    
57. summary(mammalsleep$pdr)
    
58. modl <- glm(pdr ~ log(body)+log(brain)+log(lifespan)+log(gestation)+predation+exposure+danger, family=quasibinomial, mammalsleep)
    
59. drop1(modl,test="F")
    
60. modl <- glm(pdr ~ log(body)+log(lifespan)+danger, family=quasibinomial, mammalsleep )
    
61. summary(modl)
    
62. ll <- row.names(na.omit(mammalsleep[,c(1,6,10,11)]))
    
63. halfnorm(cooks.distance(modl),labs=ll)
    
64. plot(predict(modl),residuals(modl,type="pearson"),xlab="Linear Predictor", ylab="Pearson Residuals")

复制代码

1. library(faraway)
   
2. data(pulp)
   
3. op <- options(contrasts=c("contr.sum", "contr.poly"))
   
4. lmod <- aov(bright ~ operator, pulp)
   
5. summary(lmod)
   
6. coef(lmod)
   
7. options(op)
   
8. (0.447-0.106)/5
   
9. library(lme4)
   
10. mmod <- lmer(bright ~ 1+(1|operator), pulp)
    
11. summary(mmod)
    
12. smod <- lmer(bright ~ 1+(1|operator), pulp, method="ML")
    
13. summary(smod)
    
14. nullmod <- lm(bright ~ 1, pulp)
    
15. as.numeric(2*(logLik(smod)-logLik(nullmod)))
    
16. pchisq(2.5684,1,lower=FALSE)
    
17. y <- simulate(nullmod)
    
18. lrstat <- numeric(1000)
    
19. for(i in 1:1000){
    
20. y <- unlist(simulate(nullmod))
    
21. bnull <- lm(y ~ 1)
    
22. balt <- lmer(y ~ 1 + (1|operator),pulp,method="ML")
    
23. lrstat[i] <- as.numeric(2*(logLik(balt)-logLik(bnull)))
    
24. }
    
25. mean(lrstat < 0.00001)
    
26. mean(lrstat > 2.5684)
    
27. sqrt(0.02*0.98/1000)
    
28. ranef(mmod)$operator
    
29. (cc <- model.tables(lmod))
    
30. cc[[1]]$operator/ranef(mmod)$operator
    
31. fixef(mmod)+ranef(mmod)$operator
    
32. qqnorm(resid(mmod),main="")
    
33. plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals")
    
34. abline(0,0)
    
35. data(penicillin)
    
36. summary(penicillin)
    
37. op <- options(contrasts=c("contr.sum", "contr.poly"))
    
38. lmod <- aov(yield ~ blend + treat, penicillin)
    
39. summary(lmod)
    
40. coef(lmod)
    
41. mmod <- lmer(yield ~ treat + (1|blend), penicillin)
    
42. summary(mmod)
    
43. options(op)
    
44. ranef(mmod)$blend
    
45. anova(mmod)
    
46. amod <- lmer(yield ~ treat + (1|blend), penicillin, method="ML")
    
47. nmod <- lmer(yield ~ 1 + (1|blend), penicillin, method="ML")
    
48. anova(amod,nmod)
    
49. lrstat <- numeric(1000)
    
50. for(i in 1:1000){
    
51. ryield <- unlist(simulate(nmod))
    
52. nmodr <- lmer(ryield ~ 1 + (1|blend), penicillin, method="ML")
    
53. amodr <- lmer(ryield ~ treat + (1|blend), penicillin, method="ML")
    
54. lrstat[i] <- 2*(logLik(amodr)-logLik(nmodr))
    
55. }
    
56. plot(qchisq((1:1000)/1001,3),sort(lrstat),xlab=expression(chi[3]^2),ylab="Simulated LRT")
    
57. abline(0,1)
    
58. mean(lrstat > 4.05)
    
59. rmod <- lmer(yield ~ treat + (1|blend), penicillin)
    
60. nlmod <- lm(yield ~ treat, penicillin)
    
61. 2*(logLik(rmod)-logLik(nlmod,REML=TRUE))
    
62. lrstatf <- numeric(1000)
    
63. for(i in 1:1000){
    
64. ryield <-  unlist(simulate(nlmod))
    
65. nlmodr <- lm(ryield ~ treat, penicillin)
    
66. rmodr <- lmer(ryield ~ treat + (1|blend), penicillin)
    
67. lrstatf[i] <- 2*(logLik(rmodr)-logLik(nlmodr,REML=TRUE))
    
68. }
    
69. mean(lrstatf < 0.00001)
    
70. cs <- lrstatf[lrstatf > 0.00001]
    
71. ncs <- length(cs)
    
72. plot(qchisq((1:ncs)/(ncs+1),1),sort(cs),xlab=expression(chi[1]^2),ylab="Simulated LRT")
    
73. abline(0,1)
    
74. mean(lrstatf > 2.7629)
    
75. data(irrigation)
    
76. summary(irrigation)
    
77. lmod <- lmer(yield ~ irrigation * variety + (1|field) +(1|field:variety),data=irrigation)
    
78. logLik(lmod)
    
79. lmodr <- lmer(yield ~ irrigation * variety + (1|field),data=irrigation)
    
80. logLik(lmodr)
    
81. summary(lmodr)
    
82. anova(lmodr)
    
83. plot(fitted(lmodr),resid(lmodr),xlab="Fitted",ylab="Residuals")
    
84. qqnorm(resid(lmodr),main="")
    
85. mod <- lm(yield ~ irrigation * variety, data=irrigation)
    
86. anova(mod)
    
87. data(eggs)
    
88. summary(eggs)
    
89. cmod <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs)
    
90. summary(cmod)
    
91. cmodr <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs)
    
92. anova(cmod,cmodr)
    
93. VarCorr(cmodr)
    
94. lrstat <- numeric(1000)
    
95. for(i in 1:1000){
    
96. rFat <- unlist(simulate(cmodr))
    
97. nmod <- lmer(rFat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs)
    
98. amod <- lmer(rFat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs)
    
99. lrstat[i] <- 2*(logLik(amod)-logLik(nmod))
    
100. }
     
101. mean(lrstat < 0.00001)
     
102. 2*(logLik(cmod)-logLik(cmodr))
     
103. mean(lrstat > 1.6034)
     
104. data(abrasion)
     
105. matrix(abrasion$material,4,4)
     
106. lmod <- aov(wear ~ material + run + position, abrasion)
     
107. summary(lmod)
     
108. mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion)
     
109. anova(mmod)
     
110. summary(mmod)
     
111. data(jsp)
     
112. jspr <- jsp[jsp$year==2,]
     
113. plot(jitter(math)~jitter(raven),data=jspr,xlab="Raven score",ylab="Math score")
     
114. boxplot(math ~ social, data=jspr,xlab="Social class",ylab="Math score")
     
115. glin <- lm(math ~ raven*gender*social,jspr)
     
116. anova(glin)
     
117. glin <- lm(math ~ raven*social,jspr)
     
118. anova(glin)
     
119. glin <- lm(math ~ raven+social,jspr)
     
120. summary(glin)
     
121. table(jspr$school)
     
122. mmod <- lmer(math ~ raven*social*gender+(1|school)+(1|school:class),data=jspr)
     
123. anova(mmod)
     
124. jspr$craven <- jspr$raven-mean(jspr$raven)
     
125. mmod <- lmer(math ~ craven*social+(1|school)+(1|school:class),jspr)
     
126. summary(mmod)
     
127. qqnorm(resid(mmod),main="")
     
128. plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals")
     
129. qqnorm(ranef(mmod)$school[,1],main="School effects")
     
130. qqnorm(ranef(mmod)$"school:class"[,1],main="Class effects")
     
131. adjscores <- ranef(mmod)$school[,1]
     
132. rawscores <- coef(lm(math ~ school-1,jspr))
     
133. rawscores <- rawscores-mean(rawscores)
     
134. plot(rawscores,adjscores)
     
135. sint <- c(9,14,29)
     
136. text(rawscores[sint],adjscores[sint]+0.2,c("9","15","30"))
     
137. schraven <- lm(raven ~ school, jspr)$fit
     
138. mmodc <- lmer(math ~ craven*social+schraven*social+(1|school)+(1|school:class),jspr)
     
139. anova(mmodc)

复制代码

1. library(faraway)
   
2. data(psid)
   
3. head(psid)
   
4. library(lattice)
   
5. xyplot(income ~ year | person, psid, type="l", subset=(person < 21),strip=FALSE)
   
6. xyplot(log(income+100) ~ year | sex, psid, type="l")
   
7. lmod <- lm(log(income) ~ I(year-78), subset=(person==1), psid)
   
8. coef(lmod)
   
9. slopes <- numeric(85);intercepts <- numeric(85)
   
10. for(i in 1:85){
    
11. lmod <- lm(log(income) ~ I(year-78), subset=(person==i), psid)
    
12. intercepts[i] <- coef(lmod)[1]
    
13. slopes[i] <- coef(lmod)[2]
    
14. }
    
15. plot(intercepts,slopes,xlab="Intercept",ylab="Slope")
    
16. psex <- psid$sex[match(1:85,psid$person)]
    
17. boxplot(split(slopes,psex))
    
18. t.test(slopes[psex=="M"],slopes[psex=="F"])
    
19. t.test(intercepts[psex=="M"],intercepts[psex=="F"])
    
20. library(lme4)
    
21. psid$cyear <- psid$year-78
    
22. mmod <- lmer(log(income) ~ cyear*sex +age+educ+(cyear|person),psid)
    
23. summary(mmod)
    
24. qqmath(~resid(mmod) | sex, psid)
    
25. xyplot(resid(mmod) ~ fitted(mmod) | cut(educ,c(0,8.5,12.5,20)), psid, layout=c(3,1),xlab="Fitted",ylab="Residuals")
    
26. data(vision)
    
27. vision$npower <- rep(1:4,14)
    
28. xyplot(acuity~npower|subject,vision,type="l",groups=eye,lty=1:2,layout=c(4,2))
    
29. mmod <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision)
    
30. summary(mmod)
    
31. 4.64^2/(4.64^2+3.21^2+4.07^2)
    
32. (4.64^2+3.21^2)/(4.64^2+3.21^2+4.07^2)
    
33. anova(mmod)
    
34. mmodr <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision,subset=-43)
    
35. anova(mmodr)
    
36. summary(mmodr)
    
37. op <- options(contrasts=c("contr.helmert", "contr.poly"))
    
38. mmodr <- lmer(acuity~power+(1|subject)+(1|subject:eye),vision,subset=-43)
    
39. summary(mmodr)
    
40. options(op)
    
41. contr.helmert(4)
    
42. plot(resid(mmodr) ~ fitted(mmodr),xlab="Fitted",ylab="Residuals")
    
43. abline(h=0)
    
44. qqnorm(ranef(mmodr)$"subject:eye"[[1]],main="")
    
45. data(jsp)
    
46. jspr <- jsp[jsp$year==2,]
    
47. mjspr <- data.frame(rbind(jspr[,1:6],jspr[,1:6]),subject=factor(rep(c("english","math"),c(953,953))),score=c(jspr$english/100,jspr$math/40))
    
48. xyplot(score ~ raven| subject*gender, mjspr)
    
49. mjspr$craven <- mjspr$raven-mean(mjspr$raven)
    
50. mmod <- lmer(score ~ subject*gender+craven*subject+social+(1|school)+(1|school:class)+(1|school:class:id),mjspr)
    
51. anova(mmod)
    
52. summary(mmod)
    
53. 0.1013^2/(0.1013^2+0.1166^2)
    
54. xyplot(residuals(mmod) ~ fitted(mmod)|subject,mjspr,pch=".",xlab="Fitted",ylab="Residuals")

复制代码