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

1. library(faraway)
   
2. data(ctsib)
   
3. ctsib$stable <- ifelse(ctsib$CTSIB==1,1,0)
   
4. summary(ctsib)
   
5. gf <- glm(stable ~ Sex+Age+Height+Weight+Surface+Vision,binomial,data=ctsib)
   
6. summary(gf)
   
7. gfs <- glm(stable ~ Sex+Age+Height+Weight+Surface+Vision+factor(Subject),binomial,data=ctsib)
   
8. anova(gf,gfs,test="Chi")
   
9. library(brlr)
   
10. modbr <- brlr(stable ~ Sex+Age+Height+Weight+Surface+Vision+factor(Subject), data=ctsib)
    
11. summary(modbr)
    
12. library(MASS)
    
13. gg <- glmmPQL(stable ~ Sex+Age+Height+Weight+Surface+Vision,random=~1|Subject,  family=binomial,data=ctsib)
    
14. summary(gg)
    
15. sd(coef(modbr)[9:43])
    
16. data(ctsib)
    
17. ctsib$stable <- ifelse(ctsib$CTSIB==1,1,0)
    
18. library(gee)
    
19. gg <- gee(stable ~ Sex+Age+Height+Weight+Surface+Vision,id=Subject,family=binomial,data=ctsib,corstr="exchangeable",scale.fix=TRUE)
    
20. summary(gg)
    
21. data(epilepsy)
    
22. epilepsy[1:10,]
    
23. with(epilepsy,by(seizures/timeadj,list(treat,expind),mean))
    
24. y <- matrix(epilepsy$seizures,nrow=5)
    
25. matplot(1:4,sqrt(y[-1,]),type="l",lty=epilepsy$treat[5*(1:59)]+1,xlab="Period",ylab="Sqrt(Seizures)")
    
26. my <- apply(y[-1,],2,mean)/2
    
27. plot(sqrt(epilepsy$seizures[epilepsy$expind == 0]/8),sqrt(my),pch=epilepsy$treat[5*(1:59)]+2,xlab="sqrt(Baseline seizures)",ylab="sqrt(Experiment seizures)")
    
28. abline(0,1)
    
29. g <- gee(seizures ~offset(log(timeadj))+expind+treat+I(expind*treat),  id,family=poisson,corstr="AR-M",Mv=1,data=epilepsy,subset=(id!=49))
    
30. summary(g)

复制代码

1. library(faraway)
   
2. data(ozone)
   
3. olm <- lm(O3 ~ temp + ibh + ibt, ozone)
   
4. summary(olm)
   
5. library(gam)
   
6. amgam <- gam(O3 ~ lo(temp) + lo(ibh) + lo(ibt), data=ozone)
   
7. summary(amgam)
   
8. 1-5935.1/21115
   
9. amgamr <- gam(O3 ~ lo(temp) + lo(ibh) , data=ozone)
   
10. anova(amgamr,amgam,test="F")
    
11. plot(amgam,residuals=TRUE,se=TRUE,pch=".")
    
12. library(mgcv)
    
13. ammgcv <- gam(O3 ~ s(temp)+s(ibh)+s(ibt),data=ozone)
    
14. plot(ammgcv)
    
15. am1 <- gam(O3 ~ s(temp)+s(ibh),data=ozone)
    
16. am2 <- gam(O3 ~ temp+s(ibh),data=ozone)
    
17. anova(am2,am1,test="F")
    
18. amint <- gam(O3 ~ s(temp,ibh)+s(ibt),data=ozone)
    
19. summary(amint)
    
20. anova(ammgcv,amint,test="F")
    
21. plot(amint)
    
22. vis.gam(amint,theta=-45,color="gray")
    
23. rhs <- function(x,c) ifelse(x > c, x-c, 0)
    
24. lhs <- function(x,c) ifelse(x < c, c-x, 0)
    
25. olm2 <- lm(O3 ~ rhs(temp,60)+lhs(temp,60)+rhs(ibh,1000)+lhs(ibh,1000),ozone)
    
26. summary(olm2)
    
27. predict(ammgcv,data.frame(temp=60,ibh=2000,ibt=100),se=T)
    
28. predict(ammgcv,data.frame(temp=120,ibh=2000,ibt=100),se=T)
    
29. plot(predict(ammgcv),residuals(ammgcv),xlab="Predicted",ylab="Residuals")
    
30. qqnorm(residuals(ammgcv),main="")
    
31. amred <- gam(O3 ~ s(vh)+s(wind)+s(humidity)+s(temp)+s(dpg)+s(vis)+s(doy),data=ozone)
    
32. summary(amred)
    
33. alm <- lm(O3 ~ vis+doy+ibt+humidity+temp,data=ozone)
    
34. gammgcv <- gam(O3 ~ s(temp)+s(ibh)+s(ibt),family=poisson,scale=-1,data=ozone)
    
35. summary(gammgcv)
    
36. plot(gammgcv,residuals=TRUE)
    
37. x <- ozone[,c("temp","ibh","ibt")]
    
38. library(acepack)
    
39. acefit <- ace(x,ozone$O3)
    
40. summary(lm(acefit$ty ~ acefit$tx))
    
41. plot(ozone$O3,acefit$ty,xlab="O3",ylab=expression(theta(O3)))
    
42. plot(x[,1],acefit$tx[,1],xlab="temp",ylab="f(temp)")
    
43. plot(x[,2],acefit$tx[,2],xlab="ibh",ylab="f(ibh)")
    
44. plot(x[,3],acefit$tx[,3],xlab="ibt",ylab="f(ibt)")
    
45. x <- ozone[,-1]
    
46. acefit <- ace(x,ozone$O3)
    
47. summary(lm(acefit$ty ~ acefit$tx))
    
48. y <- cbind(ozone$O3,ozone$O3^2,sqrt(ozone$O3))
    
49. x <- ozone[,c("temp","ibh","ibt")]
    
50. cancor(x,y)
    
51. avasfit <- avas(x,ozone$O3)
    
52. plot(ozone$O3,avasfit$ty,xlab="O3",ylab=expression(theta(O3)))
    
53. plot(x[,1],avasfit$tx[,1],xlab="temp",ylab="f(temp)")
    
54. plot(x[,2],avasfit$tx[,2],xlab="ibh",ylab="f(ibh)")
    
55. plot(x[,3],avasfit$tx[,3],xlab="ibt",ylab="f(ibt)")
    
56. i <- order(ozone$O3)
    
57. plot(ozone$O3[i],avasfit$ty[i],type="l",xlab="O3",ylab=expression(theta(O3)))
    
58. gs <- lm(avasfit$ty[i] ~ sqrt(ozone$O3[i]))
    
59. lines(ozone$O3[i],gs$fit,lty=2)
    
60. gl <- lm(avasfit$ty[i] ~ log(ozone$O3[i]))
    
61. lines(ozone$O3[i],gl$fit,lty=5)
    
62. lmod <- lm(avasfit$ty ~ avasfit$tx)
    
63. summary(lmod)
    
64. plot(predict(lmod),residuals(lmod),xlab="Fitted",ylab="Residuals")
    
65. data(epilepsy)
    
66. egamm <- gamm(seizures ~ treat*expind+s(age),family=poisson,random=list(id=~1),data=epilepsy,subset=(id!=49))
    
67. summary(egamm$gam)
    
68. library(mda)
    
69. data(ozone)
    
70. a <- mars(ozone[,-1],ozone[,1])
    
71. summary(lm(ozone[,1] ~ a$x-1))
    
72. a <- mars(ozone[,-1],ozone[,1],nk=7)
    
73. summary(lm(ozone[,1] ~ a$x-1))
    
74. a <- mars(ozone[,-1],ozone[,1],nk=10,degree=2)
    
75. summary(lm(ozone[,1] ~ a$x-1))
    
76. a$factor[a$selected.terms,]
    
77. plot(ozone[,6],a$x[,4]*a$coef[4],xlab="ibh",ylab="Contribution of ibh")
    
78. plot(ozone[,10],a$x[,7]*a$coef[7]+a$x[,6]*a$coef[6],xlab="Day",ylab="Contribution of day")
    
79. humidity <- seq(10,100,len=20)
    
80. temp <- seq(20,100,len=20)
    
81. medians <- apply(ozone,2,median)
    
82. pdf <- matrix(medians,nrow=400,ncol=10,byrow=T)
    
83. pdf[,4] <- rep(humidity,20)
    
84. pdf[,5] <- rep(temp,rep(20,20))
    
85. pdf <- as.data.frame(pdf)
    
86. names(pdf) <- names(medians)
    
87. z <- predict.mars(a,pdf[,-1])
    
88. zm <- matrix(z,ncol=20,nrow=20)
    
89. contour(humidity,temp,zm,xlab="Humidity",ylab="Temperature")
    
90. persp(humidity,temp,zm,xlab="Humidity",ylab="Temperature",zlab="Ozone",theta=-30)
    
91. qqnorm(a$res,main="")
    
92. plot(a$fit,a$res,xlab="Fitted",ylab="Residuals")

复制代码

1. library(faraway)
   
2. data(ozone)
   
3. summary(ozone)
   
4. pairs(ozone,pch=".")
   
5. library(rpart)
   
6. (roz <- rpart(O3 ~ .,ozone))
   
7. plot(roz,margin=.10)
   
8. text(roz)
   
9. plot(roz,compress=T,uniform=T,branch=0.4,margin=.10)
   
10. text(roz)
    
11. plot(predict(roz),residuals(roz),xlab="Fitted",ylab="Residuals")
    
12. qqnorm(residuals(roz))
    
13. (x0 <- apply(ozone[,-1],2,median))
    
14. predict(roz,data.frame(t(x0)))
    
15. roze <- rpart(O3 ~ .,ozone,cp=0.001)
    
16. printcp(roze)
    
17. rozr <- prune.rpart(roze,0.01091)
    
18. plotcp(roz)
    
19. post(roz,filename="")
    
20. rozr <- prune.rpart(roz,0.0154)
    
21. 1-sum(residuals(rozr)^2)/sum((ozone$O3-mean(ozone$O3))^2)
    
22. data(kanga)
    
23. x0 <- c(1115,NA,748,182,NA,NA,178,311,756,226,NA,NA,NA,48,1009,NA,204,593)
    
24. kanga <- kanga[,c(T,F,!is.na(x0))]
    
25. kanga[1:2,]
    
26. apply(kanga,2,function(x) sum(is.na(x)))
    
27. round(cor(kanga[,-1],use="pairwise.complete.obs")[,c(3,9)],2)
    
28. newko <- na.omit(kanga[,-c(4,10)])
    
29. dim(newko)
    
30. dim(na.omit(kanga))
    
31. plot(foramina.length ~ zygomatic.width,data=newko,pch=substring(species,1,1))
    
32. kt <- rpart(species ~ ., data=newko,cp=0.001)
    
33. printcp(kt)
    
34. ktp <- prune(kt,cp=0.0211)
    
35. ktp
    
36. plot(ktp,compress=T,uniform=T,branch=0.4,margin=0.1)
    
37. text(ktp)
    
38. (tt <- table(actual=newko$species,predicted=predict(ktp,type="class")))
    
39. 1-sum(diag(tt))/sum(tt)
    
40. pck <- princomp(newko[,-1])
    
41. pcdf <- data.frame(species=newko$species,pck$scores)
    
42. kt <- rpart(species ~ ., pcdf,cp=0.001)
    
43. printcp(kt)
    
44. nx0 <- x0[!is.na(x0)]
    
45. nx0 <- nx0[-c(3,9)]
    
46. nx0 <- (nx0-pck$center)/pck$scale
    
47. nx0 %*% pck$loadings
    
48. ktp <- prune.rpart(kt,0.0421)
    
49. ktp
    
50. (tt <- table(newko$species,predict(ktp,type="class")))
    
51. 1-sum(diag(tt))/sum(tt)

复制代码

1. library(faraway)
   
2. library(nnet)
   
3. data(ozone)
   
4. nnmdl <- nnet(O3 ~ temp + ibh + ibt, ozone, size=2, linout=T)
   
5. sx <- scale(ozone)
   
6. bestrss <- 10000
   
7. for(i in 1:100){
   
8. nnmdl <- nnet(O3 ~ temp + ibh + ibt, sx, size=2, linout=T, trace=F)
   
9. cat(nnmdl$value,"\n")
   
10. if(nnmdl$value < bestrss){
    
11. bestnn <- nnmdl
    
12. bestrss <- nnmdl$value
    
13. }}
    
14. bestnn$value
    
15. summary(bestnn)
    
16. 1-88.03/sum((sx[,1]-mean(sx[,1]))^2)
    
17. ozmeans <- attributes(sx)[        DISCUZ_CODE_0        ]quot;scaled:center"
    
18. ozscales <- attributes(sx)[        DISCUZ_CODE_0        ]quot;scaled:scale"
    
19. xx <- expand.grid(temp=seq(-3,3,0.1),ibh=0,ibt=0)
    
20. plot(xx$temp*ozscales['temp']+ozmeans['temp'],predict(bestnn,new=xx)*ozscales['O3']+ozmeans['O3'],xlab="Temp",ylab="O3")
    
21. xx <- expand.grid(temp=0,ibh=seq(-3,3,0.1),ibt=0)
    
22. plot(xx$ibh*ozscales['ibh']+ozmeans['ibh'],predict(bestnn,new=xx)*ozscales['O3']+ozmeans['O3'],xlab="IBH",ylab="O3")
    
23. xx <- expand.grid(temp=0,ibh=0,ibt=seq(-3,3,0.1))
    
24. plot(xx$ibt*ozscales['ibt']+ozmeans['ibt'],predict(bestnn,new=xx)*ozscales['O3']+ozmeans['O3'],xlab="IBT",ylab="O3")
    
25. bestrss <- 10000
    
26. for(i in 1:100){
    
27. nnmdl <- nnet(O3 ~ temp + ibh + ibt, sx, size=2,linout=T,decay=0.001,trace=F)
    
28. cat(nnmdl$value,"\n")
    
29. if(nnmdl$value < bestrss){
    
30. bestnn <- nnmdl
    
31. bestrss <- nnmdl$value
    
32. }}
    
33. bestnn$value
    
34. xx <- expand.grid(temp=seq(-3,3,0.1),ibh=0,ibt=0)
    
35. plot(xx$temp*ozscales['temp']+ozmeans['temp'],  predict(bestnn,new=xx)*ozscales['O3']+ozmeans['O3'],xlab="Temp",ylab="O3")
    
36. xx <- expand.grid(temp=0,ibh=seq(-3,3,0.1),ibt=0)
    
37. plot(xx$ibh*ozscales['ibh']+ozmeans['ibh'],  predict(bestnn,new=xx)*ozscales['O3']+ozmeans['O3'],xlab="IBH",ylab="O3")
    
38. xx <- expand.grid(temp=0,ibh=0,ibt=seq(-3,3,0.1))
    
39. plot(xx$ibt*ozscales['ibt']+ozmeans['ibt'],  predict(bestnn,new=xx)*ozscales['O3']+ozmeans['O3'],xlab="IBT",ylab="O3")
    
40. bestrss <- 10000
    
41. for(i in 1:100){
    
42. nnmdl <- nnet(O3 ~ ., sx, size=4, linout=T,trace=F)
    
43. cat(nnmdl$value,"\n")
    
44. if(nnmdl$value < bestrss){
    
45. bestnn <- nnmdl
    
46. bestrss <- nnmdl$value
    
47. }}
    
48. 1-bestnn$value/sum((sx[,1]-mean(sx[,1]))^2)

复制代码

1. library(faraway)
   
2. loglik <- function(p,y,n) lchoose(n,y) + y*log(p) + (n-y)*log(1-p)
   
3. nloglik <- function(p,y,n) loglik(p,y,n) - loglik(y/n,y,n)
   
4. pr <- seq(0.05,0.95,by=0.01)
   
5. matplot(pr,cbind(nloglik(pr,10,25),nloglik(pr,20,50)),type="l",xlab="p",ylab="log-likelihood")
   
6. f <- function(x) -loglik(x,10,25)
   
7. mm <- nlm(f,0.5,hessian=T)
   
8. mm$estimate
   
9. c(1/mm$hessian,0.4*(1-0.4)/25)
   
10. pr <- seq(0.25,0.55,by=0.01)
    
11. plot(pr,nloglik(pr,40,100),type="l",xlab="p",ylab="log-likelihood")
    
12. abline(h=-qchisq(0.95,1)/2)
    
13. g <- function(x) nloglik(x,40,100)+qchisq(0.95,1)/2
    
14. uniroot(g,c(0.45,0.55))$root
    
15. uniroot(g,c(0.25,0.35))$root
    
16. abline(v=c(0.49765,0.30743))
    
17. se <- sqrt(0.4*(1-0.4)/100)
    
18. cv <- qnorm(0.975)
    
19. c(0.4-cv*se,0.4+cv*se)
    
20. (lrstat <- 2*(loglik(0.4,40,100)-loglik(0.5,40,100)))
    
21. pchisq(lrstat,1,lower=F)
    
22. (z <- (0.5-0.4)/se)
    
23. 2*pnorm(z,lower=F)
    
24. (sc <- 40/0.5-(100-40)/(1-0.5))
    
25. (obsinf <- 40/0.5^2+(100-40)/(1-0.5)^2)
    
26. (score.test <- 40*40/400)
    
27. pchisq(4,1,lower=F)

复制代码

怎么下载呢  楼主？