# These self contained codes illustrate most of the estimation procedures for # Kostov, P. (2010) Model boosting for spatial weighting matrix selection in spatial lag models, # Environment and Planning B: Planning and Design, 37(3), 533-549. # the R packages spdep and mboost are required # The original implementation used older version of mboost, so it may need some changes # to work with more up-to-date versions of the package # I have similar codes for robust, as well as quantile regression estimation along the same lines #load spdep package library (spdep) data(boston) attach (boston.c) my=log(MEDV) mx=cbind(CRIM , ZN , INDUS , CHAS , I(NOX^2) , I(RM^2) , AGE , log(DIS) , log(RAD) , TAX , PTRATIO , B , log(LSTAT)) dimnames(mx)[[2]] <- c("CRIM", "ZN" , "INDUS" , "CHAS" , "NOX^2" , "RM^2" , "AGE" , "log(DIS)", "log(RAD)" , "TAX" , "PTRATIO", "B", "log(LSTAT)" ) coords=cbind(LON,LAT) detach(boston.c) ############################################################################### get.list=function(coords,k,ws) { my.nb=knn2nb(knearneigh(coords, k=k),sym=T ) dlist <- nbdists(my.nb, coords) dlist <- lapply(dlist, ws) wlist<-nb2listw(my.nb, glist=dlist, style="W", zero.policy=T) wlist } create.instr=function(my,mx,wlist){ ly= lag.listw(wlist,my) res=matrix(data = NA, nrow = nrow(mx), ncol = ncol(mx)) for (i in 1:ncol(mx)) { res[ ,i]= lag.listw(wlist,mx[,i]) } instr=lm(ly~ res)$fitted.values instr } ############################################################################### degs=seq(from = 0.4, to = 4,by=0.1) funs=as.list(rep(NA,length(degs))) for(i in 1:length(degs)) { funs[[i]]=function(x) 1/(x^degs[i]) } res=matrix(NA,length(my),length(funs)) X=mx #set the initial such matrix #put a smalll number n=51 for(J in 1:(n-1)) { for(i in 1:length(funs)) { a1 <-funs[[i]] a<- get.list(coords,k=J,a1) b<-create.instr(my,mx,a) names(b)=NULL res[,i]=b } mn=paste("n",J,"w",degs,sep="") res=data.frame(res) names(res)= mn res=as.matrix(res) X=cbind(X,res) } ############################################################################### n =nrow(X) save(X,my,n, file="e:/boston.w.Rdata") ############################################################################### ############################################################################### #load("e:/boston.w.Rdata") #prepare folds here n =nrow(X) k = 10 ntest <- floor(n / k) cv10f <- matrix(c(rep(c(rep(0, ntest), rep(1, n)), k - 1), rep(0, n * k - (k - 1) * (n + ntest))), nrow = n) k = 8 ntest <- floor(n / k) cv8f <- matrix(c(rep(c(rep(0, ntest), rep(1, n)), k - 1), rep(0, n * k - (k - 1) * (n + ntest))), nrow = n) set.seed(290875) bs25 <- rmultinom(25, n, rep(1, n)/n) bs100 <- rmultinom(100, n, rep(1, n)/n) library(mboost) ######################################################################### m1=glmboost(X,my, control = boost_control(mstop = 10000,nu=0.3)) aic <- AIC(m1,method="corrected") aic mbest=m1[mstop(aic)] names(coef(mbest)[abs(coef(mbest)) > 0]) gMDL <- AIC(m1,method="gMDL") gMDL mbest=m1[mstop(gMDL)] names(coef(mbest)[abs(coef(mbest)) > 0]) sel.names=names(coef(mbest)[abs(coef(mbest)) > 0]) cvm <- cvrisk(m1, folds = cv10f) mstop(cvm) mbest= m1[mstop(cvm)] names(coef(mbest)[abs(coef(mbest)) > 0]) cvm <- cvrisk(m1, folds = cv8f) mstop(cvm) mbest= m1[mstop(cvm)] names(coef(mbest)[abs(coef(mbest)) > 0]) ### 25 bootstrap iterations cvm <- cvrisk(m1, folds = bs25) mstop(cvm) mbest= m1[mstop(cvm)] names(coef(mbest)[abs(coef(mbest)) > 0]) ### 100 bootstrap iterations cvm <- cvrisk(m1, folds = bs100) mstop(cvm) mbest= m1[mstop(cvm)] names(coef(mbest)[abs(coef(mbest)) > 0]) ##############################################################################