# Computer code for example 3.5
# HINT:
# The program assumes that all files are stored in
# C:\monographregression\computercode
# and corresponding subfolders.
# Change directory if code is located elsewhere.
# change working directory
setwd("c:/monographregression/computercode")
library(foreign)
# read data
rent99 <- read.dta(file="data/stata/rent99.dta")
# estimate the model rentsqm = beta_0 + beta_1*(1/area) + beta_2*yearc + epsilon
# note that the results printed in the book are obatined with STATA. Due to rounding errors the results obtained with R are slightly different
mod1 <- lm(rentsqm ~ I(1/area) + yearc, data = rent99)
summary(mod1)
# estimate the model rentsqm = beta_0 + beta_1*(1/area) + beta_2*yearc + beta_3*yearc^2 + beta_4*yearc^3 + epsilon
# note that the results printed in the book are obatined with STATA. Due to rounding errors the results obtained with R are slightly different
mod2 <- lm(rentsqm ~ I(1/area) + yearc + I(yearc^2) + I(yearc^3), data = rent99)
summary(mod2)
# estimate the model rentsqm = beta_0 + beta_1*(1/area) + beta_2*yearc + beta_3*yearc^2 + beta_4*yearc^3 + epsilon
# where the value 1900 is substracted form yearc
# note that the results printed in the book are obatined with STATA. Due to rounding errors the results obtained with R are slightly different
mod3 <- lm(rentsqm ~ I(1/area) + I(yearc-1900) + I((yearc-1900)^2) + I((yearc-1900)^3), data = rent99)
summary(mod3)
# plot the effect of yearc
curve(coef(mod2)[3]*x + coef(mod2)[4]*x^2 +coef(mod2)[5]*x^3 , 1915, 2000,
main = "a) effect of year of construction, coded from 18 to 98",
xlab = "year of construction", ylab = "effect of year of construction")
curve(coef(mod3)[3]*x + coef(mod3)[4]*x^2 +coef(mod3)[5]*x^3 , 15, 100,
main = "b) effect of year of construction, coded from 18 to 98",
xlab = "year of construction", ylab = "effect of year of construction")
# estimate the model rentsqm = beta_0 + beta_1*(1/area) + beta_2*yearc + beta_3*yearc^2 + beta_4*yearc^3 + epsilon
# using 1/area centered and orthogonal polynomials for yearc
# note that the results printed in the book are obatined with STATA. In contrast to R STATA uses orthoNORMAL polynomials. Therefore
# the estimated coefficients are different to the printed results.
mod4 <- lm(rentsqm ~ I(1/area - mean(1/area)) + poly(yearc,3) , data = rent99)
summary(mod4)
# plot estimated curves including partial residuals
pred <- predict(mod4, type="terms")
resid <- residuals(mod4, type="partial")
plot(rent99$yearc,resid[,2],ylab="effect year of construction/partial residuals",xlab="year of construction")
help <- data.frame(yearc=rent99$yearc,fyear = pred[,2])
o<-order(help$yearc)
lines(help[o,1],help[o,2])