# NOT RUN {
require(sp)
require(spdep)
require(rms)
data(meuse)
coordinates(meuse) <- ~x+y
meuse@data <- data.frame(DepVar=rbinom(dim(meuse)[1], 1, 0.5),
meuse@data)
#### Logistic model
lmodel <- logistic.regression(meuse, y='DepVar',
x=c('dist','cadmium','copper'))
lmodel$model
lmodel$diagTable
lmodel$coefTable
#### Logistic model with factorial variable
lmodel <- logistic.regression(meuse, y='DepVar',
x=c('dist','cadmium','copper', 'soil'))
lmodel$model
lmodel$diagTable
lmodel$coefTable
### Auto-logistic model using 'autocov_dist' in 'spdep' package
lmodel <- logistic.regression(meuse, y='DepVar',
x=c('dist','cadmium','copper'), autologistic=TRUE,
coords=coordinates(meuse), bw=5000)
lmodel$model
lmodel$diagTable
lmodel$coefTable
est <- predict(lmodel$model, type='fitted.ind')
#### Add residuals, standardized residuals and estimated probabilities
VarNames <- rownames(lmodel$model$var)[-1]
meuse@data$AutoCov <- lmodel$AutoCov
meuse@data <- data.frame(meuse@data, Residual=lmodel$Residuals[,1],
StdResid=lmodel$Residuals[,2], Probs=predict(lmodel$model,
meuse@data[,VarNames],type='fitted') )
#### Plot fit and probabilities
resid(lmodel$model, "partial", pl="loess")
# plot residuals
resid(lmodel$model, "partial", pl=TRUE)
# global test of goodness of fit
resid(lmodel$model, "gof")
# Approx. leave-out linear predictors
lp1 <- resid(lmodel$model, "lp1")
# Approx leave-out-1 deviance
-2 * sum(meuse@data$DepVar * lp1 + log(1-plogis(lp1)))
# plot estimated probabilities at points
spplot(meuse, c('Probs'))
# }
Run the code above in your browser using DataLab