I am trying to graph 2 rates together. Codes are as below. None of my code produces error messages so I really have no clue why this line matshade(np$c_year, ci.pred(miapF, np), col=cl[i])
does not appear. Any help would be greatly appreciated. I am new to R and my apologies if this is too simple question.
dataex:
' c_year age_group noofnewcases popno
2005 20 112 1116
2005 30 119 1184
2005 40 257 1176
2006 20 121 1141
2006 30 138 1161
2006 40 268 1164
2007 20 203 1167
2007 30 179 1145
2007 40 258 1151'
' library (Epi)
require(mgcv)
yrF <-transform (yrF, A=age_group+1/2, P=c_year+1/2)
gam.d <- gam(
cbind(noofnewcases, popno) ~ s(c_year, by = age_group),
family = poisreg,
data = yrF
)
ag <- seq(20, 40, 10)
n20 <- data.frame(c_year = 2005:2007, age_group = 20)
n30 <- data.frame(c_year = 2005:2007, age_group = 30)
n40 <- data.frame(c_year = 2005:2007, age_group = 40)
miapF <- gam(
cbind(n1, popno / 10^5) ~ s(A, c_year, k = 15),
family = poisreg,
data = yrF
)
ag <- seq(20, 40, 10)
cl <- gray(seq(.7, .1, length = 3))
np <- data.frame(c_year = seq(2005, 2007, 0.5))
plot(
NA , log = "y", ylim = c(15, 1500), xlim = c(2005, 2007),
xlab = "Year", ylab = "inc (per 100,000 Y)"
)
for(i in 1:2) {
np$A <- ag[i]
matshade(np$c_year, ci.pred(miapF, np), col = cl[i])
matshade(
n20$c_year,
cbind(
ci.pred(gam.d, n20) * 100000,
ci.pred(gam.d, n30) * 100000,
ci.pred(gam.d, n40) * 100000
),
plot = TRUE,
col = c("red", "orange", "yellow"),
alpha = 0,
lty = "44",
lend = "butt"
)
}'