Poission_ML

Poission_ML

setwd("/Users/subasishdas1/Desktop/TRB2021/Book/Chap10")
load("annual.RData")

##annual = read.table("http://myweb.fsu.edu/jelsner/data/AnnualData.txt", header = TRUE)
dat = subset(annual, Year >= 1950)
df = data.frame(H = dat$B.1, SOI = dat$soi, SST = dat$sst)
head(df)

##    H        SOI           SST
## 1 11  3.4333333 -0.0007373333
## 2  8 -3.4666667  0.2932626667
## 3  6 -0.5666667  0.3989293333
## 4  6 -3.6000000  0.2565960000
## 5  8  1.0666667  0.0112626667
## 6  9  4.2333333  0.2522626667

require(tree)

## Loading required package: tree

## Loading required package: tree
tr = tree(H ~ SOI + SST, data = df)
plot(tr)
text(tr)

tr2 = tree(H ~ SOI + SST, data = df[-61, ])

require(randomForest)

## Loading required package: randomForest

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## Loading required package: randomForest
## randomForest 4.6-7
## Type rfNews() to see new features/changes/bug fixes.
rf = randomForest(H ~ SOI + SST, data = df)

rf = randomForest(H ~ SOI + SST, data = df[-61, ])


n = length(df$H)
rfx = numeric(n)
prx = numeric(n)
for (i in 1:n) {
  rfm = randomForest(H ~ SOI + SST, data = df[-i, ])
  prm = glm(H ~ SOI + SST, data = df[-i, ], family = "poisson")
  new = df[i, ]
  rfx[i] = predict(rfm, newdata = new)
  prx[i] = predict(prm, newdata = new, type = "response")
}

newdat = expand.grid(SST = seq(-0.5, 0.7, 0.01), SOI = seq(-5, 4, 0.1))
z1 = predict(rf, newdata = newdat)
prm = glm(H ~ SOI + SST, data = df, family = "poisson")
z2 = predict(prm, newdata = newdat, type = "response")
newdat$Hrf = z1
newdat$Hpr = z2
require(lattice)

## Loading required package: lattice

require(colorRamps)

## Loading required package: colorRamps

## Loading required package: colorRamps
require(grid)

## Loading required package: grid

## Loading required package: grid
cr = blue2green(20)
p1 = levelplot(Hrf ~ SST + SOI, data = newdat, at = seq(2, 12, 2), scales = list(tck = 0.5, 
                                                                                 alternating = 2, cex = 0.7), xlab.top = textGrob("SST [C]", gp = gpar(cex = 0.8)), 
               ylab.right = textGrob("SOI [s.d.]", gp = gpar(cex = 0.8), rot = 90), xlab = "", 
               ylab = "", col.regions = cr, colorkey = list(space = "bottom", cex = 0.7), 
               sub = list("Hurricane count", cex = 0.9, font = 1))
p2 = levelplot(Hpr ~ SST + SOI, data = newdat, at = seq(2, 12, 2), scales = list(tck = 0.5, 
                                                                                 alternating = 2, cex = 0.7), xlab.top = textGrob("SST [C]", gp = gpar(cex = 0.8)), 
               ylab.right = textGrob("SOI [s.d.]", gp = gpar(cex = 0.8), rot = 90), xlab = "", 
               ylab = "", col.regions = cr, colorkey = list(space = "bottom", cex = 0.7), 
               sub = list("Hurricane count", cex = 0.9, font = 1))
p1 = update(p1, main = textGrob("a", x = unit(0.05, "npc")))
p2 = update(p2, main = textGrob("b", x = unit(0.05, "npc")))
print(p1, split = c(1, 1, 2, 1), more = TRUE)
print(p2, split = c(2, 1, 2, 1), more = FALSE)