ISM 602 Week 6 Outline

• Review
  • assumptions of linear regression
  • how to interpret regression results of summary() and plot
• Required Packages and Data: MASS, caret, homes.csv
• Three Unifying Features of Generalized Linear Models (slides)
• Which model works for predicting which business variables (slides)
• Assumptions of gamma regression: gamma response, independent cases, constant dispersion, and no collinearity 
• Comparisons of OLS and MLE
• Concept of MLE:
  • computing likelihood and log likelihood of gamma variable observations y = c(10, 15, 3, 19) or y = homes$totalvalue
  • comparing which model parameters improve the likelihood alpha=1, beta=1, alpha=1, beta = 2, ...
• Normality Check and Transformation:
  • check: hist, qqnorm, qqline
  • Use boxcox() in MASS to find lambda: boxcox(lm(totalvalue ~ finsqft, homes))
  • box-cox transformation and inverse box-cox transformation (create custom functions)

BoxCox= function(z, lambda) {
if (lambda == 0) {
    return(log(z))
} else {
    return((z^lambda - 1)/ lambda)
}
}

invBoxCox= function(z, lambda) {
if (lambda == 0) {
    return(exp(z))
} else {
    return((lambda * z + 1)^(1 / lambda))
}
}

• Build prediction models
  
  • rescale total value and finished square feet by dividing by the max value
  • transform scaled total value using box-cox: homes$totalValue.bc = BoxCox(homes$totalvalue, lambda)
  • split dataset into train and test sets
  • m1 = lm(scaledtotalvalue ~ scaledsqft, train) 
  • m2 = lm(log(totalvalue) ~ scaledsqft, train)
  • m3 = lm(totalValue.bc ~ scaledsqft, train)
  • m4 = glm(totalvalue ~ scaledsqft, train, family = Gamma(link="log"))

• Prediction of rescaled responses: 

        m1 = lm(scaledtotalvalue ~ scaledsqft, train)
        prediction1 = predict(m1, test)
        RMSE(prediction1, test$scaledtotalvalue)

        #against acctual value
        prediction1.ac = prediction1 *maxvalue
        plot(homes$finsqft, homes$totalvalue)
        curve(predict(m1, data.frame(scaledsqft=x/maxsqft))*maxvalue, col="red", add=TRUE) 

• Prediction of log-transformed responses:

m2 = lm(log(totalvalue) ~ scaledsqft, train)
summary(m2)
plot(m2)

#prediction must consider Duan's smearing factor or the predicted value is median, not mean
smear = mean(exp(residuals(m2)))
prediction2 = exp(predict(m2, test)) * smear
RMSE(prediction2, test$totalvalue) #in different unit

plot(homes$finsqft, homes$totalvalue)
curve(exp(predict(m2, data.frame(scaledsqft = x/maxsqft)))*smear, col="red", add=TRUE)

• Prediction of box-cox-transformed response using linear regression:

m3 = lm(totalvalue.bc ~ scaledsqft, train)
summary(m3)
plot(m3)
predicion3 = predict(m3, test)
RMSE(prediction3, test$totalvalue.bc)

# Duan's smearing factor is more complicated
#for each predcited value z, you apply inverse boxcox to z + residual of each training case and compute the average
#predict total value for 2000 sqft home

prediction.2000 = predict(m3, data.frame(scaledsqft = 2000/maxsqft))
#smearing factored prediction
prediction.smear = mean(invBoxCox(prediction.2000 + residuals(m3), lambda))

#get all predicted values for test data set
prediction_plus_e = outer(predict(m3, test), residuals(m3), "+")
prediction_plus_e = invBoxCox(prediction_plus_e, lambda)
prediction3 = rowMeans(prediction_plus_e)
RMSE(prediction3, test$otalvalue)

#create a function to smear box-cox transformed predictions

predict.smear = function(x) {
    prediction_plus_e = outer(predict(m3, data.frame(scaledsqft = x/maxsqft)), residuals(m3), "+")
    prediction_plus_e = invBoxCox(prediction_plus_e, lambda)
    prediction = rowMeans(prediction_plus_e)
    return(prediction)
}

plot(homes$finsqft, homes$totalvalue)
curve(predict.smear(x), col="red", add=TRUE)

• Prediction of un-engineered responses using gamma regression:

m4 = glm(totalvalue ~ finsqft, family = Gamma(link = "log"), data = train)
summary(m4)

prediction4 = predict(m4, test, type="response")
RMSE(prediction, test$totalvalue)

plot(homes$finsqft, homes$totalvalue)
curve(predict(m4, data.frame(finsqft=x), type="response"), col="red", add=TRUE)

• Interpret gamma regression result:
  • regression line: log(y) = a + bx and t-tests results
  • exp(b) - 1 is the percentage of y increase for each unit increase of x if b > 0
  • deviance, residual deviance, deviance residuals
• goodness of model: 
  • deviance (= residual deviance): D = - 2 × log-likelihood of fitted model + 2 × log-likelihood of saturated model --- the two times of the likelihood difference between your model and the saturated model (similar to residual standard error in linear regression)
  • deviance residuals: the contribution of each case to the the residual deviance: sum of squared deviance residuals = residual deviance (similar to residuals in linear regression)
  • R-squared, chi-squared, AIC = -2 × log-likelihood of fitted model + 2k, where k is the number of estimated parameters
• Homework:
  • Reading: Lecture Note 13 -- Gamma Regressions
  • Writing: 
    • Multiple Choice Questions (on ecourse.org): due before the next week
    • Hands-on: Q10 of Lecture Note 13: due before the next week