{"word_count": 535, "character_count": 5710}
Visualization is a key to success in regression analysis. This is one of the (many) reasons I am also suspicious when I read an article with a quantitative (econometric) analysis without any graph. Consider for instance the following dataset, obtained from http://data.worldbank.org/, with, for each country, the GDP per capita (in some common currency) and the infant mortality rate (deaths before the age of 5),&#62; library(gdata)&#62; XLS1=read.xls("http://api.worldbank.org/datafiles/NY.GDP.PCAP.PP.CD_Indicator_MetaData_en_EXCEL.xls", sheet = 1)&#62; data1=XLS1[-(1:28),c("Country.Name","Country.Code","X2010")]&#62; names(data1)[3]="GDP"&#62; XLS2=read.xls("http://api.worldbank.org/datafiles/SH.DYN.MORT_Indicator_MetaData_en_EXCEL.xls", sheet = 1)&#62; data2=XLS2[-(1:28),c("Country.Code","X2010")]&#62; names(data2)[2]="MORTALITY"&#62; data=merge(data1,data2)&#62; head(data)Country.Code         Country.Name       GDP MORTALITY1          ABW                Aruba        NA        NA2          AFG          Afghanistan  1207.278     149.23          AGO               Angola  6119.930     160.54          ALB              Albania  8817.009      18.45          AND              Andorra        NA       3.86          ARE United Arab Emirates 47215.315       7.1If we estimate a simple linear regression -   - we get&#62; regBB=lm(MORTALITY~GDP,data=data)&#62; summary(regBB)Call:lm(formula = MORTALITY ~ GDP, data = data)Residuals:Min     1Q Median     3Q    Max-45.24 -29.58 -12.12  16.19 115.83Coefficients:Estimate Std. Error t value Pr(&#62;|t|)(Intercept) 67.1008781  4.1577411  16.139  &#60; 2e-16 ***GDP         -0.0017887  0.0002161  -8.278 3.83e-14 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 39.99 on 167 degrees of freedom(47 observations deleted due to missingness)Multiple R-squared: 0.2909,Adjusted R-squared: 0.2867F-statistic: 68.53 on 1 and 167 DF,  p-value: 3.834e-14We can look at the scatter plot, including the linear regression line, and some confidence bounds,&#62; plot(data$GDP,data$MORTALITY,xlab="GDP per capita",+ ylab="Mortality rate (under 5)",cex=.5)&#62; text(data$GDP,data$MORTALITY,data$Country.Name,pos=3)&#62; x=seq(-10000,100000,length=101)&#62; y=predict(regBB,newdata=data.frame(GDP=x),+ interval="prediction",level = 0.9)&#62; lines(x,y[,1],col="red")&#62; lines(x,y[,2],col="red",lty=2)&#62; lines(x,y[,3],col="red",lty=2)We should be able to do a better job here. For instance, if we look at the Box-Cox profile likelihood,&#62; boxcox(regBB)it looks like taking the logarithm of the mortality rate should be better, i.e.  or :&#62; regLB=lm(log(MORTALITY)~GDP,data=data)&#62; summary(regLB)Call:lm(formula = log(MORTALITY) ~ GDP, data = data)Residuals:Min      1Q  Median      3Q     Max-1.3035 -0.5837 -0.1138  0.5597  3.0583Coefficients:Estimate Std. Error t value Pr(&#62;|t|)(Intercept)  3.989e+00  7.970e-02   50.05   &#60;2e-16 ***GDP         -6.487e-05  4.142e-06  -15.66   &#60;2e-16 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.7666 on 167 degrees of freedom(47 observations deleted due to missingness)Multiple R-squared: 0.5949,Adjusted R-squared: 0.5925F-statistic: 245.3 on 1 and 167 DF,  p-value: &#60; 2.2e-16&#62; plot(data$GDP,data$MORTALITY,xlab="GDP per capita",+ ylab="Mortality rate (under 5) log scale",cex=.5,log="y")&#62; text(data$GDP,data$MORTALITY,data$Country.Name)&#62; x=seq(300,100000,length=101)&#62; y=exp(predict(regLB,newdata=data.frame(GDP=x)))*+ exp(summary(regLB)$sigma^2/2)&#62; lines(x,y,col="red")&#62; y=qlnorm(.95, meanlog=predict(regLB,newdata=data.frame(GDP=x)),+ sdlog=summary(regLB)$sigma^2)&#62; lines(x,y,col="red",lty=2)&#62; y=qlnorm(.05, meanlog=predict(regLB,newdata=data.frame(GDP=x)),+ sdlog=summary(regLB)$sigma^2)&#62; lines(x,y,col="red",lty=2)on the log scale or&#62; plot(data$GDP,data$MORTALITY,xlab="GDP per capita",+ ylab="Mortality rate (under 5) log scale",cex=.5)on the standard scale. Here we use quantiles of the log-normal distribution to derive confidence intervals.But why shouldn't we take also the logarithm of the GDP ? We can fit a model  or equivalently .&#62; regLL=lm(log(MORTALITY)~log(GDP),data=data)&#62; summary(regLL)Call:lm(formula = log(MORTALITY) ~ log(GDP), data = data)Residuals:Min       1Q   Median       3Q      Max-1.13200 -0.38326 -0.07127  0.26610  3.02212Coefficients:Estimate Std. Error t value Pr(&#62;|t|)(Intercept) 10.50192    0.31556   33.28   &#60;2e-16 ***log(GDP)    -0.83496    0.03548  -23.54   &#60;2e-16 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.5797 on 167 degrees of freedom(47 observations deleted due to missingness)Multiple R-squared: 0.7684,Adjusted R-squared: 0.767F-statistic:   554 on 1 and 167 DF,  p-value: &#60; 2.2e-16&#62; plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",+ ylab="Mortality rate (under 5)",cex=.5,log="xy")&#62; text(data$GDP,data$MORTALITY,data$Country.Name)&#62; x=exp(seq(1,12,by=.1))&#62; y=exp(predict(regLL,newdata=data.frame(GDP=x)))*+ exp(summary(regLL)$sigma^2/2)&#62; lines(x,y,col="red")&#62; y=qlnorm(.95, meanlog=predict(regLL,newdata=data.frame(GDP=x)),+ sdlog=summary(regLL)$sigma^2)&#62; lines(x,y,col="red",lty=2)&#62; y=qlnorm(.05, meanlog=predict(regLL,newdata=data.frame(GDP=x)),+ sdlog=summary(regLL)$sigma^2)&#62; lines(x,y,col="red",lty=2)on the log scales or&#62; plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",+ ylab="Mortality rate (under 5)",cex=.5)on the standard scale. If we compare the last two predictions, we havewith in blue is the log model, and in red is the log-log model (I did not include the first one for obvious reasons).