Data Analysis

Final Project of 805Jing Zhao1. DescriptionBefore I come to America from China, many people have told me that the public security is poor in some places of America. So I noticed the data set in my SAS book of my undergraduate study. (The name of the book is “The Methods of Data Analysis and The System of SAS”). And the data set is the ratio of commitment of 4 kinds of crime in 50 states of America. I have checked the reference of the book and noticed that the data is from /about-us/cjis/ucr/ucr/ which is the website of FBI. I want to analysis the relationship between those three kinds of crime which are rape crime, assault crime and robbery with the murder crime. Since the murder crime is more serious and more evil than those three. And I want to know which of the three are more important to the commitment of murder. And if any of those three crimes will cause more murder when the ratio of them increase. So I use the data set.I set the ratio of rape crime as x1, the ratio of robbery as x2 and the ratio of assault crime as x3 these three predictors. I also set murder as y. There will be 50 observationsin the data set.2. Data SetThe chance of commitments of four kinds of crime in 50 states of America in every 100,000 people.3. Analysis the data set(1). Select the modelWe want to know the relationship between murder crime and other three kinds of crime. So we let y equals the murder ratio, and x1 is the ratio of rape crime, and x2 is the ratio of robbery, and x3 is the ratio of assault crime.>mydata=read.table("C:/Users/jacky/Desktop/crime.txt")> y=mydata[,1] #murder>x1=mydata[,2] #rape crime>x2=mydata[,3] #robbery>x3=mydata[,4] #assaultSince we have 3 predictors so we will have 3217-=kinds of regressions.> LM1=lm(y~x1)> LM2=lm(y~x2)> LM3=lm(y~x3)> LM4=lm(y~x1+x2)> LM5=lm(y~x1+x3)> LM6=lm(y~x2+x3)> LM7=lm(y~x1+x2+x3)These are all possible models. And then we want to check which is the best and most simple one. We use the AIC method.> n=length(y)>AIC(LM1,k=2)[1] 259.6966>AIC(LM2,k=2)[1] 268.7988>AIC(LM3,k=2)[1] 254.836>AIC(LM4,k=2)[1] 259.6864>AIC(LM5,k=2)[1] 253.9535>AIC(LM6,k=2)[1] 254.9248>AIC(LM7,k=2)[1] 255.1303Now, we can see that LM5( That is y~x1+x3) is the best, since AIC(LM5,k=2)=253.9535 is the smallest AIC value.Then we use 2R and 2R to check the result.adj>summary(LM1)Call:lm(formula = y ~ x1)Residuals:Min 1Q Median 3Q Max-6.1400 -1.8733 -0.4707 1.7497 8.1813Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 1.88378 1.15443 1.632 0.109x1 0.21607 0.04145 5.213 3.89e-06 ***---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 3.122 on 48 degrees of freedomMultiple R-squared: 0.3615, Adjusted R-squared: 0.3482F-statistic: 27.17 on 1 and 48 DF, p-value: 3.892e-06>summary(LM2)Call:lm(formula = y ~ x2)Residuals:Min 1Q Median 3Q Max-5.2968 -2.9440 -0.2921 2.2146 8.0922Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 4.816907 0.839399 5.739 6.27e-07 ***x2 0.021171 0.005529 3.829 0.000373 *** ---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 3.419 on 48 degrees of freedom Multiple R-squared: 0.234, Adjusted R-squared: 0.218 F-statistic: 14.66 on 1 and 48 DF, p-value: 0.0003727>summary(LM3)Call:lm(formula = y ~ x3)Residuals:Min 1Q Median 3Q Max-4.8518 -2.2835 -0.5756 1.5642 7.4113Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 2.158405 0.989237 2.182 0.0340 *x3 0.025015 0.004238 5.903 3.52e-07 *** ---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.974 on 48 degrees of freedom Multiple R-squared: 0.4206, Adjusted R-squared: 0.4085 F-statistic: 34.85 on 1 and 48 DF, p-value: 3.524e-07>summary(LM4)Call:lm(formula = y ~ x1 + x2)Residuals:Min 1Q Median 3Q Max-5.405 -1.775 -0.569 1.429 8.428Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 1.892178 1.143444 1.655 0.10463x1 0.174204 0.050935 3.420 0.00130 ** x2 0.008613 0.006203 1.389 0.17151---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 3.092 on 47 degrees of freedom Multiple R-squared: 0.3866, Adjusted R-squared: 0.3605 F-statistic: 14.81 on 2 and 47 DF, p-value: 1.027e-05>summary(LM5)Call:lm(formula = y ~ x1 + x3)Residuals:Min 1Q Median 3Q Max-5.221 -2.009 -0.512 1.771 7.832Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 1.296838 1.099797 1.179 0.24427x1 0.096301 0.057662 1.670 0.10155x3 0.017364 0.006189 2.806 0.00728 ** ---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.92 on 47 degrees of freedom Multiple R-squared: 0.4531, Adjusted R-squared: 0.4298 F-statistic: 19.47 on 2 and 47 DF, p-value: 6.939e-07>summary(LM6)Call:lm(formula = y ~ x2 + x3)Residuals:Min 1Q Median 3Q Max-5.1240 -1.8898 -0.5108 1.8066 7.7803Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 2.000245 0.987719 2.025 0.048556 *x2 0.007769 0.005741 1.353 0.182448x3 0.021201 0.005059 4.191 0.000122 *** ---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.948 on 47 degrees of freedomMultiple R-squared: 0.4423, Adjusted R-squared: 0.4186 F-statistic: 18.64 on 2 and 47 DF, p-value: 1.095e-06>summary(LM7) Call:lm(formula = y ~ x1 + x2 + x3)Residuals:Min 1Q Median 3Q Max -4.8858 -1.8890 -0.4094 1.8296 8.0070Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 1.343193 1.103849 1.217 0.2299 x1 0.079151 0.061048 1.297 0.2013 x2 0.005260 0.006019 0.874 0.3867 x3 0.016144 0.006359 2.539 0.0146 * ---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 2.927 on 46 degrees of freedomMultiple R-squared: 0.462, Adjusted R-squared: 0.4269 F-statistic: 13.17 on 3 and 46 DF, p-value: 2.459e-06Now we can see that in the same size of predictor the 2R of LM5 is the largest and it equals 0.4531. And in the different size of predictor the 2adj R of LM5 is also the largest and it equals 0.4298. Also we know the model is131.2968380.096301*0.017364*Y X X =++(2). Use Hypothesis test to check the modelWe set 0212:0 :0H vs H ββ=≠. And the significant level 0.01α=>anova(LM7)Analysis of Variance TableResponse: yDf Sum Sq Mean Sq F value Pr(>F)x1 1 264.83 264.83 30.9064 1.318e-06 *** x2 1 18.44 18.44 2.1515 0.14924 x3 1 55.22 55.22 6.4450 0.01457 * Residuals 46 394.16 8.57 ---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >qf(1-0.01,p-1,n-p) [1] 5.087373Now we can see that since F(=2.1515)<qf(1-0.01,p-1,n-p)=5.087373 so we include that we fail to reject H0. That is 20β=.(3). Check the multicollinearity of the modelNow we want to check the multicollinearity of the model131.2968380.096301*0.017364*Y X X =++.>source("/~cspark/805/R/VIF.R") >vif(LM5)x1 x32.212312 2.212312We can see that since max(VIF1,VIF2)<10. So there is no multicollinearity in the model.(4). Check the Boxcox transformation.>source("/~cspark/805/R/bxcx.R")>source("/~cspark/805/R/BoxCox.R") #Box cox trans> y1=fitted(LM5)> e1=y-y1> plot(y,e1)> lam=seq(-1,1,0.1)> SSE=BoxCox(y~x1+x3,lambda=lam,SSE=TRUE)>cbind(lam,SSE)lam SSE[1,] -1.0 1726.2062[2,] -0.9 1412.3691[3,] -0.8 1170.1879[4,] -0.7 982.4830[5,] -0.6 836.4204[6,] -0.5 722.3825[7,] -0.4 633.1395[8,] -0.3 563.2405[9,] -0.2 508.5651[10,] -0.1 465.9918[11,] 0.0 433.1527[12,] 0.1 408.2517[13,] 0.2 389.9292[14,] 0.3 377.1613[15,] 0.4 369.1853[16,] 0.5 365.4434[17,] 0.6 365.5417[18,] 0.7 369.2195[19,] 0.8 376.3275[20,] 0.9 386.8112[21,] 1.0 400.7009Since when lam=0.5 then the SSE is the smallest one which equals 365.4434. So lam should be 0.5. And then check if lambda=0.5 is better.>yprime=sqrt(y)>LM.prime=lm(yprime~x1+x3) > plot(x1+x3,yprime)>plot(lam,SSE)1002003004005001.01.52.02.53.03.54.x1 + x3y p r i m e-1.0-0.50.00.5 1.0400600800100012001400160lamS S E>qqnorm(e1)Now from the pictures above we know that when lambda=0.5 is better.-2-1012-4-202468Normal Q-Q PlotTheoretical QuantilesS a m p l e Q u a n t i l e s(5). Matrix Plot>colnames(mydata)=c("Murder","RapeCrime","Robbery","Assault") >pairs(mydata, cex=0.5, pch=1)>cor(mydata)Murder Rape Crime Robbery Assault Murder 1.0000000 0.6012205 0.4837076 0.6485505 Rape Crime 0.6012205 1.0000000 0.5918793 0.7402595 Robbery 0.4837076 0.5918793 1.0000000 0.5570782 Assault 0.6485505 0.7402595 0.5570782 1.0000000We confirm the model again from the Matrix plotMurder102030405010020030040050051015102030405Rape CrimeRobbery010030051015100200300400500100300Assault(6). Check the constancy of error varianceWe assume that0:H the error variance is constant vs1:H the error variance is not constant>source("/~cspark/805/R/Breusch-Pagan.R") #constancy of error variance>>BP.test(y~x1+x3)$test.stat[1] 0.8323462$df[1] 2$p.value[1] 0.6595661>qchisq(1-0.01,2)[1] 9.21034Since 0.8323462<qchisq(1-0.01,2)=9.21034, so we reject H1. That is the error variance is constant.(7). Check the independence of the predictorsNow we want to check the independence of x1 and x3.> LM8=lm(y~x3+x1)>anova(LM8)Analysis of Variance TableResponse: yDf Sum Sq Mean Sq F value Pr(>F)x3 1 308.16 308.16 36.1458 2.58e-07 ***x1 1 23.78 23.78 2.7892 0.1015Residuals 47 400.70 8.53---Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1>anova(LM3) # We can see that SSR(x3)=SSR(x3|x1) So independent. Analysis of Variance TableResponse: yDf Sum Sq Mean Sq F value Pr(>F)x3 1 308.16 308.16 34.847 3.524e-07 *** Residuals 48 424.48 8.84 ---Signif.codes : 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >Since SSR(x3)=SSR(x3|x1), so the two predictors are independent.(8). Check the linearity of y and residual.> y1=fitted(LM5) > e1=y-y1 > plot(y,e1)51015-4-202468ye 1(9). Check if the residual follows the norm distribution> e2=resid(LM5) > c1=sort(e2)>LM.sum=summary(LM5) > s=LM.sum[["sigma"]] > k=1:n> c2=(k-0.375)/(n+0.25) > c3=qnorm(c2) > c4=s*c3 >plot(c4,c1) >lines(c4,c4)>cor=cor(c4,c1) >cor[1] 0.9860581We can see that the error terms almost normally distribute.-6-4-20246-4-202468c4c 1(10). Check if x1 and x3 are independent with residual.>plot(x1,e1)> plot(x3,e1)We can see that x1 or x3 are independent with residual.1020304050-4-202468x1e 1100200300400500-4-202468x3e 1(11). Check the median of the residuals and the maximum and minimum of the residuals.>boxplot(e1) #check the median of the residuals and the maximum and minimum of the residuals.We can obtain the median of the residuals and the maximum and minimum of the residuals from the picture.-4-202468。