## -----------------------------------------------------------------------------
# independent variable
ind <- c(87, 95, 72, 37, 44, 24, 
         40, 55, 55, 38, 88, 34,
         41, 30, 26, 35, 38, 24,
         14, 56, 37, 34,  8, 18, 
         49, 44, 51, 67, 17, 37, 
         55, 25, 33, 32, 59, 54)

# dependent variable
dep <- c(72, 75, 85, 29, 58, 30, 
         50, 60, 49, 46, 84, 23, 
         21, 46, 22, 42, 45, 14,
         19, 36, 48, 23,  8, 29, 
         38, 47, 52, 52, 22, 48, 
         58, 40, 46, 38, 35, 55)


## ----ch2-1--------------------------------------------------------------------
plot(ind, dep)


## -----------------------------------------------------------------------------
m <- glm(dep ~ ind)


## -----------------------------------------------------------------------------
m


## -----------------------------------------------------------------------------
s <- summary(m)
s


## ----ch2-2--------------------------------------------------------------------
plot(ind, dep)
abline(m)


## ----ch2-3--------------------------------------------------------------------
plot(ind, dep, pch=18, xlim=c(0,100), ylim=c(0,100), 
    axes=FALSE, xlab='', ylab='', yaxs="i", xaxs="i")
axis(1, at=(0:5)*20)
axis(2, at=(0:5)*20, las=1)

# create regression formula
f <- paste0('y = ', round(s$coefficients[2], 4), 'x + ', round(s$coefficients[1], 4))
# add the text in variable f to the plot 
text(0, 96, f, pos=4)
# compute r-squared
R2 <- cor(dep, predict(m))^2

# set up the expression (a bit complex, this)
r2 <- bquote(italic(R)^2 == .(round(R2, 4)))
# and add it to the plot
text(0, 85, r2, pos=4)

# compute regression line
# create a data.frame with the range (minimum and maximum) of values of ind
px <- data.frame(ind = range(ind))
# use the regression model to get predicted value for dep at these two extremes
py <- predict(m, px)
# combine the min and max values and the predicted values
ln <- cbind(px, py)
# add to the plot as a line
lines(ln, lwd=2)


## -----------------------------------------------------------------------------
mi <- matrix(ind, ncol=6, nrow=6, byrow=TRUE)
md <- matrix(dep, ncol=6, nrow=6, byrow=TRUE)


## ----ch2-4a-------------------------------------------------------------------
# load package
library(terra)

# turn matrices into SpatRaster objects
ri <- rast(mi)
rd <- rast(md)


## ----ch2-4b-------------------------------------------------------------------
ri
plot(ri, legend=FALSE)
text(ri)


## ----ch2-5a-------------------------------------------------------------------
ai1 <- aggregate(ri, c(2, 1), fun=mean)
ad1 <- aggregate(rd, c(2, 1), fun=mean)


## ----ch2-5b-------------------------------------------------------------------
as.matrix(ai1)
plot(ai1)
text(ai1, digits=1)


## ----ch2-6, fig.width=9, fig.height=4-----------------------------------------
s1 <- c(ai1, ad1)
names(s1) <- c("ind", "dep")
s1
plot(s1)


## -----------------------------------------------------------------------------
d1 <- as.data.frame(s1)
head(d1)


## -----------------------------------------------------------------------------
ma1 <- glm(dep~ind, data=d1)


## ----ch2-7--------------------------------------------------------------------
ai2 <- aggregate(ri, c(1, 2), fun=mean)
ad2 <- aggregate(rd, c(1, 2), fun=mean)
plot(ai2)
text(ai2, digits=1, srt=90)

s2 <- c(ai2, ad2)
names(s2) <- c('ind', 'dep')
# coerce to data.frame
d2 <- as.data.frame(s2)
ma2 <- glm(dep ~ ind, data=d2)


## -----------------------------------------------------------------------------
m$coefficients 
ma1$coefficients 
ma2$coefficients 


## ----include=FALSE------------------------------------------------------------
plotMAUP <- function(r1, r2, title="") {
  # get and plot the raster cell borders
	b1 <- as.polygons(r1, dissolve=FALSE) 
    plot(b1, axes=FALSE)
  # plot the raster cell values
	text(r1, cex=1.75)
	
  # add the title	
	if (length(title) == 1) {
		text(1.2, 1.1, title, xpd =NA, cex=2)
	} else {
		mtext(paste0('      ', title[1]), xpd =NA, cex=1.25)
	}
  # same for r2
	plot(as.polygons(r2, dissolve=FALSE), axes=FALSE)
	text(r2, cex=1.75)

  # add the title	
	if (length(title) > 1) {
		mtext(paste0(title[2], '      '), xpd =NA, cex=1.25)
	}

	
  # scatter plot
	i <- as.vector(values(r1))
	d <- as.vector(values(r2))
	plot(i, d, pch=18, xlim=c(0,100), ylim=c(0,100),
			axes=FALSE, ylab='', xlab='', yaxs="i", xaxs="i")
	axis(1, at=seq(0, 100, by=20), cex.axis=1.5)
	axis(2, at=seq(0, 100, by=20), cex.axis=1.5, las=1)

  # fit regression model
    m <- glm(d~i)
	s <- summary(m)
  # regression formula with coefficients
    f <- paste0('y = ', round(s$coefficients[2], 4), 'x + ', round(s$coefficients[1], 3))
	text(0, 95, f, pos=4, cex=1.75)
  # R-squared
    Rsq <- cor(d, predict(m))^2
	Rsq <- bquote(italic(R)^2 == .(round(Rsq, 4)))
	text(0, 85, Rsq, pos=4, cex=1.75)
	# regression line
    px <- data.frame(i = range(i))
	py <- predict(m, px)
	lines(cbind(px, py), lwd=2)
}



## ----figmaup, fig.cap='Figure 2.1   An illustration of MAUP', fig.width=10, fig.height=9----
# plotting parameters
par(mfrow=c(3,3), mai=c(0.25,0.15,0.25,0.15))

# Now call plotMAUP 3 times
plotMAUP(ri, rd, title=c('Independent variable', 'Dependent variable'))
# aggregation scheme 1
plotMAUP(ai1, ad1, title='Aggregation scheme 1')
# aggregation scheme 2
plotMAUP(ai2, ad2, title='Aggregation scheme 2')


## -----------------------------------------------------------------------------
A <- c(40, 43)
B <- c(1, 101)
C <- c(54, 111)
D <- c(104, 65)
E <- c(60, 22)
F <- c(20, 2)
pts <- rbind(A,B,C,D,E,F)
head(pts)


## ----points-------------------------------------------------------------------
plot(pts, xlim=c(0,120), ylim=c(0,120), pch=20, cex=2, col='red', xlab='X', ylab='Y', las=1)
text(pts+5, LETTERS[1:6])


## -----------------------------------------------------------------------------
dis <- dist(pts)
dis
D <- as.matrix(dis)
round(D)


## -----------------------------------------------------------------------------
a <-  D < 50
a


## -----------------------------------------------------------------------------
diag(a) <- NA
adj50 <- a * 1
adj50


## -----------------------------------------------------------------------------
cols <- apply(D, 1, order)
# we need to transpose the result
cols <- t(cols)


## -----------------------------------------------------------------------------
cols <- cols[, 2:4]
cols


## -----------------------------------------------------------------------------
rowcols <- cbind(rep(1:6, each=3), as.vector(t(cols)))
head(rowcols)


## -----------------------------------------------------------------------------
Ak3 <- adj50 * 0
Ak3[rowcols] <- 1
Ak3


## -----------------------------------------------------------------------------
W <- 1 / D
round(W, 4)


## -----------------------------------------------------------------------------
W[!is.finite(W)] <- NA


## -----------------------------------------------------------------------------
rtot <- rowSums(W, na.rm=TRUE)
# this is equivalent to
# rtot <- apply(W, 1, sum, na.rm=TRUE)
rtot


## -----------------------------------------------------------------------------
W <- W / rtot
rowSums(W, na.rm=TRUE)


## -----------------------------------------------------------------------------
colSums(W, na.rm=TRUE)


## -----------------------------------------------------------------------------
v <- voronoi(vect(pts))


## ----ch2-20, fig.width=6, fig.height=6----------------------------------------
par(mai=rep(0,4))
plot(v, lwd=4, border='gray', col=rainbow(6))
points(pts, pch=20, cex=2)
text(pts+5, toupper(letters[1:6]), cex=1.5)


## -----------------------------------------------------------------------------
class(v)
v