Obs-plots/obs plots areal data.R at master · wcwetzel/Obs-plots · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
### Analyzing obs plots with plants as areal data
## 20 Feb 2012

# load packages
library(spdep)
library(bbmle)
library(sp)
library(pgirmess)

# load data
d2 = read.csv(
	'~/Documents/DATA/2010 DATA/FIELD/plot 2/plot2-2010-2011.csv')
d2 = d2[!is.na(d2$galls2010),]
d2$mgalls = with(d2, (galls2011 + galls2010)/2 )
d2$vol = with(d2, d1/100 * d2/100 * h/100) # cubic meters
d2$density = d2$galls2011/d2$vol
d2$change = d2$galls2011 - d2$galls2010

# turn x-y coordinates into a square distance matrix
dists = as.matrix(dist(cbind(d2$x,d2$y), diag=TRUE, upper=TRUE))

# now convert distance matrix into weights matrix
# with negative exp decline in weight with distance
curve(exp(-x * 1), 0, max(dists))
abline(h=0.05)
abline(h=0.1)
abline(h=0.25)
nexpW = exp(-dists)

# mat2listw converts a square spatial weights matrix into
# a spdep weights list object
nelistW = mat2listw(nexpW, row.names=d2$tag)

# now try a neg exp weight with a cut off at 0.15
# with negative exp decline in weight with distance
nexpcoW = nexpW
nexpcoW[nexpcoW < 0.25] = 0

# mat2listw converts a square spatial weights matrix into
# a spdep weights list object
necolistW = mat2listw(nexpcoW, row.names=d2$tag)


# some things we can do with spdep spatial weights objects
summary(nelistW)
summary(necolistW)

coordinates(d2) = cbind(d2$x, d2$y)

plot(nelistW, coordinates(d2))
plot(necolistW, coordinates(d2))


# OK now Moran's I
uncorr = rpois(nrow(d2), lambda=mean(d2$mgalls))
mt.ne = moran.test(d2$galls2011, listw=nelistW)
mt.neco = moran.test(d2$galls2011, listw=necolistW)
moran.test(d2$mgalls, listw=necolistW)
moran.test(d2$galls2010, listw=necolistW)
moran.test(uncorr, listw=necolistW)
# it's not working. i'm getting very sig autocorr for
# a random distribution of galls. crazy!

moran.mc(uncorr, necolistW, nsim=4999)
moran.mc(d2$galls2011, necolistW, nsim=4999)
moran.mc(d2$galls2010, necolistW, nsim=4999)
moran.mc(d2$mgalls, necolistW, nsim=4999)
mcm.neco = moran.mc(d2$change, necolistW, nsim=4999)
mcm.ne = moran.mc(d2$change, nelistW, nsim=4999)


hikingboots = mcm.ne$res[-length(mcm.ne$res)]
bootdens = density(hikingboots)
hist(hikingboots, prob=TRUE)
lines(bootdens)
abline(v=mcm.ne$statistic,lwd=2,col=4)


# moran scatter plot:
# this shows how similar a point is to its neighbors
# x-axis is d2$change
# y-axis is a weighted average of neighbors
moran.plot(d2$change, nelistW)
# there seems to be one far out influential point
# but even with that point removed, there is a strong trend
# for points with a bigger change
# to have neighbors with a bigger change


# what about using residuals instead of raw change?

model = glm(galls2011 ~ galls2010, family=poisson, data=d2)
res = resid(model)
mcm.ne.res = moran.mc(res, nelistW, nsim=4999)
moran.plot(res, nelistW)


hikingboots2 = mcm.ne.res$res[-length(mcm.ne.res$res)]
bootdens2 = density(hikingboots2)
hist(hikingboots2, prob=TRUE)
lines(bootdens2)
abline(v=mcm.ne.res$statistic,lwd=2,col=4)

## look for non-stationarity ##
plot(galls2011 ~ x, data=d2)
plot(galls2010 ~ x, data=d2)
plot(change ~ x, data=d2)

plot(res ~ x, data=d2)
summary(lm(res~ d2$x))
abline(lm(res~ d2$x))

plot(res ~ y, data=d2)
summary(lm(res~ d2$y))
abline(lm(res~ d2$y))
# It looks as though change in galls increases with x coord!

# residuals from a model that has x coords in it too
model2 = glm(galls2011 ~ galls2010 + x, family=poisson, data=d2)
res2 = resid(model2)
mcm.ne.res2 = moran.mc(res2, nelistW, nsim=4999)
moran.plot(res2, nelistW)


# spatial correlograms
#neighbors = listw2sn(nelistW)
neighbors = dnearneigh(coordinates(d2), d1=0, d2=1.5, row.names=d2$tag)
summary(neighbors)
plot(neighbors, coordinates(d2))
cor1 = sp.correlogram(neighbors, var=d2$change, order=8)
print(cor1, p.adj.method='holm')
plot(cor1)


# pgirmess correlog
cor2 = correlog(coordinates(d2), d2$change)
plot(cor2)

cor3 = correlog(coordinates(d2), d2$galls2011)
plot(cor3)