# 
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

library(alr3)
data(UN1)
attach(UN1)

# 2.6.1: (1.3.3 -> taking logarithms makes a linear regression makes much more reasonble in this case:)
#
LPPgdb     = log(PPgdp,base=10)
LFertility = log(Fertility,base=10)

plot( LPPgdb, LFertility )

m <- lm( LFertility ~ LPPgdb ) 

summary(m)
anova(m)

# 2.6.2:
#
postscript("../../WriteUp/Graphics/Chapter2/prob_6_UN.eps", onefile=FALSE, horizontal=FALSE)
plot( LPPgdb, LFertility, xlab="LPPgdb", ylab="LFertility" ) 
abline( m ) 
dev.off()

# 2.6.3: test for significance of \beta_1 \ne 0: 
#
t = -0.22116 / 0.01737
n = length(LPPgdb)
pt(t,n-2)

RSS       = sum( m$residuals^2 )
n         = length(m$residuals)
sigma_hat = sqrt( RSS / (n-2) ) 
SYY       = sum( ( LFertility - mean(LFertility) )^2 ) 

# compute the F statistic: 
F = ( ( SYY - RSS )/1 ) / ( sigma_hat^2 ) 

1 - pf(F,1,n-2)

# 2.6.4: summary(m) gives
#
# R^2 = 0.4591
#

# 2.6.5: WWX: How do?
#

# 2.6.6:
#
l10f <- predict( m, newdata=data.frame( LPPgdb=log(1000,base=10) ), interval="prediction", level=0.95 )

# map these endpoint variables into the true domain of interest:
#
c( 10^( l10f[2] ), 10^( l10f[1] ), 10^( l10f[3] ) ) 

# 2.6.7:
#
# find the location of maximum fertility (Niger):
#
max_fert <- max(UN1$Fertility)
inds <- UN1$Fertility == max_fert
UN1$Locality[inds]

# find the location of minimum fertility (Hong.Kong):
#
min_fert <- min(UN1$Fertility)
inds <- UN1$Fertility == min_fert
UN1$Locality[inds]

# find the location of the two smallest residuals ... 
#
# indices: 6, 181
#
sr <- sort( m$residuals )
sr[1:2]

UN1$Locality[6]   # Armenia, 
UN1$Locality[181] # Ukraine

# find the location of the two largest residuals ...  
#
# indices: 55, 129 
#
sr <- sort( m$residuals, decreasing=TRUE )
sr[1:2]

UN1$Locality[55]  # Equatorial.Guinea,
UN1$Locality[129] # Oman