Issuu

SOLUTIONS MANUAL FOR Korosteleva, O. (2018). Advanced Regression Models with SAS and R, CRC Press By OLGA KOROSTELEVA Department of Mathematics and Statistics California State University, Long Beach

TABLE OF CONTENTS CHAPTER 1 ……………………………………………………………………………………. 3 CHAPTER 2 ……………………………………………………………………………………. 24 CHAPTER 3 ……………………………………………………………………………………. 58 CHAPTER 4 ……………………………………………………………………………………. 92 CHAPTER 5 ……………………………………………………………………………………. 131 CHAPTER 6 ……………………………………………………………………………………. 163 CHAPTER 7 ……………………………………………………………………………………. 187 CHAPTER 8 ……………………………………………………………………………………. 218 CHAPTER 9 ……………………………………………………………………………………. 284 CHAPTER 10 …………………………………………………………………………………. 315

CHAPTER 1 EXERCISE 1.1. Show that the normal distribution belongs to the exponential family of distributions. 𝑓(𝑦, 𝜇, 𝜎 ) = √

exp −

(

)

= exp − ln(2𝜋𝜎 ) −

(𝑦 − 2𝑦𝜇 + 𝜇 ) . Let 𝜃 = 𝜇

and 𝜙 = 𝜎 . Then, we can write 𝑓(𝑦, 𝜃, 𝜙) = exp − ln(2𝜋𝜙 ) −

= exp

− ln(2𝜋𝜙) −

= exp

( )

(𝑦 − 2𝑦𝜃 + 𝜃 )

+ ℎ(𝑦, 𝜙) where 𝑐(𝜃) =

, and

1 𝑦 ℎ(𝑦, 𝜙) = − ln(2𝜋𝜙) − . 2 2𝜙

EXERCISE 1.2. (a) Verify normality of the response variable, then fit the linear regression model to the data. State the fitted model. Give estimates for all parameters. In SAS: data weightloss; input drug$ age gender$ EWL @@; cards; A 49 F 14.2 A 54 M 25.4 A 37 F 14.1 A 34 F 15.9 A 51 F 17.4 A 54 F 22.8 A 44 M 8.4 A 56 M 11.2 A 44 M 17.3 B 51 M 21.9 B 44 F 23.6 B 53 F 23.8 B 26 M 14.1 B 56 F 24.6 B 28 F 17.8 B 52 F 15.7 B 54 F 23.7 ;

A 43 F 20.0 A 45 F 16.7 A 47 M 20.5 B 55 M 7.4 B 34 M 27.8

A 57 A 36 A 44 B 30 B 43

M 11.7 A 48 M 16.6 M 12.7 A 57 M 15.0 F 6.7 B 52 F 29.4 F 23.1 B 47 M 16.8 M 10.6 B 55 M 26.8

/*running normality check*/ proc univariate; var EWL; histogram/normal; run;

Goodness-of-Fit Tests for Normal Distribution Test Statistic p Value Kolmogorov-Smirnov D 0.10216310 Pr > D >0.150 Cramer-von Mises W-Sq 0.05103595 Pr > W-Sq >0.250 Anderson-Darling A-Sq 0.28788730 Pr > A-Sq >0.250

Based on the large p-values of the normality tests and the histogram, we can conclude that the response variable follows a normal distribution. /*fitting general linear model*/ proc genmod; class drug(ref="A") gender; model EWL = drug age gender / dist=normal link=identity; run; Log Likelihood -98.4395 Analysis Of Maximum Likelihood Parameter Estimates Parameter DF Estimate Standard Wald 95% Confidence Wald Chi- Pr > ChiSq Error Limits Square Intercept 1 9.2146 5.3301 -1.2322 19.6614 2.99 0.0838 drug B 1 4.8103 1.8697 1.1456 8.4749 6.62 0.0101 drug A 0 0.0000 0.0000 0.0000 0.0000 . . age 1 0.1102 0.1067 -0.0988 0.3192 1.07 0.3015 gender F 1 2.7235 1.8664 -0.9346 6.3815 2.13 0.1445 gender M 0 0.0000 0.0000 0.0000 0.0000 . . Scale 1 5.2451 0.6556 4.1054 6.7012

The fitted model is 𝐸 (𝐸𝑊𝐿) = 9.2146 + 4.8103 ∙ 𝑑𝑟𝑢𝑔𝐵 + 0.1102 ∙ 𝑎𝑔𝑒 + 2.7235 ∙ 𝑓𝑒𝑚𝑎𝑙𝑒 , and 𝜎 = 5.2451. In R: weightloss.data<- read.csv(file="C:/./Exercise1.2Data.csv", header = TRUE, sep = ",") #running normality check library(rcompanion) plotNormalHistogram(weightloss.data$EWL)

shapiro.test(weightloss.data$EWL) Shapiro-Wilk normality test W = 0.97424, p-value = 0.6234 #specifying reference levels drug.rel<- relevel(weightloss.data$drug, ref="A") gender.rel<- relevel(weightloss.data$gender, ref="M") #fitting general linear model summary(fitted.model<- glm(EWL ~ drug.rel + age + gender.rel, data = weightloss.data, family=gaussian(link=identity))) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.2146 5.6981 1.617 0.1171 drug.relB 4.8103 1.9988 2.407 0.0229 age 0.1102 0.1140 0.967 0.3420 gender.relF 2.7235 1.9952 1.365 0.1831 #outputting estimated sigma sigma(fitted.model) 5.607257

(b) Which regression coefficients turn out to be significant at the 5%? Discuss goodness of fit of the model. Drug B is the only significant predictor in the model at the 5% significance level since the corresponding p-value is the only one under 0.05. In SAS: /*checking model fit*/ proc genmod; model EWL = / dist=normal link=identity; run; Log Likelihood -102.6326 data deviance_test; deviance = -2*(-102.6326 - (-98.4395)); pvalue = 1 - probchi(deviance,3); run; proc print noobs; run; deviance 8.3862

pvalue 0.038669

The p-value for the deviance test is less than 0.05, indicating a good fit of the model. The R code and output are: #checking model fit null.model<- glm(EWL ~ 1, data=weightloss.data, family=gaussian(link=identity)) print(deviance<- -2*(logLik(null.model)-logLik(fitted.model)))

8.386158 print(p.value<- pchisq(deviance, df=3, lower.tail=FALSE)) 0.03867005

(c) Is one of the drugs more efficient for weight loss than the other? Interpret all estimated significant coefficients. The estimated average EWL for subjects taking drug B is 4.8103 percent higher than that for subjects taking drug A, keeping all the other predictors fixed. It means that drug B is more efficient than drug A. (d) According to the model, what is the predicted percent decrease in excess body weight for a 35year old male who is taking drug A? The predicted percent decrease in excess body weight for a 35-year old male who is taking drug A is computed by hand as: 𝐸𝑊𝐿 = 9.2146 + 0.1102 ∙ 35 = 13.0716. In SAS: /*using fitted model for prediction*/ data predict; input drug$ age gender$; cards; A 35 M ; data weightloss; set weightloss predict; run; proc genmod; class drug gender; model EWL = drug age gender / dist=normal link=identity; output out=outdata p=pEWL; run; proc print data=outdata (firstobs=33) noobs; var pEWL; run; pEWL 13.0718

In R: #using fitted model for prediction print(predict(fitted.model, data.frame(drug.rel="A", age=35, gender.rel="M"))) 13.7178

EXERCISE 1.3. (a) Reduce the car price by the factor of 1000. Check that the distribution of the price is normal. Fit a general linear regression model to predict the price of a car. Write down the fitted model, specifying all estimated parameters. In SAS: data carsales; input bodystyle$ 1-9 country$ hwy doors leather$ price @@; priceK=price/1000; cards; coupe USA 26 4 no 17445 coupe USA 40 4 no 23500 coupe USA 35 2 no 19600 coupe Germany 37 4 no 23400 coupe Germany 25 4 no 24100 coupe Germany 24 2 no 12400 coupe Japan 26 2 no 13300 coupe Japan 27 4 no 15550 coupe Japan 20 4 yes 29345 hatchback USA 30 2 no 12540 hatchback USA 39 4 no 17595 hatchback USA 38 2 no 17300 hatchback Germany 38 4 no 17800 hatchback Germany 32 4 no 22500 hatchback Germany 34 4 no 20300 hatchback Japan 38 4 yes 27300 hatchback Japan 38 2 yes 23300 hatchback Japan 38 2 yes 29300 sedan USA 29 4 no 32000 sedan USA 25 2 yes 34200 sedan USA 33 4 yes 33395 sedan Germany 40 4 no 22850 sedan Germany 23 2 yes 36000 sedan Germany 25 4 no 19900 sedan Japan 40 4 yes 36700 sedan Japan 35 4 yes 31600 sedan Japan 37 4 no 24600 run; /*running normality check*/ proc univariate; var priceK; histogram/normal; run;

Goodness-of-Fit Tests for Normal Distribution Test Statistic p Value Kolmogorov-Smirnov D 0.11287889 Pr > D >0.150 Cramer-von Mises W-Sq 0.05867848 Pr > W-Sq >0.250 Anderson-Darling A-Sq 0.37263698 Pr > A-Sq >0.250

P-values for the normality tests are all in excess of 0.05, indicating that normality holds. The histogram also displays a distribution close to bell-shaped. 7

/*fitting general linear model*/ proc genmod; class bodystyle(ref="hatchback") country(ref="Japan") leather(ref="no"); model priceK=bodystyle country hwy doors leather/dist=normal link=identity; run; Log Likelihood -67.2613 Analysis Of Maximum Likelihood Parameter Estimates Parameter DF Estimate Standard Wald 95% Confidence Wald Chi- Pr > ChiSq Error Limits Square Intercept 1 5.1353 4.6900 -4.0570 14.3276 1.20 0.2735 bodystyle coupe 1 2.2698 1.6836 -1.0301 5.5696 1.82 0.1776 bodystyle sedan 1 6.4107 1.5477 3.3772 9.4441 17.16 <.0001 bodystyle hatchback 0 0.0000 0.0000 0.0000 0.0000 . . country Germany 1 3.1959 1.6859 -0.1085 6.5002 3.59 0.0580 country USA 1 3.2128 1.5780 0.1199 6.3058 4.15 0.0418 country Japan 0 0.0000 0.0000 0.0000 0.0000 . . hwy 1 0.1305 0.1117 -0.0884 0.3494 1.36 0.2427 doors 1 1.5554 0.6630 0.2560 2.8549 5.50 0.0190 leather yes 1 12.1757 1.6217 8.9972 15.3541 56.37 <.0001 leather no 0 0.0000 0.0000 0.0000 0.0000 . . Scale 1 2.9219 0.3976 2.2378 3.8150

The fitted model is 𝐸 (𝑝𝑟𝑖𝑐𝑒𝐾) = 5.1353 + 2.2698 ∙ 𝑐𝑜𝑢𝑝𝑒 + 6.4107 ∙ 𝑠𝑒𝑑𝑎𝑛 + 3.1959 ∙ 𝐺𝑒𝑟𝑚𝑎𝑛𝑦 + 3.2128 ∙ 𝑈𝑆𝐴 + 0.1305 ∙ ℎ𝑤𝑦 + 1.5554 ∙ 𝑑𝑜𝑜𝑟𝑠 + 12.1757 ∙ 𝑙𝑒𝑎𝑡ℎ𝑒𝑟, and 𝜎 = 2.9219. In R: carsales.data<- read.csv(file="C:/./Exercise1.3Data.csv",header=TRUE, sep=",") #rescaling price priceK<- carsales.data$price/1000 #running normality check library(rcompanion) plotNormalHistogram(priceK)

shapiro.test(priceK) Shapiro-Wilk normality test W = 0.95482, p-value = 0.28

#specifying reference levels bodystyle.rel<- relevel(carsales.data$bodystyle, ref="hatchback") country.rel<- relevel(carsales.data$country, ref="Japan") leather.rel<- relevel(carsales.data$leather, ref="no") #fitting general linear model summary(fitted.model<- glm(priceK ~ bodystyle.rel + country.rel + hwy + doors + leather.rel, data=carsales.data, family=gaussian(link=identity))) Coefficients: (Intercept) bodystyle.relcoupe bodystyle.relsedan country.relGermany country.relUSA hwy doors leather.relyes

Estimate Std. Error t value Pr(>|t|) 5.1353 5.5909 0.919 0.36986 2.2698 2.0070 1.131 0.27216 6.4107 1.8450 3.475 0.00254 3.1959 2.0098 1.590 0.12829 3.2128 1.8812 1.708 0.10394 0.1305 0.1332 0.980 0.33937 1.5554 0.7904 1.968 0.06384 12.1757 1.9332 6.298 4.79e-06

#outputting estimated sigma sigma(fitted.model) 3.483088

(b) How good is the model fit? Discuss significance of the regression coefficients. The p-value in the deviance test is way below 0.05, indicating a good model fit. Significant variables are sedan body style and leather interior. In SAS: /*checking model fit*/ proc genmod; model priceK = / dist=normal link=identity; run; Log Likelihood -91.1942 data deviance_test; deviance = -2*(-91.1942 - (-67.2613)); pvalue = 1 - probchi(deviance,7); run; proc print noobs; run; deviance 47.8658

pvalue 3.7823E-8

In R: #checking model fit null.model<- glm(priceK ~ 1, data=carsales.data, family=gaussian(link=identity)) print(deviance<- -2*(logLik(null.model)-logLik(fitted.model))) 47.86586 print(p.value<- pchisq(deviance, df=7, lower.tail = FALSE))

3.78218e-08

(c) Interpret the estimates of those regression coefficients that differ significantly from zero. As estimated, sedan costs on average $6,410.70 more than a hatchback, under all other equal conditions. The estimated average price of a car with leather interior is $12,175.70 larger compared to a car without leather interior. (d) What is the predicted price of a sedan made in USA that has 4 doors, leather seats, and runs 30 mpg on highway? The predicted price of a sedan that is made in USA, has 4 doors, leather seats, and runs 30 mpg on highway is calculated as: 𝑝𝑟𝑖𝑐𝑒 = $1,000(5.1353 + 6.4107 + 3.2128 + 0.1305 ∙ 30 + 1.5554 ∙ 4 + 12.1757) = $37,071.10. In SAS: /*using fitted model for prediction*/ data predict; input bodystyle$ country$ hwy doors leather$; cards; sedan USA 30 4 yes ; data carsales; set carsales predict; run; proc genmod; class bodystyle country leather; model priceK = bodystyle country hwy doors leather / dist=normal link=identity; output out=outdata p=ppriceK; run; data final_prediction; set outdata; pprice=ppriceK*1000; run; proc print data=final_prediction (firstobs=28) noobs; var pprice; run; pprice 37071.14

In R: #using fitted model for prediction prediction<- (predict(fitted.model, data.frame(bodystyle.rel="sedan", country.rel ="USA", hwy=30, doors=4, leather.rel="yes"))) print(prediction*1000) 37071.14

EXERCISE 1.4. (a) Show normality of the distribution of the number of hours of sleep per night. Regress the number of hours of sleep on all the given factors. Write explicitly what the fitted model is. In SAS: data sleep; input age gender$ quiettime nchildren stresslevel jobstatus$ nactivities pastvac sleephours @@; cards; 62 F 60 1 5 unempl 1 15 7.7 28 F 15 1 6 unempl 5 11 5.3 50 M 15 0 5 unempl 1 19 6.4 36 M 60 1 6 full 1 21 7.7 56 F 50 0 3 part 4 5 7.6 48 M 180 0 5 full 0 6 6.4 55 M 40 0 8 full 8 23 7.0 26 F 80 0 7 student 9 8 8.3 44 M 180 1 3 part 6 20 9.6 49 F 5 0 7 unempl 5 15 5.5 29 M 60 2 5 student 5 7 7.7 56 M 10 1 4 unempl 4 17 5.7 46 F 40 1 7 part 3 3 7.4 41 F 5 2 6 full 9 10 6.2 22 M 15 0 8 full 4 3 6.3 36 F 45 2 5 part 8 14 7.5 54 F 120 1 8 part 7 10 8.5 42 F 60 3 1 full 9 11 6.3 58 F 5 1 7 full 1 17 5.3 33 M 100 2 1 full 9 5 8.3 50 F 2 2 6 full 3 12 5.1 59 M 30 2 5 full 2 6 6.9 32 M 30 1 8 full 5 9 6.9 50 M 60 2 8 part 8 13 8.0 56 F 10 0 3 unempl 7 7 6.1 42 F 240 0 1 part 8 21 8.8 58 F 10 2 7 full 9 4 6.2 57 F 15 1 6 full 2 16 6.3 30 F 30 0 2 full 8 9 8.3 54 M 20 2 8 full 6 7 6.5 57 M 45 2 4 full 7 18 7.5 45 F 120 0 9 part 2 13 6.6 33 F 40 1 6 unempl 9 24 7.0 56 F 120 0 5 part 2 20 8.7 59 F 60 2 9 part 4 19 8.1 41 M 60 2 3 student 2 3 7.5 62 M 40 0 1 unempl 0 2 8.6 29 M 15 1 7 unempl 3 20 6.3 34 F 30 0 7 unempl 9 0 6.6 32 F 20 3 7 unempl 2 8 7.8 46 F 20 2 3 unempl 9 18 7.9 45 M 60 0 2 unempl 0 22 9.0 23 M 45 0 6 part 4 12 7.6 38 M 60 4 5 full 3 5 7.8 45 M 30 0 5 unempl 9 7 6.8 63 F 40 0 6 unempl 5 5 7.3 27 F 120 0 4 student 1 16 7.3 30 F 45 0 7 part 8 10 7.7 34 F 5 3 6 full 0 4 6.0 62 M 10 0 10 part 8 11 6.0 ; /*running normality check*/ proc univariate; var sleephours; histogram/normal; run;