# Start the library library(GLDreg) ?GLDreg # Engel dataset example data(engel) set.seed(100) engel.fit.all<-GLD.lm.full(foodexp~income,data=engel,param="fmkl",fun=fun.RMFMKL.ml.m,summary.plot=F) result<-GLD.quantreg(seq(0.1,.9,length=9),engel.fit.all,intercept="fixed") sapply(1:9,function(i) sum(engel\$foodexp-cbind(1,engel\$income)%*%(result[1:2,i])<0)/nrow(engel)*100) result.emp<-GLD.quantreg(seq(0.1,.9,length=9),engel.fit.all,intercept="fixed",emp=T) sapply(1:9,function(i) sum(engel\$foodexp-cbind(1,engel\$income)%*%(result.emp[1:2,i])<0)/nrow(engel)*100) fun.plot.q(x=engel\$income,y=engel\$foodexp,fit=engel.fit.all[[1]],result,xlab="Income",ylab="Food Expense") summaryGraphics.gld.lm(engel.fit.all) library(quantreg) library(MASS) par(mfrow=c(1,1)) plot(foodexp~income,data=engel,xlab="income",ylab="food expense",main="Belgian Engel Dataset") abline(lm(foodexp~income,data=engel),lty=2,lwd=3) abline(rlm(foodexp~income,data=engel),lty=4,lwd=3) abline(rq(foodexp~income,data=engel),lty=3,lwd=2) abline(engel.fit.all[[1]][[3]][1:2],lty=1,lwd=3) legend("bottomright",c("Standard Regression","Robust Regression","Quantile Regression","GLD Regression"), lty=c(2,4,3,1),lwd=c(3,3,2,3)) # mcycle data example library(MASS) library(splines) cutoff<-15 mcycle.part1<-mcycle[mcycle\$times=cutoff,] mcycle.p1.fit<-GLD.lm.full(accel~bs(times,df=8),data=mcycle.part1,param="rs",fun=fun.RPRS.ml.m,summary.plot=F,n=100) mcycle.p2.fit<-GLD.lm.full(accel~bs(times,df=15),data=mcycle.part2,param="fmkl",fun=fun.RMFMKL.ml.m,summary.plot=F,n=100) mcycle.p1.fit.quant<-GLD.quantreg(seq(0.1,.9,length=9),mcycle.p1.fit,slope="fixed") mcycle.p2.fit.quant<-GLD.quantreg(seq(0.1,.9,length=9),mcycle.p2.fit,slope="fixed") plot(mcycle,ylim=c(-150,75)) sapply(1:9,function(i) lines(mcycle\$times,c( cbind(1,bs(mcycle[which(mcycle\$times=cutoff),]\$times,df=15))%*%mcycle.p2.fit.quant[1:16,i]))) # crime data example ### RS GLD ### require(foreign) require(MASS) require(GLDreg) require(quantreg) cdata <- read.dta("http://www.ats.ucla.edu/stat/data/crime.dta") mcdata<-subset(cdata,select=c("crime","single")) mcdata[51,1]<-10 mcdata[25,2]<-22 par(mfrow=c(2,2)) plot(crime~single,data=mcdata,ylab="Violent crimes per 100,000 people", xlab="Percentage of population that are single parents",main="GLD Regression fails due to improper initial value selection (A)") abline(rlm(crime~single,data=mcdata),lty=2,lwd=3) abline(lm(crime~single,data=mcdata),lty=4,lwd=3) abline(rq(crime~single,data=mcdata),lty=3,lwd=2) abline(GLD.lm(crime~single,data=mcdata, param="rs",fun=fun.RPRS.ml.m,diagnostics=FALSE)[[3]][1:2],lty=1,lwd=3) legend("topright",c("Robust Regression","Standard Regression","Quantile Regression","GLD Regression"), lty=c(2,4,3,1),lwd=c(3,3,2,3),bg="white") init1<-rq(crime~single,data=mcdata)\$coeff init2<-rlm(crime~single,data=mcdata)\$coeff init3<-GLD.lm(crime~single,data=rbind(mcdata[1:24,],mcdata[26:50,]), param="rs",fun=fun.RPRS.ml.m,diagnostics=FALSE)[[3]][1:2] init4<-GLD.lm(crime~single,data=rbind(mcdata[1:24,],mcdata[26:50,],mcdata), param="rs",fun=fun.RPRS.ml.m,diagnostics=FALSE)[[3]][1:2] plot(crime~single,data=mcdata,ylab="Violent crimes per 100,000 people", xlab="Percentage of population that are single parents",main="GLD Regression with quantile regression coefficients as initial values (B)") abline(rlm(crime~single,data=mcdata),lty=2,lwd=3) abline(lm(crime~single,data=mcdata),lty=4,lwd=3) abline(rq(crime~single,data=mcdata),lty=3,lwd=2) abline(GLD.lm(crime~single,data=mcdata, param="rs",fun=fun.RPRS.ml.m,diagnostics=FALSE,init=init1)[[3]][1:2],lty=1,lwd=3) plot(crime~single,data=mcdata,ylab="Violent crimes per 100,000 people", xlab="Percentage of population that are single parents",main="GLD Regression with robust regression coefficients as initial values (C)") abline(rlm(crime~single,data=mcdata),lty=2,lwd=3) abline(lm(crime~single,data=mcdata),lty=4,lwd=3) abline(rq(crime~single,data=mcdata),lty=3,lwd=2) abline(GLD.lm(crime~single,data=mcdata, param="rs",fun=fun.RPRS.ml.m,diagnostics=FALSE,init=init2)[[3]][1:2],lty=1,lwd=3) plot(crime~single,data=mcdata,ylab="Violent crimes per 100,000 people", xlab="Percentage of population that are single parents",main="GLD Regression using modified data fitted by GLD regression as initial values (D)") abline(rlm(crime~single,data=mcdata),lty=2,lwd=3) abline(lm(crime~single,data=mcdata),lty=4,lwd=3) abline(rq(crime~single,data=mcdata),lty=3,lwd=2) abline(GLD.lm(crime~single,data=mcdata, param="rs",fun=fun.RPRS.ml.m,diagnostics=FALSE,init=init3)[[3]][1:2],lty=1,lwd=3) # Occupation example library(MASS) library(car) data(Duncan) job.fit.full<-GLD.lm.full(income~education+prestige, data=Duncan,param="fkml",fun=fun.RMFMKL.lm, init=rlm(income~education+prestige,data=Duncan)\$coeff, summary.plot=F) summaryGraphics.gld.lm(job.fit.full) plot(job.fit.full[[1]]\$Fitted,job.fit.full[[1]]\$y,xlab="Fitted",ylab="Actual") abline(0,1)