# House data set #required R libraries spdep, Matrix # select systematic sample of size ns=1015 or ns=5072 ?? #choice=1: ns=1015 #choice=2: ns=5072 choice<-1 model<-"CAR" library(spdep) data(house) hform <- formula(log(price) ~ age + I(age^2) + I(age^3) + + log(lotsize) + rooms + log(TLA) + beds + syear) hmm0 <- model.matrix(hform, data = house) hlw <- nb2listw(LO_nb,style="B") #binary W<- as(hlw, "CsparseMatrix") W<-forceSymmetric(W) # as this is symmetrix hlw <- mat2listw(W,style="B") #make sure it is symmetric and binary X0<-hmm0 n<-dim(W)[1] Y<-log(house$price) n <- length(Y) isSymmetric(W) d1<-Diagonal(x=c(1/sqrt(rowSums(W)))) W0<-d1%*%W%*%d1 # normalise such that 1/min(eigenvalues)<=rho <=1 hsn0 <- listw2sn(hlw) hsn0[,3]<-W0@x hlw0<- sn2listw(hsn0) hlw0$style <- "B" if(choice==2){ s<- seq(1,n,by=5) # sample }else{ s<- seq(1,n,by=25) #sample } ns<-length(s) u<-setdiff(1:n,s) # non-sample units # full data set, all n=25357 observations # using spdep to fit full data set without missing observations rho.hat.se<-omega.hat.se<-rho.hat <- omega.hat<- times <-times<- iter<- L <- rep(NA,6) beta.hat.se<-beta.hat<-matrix(NA,6,13) st<-system.time( res.CAR <- spautolm(hform, data=house, listw=hlw0, interval = c(0,1), family="CAR", method="Matrix") ) summary(res.CAR) times[6]<-st[1] rho.hat[6] <- res.CAR$lambda # lambda 0.900805 rho.hat.se[6] <-res.CAR$lambda.se omega.hat[6]<-res.CAR$fit$s2 beta.hat[6,]<-res.CAR$fit$coefficients beta.hat.se[6,]<-sqrt(diag((res.CAR$fit$imat))) L[6]<- res.CAR$LL res.CAR$timings source("ARmodels.sparse.r") imax<-1e4 # extract sample information and re-order X and W, such that first ns units refer to sample (observed) # and remaing nu=n-ns units refer to non-sample (missing) Ys<-Y[s] X1<-X0[c(s,u),] W1<-W0[c(s,u),c(s,u)] # CAR tol<-1e-8 REML<-F st<-system.time({ lag.exact<-try(AR.pop.exact.sparse(Ys,W1,X1,model="CAR",se=F,approx.se=F,tol=tol,approx=F,REML=REML,pos.neg.sep=FALSE)) }) lag.exact$L # log-lik # estimates rho.hat[1] <- lag.exact$rho; # rho beta.hat[1,] <- lag.exact$beta; #beta omega.hat[1]<- lag.exact$sigma2; #sigma^2 times[1]<- st[1] # only use Running INLA, without pre and post-processing iter[1]<- 1 L[1]<-AR.exact.log.lik.fct(par=rho.hat[1],Ys,W1,X1,WplusWt[c(s,u),c(s,u)],WtW[c(s,u),c(s,u)],ARmod=3,pr=TRUE,REM=REML) #about 7 seconds hilf.exact<-try(AR.exact.Info(Ys,W1,X1,beta=lag.exact$beta,rho=lag.exact$rho,sigma2=lag.exact$sigma2,model=model,approx=F)) # standard errors based on Expected Info matrix beta.hat.se[1,]<-sqrt(diag(hilf.exact$Cov.beta)) #beta.se rho.hat.se[1]<-sqrt(hilf.exact$Cov.rho) #rho.se omega.hat.se[1] <- sqrt(hilf.exact$Cov.sigma2) #omega.se m<-1 t<-system.time({ r<-try(EM.algorithm.CAR(Ys,W=W1,X=X1,tol=tol,rho0=0,method=m,iter.max=imax,omega0=NULL,beta0=NULL,se=F) )}) rho.hat[m+1] <- r$rho # rho omega.hat[m+1] <- r$sigma2 # omega beta.hat[m+1,] <- r$beta #beta times[m+1]<- t[1] # only use Running INLA, without pre and post-processing iter[m+1]<- r$iter L[m+1]<-AR.exact.log.lik.fct(par=rho.hat[m+1],Ys,W1,X1,NULL,NULL,ARmod=ARmod,pr=TRUE,REM=REML) hilf.exact<-try(AR.exact.Info(Ys,W1,X1,beta=beta.hat[m+1,],rho=rho.hat[m+1],sigma2=omega.hat[m+1],model=model,approx=F)) # standard errors based on Expected Info matrix beta.hat.se[m+1,]<-sqrt(diag(hilf.exact$Cov.beta)) #beta.se rho.hat.se[m+1]<-sqrt(hilf.exact$Cov.rho) #rho.se omega.hat.se[m+1] <- sqrt(hilf.exact$Cov.sigma2) #omega.se # compare with INLA for CAR # to install INLA use: install.packages("INLA", repos="https://inla.r-inla-download.org/R/stable") library(INLA) zero.variance <- list(prec=list(initial = 25, fixed=TRUE)) # necessary for fixed variance? # Number of covariates (manually set for convenience) n_covariates <- 12 #Variance-covarinace matrix for beta coeffients' prior betaprec<-.0001 Q.beta = Diagonal(n=n_covariates, x=1) Q.beta = betaprec*Q.beta # Define the parameters for the SLM #args.slm = list( # rho.min = 0 , # rho.max = 1, # W = W0, # X = matrix(0, n, 0), # Q.beta = matrix(1,0,0) #) #Definition of priors for precision of the error effect and the spatial autocorrelation parameter (in the (0,1) interval). hyper.generic = list( prec = list( prior = "loggamma", param = c(0.01, 0.01)), rho = list(initial=0, prior = "logitbeta", param = rep(1,13)) ) ## Pretend we observe the first 50 house prices A <- sparseMatrix(i = 1:ns, j = 1:ns, x = 1, dims = c(ns,n)) ## Create a data stack containing the data and covariates stk <- inla.stack(data = list(Ys= Ys), A = list(A,1,1), effects=list(idx = 1:n, Intercept = rep(1,ns), X=X[1:ns,2:13]),tag='est') hilf<-eigen(W1,only.values=TRUE) car <- inla(Ys ~ -1 + Intercept + X + f(idx, model = "generic1", Cmatrix = W1), #hyper = hyper.generic), data=inla.stack.data(stk), family = "gaussian", control.family = list(hyper = zero.variance), control.predictor = list(A = inla.stack.A(stk))) summary(car) rho.hat[5] <- car$summary.hyperpar$mean[2] # rho omega.hat[5] <- 1/car$summary.hyperpar$mean[1] # omega beta.hat[5,] <- car$summary.fixed[,1] #beta times[5]<- car$cpu.used[2] # only use Running INLA, without pre and post-processing iter[5]<- NA L[5]<-AR.exact.log.lik.fct(par=rho.hat[5],Ys,W1,X1,WplusWt[c(s,r),c(s,r)],WtW[c(s,r),c(s,r)],ARmod=ARmod,pr=TRUE,REM=REML) beta.hat.se[5,]<- car$summary.fixed[,2] rho.hat.se[5]<- car$summary.hyperpar$sd[2] omega.hat.se[5]<- (1/car$summary.hyperpar$mean[1])^2*car$summary.hyperpar$sd[1] save(rho.hat,omega.hat,beta.hat,times,iter,L,beta.hat.se,rho.hat.se,omega.hat.se,file="Example.CAR.RData") # OUTPUT data example # method.names<-c("MML","EM","Martin84","GBH89","INLA","FULL") method.names <- c("Intercept","$age$","$age^2$","$age^3$","$\\log(lotsize)$","rooms","$\\log(LTA)$","beds", "syear1994","syear1995","syear1996","syear1997","syear1998","$\\rho$","$\\omega$") hilf<-cbind(t(beta.hat),t(beta.hat.se)) hilf<-rbind(hilf,c(rho.hat,rho.hat.se)) hilf<-rbind(hilf,c(omega.hat,omega.hat.se)) hilf<-hilf[,c(1,7,2,8,3,9,4,10,6,12,5,11)] output1(hilf,method.names,nsmall=3,digits=1) hilf<-rbind(L,iter,times) hilf<-hilf[,c(1:4,6,5)] output(hilf,NULL,nsmall=4)