[R] Using predict() After Adding a Factor to a glm.nb() Model

[077315q at acadiau.ca]

# Hello,

# I have a data set that looks something like the following:

site<-c(rep('a',5),rep('b',2),rep('c',4),rep('d',11))
year<-c(1980, 1981, 1982, 1993, 1995, 1980, 1983, 1981, 1993, 1995,  
1999, c(1980:1990))
count<-c(60,35,36,12,8,112,98,20,13,15,15,65,43,49,51,34,33,33,33,40,11,0)

data<-data.frame(site, year, count)

# > site year count
# 1     a 1980    60
# 2     a 1981    35
# 3     a 1982    36
# 4     a 1993    12
# .
# .
# .

# I ran a negative binomial glm using:

library(MASS)
model_a<-glm.nb(count~year, data=data)

# then extracted predicted values using:

py<-data.frame(seq(from=min(data$year), to=max(data$year), by=1))
p1<-predict(model_a, newdata=py, se.fit=TRUE, type='response')
f1<-p1$fit

plot(count~year, data=data)
lines(py$year, f1)

# Works perfectly.

# Now, I want to add site as a factor:
model_b<-glm.nb(count~year+as.factor(site), data=data)

# I have tried a couple different ways of predicting the values based  
on model_b, but am having trouble.

py<-data.frame(seq(from=min(data$year), to=max(data$year), by=1))
p1<-predict(model_b, newdata=py, se.fit=TRUE, type='response')

# gives:
# >Error in model.frame.default(Terms, newdata, na.action = na.action,  
xlev = object$xlevels) :
#   variable lengths differ (found for 'as.factor(site)')
# In addition: Warning message:
# 'newdata' had 20 rows but variable(s) found have 22 rows

py<-expand.grid(data$site, data$year)
names(py)<-c('site','year')
p1<-predict(model_b, newdata=py)

# This works, but results in 484 values, and I can't plot a line over  
my points.

# There is probably a simple solution, but I'm having trouble wrapping  
my mind around it. Mind you, this is also a last
# minute change to my thesis, and I haven't slept in about three days.

# Any suggestions? I would be extremely grateful...

# Banging my head against the wall,
# A stressed out grad student