home *** CD-ROM | disk | FTP | other *** search
/ Big Green CD 8 / BGCD_8_Dev.iso / NEXTSTEP / UNIX / Educational / R-0.49-MI / R-0.49-I / demos / models / lm+glm < prev   
Encoding:
Text File  |  1996-11-24  |  4.0 KB  |  130 lines
  1. ###-*- R -*-
  2. ### Examples from: "An Introduction to Statistical Modelling"
  3. ### By Annette Dobson
  4.  
  5.  
  6. ## Plant Weight Data (Page 9)
  7. ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
  8. trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
  9. group <- gl(2,10,20,labels=c("Ctl","Trt"))
  10. weight <- c(ctl,trt)
  11. anova(lm(weight~group))
  12. summary(lm(weight~group-1))
  13.  
  14.  
  15. ## Birth Weight Data (Page 14)
  16. tmp <- matrix(c(
  17. 40,2968,38,2795,40,3163,35,2925,36,2625,37,2847,
  18. 41,3292,40,3473,37,2628,38,3176,40,3421,38,2975,
  19. 40,3317,36,2729,40,2935,38,2754,42,3210,39,2817,
  20. 40,3126,37,2539,36,2412,38,2991,39,2875,40,3231
  21. ),nc=2,byrow=T)
  22. age <- tmp[,1]
  23. birthw <- tmp[,2]
  24. sex <- gl(2,12,24,labels=c("M","F"))
  25. plot(age,birthw,col=as.integer(sex))
  26. summary(lm(birthw ~ sex + age), cor=T)
  27. summary(lm(birthw ~ sex + age -1), cor=T)
  28. summary(lm(birthw ~ sex + sex:age -1), cor=T)
  29. summary(glm(birthw ~ sex + age, family=gaussian()))
  30. summary(z <- glm(birthw ~ sex + age - 1, family=gaussian()))
  31.  
  32.  
  33. ## Poisson Regression Data (Page 42)
  34. x <- c(-1,-1,0,0,0,0,1,1,1)
  35. y <- c(2,3,6,7,8,9,10,12,15)
  36. summary(glm(y~x,family=poisson(link="identity")))
  37.  
  38.  
  39. ## Calorie Data (Page 45)
  40. tmp <- matrix(c(
  41. 33,33,100,14, 40,47, 92,15, 37,49,135,18, 27,35,144,12,
  42. 30,46,140,15, 43,52,101,15, 34,62, 95,14, 48,23,101,17,
  43. 30,32, 98,15, 38,42,105,14, 50,31,108,17, 51,61, 85,19,
  44. 30,63,130,19, 36,40,127,20, 41,50,109,15, 42,64,107,16,
  45. 46,56,117,18, 24,61,100,13, 35,48,118,18, 37,28,102,14),nc=4,byrow=T)
  46. carb <- tmp[,1]
  47. age <- tmp[,2]
  48. wgt <- tmp[,3]
  49. prot <- tmp[,4]
  50. summary(lm(carb~age+wgt+prot))
  51.  
  52.  
  53. ## Extended Plant Data (Page 59)
  54. ctl <-  c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
  55. trtA <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
  56. trtB <- c(6.31,5.12,5.54,5.50,5.37,5.29,4.92,6.15,5.80,5.26)
  57. group <- gl(3,10,30,labels=c("Ctl","A","B"))
  58. weight <- c(ctl,trtA,trtB)
  59. anova(lm(weight~group))
  60. summary(lm(weight~group))
  61.  
  62.  
  63. ## Fictitious Anova Data (Page 64)
  64. y <- c(6.8,6.6,5.3,6.1,7.5,7.4,7.2,6.5,7.8,9.1,8.8,9.1)
  65. a <- gl(3,4,12)
  66. b <- gl(2,2,12)
  67. anova(z <- lm(y~a*b))
  68.  
  69.  
  70. ## Achievement Scores (Page 70)
  71. y <- c(6,4,5,3,4,3,6,8,9,7,9,8,5,7,6,7,7,7,8,5,7)
  72. x <- c(3,1,3,1,2,1,4,4,5,5,4,3,1,2,3,2,2,3,4,1,4)
  73. m <- gl(3,7,21)
  74. print(anova(z <- lm(y~x+m)))
  75.  
  76.  
  77. ## Beetle Data (Page 78)
  78. tmp <- matrix(c(
  79. 1.6907,59, 6, 1.7242,60,13, 1.7552,62,18, 1.7842,56,28,
  80. 1.8113,63,52, 1.8369,59,53, 1.8610,62,61, 1.8839,60,60
  81. ),nc=3,byrow=T)
  82. dose <- tmp[,1]
  83. n <- tmp[,2]
  84. x <- tmp[,3]
  85. dead <- cbind(x, n-x)
  86. summary(z <- glm(dead~dose,family=binomial(link=logit)))
  87. summary(z <- glm(dead~dose,family=binomial(link=probit)))
  88. summary(z <- glm(dead~dose,family=binomial(link=cloglog)))
  89.  
  90.  
  91. ## Anther Data (Page 84)
  92. ## Note that the proportions below are not exactly
  93. ## in accord with the sample sizes quoted below.
  94. ## In particular, the value 0.555 does not seem sensible.
  95. n <- c(102,99,108,76,81,90)
  96. p <- c(0.539,0.525,0.528,0.724,0.617,0.555)
  97. # x <- round(n*p)
  98. x <- n*p
  99. y <- cbind(x,n-x)
  100. f <- rep(c(40,150,350),2)
  101. g <- gl(2,3,6)
  102. summary(glm(y~g*f,family=binomial(link="logit")))
  103. summary(glm(y~g+f,family=binomial(link="logit")))
  104. summary(glm(y~f,family=binomial(link="logit")))
  105.  
  106.  
  107. ## Tumour Data (Page 92)
  108. counts <- c(22,2,10,16,54,115,19,33,73,11,17,28)
  109. type <- gl(4,3,12,labels=c("freckle","superficial","nodular","indeterminate"))
  110. site <- gl(3,1,12,labels=c("head/neck","trunk","extremities"))
  111. data.frame(counts,type,site)
  112. summary(z <- glm(counts ~ type + site,family=poisson()))
  113.  
  114.  
  115. ## Randomized Controlled Trial (Page 93)
  116. counts <- c(18,17,15,20,10,20,25,13,12)
  117. outcome <- gl(3,1,9)
  118. treatment <- gl(3,3,9)
  119. summary(z <- glm(counts ~ outcome + treatment,family=poisson()))
  120.  
  121. ## Peptic Ulcers and Blood Groups
  122. counts <- c(579,4219,911,4578,246,3775,361,4532,291,5261,396,6598)
  123. group <- gl(2,1,12,labels=c("cases","controls"))
  124. blood <- gl(2,2,12,labels=c("A","O"))
  125. city <- gl(3,4,12,labels=c("London","Manchester","Newcastle"))
  126. cbind(codes(group),codes(blood),codes(city),counts)
  127. summary(z1 <- glm(counts ~ group*city+blood+group:blood, family=poisson()))
  128. summary(z2 <- glm(counts ~ group*city+blood, family=poisson()))
  129. anova(z2,z1)
  130.