﻿### Soc 6708 ### Advanced Statistics ### Robert Andersen, University of Toronto ### Wednesday, May 16, 2012 ###################### ## LINEAR MODELS ###################### ############################################ # 1. Fitting Linear Models ############################################ library(car) data(Prestige) attach(Prestige) # Fitting linear model mod1<-lm(prestige~income+education, data=Prestige) summary(mod1) ## model summary anova(mod1) ## Type I F-test Anova(mod1) ## Type II F-tests (compare to 'anova') ##From 'car' library # Adding categorical predictor 'type' of occupation mod2<-lm(prestige~income+education+type, data=Prestige) anova(mod1, mod2) ## Analysis of variance ## Check the dataset to learn why there is an error mod1b<-lm(prestige~income+education, data=na.omit(Prestige)) #'na.omit' removes missing observations mod2b<-lm(prestige~income+education+type, data=na.omit(Prestige)) # Should 'type' be included in the model? summary(mod2b) anova(mod1b, mod2b) Anova(mod2b) ## compare to 'anova (mod1b, mod2b)' AIC(mod2b) ################################################ # 2. Dealing with categorical predictors ################################################ #checking reference category for 'type' contrasts(type) # reference coded '0' # alternatively: levels(type) # first category is reference # changing the reference category to 'prof' contrasts(type)<-contr.treatment(levels(type), base=2) # second row contrasts(type) # alternatively: type2<-relevel(type, 'prof') contrasts(type2) ########################################################## # 3. Quasi-variances ######################################################### library(qvcalc) qvtype<-qvcalc(mod2b, "type") summary(qvtype) ## Is there a problem with the quasi-variances'? ## see '?qvcalc' for details ?qvcalc ## Refitting model wihout 'income': mod3<-lm(prestige~type, data=na.omit(Prestige)) qvtype.mod3<-qvcalc(mod3, "type") plot(qvtype.mod3) # plotting 95% confidence intervals for each category # based on quasi-variances #Making function for plot of regular standard errors: plot.se<-function (x, width = 2, ylab = "estimate", main = "Intervals based on standard errors", ...) { frame <- x\$qvframe if (is.na(frame\$estimate[1])) stop("x has no parameter estimates to plot") if (any(is.nan(frame\$quasiSE))) levels <- factor(row.names(frame), levels = row.names(frame)) xvalues <- seq(along = levels) est <- frame\$estimate se <- frame\$SE tops <- est + width * se tails <- est - width * se range <- max(tops) - min(tails) ylim <- c(min(tails) - range/10, max(tops) + range/10) plot(levels, frame\$estimate, border = "transparent", ylim = ylim, xlab = x\$factorname, ylab = ylab, main = main, ...) points(frame\$estimate, ...) segments(xvalues, tails, xvalues, tops) } #Plotting graphs of SE ad Quasi SEs side-by-side split.screen(figs=c(1,2)) # an alternative to 'par(mfrow=c(1,2))' screen(1) plot.se(qvtype.mod3, cex=1.5, main="Regular SEs") screen(2) plot(qvtype.mod3, cex=1.5, main="Quasi SEs") close.screen(all=TRUE) ########################################## # 4. Simple nonlinear relationships ########################################## # Fitting model with quadratic effect for income # creating income-squared term: Prestige\$income2<-Prestige\$income^2 Prestige ## Fitting model mod.poly<-lm(prestige~education+income+income2, data=Prestige) summary(mod.poly) ## Compare model 'mod.poly' to model with orthogonal polynomial: mod.poly2<-lm(prestige~education+poly(income,2), data=Prestige) summary(mod.poly) summary(mod.poly2) # How and why do the models differ? # How and why are they similar? mod.inc2<-lm(prestige~education+income2, data=Prestige) summary(mod.inc2) #Does this model differ from the quadratic polynomial? summary(mod.poly) #################################### # 5. Effect displays #################################### library(effects) ?effect # income effect for quadratic polynomial: inc.eff<-effect("poly(income,2)", mod.poly2) inc.eff # see what is in the 'effect' object summary(inc.eff) plot(inc.eff, main="Income effect", xlab="Income", ylab="Fitted Prestige score") ######Creating my own plot from the effects # a) Creating data.frame from effects output: inc.eff2<-data.frame(inc.eff) inc.eff2 # b) Extracting the fitted values and 95% confidence bands income.eff<-inc.eff2[1:nrow(inc.eff2), 1] fit.prestige<-inc.eff2[1:nrow(inc.eff2), 2] lower.prestige<-inc.eff2[1:nrow(inc.eff2), 4] upper.prestige<-inc.eff2[1:nrow(inc.eff2), 5] # c) making the plot plot(income.eff, fit.prestige, main="Effect display for income", xlab="Income", ylab="Fitted prestige", type="l", lwd=2, ylim=c(20, 75)) lines(income.eff, lower.prestige, lty=3) lines(income.eff, upper.prestige, lty=3) ### Creating an Effect Display using the 'predict' function ?predict income<-seq(0,26000, length=10) length(income) education<-seq(0,16, length=10) length(education) new.data<-data.frame(income, education) fitted.data<-data.frame(income, education, predict(mod.poly2, new.data, se.fit=TRUE)) fitted.data plot(fitted.data\$income, fitted.data\$fit, type="l", main="Fitted Prestige by Income", xlab="Income", ylab="Fitted prestige") ############################################# # 6. Standardization and Relative importance ############################################# #standardizing 'income' and 'education' using scale function: mod1<-lm(prestige~scale(income)+scale(education)+type, data=Prestige) summary(mod1) # relative importance of several terms: mod1b<-lm(prestige~income+education+type, data=Prestige) summary(mod1b) library(relimp) # comparing 'income' (2nd terms in summary) # and occupation type (4th and 5th terms in summary) relimp(mod1b, set1=2, set2=c(4,5)) ############################################## # 7. Interactions ############################################## mod2a<-lm(prestige~income*type, data=Prestige) ##alternatively mod2b<-lm(prestige~income+type+income:type, data=Prestige) Anova(mod2b) library(effects) plot(allEffects(mod2b)) ## type: ?AllEffects for help