Respuestas Laboratorio 2

Exercise 1

library(spatstat)
plot(japanesepines)

intensity(japanesepines)
## [1] 65
D <- density(japanesepines, sigma = 0.1)
plot(D, main = "")
plot(D, main = "")
plot(japanesepines, add = TRUE, cols = "white", cex = 0.5, pch = 16)

plot(D, main = "", ribbon = FALSE)
plot(japanesepines, add = TRUE, cols = "white", cex = 0.5, pch = 16)

persp(D, theta=70, phi=25, shade=0.4)

Al buscar en la documentación, theta da la dirección azimutal y phi la colatitud, la sombra controla el sombreado de las facetas de la superficie.

range(D)
## [1]  17.47221 157.95229
round(integral(D))
## [1] 64

Exercise 2

elev <- bei.extra$elev
plot(elev)
plot(bei, add = TRUE, col = "black")

rh <- rhohat(bei, elev)
plot(rh)

prh <- predict(rh)
plot(prh, main = "")
plot(bei, add = TRUE, cols = "white", cex = .2, pch = 16)

dbei <- density(bei, sigma = bw.scott)
plot(dbei, main = "")
plot(bei, add = TRUE, cols = "white", cex = .2, pch = 16)

l <- solist(prh, dbei)
plot(l, equal.ribbon = TRUE, main = "",
     main.panel = c("rhohat prediction", "kernel smoothing"))

Exercise 3

replicate(3, plot(rpoispp(lambda = 100), main = ""))

## [[1]]
## Symbol map with no parameters
## 
## [[2]]
## Symbol map with constant values
## cols: #000000FE
## 
## [[3]]
## Symbol map with constant values
## cols: #000000FF
plot(rpoispp(lambda = 1.5), main = "")

Se esperan 1.5 puntos por cada realización.

Exercise 4

ppm(japanesepines~1)
## Stationary Poisson process
## Intensity: 65
##             Estimate      S.E.  CI95.lo  CI95.hi Ztest     Zval
## log(lambda) 4.174387 0.1240347 3.931284 4.417491   *** 33.65499
m.jp <- ppm(japanesepines ~ 1)
print(m.jp)
## Stationary Poisson process
## Intensity: 65
##             Estimate      S.E.  CI95.lo  CI95.hi Ztest     Zval
## log(lambda) 4.174387 0.1240347 3.931284 4.417491   *** 33.65499
unname(exp(coef(m.jp)))
## [1] 65
intensity(japanesepines)
## [1] 65

Los valores coinciden, con la función intensity y el logaritmo del lambda del coeficiente.

Exercise 5

jp.dens <- density(japanesepines, sigma = bw.scott)
plot(jp.dens)
plot(japanesepines, col = "white", cex = .4, pch = 16, add = TRUE)

jp.m <- ppm(japanesepines ~ x + y)
jp.m2 <- ppm(japanesepines ~ polynom(x, y, 2) )

coef(jp.m)
## (Intercept)           x           y 
##   4.0670790  -0.2349641   0.4296171
coef(jp.m2)
## (Intercept)           x           y      I(x^2)    I(x * y)      I(y^2) 
##   4.0645501   1.1436854  -1.5613621  -0.7490094  -1.2009245   2.5061569
par(mar=rep(0,4))
plot(predict(jp.m), main = "")

plot(predict(jp.m, se=TRUE)$se, main = "")

plot(predict(jp.m2), main = "")

plot(predict(jp.m2, se=TRUE)$se, main = "")

anova(jp.m, jp.m2)
## Analysis of Deviance Table
## 
## Model 1: ~x + y   Poisson
## Model 2: ~x + y + I(x^2) + I(x * y) + I(y^2)      Poisson
##   Npar Df Deviance
## 1    3            
## 2    6  3   3.3851
par(mar=rep(0.5,4))
plot(simulate(jp.m2, nsim=10), main = "")
## Generating 10 simulated patterns ...1, 2, 3, 4, 5, 6, 7, 8, 9,  10.

Exercise 6

fit0 <- ppm(japanesepines ~ 1)
fit1 <- update(fit0, . ~ x)
fit1
## Nonstationary Poisson process
## 
## Log intensity:  ~x
## 
## Fitted trend coefficients:
## (Intercept)           x 
##   4.2895587  -0.2349362 
## 
##               Estimate      S.E.   CI95.lo   CI95.hi Ztest       Zval
## (Intercept)  4.2895587 0.2411952  3.816825 4.7622926   *** 17.7845936
## x           -0.2349362 0.4305416 -1.078782 0.6089098       -0.5456759

Al agregar la covariable x esta resulta no significativa en el modelo.

fit2 <- update(fit1, . ~ . + y)
fit2
## Nonstationary Poisson process
## 
## Log intensity:  ~x + y
## 
## Fitted trend coefficients:
## (Intercept)           x           y 
##   4.0670790  -0.2349641   0.4296171 
## 
##               Estimate      S.E.    CI95.lo   CI95.hi Ztest       Zval
## (Intercept)  4.0670790 0.3341802  3.4120978 4.7220602   *** 12.1703167
## x           -0.2349641 0.4305456 -1.0788181 0.6088898       -0.5457357
## y            0.4296171 0.4318102 -0.4167154 1.2759495        0.9949211

Asi mismo, en este último modelo, las covariables x y y tampoco resultan significativas.

step(fit2)
## Start:  AIC=-407.96
## ~x + y
## 
##        Df     AIC
## - x     1 -409.66
## - y     1 -408.97
## <none>    -407.96
## 
## Step:  AIC=-409.66
## ~y
## 
##        Df     AIC
## - y     1 -410.67
## <none>    -409.66
## 
## Step:  AIC=-410.67
## ~1
## Stationary Poisson process
## Intensity: 65
##             Estimate      S.E.  CI95.lo  CI95.hi Ztest     Zval
## log(lambda) 4.174387 0.1240347 3.931284 4.417491   *** 33.65499

Finalmente, al eliminar x da como resultado el menor AIC, en el segundo paso, eliminar y da como resultado un AIC más bajo que no eliminar nada, quedando el modelo nulo.

Exercise 7

bei.m <- ppm(bei ~ elev + grad, data = bei.extra)
coef(bei.m)
## (Intercept)        elev        grad 
## -8.56355220  0.02143995  5.84646680
plot(predict(bei.m), main = "")
plot(bei, cex = 0.3, pch = 16, cols = "white", add = TRUE)

vcov(bei.m)
##               (Intercept)          elev          grad
## (Intercept)  0.1163586583 -7.774771e-04 -0.0354792767
## elev        -0.0007774771  5.234331e-06  0.0001992266
## grad        -0.0354792767  1.992266e-04  0.0654239289
std.err <- predict(bei.m, se = TRUE)$se
plot(std.err, main = "")

Exercise 8

coef(ppm1 <- ppm(japanesepines ~ 1))
## log(lambda) 
##    4.174387
coef(ppm2 <- ppm(japanesepines ~ x))
## (Intercept)           x 
##   4.2895587  -0.2349362
coef(ppm3 <- ppm(japanesepines ~ sin(x)))
## (Intercept)      sin(x) 
##   4.2915935  -0.2594537
coef(ppm4 <- ppm(japanesepines ~ x + y))
## (Intercept)           x           y 
##   4.0670790  -0.2349641   0.4296171
coef(ppm5 <- ppm(japanesepines ~ polynom(x, y, 2)))
## (Intercept)           x           y      I(x^2)    I(x * y)      I(y^2) 
##   4.0645501   1.1436854  -1.5613621  -0.7490094  -1.2009245   2.5061569
coef(ppm6 <- ppm(japanesepines ~ factor(x < 0.4)))
##         (Intercept) factor(x < 0.4)TRUE 
##           4.1048159           0.1632665

Exercise 9

plot(predict(ppm1), main = "")

plot(predict(ppm2), main = "")

plot(predict(ppm3), main = "")

plot(predict(ppm4), main = "")

plot(predict(ppm5), main = "")

plot(predict(ppm6), main = "")

Exercise 10

plot(split(hamster), main = "")

plot(density(split(hamster)), main = "")

plot(relrisk(hamster, hmax = 1, relative = TRUE, control = "dividing"))

Exercise 11

plot(ants)

fit1 <- ppm(ants ~ marks)
exp(coef(fit1)[1])
##  (Intercept) 
## 6.762949e-05
exp(coef(fit1)[1] + coef(fit1)[2])
##  (Intercept) 
## 0.0001585795

Este último resultado coincide con el siguiente a través de summary:

summary(ants)
## Marked planar point pattern:  97 points
## Average intensity 0.0002261486 points per square unit (one unit = 0.5 feet)
## 
## Coordinates are integers
## i.e. rounded to the nearest unit (one unit = 0.5 feet)
## 
## Multitype:
##             frequency proportion    intensity
## Cataglyphis        29  0.2989691 6.761144e-05
## Messor             68  0.7010309 1.585372e-04
## 
## Window: polygonal boundary
## single connected closed polygon with 11 vertices
## enclosing rectangle: [-25, 803] x [-49, 717] units
##                      (828 x 766 units)
## Window area = 428922 square units
## Unit of length: 0.5 feet
## Fraction of frame area: 0.676
fit2 <- ppm(ants ~ marks + x)
(co <- coef(fit2))
##   (Intercept)   marksMessor             x 
## -9.5243832518  0.8522118655 -0.0002041438

Expresiones: \(\lambda( (x,y) ) = \exp(-9.5243833 + -2.0414381\times 10^{-4} \cdot x)\) \(\lambda( (x,y) ) = \exp(-9.5243833 + 0.8522119 + -2.0414381\times 10^{-4} \cdot x)\)

fit3 <- ppm(ants ~ marks * x)
(co <- coef(fit3))
##   (Intercept)   marksMessor             x marksMessor:x 
## -9.605698e+00  9.676854e-01  1.107981e-05 -3.071343e-04

Expresiones: \(\lambda( (x,y) ) = \exp(-9.605698 + 1.1079805\times 10^{-5} \cdot x)\) \(\lambda( (x,y) ) = \exp(-9.605698 + 0.9676854 + (1.1079805\times 10^{-5} + 0.9676854) \cdot x)\)

pred <- c(predict(fit1), predict(fit2), predict(fit3))
plot(as.solist(pred), ncols = 2, main = "")

Para el modelo aditivo el efecto de la coordenada x es el mismo para ambos tipos de hormigas, mientras que el efecto de x difiere en el modelo multiplicativo.

Exercise 12

fs <- function(x,y) {
  ends <- ants.extra$fieldscrub
  angle <- atan(diff(ends$y)/diff(ends$x))
  normal <- angle + pi/2
  project <- (x - ends$x[1]) * cos(normal) + (y - ends$y[1]) * sin(normal)
  factor(ifelse(project > 0, "scrub", "field"))
}
fit1 <- ppm(ants ~ marks + side, data = list(side=fs))
fit2 <- ppm(ants ~ marks * side, data = list(side=fs))
fit1
## Nonstationary multitype Poisson process
## 
## Possible marks: 'Cataglyphis' and 'Messor'
## 
## Log intensity:  ~marks + side
## 
## Fitted trend coefficients:
## (Intercept) marksMessor   sidescrub 
## -9.57189990  0.85221187 -0.07380742 
## 
##                Estimate      S.E.    CI95.lo    CI95.hi Ztest        Zval
## (Intercept) -9.57189990 0.2027872 -9.9693554 -9.1744444   *** -47.2017058
## marksMessor  0.85221187 0.2217851  0.4175210  1.2869027   ***   3.8425114
## sidescrub   -0.07380742 0.2080023 -0.4814844  0.3338695        -0.3548395
fit2
## Nonstationary multitype Poisson process
## 
## Possible marks: 'Cataglyphis' and 'Messor'
## 
## Log intensity:  ~marks * side
## 
## Fitted trend coefficients:
##           (Intercept)           marksMessor             sidescrub 
##            -9.3509797             0.5198755            -0.7789884 
## marksMessor:sidescrub 
##             0.9682016 
## 
##                         Estimate      S.E.      CI95.lo     CI95.hi Ztest
## (Intercept)           -9.3509797 0.2132007 -9.768845473 -8.93311402   ***
## marksMessor            0.5198755 0.2692240 -0.007793922  1.04754484      
## sidescrub             -0.7789884 0.4339489 -1.629512751  0.07153586      
## marksMessor:sidescrub  0.9682016 0.4975910 -0.007058797  1.94346199      
##                             Zval
## (Intercept)           -43.859983
## marksMessor             1.931014
## sidescrub              -1.795115
## marksMessor:sidescrub   1.945778

En el primer modelo la intensidad ajustada es menor en el srub que en el field.

En el segundo modelo, la intensidad ajustada de Cataglyphis es menor en el scrub que la intensidad de Cataglyphis en el campo, mientras que para Messor es al revés. Cuando se permite el efecto diferente entre los tipos de hormigas, la covariable matorral/campo es significativa.