LA Pollution-Mortality Study: 1970-1979, 508 observations, 6-day spacing. Weekly FILTERED data. The data were lowpass filtered, filtering out frequencies above 0.1 cycles per day. Mortality: (1) Mrt: Total Mortality (2) Rsp: Respiratory Mortality (3) Crd: Cardiovascular Mortality Weather: (4) Tmp: Temperature (5) Hum: Relative Humidity Pollution: (6) Crb: Carbon Monoxide (7) Slf: Sulfur Dioxideglm.LAshumway (8) Nit: Nitrogen Dioxide (9) Hdr: Hydrocarbons (10) Ozn: Ozone (11) Par: Particulates A <- matrix(scan("LAshumway"), byrow=T,ncol=11, dimnames=list(NULL,c('Mrt','Rsp','Crd','Tmp','Hum','Crb', 'Slf','Nit','Hdr','Ozn','Par'))) ###OK!!! LA <- data.frame(A) ##Adding Mrt(t-1)=y(t-1), Mrt(t-2)=y(t-2) ##Autoregress on y(t-1),y(t-2)!!! y1 <- c(169,LA$Mrt)[-509] y2 <- c(169,169,LA$Mrt)[-c(509,510)] LA.y1y2 <- cbind(LA,y1,y2) T1 <- c(mean(LA.y1y2$Tmp),LA.y1y2$Tmp)[-509] T2 <- c(mean(LA.y1y2$Tmp),mean(LA.y1y2$Tmp),LA.y1y2$Tmp)[-c(509,510)] LA1 <- cbind(LA.y1y2,T1,T2) y1y2T1LogC.Poisson <- glm(formula = Mrt ~ y1 + y2 + T1 + log(Crb), family = poisson, data = LA1) > is.data.frame(LA1) [1] TRUE Ridge Regression =================== library(MASS) ##Need this library help(lm.ridge) lm.ridge(formula, data, subset, na.action, lambda = 0, model = FALSE, x = FALSE, y = FALSE, contrasts = NULL, ...) lambda: A scalar or vector of ridge constants. coef: matrix of coefficients, one row for each value of 'lambda'. Note that these are not on the original scale and are for use by the 'coef' method. NOTE: Coefficients not on the original scale!!! ###Compare multiple regression with ridge reg: ##Nultp Reg -------------- y1y2T1LogC.LM <- lm(formula = Mrt ~ y1+y2+T1+log(Crb), data = LA1) > y1y2T1LogC.LM Call: lm(formula = Mrt ~ y1 + y2 + T1 + log(Crb), data = LA1) Coefficients: (Intercept) y1 y2 T1 log(Crb) 58.5753 0.3374 0.3240 -0.2242 7.7998 ##Ridge Reg ------------ y1y2T1LogC.Ridge <- lm.ridge(formula = Mrt ~ y1+y2+T1+log(Crb), lambda=seq(0,0.1,0.001), data = LA1) > names(y1y2T1LogC.Ridge) [1] "coef" "scales" "Inter" "lambda" "ym" "xm" "GCV" "kHKB" [9] "kLW" > y1y2T1LogC.Ridge$lambda [1] 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 [13] 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.020 0.021 0.022 0.023 [25] 0.024 0.025 0.026 0.027 0.028 0.029 0.030 0.031 0.032 0.033 0.034 0.035 [37] 0.036 0.037 0.038 0.039 0.040 0.041 0.042 0.043 0.044 0.045 0.046 0.047 [49] 0.048 0.049 0.050 0.051 0.052 0.053 0.054 0.055 0.056 0.057 0.058 0.059 [61] 0.060 0.061 0.062 0.063 0.064 0.065 0.066 0.067 0.068 0.069 0.070 0.071 [73] 0.072 0.073 0.074 0.075 0.076 0.077 0.078 0.079 0.080 0.081 0.082 0.083 [85] 0.084 0.085 0.086 0.087 0.088 0.089 0.090 0.091 0.092 0.093 0.094 0.095 [97] 0.096 0.097 0.098 0.099 0.100 > y1y2T1LogC.Ridge lambda y1 y2 T1 log(Crb) 0.000 58.57532 0.3373620 0.3239660 -0.2241764 7.799826 0.001 58.57546 0.3373616 0.3239657 -0.2241765 7.799818 0.002 58.57560 0.3373612 0.3239654 -0.2241767 7.799811 0.003 58.57574 0.3373609 0.3239651 -0.2241768 7.799803 0.004 58.57588 0.3373605 0.3239648 -0.2241770 7.799796 0.005 58.57602 0.3373601 0.3239644 -0.2241771 7.799788 0.006 58.57616 0.3373597 0.3239641 -0.2241773 7.799781 0.007 58.57631 0.3373594 0.3239638 -0.2241774 7.799774 0.008 58.57645 0.3373590 0.3239635 -0.2241775 7.799766 0.009 58.57659 0.3373586 0.3239632 -0.2241777 7.799759 0.010 58.57673 0.3373583 0.3239629 -0.2241778 7.799751 0.011 58.57687 0.3373579 0.3239625 -0.2241780 7.799744 0.012 58.57701 0.3373575 0.3239622 -0.2241781 7.799736 0.013 58.57716 0.3373572 0.3239619 -0.2241783 7.799729 0.014 58.57730 0.3373568 0.3239616 -0.2241784 7.799721 0.015 58.57744 0.3373564 0.3239613 -0.2241786 7.799714 0.016 58.57758 0.3373561 0.3239610 -0.2241787 7.799706 0.017 58.57772 0.3373557 0.3239606 -0.2241789 7.799699 0.018 58.57786 0.3373553 0.3239603 -0.2241790 7.799692 0.019 58.57800 0.3373550 0.3239600 -0.2241792 7.799684 0.020 58.57815 0.3373546 0.3239597 -0.2241793 7.799677 0.021 58.57829 0.3373542 0.3239594 -0.2241795 7.799669 0.022 58.57843 0.3373538 0.3239591 -0.2241796 7.799662 0.023 58.57857 0.3373535 0.3239587 -0.2241798 7.799654 0.024 58.57871 0.3373531 0.3239584 -0.2241799 7.799647 0.025 58.57885 0.3373527 0.3239581 -0.2241800 7.799639 0.026 58.57899 0.3373524 0.3239578 -0.2241802 7.799632 0.027 58.57914 0.3373520 0.3239575 -0.2241803 7.799624 0.028 58.57928 0.3373516 0.3239572 -0.2241805 7.799617 0.029 58.57942 0.3373513 0.3239568 -0.2241806 7.799610 0.030 58.57956 0.3373509 0.3239565 -0.2241808 7.799602 0.031 58.57970 0.3373505 0.3239562 -0.2241809 7.799595 0.032 58.57984 0.3373502 0.3239559 -0.2241811 7.799587 0.033 58.57998 0.3373498 0.3239556 -0.2241812 7.799580 0.034 58.58013 0.3373494 0.3239553 -0.2241814 7.799572 0.035 58.58027 0.3373491 0.3239549 -0.2241815 7.799565 0.036 58.58041 0.3373487 0.3239546 -0.2241817 7.799557 0.037 58.58055 0.3373483 0.3239543 -0.2241818 7.799550 0.038 58.58069 0.3373480 0.3239540 -0.2241820 7.799542 0.039 58.58083 0.3373476 0.3239537 -0.2241821 7.799535 0.040 58.58098 0.3373472 0.3239533 -0.2241823 7.799528 0.041 58.58112 0.3373468 0.3239530 -0.2241824 7.799520 0.042 58.58126 0.3373465 0.3239527 -0.2241825 7.799513 0.043 58.58140 0.3373461 0.3239524 -0.2241827 7.799505 0.044 58.58154 0.3373457 0.3239521 -0.2241828 7.799498 0.045 58.58168 0.3373454 0.3239518 -0.2241830 7.799490 0.046 58.58182 0.3373450 0.3239514 -0.2241831 7.799483 0.047 58.58197 0.3373446 0.3239511 -0.2241833 7.799475 0.048 58.58211 0.3373443 0.3239508 -0.2241834 7.799468 0.049 58.58225 0.3373439 0.3239505 -0.2241836 7.799460 0.050 58.58239 0.3373435 0.3239502 -0.2241837 7.799453 0.051 58.58253 0.3373432 0.3239499 -0.2241839 7.799446 0.052 58.58267 0.3373428 0.3239495 -0.2241840 7.799438 0.053 58.58281 0.3373424 0.3239492 -0.2241842 7.799431 0.054 58.58296 0.3373421 0.3239489 -0.2241843 7.799423 0.055 58.58310 0.3373417 0.3239486 -0.2241845 7.799416 0.056 58.58324 0.3373413 0.3239483 -0.2241846 7.799408 0.057 58.58338 0.3373410 0.3239480 -0.2241848 7.799401 0.058 58.58352 0.3373406 0.3239476 -0.2241849 7.799393 0.059 58.58366 0.3373402 0.3239473 -0.2241850 7.799386 0.060 58.58380 0.3373398 0.3239470 -0.2241852 7.799378 0.061 58.58395 0.3373395 0.3239467 -0.2241853 7.799371 0.062 58.58409 0.3373391 0.3239464 -0.2241855 7.799364 0.063 58.58423 0.3373387 0.3239461 -0.2241856 7.799356 0.064 58.58437 0.3373384 0.3239457 -0.2241858 7.799349 0.065 58.58451 0.3373380 0.3239454 -0.2241859 7.799341 0.066 58.58465 0.3373376 0.3239451 -0.2241861 7.799334 0.067 58.58479 0.3373373 0.3239448 -0.2241862 7.799326 0.068 58.58494 0.3373369 0.3239445 -0.2241864 7.799319 0.069 58.58508 0.3373365 0.3239442 -0.2241865 7.799311 0.070 58.58522 0.3373362 0.3239438 -0.2241867 7.799304 0.071 58.58536 0.3373358 0.3239435 -0.2241868 7.799296 0.072 58.58550 0.3373354 0.3239432 -0.2241870 7.799289 0.073 58.58564 0.3373351 0.3239429 -0.2241871 7.799282 0.074 58.58578 0.3373347 0.3239426 -0.2241873 7.799274 0.075 58.58593 0.3373343 0.3239422 -0.2241874 7.799267 0.076 58.58607 0.3373340 0.3239419 -0.2241875 7.799259 0.077 58.58621 0.3373336 0.3239416 -0.2241877 7.799252 0.078 58.58635 0.3373332 0.3239413 -0.2241878 7.799244 0.079 58.58649 0.3373328 0.3239410 -0.2241880 7.799237 0.080 58.58663 0.3373325 0.3239407 -0.2241881 7.799229 0.081 58.58678 0.3373321 0.3239403 -0.2241883 7.799222 0.082 58.58692 0.3373317 0.3239400 -0.2241884 7.799214 0.083 58.58706 0.3373314 0.3239397 -0.2241886 7.799207 0.084 58.58720 0.3373310 0.3239394 -0.2241887 7.799200 0.085 58.58734 0.3373306 0.3239391 -0.2241889 7.799192 0.086 58.58748 0.3373303 0.3239388 -0.2241890 7.799185 0.087 58.58762 0.3373299 0.3239384 -0.2241892 7.799177 0.088 58.58777 0.3373295 0.3239381 -0.2241893 7.799170 0.089 58.58791 0.3373292 0.3239378 -0.2241895 7.799162 0.090 58.58805 0.3373288 0.3239375 -0.2241896 7.799155 0.091 58.58819 0.3373284 0.3239372 -0.2241898 7.799147 0.092 58.58833 0.3373281 0.3239369 -0.2241899 7.799140 0.093 58.58847 0.3373277 0.3239365 -0.2241900 7.799132 0.094 58.58861 0.3373273 0.3239362 -0.2241902 7.799125 0.095 58.58876 0.3373270 0.3239359 -0.2241903 7.799118 0.096 58.58890 0.3373266 0.3239356 -0.2241905 7.799110 0.097 58.58904 0.3373262 0.3239353 -0.2241906 7.799103 0.098 58.58918 0.3373258 0.3239350 -0.2241908 7.799095 0.099 58.58932 0.3373255 0.3239346 -0.2241909 7.799088 0.100 58.58946 0.3373251 0.3239343 -0.2241911 7.799080 > y1y2T1LogC.Ridge$lambda [1] 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 [13] 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.020 0.021 0.022 0.023 [25] 0.024 0.025 0.026 0.027 0.028 0.029 0.030 0.031 0.032 0.033 0.034 0.035 [37] 0.036 0.037 0.038 0.039 0.040 0.041 0.042 0.043 0.044 0.045 0.046 0.047 [49] 0.048 0.049 0.050 0.051 0.052 0.053 0.054 0.055 0.056 0.057 0.058 0.059 [61] 0.060 0.061 0.062 0.063 0.064 0.065 0.066 0.067 0.068 0.069 0.070 0.071 [73] 0.072 0.073 0.074 0.075 0.076 0.077 0.078 0.079 0.080 0.081 0.082 0.083 [85] 0.084 0.085 0.086 0.087 0.088 0.089 0.090 0.091 0.092 0.093 0.094 0.095 [97] 0.096 0.097 0.098 0.099 0.100 To get "optimal k=lambda" use: select(y1y2T1LogC.Ridge <- lm.ridge(formula = Mrt ~ y1+y2+T1+log(Crb), data = LA1, lambda = seq(0,0.1,0.001))) modified HKB estimator is 2.018087 (HKB: Ridge k of Hoerl, Kennard, and Baldwin) modified L-W estimator is 0.8647981 (L-W: Ridge k of Lawless and Wang) smallest value of GCV at 0.1 (GCV: Generalized cross validaion. Another est of k) Check HKB: k=lambda=2.018087 lm.ridge(formula = Mrt ~ y1+y2+T1+log(Crb), lambda=2.018087, data = LA1) y1 y2 T1 log(Crb) 58.8597740 0.3366203 0.3233278 -0.2244688 7.7848219 Check L-W: k=lambda=0.8647981 lm.ridge(formula = Mrt ~ y1+y2+T1+log(Crb), lambda=0.8647981, data = LA1) y1 y2 T1 log(Crb) 58.6974782 0.3370436 0.3236921 -0.2243028 7.7933860 Check GCV: k=lambda=0.1 lm.ridge(formula = Mrt ~ y1+y2+T1+log(Crb), lambda=0.1, data = LA1) y1 y2 T1 log(Crb) 58.5894625 0.3373251 0.3239343 -0.2241911 7.7990803 > y1y2T1LogC.Ridge$kHKB [1] 2.018087 ###To get coefficients not on the original scale: > y1y2T1LogC.Ridge$coef 0.000 0.001 0.002 0.003 0.004 0.005 0.006 y1 4.779927 4.779922 4.779916 4.779911 4.779906 4.779901 4.779895 y2 4.590118 4.590113 4.590109 4.590104 4.590100 4.590095 4.590091 T1 -2.018298 -2.018299 -2.018301 -2.018302 -2.018303 -2.018304 -2.018306 log(Crb) 3.499758 3.499755 3.499752 3.499748 3.499745 3.499742 3.499738 0.007 0.008 0.009 0.010 0.011 0.012 0.013 y1 4.779890 4.779885 4.779880 4.779875 4.779869 4.779864 4.779859 y2 4.590086 4.590082 4.590077 4.590073 4.590068 4.590064 4.590059 T1 -2.018307 -2.018308 -2.018310 -2.018311 -2.018312 -2.018314 -2.018315 log(Crb) 3.499735 3.499732 3.499728 3.499725 3.499722 3.499718 3.499715 0.014 0.015 0.016 0.017 0.018 0.019 0.020 y1 4.779854 4.779848 4.779843 4.779838 4.779833 4.779828 4.779822 y2 4.590055 4.590050 4.590046 4.590041 4.590037 4.590032 4.590028 T1 -2.018316 -2.018318 -2.018319 -2.018320 -2.018322 -2.018323 -2.018324 log(Crb) 3.499712 3.499708 3.499705 3.499701 3.499698 3.499695 3.499691 0.021 0.022 0.023 0.024 0.025 0.026 0.027 y1 4.779817 4.779812 4.779807 4.779801 4.779796 4.779791 4.779786 y2 4.590023 4.590019 4.590014 4.590010 4.590005 4.590001 4.589996 T1 -2.018326 -2.018327 -2.018328 -2.018330 -2.018331 -2.018332 -2.018334 log(Crb) 3.499688 3.499685 3.499681 3.499678 3.499675 3.499671 3.499668 0.028 0.029 0.030 0.031 0.032 0.033 0.034 y1 4.779781 4.779775 4.779770 4.779765 4.779760 4.779754 4.779749 y2 4.589992 4.589987 4.589983 4.589978 4.589974 4.589969 4.589965 T1 -2.018335 -2.018336 -2.018338 -2.018339 -2.018340 -2.018342 -2.018343 log(Crb) 3.499665 3.499661 3.499658 3.499655 3.499651 3.499648 3.499645 0.035 0.036 0.037 0.038 0.039 0.040 0.041 y1 4.779744 4.779739 4.779734 4.779728 4.779723 4.779718 4.779713 y2 4.589960 4.589956 4.589951 4.589947 4.589942 4.589938 4.589933 T1 -2.018344 -2.018346 -2.018347 -2.018348 -2.018350 -2.018351 -2.018352 log(Crb) 3.499641 3.499638 3.499635 3.499631 3.499628 3.499625 3.499621 0.042 0.043 0.044 0.045 0.046 0.047 0.048 y1 4.779707 4.779702 4.779697 4.779692 4.779687 4.779681 4.779676 y2 4.589929 4.589924 4.589920 4.589915 4.589911 4.589906 4.589902 T1 -2.018353 -2.018355 -2.018356 -2.018357 -2.018359 -2.018360 -2.018361 log(Crb) 3.499618 3.499615 3.499611 3.499608 3.499604 3.499601 3.499598 0.049 0.050 0.051 0.052 0.053 0.054 0.055 y1 4.779671 4.779666 4.779660 4.779655 4.779650 4.779645 4.779640 y2 4.589897 4.589893 4.589888 4.589884 4.589880 4.589875 4.589871 T1 -2.018363 -2.018364 -2.018365 -2.018367 -2.018368 -2.018369 -2.018371 log(Crb) 3.499594 3.499591 3.499588 3.499584 3.499581 3.499578 3.499574 0.056 0.057 0.058 0.059 0.060 0.061 0.062 y1 4.779634 4.779629 4.779624 4.779619 4.779613 4.779608 4.779603 y2 4.589866 4.589862 4.589857 4.589853 4.589848 4.589844 4.589839 T1 -2.018372 -2.018373 -2.018375 -2.018376 -2.018377 -2.018379 -2.018380 log(Crb) 3.499571 3.499568 3.499564 3.499561 3.499558 3.499554 3.499551 0.063 0.064 0.065 0.066 0.067 0.068 0.069 y1 4.779598 4.779593 4.779587 4.779582 4.779577 4.779572 4.779567 y2 4.589835 4.589830 4.589826 4.589821 4.589817 4.589812 4.589808 T1 -2.018381 -2.018383 -2.018384 -2.018385 -2.018387 -2.018388 -2.018389 log(Crb) 3.499548 3.499544 3.499541 3.499538 3.499534 3.499531 3.499528 0.070 0.071 0.072 0.073 0.074 0.075 0.076 y1 4.779561 4.779556 4.779551 4.779546 4.779540 4.779535 4.779530 y2 4.589803 4.589799 4.589794 4.589790 4.589785 4.589781 4.589776 T1 -2.018391 -2.018392 -2.018393 -2.018395 -2.018396 -2.018397 -2.018398 log(Crb) 3.499524 3.499521 3.499518 3.499514 3.499511 3.499507 3.499504 0.077 0.078 0.079 0.080 0.081 0.082 0.083 y1 4.779525 4.779520 4.779514 4.779509 4.779504 4.779499 4.779493 y2 4.589772 4.589767 4.589763 4.589758 4.589754 4.589749 4.589745 T1 -2.018400 -2.018401 -2.018402 -2.018404 -2.018405 -2.018406 -2.018408 log(Crb) 3.499501 3.499497 3.499494 3.499491 3.499487 3.499484 3.499481 0.084 0.085 0.086 0.087 0.088 0.089 0.090 y1 4.779488 4.779483 4.779478 4.779473 4.779467 4.779462 4.779457 y2 4.589740 4.589736 4.589731 4.589727 4.589722 4.589718 4.589713 T1 -2.018409 -2.018410 -2.018412 -2.018413 -2.018414 -2.018416 -2.018417 log(Crb) 3.499477 3.499474 3.499471 3.499467 3.499464 3.499461 3.499457 0.091 0.092 0.093 0.094 0.095 0.096 0.097 y1 4.779452 4.779446 4.779441 4.779436 4.779431 4.779426 4.779420 y2 4.589709 4.589704 4.589700 4.589695 4.589691 4.589686 4.589682 T1 -2.018418 -2.018420 -2.018421 -2.018422 -2.018424 -2.018425 -2.018426 log(Crb) 3.499454 3.499451 3.499447 3.499444 3.499441 3.499437 3.499434 0.098 0.099 0.100 y1 4.779415 4.779410 4.779405 y2 4.589677 4.589673 4.589668 T1 -2.018428 -2.018429 -2.018430 log(Crb) 3.499431 3.499427 3.499424 > Can get y1y2T1LogC.Ridge$coef graphically for I think lambda=0.01 plot(lm.ridge(Mrt ~ y1+y2+T1+log(Crb), data = LA1,lambda = seq(0,0.1,0.001))) ===================================================== Example From R: plot(lm.ridge(y ~ ., longley, lambda = seq(0,0.1,0.001))) longley # not the same as the S-PLUS dataset names(longley)[1] <- "y" lm.ridge(y ~ ., longley) plot(lm.ridge(y ~ ., longley, lambda = seq(0,0.1,0.001))) select(lm.ridge(y ~ ., longley, lambda = seq(0,0.1,0.0001)))