[R] glm.predict?
Patrick Connolly
P.Connolly at hortresearch.co.nz
Tue Jan 22 04:00:42 CET 2002
I've been attempting to calculate the predictions from a poisson glm
object, along these lines:
predict(foo.glm, type = "response")
and
predict(foo.glm, type = "response", se.fit = TRUE)
foo.glm is arrived at this way:
foo.glm <- glm(Insects ~ Dad * Mum + Location, offset = log(MM), family
= "poisson", data = model.df)
There are two levels for Dad and three each for Mum and Location but
the numbers of each level vary considerably.
The help file for predict.glm tells us that the returned value is:
If `se = FALSE', a vector or matrix of predictions. If `se =
TRUE', a list with components
fit: Predictions
se.fit: Estimated standard errors
So, I imagined I'd get the same result for predictions arrived at
these two ways:
things <- list()
things$fitA <- predict(foo.glm, type = "response")
things$fitB <- predict(foo.glm, type = "response", se.fit = TRUE)$fit
however, that is not the case
Warning messages:
1: longer object length
is not a multiple of shorter object length in: predictor + offset
2: longer object length
is not a multiple of shorter object length in: se.fit * abs(family(object)$mu.eta(fit))
Browse[1]> sapply(things,length)
fitA fitB
135 141
Browse[1]> sapply(things, function(x) length(x[is.na(x)]))
fitA fitB
0 6
So the warning message indicates something serious is happening.
Simply removing the 6 NAs does not fix the problem (though it's close
to it). Removing the offset would eliminate the problem, but it would
not give me usable answers.
Looking through the code for predict.glm, I got confused by the
relationship between type and terms and by the subtleties of the
concise coding using switch. Without any comments in the code, I
couldn't follow what was happening.
Ideas are gratefully accepted. If anyone needs the dataframe to be
able to tell what's happening, it looks like this:
dput(model.df)
structure(list(Dad = structure(c(2, 2, 1, 1, 1, 2, 2, 2, 2, 2,
2, 2, 1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 2,
2, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2,
2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1,
1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2,
2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2), .Label = c("1",
"2"), class = "factor"), Mum = structure(c(3, 3, 3, 3, 2, 1,
3, 3, 1, 1, 1, 2, 1, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 1, 1, 1, 1,
2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 1, 1, 2, 1, 1, 2, 3, 3,
3, 3, 2, 1, 3, 1, 1, 1, 2, 1, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 1,
1, 1, 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 1, 1, 2, 1, 2,
3, 3, 3, 3, 2, 1, 3, 3, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 2, 3,
3, 3, 1, 1, 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 1, 1, 2,
1, 1, 2), .Label = c("A", "B", "C"), class = "factor"), Location = structure(c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3), .Label = c("Lead1", "Lead2", "Trunk"
), class = "factor"), Insects = structure(c(2, 2, 2, 2, 1, 2,
0, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 3, 3, 31, 0, 0, 0, 38, 0, 0, 0, 0, 0, 0, 16, 0, 1, 0, 4,
1, 17, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 3, 0, 0, 1, 0, 0, 1, 15, 4, 1, 0, 0, 2, 1, 0, 0, 3, 2, 0,
11, 0, 0, 3, 4, 21, 0, 16, 6, 0, 0, 0, 0, 0, 0, 2, 0, 1, 3, 0,
1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0, 0, 8, 6, 0, 0,
0, 3, 0, 93, 1, 17, 0), .Names = c("Lead1.no1", "Lead1.no2",
"Lead1.no3", "Lead1.no4", "Lead1.no5", "Lead1.no6", "Lead1.no7",
"Lead1.no8", "Lead1.no9", "Lead1.no10", "Lead1.no11", "Lead1.no12",
"Lead1.no13", "Lead1.no14", "Lead1.no15", "Lead1.no16", "Lead1.no17",
"Lead1.no18", "Lead1.no19", "Lead1.no20", "Lead1.no21", "Lead1.no22",
"Lead1.no23", "Lead1.no25", "Lead1.no26", "Lead1.no27", "Lead1.no28",
"Lead1.no29", "Lead1.no30", "Lead1.no31", "Lead1.no32", "Lead1.no33",
"Lead1.no34", "Lead1.no35", "Lead1.no36", "Lead1.no37", "Lead1.no38",
"Lead1.no39", "Lead1.no40", "Lead1.no41", "Lead1.no42", "Lead1.no43",
"Lead1.no44", "Lead1.no45", "Lead1.no46", "Lead1.no47", "Lead2.no1",
"Lead2.no2", "Lead2.no3", "Lead2.no4", "Lead2.no5", "Lead2.no6",
"Lead2.no8", "Lead2.no9", "Lead2.no10", "Lead2.no11", "Lead2.no12",
"Lead2.no13", "Lead2.no14", "Lead2.no15", "Lead2.no16", "Lead2.no17",
"Lead2.no18", "Lead2.no19", "Lead2.no20", "Lead2.no21", "Lead2.no22",
"Lead2.no23", "Lead2.no25", "Lead2.no26", "Lead2.no27", "Lead2.no28",
"Lead2.no29", "Lead2.no30", "Lead2.no31", "Lead2.no32", "Lead2.no33",
"Lead2.no34", "Lead2.no35", "Lead2.no36", "Lead2.no37", "Lead2.no38",
"Lead2.no39", "Lead2.no40", "Lead2.no41", "Lead2.no42", "Lead2.no43",
"Lead2.no44", "Lead2.no45", "Lead2.no47", "Trunk.no1", "Trunk.no2",
"Trunk.no3", "Trunk.no4", "Trunk.no5", "Trunk.no6", "Trunk.no7",
"Trunk.no8", "Trunk.no9", "Trunk.no10", "Trunk.no11", "Trunk.no13",
"Trunk.no14", "Trunk.no15", "Trunk.no16", "Trunk.no17", "Trunk.no18",
"Trunk.no19", "Trunk.no20", "Trunk.no21", "Trunk.no22", "Trunk.no23",
"Trunk.no24", "Trunk.no25", "Trunk.no27", "Trunk.no28", "Trunk.no29",
"Trunk.no30", "Trunk.no31", "Trunk.no32", "Trunk.no33", "Trunk.no34",
"Trunk.no35", "Trunk.no36", "Trunk.no37", "Trunk.no38", "Trunk.no39",
"Trunk.no40", "Trunk.no41", "Trunk.no42", "Trunk.no43", "Trunk.no44",
"Trunk.no45", "Trunk.no46", "Trunk.no47")), MM = structure(c(123,
88, 133, 90, 153, 96, 116, 125, 89, 102, 107, 110, 112, 135,
138, 145, 94, 132, 90, 117, 130, 135, 82, 124, 158, 151, 135,
114, 126, 71, 113, 84, 97, 88, 115, 88, 120, 102, 123, 100, 92,
126, 96, 108, 130, 76, 122, 110, 205, 219, 192, 126, 153, 162,
162, 112, 180, 101, 205, 157, 239, 241, 165, 153, 146, 134, 135,
155, 113, 154, 203, 94, 108, 106, 83, 125, 118, 104, 157, 144,
140, 88, 110, 173, 67, 111, 83, 72, 76, 135, 177, 142, 177, 142,
166, 144, 164, 171, 200, 187, 169, 189, 168, 182, 208, 207, 193,
144, 178, 177, 176, 205, 153, 228, 227, 147, 173, 157, 214, 167,
140, 179, 204, 184, 151, 115, 173, 208, 135, 175, 136, 121, 189,
148, 174), .Names = c("Lead1.mm1", "Lead1.mm2", "Lead1.mm3",
"Lead1.mm4", "Lead1.mm5", "Lead1.mm6", "Lead1.mm7", "Lead1.mm8",
"Lead1.mm9", "Lead1.mm10", "Lead1.mm11", "Lead1.mm12", "Lead1.mm13",
"Lead1.mm14", "Lead1.mm15", "Lead1.mm16", "Lead1.mm17", "Lead1.mm18",
"Lead1.mm19", "Lead1.mm20", "Lead1.mm21", "Lead1.mm22", "Lead1.mm23",
"Lead1.mm25", "Lead1.mm26", "Lead1.mm27", "Lead1.mm28", "Lead1.mm29",
"Lead1.mm30", "Lead1.mm31", "Lead1.mm32", "Lead1.mm33", "Lead1.mm34",
"Lead1.mm35", "Lead1.mm36", "Lead1.mm37", "Lead1.mm38", "Lead1.mm39",
"Lead1.mm40", "Lead1.mm41", "Lead1.mm42", "Lead1.mm43", "Lead1.mm44",
"Lead1.mm45", "Lead1.mm46", "Lead1.mm47", "Lead2.mm1", "Lead2.mm2",
"Lead2.mm3", "Lead2.mm4", "Lead2.mm5", "Lead2.mm6", "Lead2.mm8",
"Lead2.mm9", "Lead2.mm10", "Lead2.mm11", "Lead2.mm12", "Lead2.mm13",
"Lead2.mm14", "Lead2.mm15", "Lead2.mm16", "Lead2.mm17", "Lead2.mm18",
"Lead2.mm19", "Lead2.mm20", "Lead2.mm21", "Lead2.mm22", "Lead2.mm23",
"Lead2.mm25", "Lead2.mm26", "Lead2.mm27", "Lead2.mm28", "Lead2.mm29",
"Lead2.mm30", "Lead2.mm31", "Lead2.mm32", "Lead2.mm33", "Lead2.mm34",
"Lead2.mm35", "Lead2.mm36", "Lead2.mm37", "Lead2.mm38", "Lead2.mm39",
"Lead2.mm40", "Lead2.mm41", "Lead2.mm42", "Lead2.mm43", "Lead2.mm44",
"Lead2.mm45", "Lead2.mm47", "Trunk.mm1", "Trunk.mm2", "Trunk.mm3",
"Trunk.mm4", "Trunk.mm5", "Trunk.mm6", "Trunk.mm7", "Trunk.mm8",
"Trunk.mm9", "Trunk.mm10", "Trunk.mm11", "Trunk.mm13", "Trunk.mm14",
"Trunk.mm15", "Trunk.mm16", "Trunk.mm17", "Trunk.mm18", "Trunk.mm19",
"Trunk.mm20", "Trunk.mm21", "Trunk.mm22", "Trunk.mm23", "Trunk.mm24",
"Trunk.mm25", "Trunk.mm27", "Trunk.mm28", "Trunk.mm29", "Trunk.mm30",
"Trunk.mm31", "Trunk.mm32", "Trunk.mm33", "Trunk.mm34", "Trunk.mm35",
"Trunk.mm36", "Trunk.mm37", "Trunk.mm38", "Trunk.mm39", "Trunk.mm40",
"Trunk.mm41", "Trunk.mm42", "Trunk.mm43", "Trunk.mm44", "Trunk.mm45",
"Trunk.mm46", "Trunk.mm47"))), .Names = c("Dad", "Mum", "Location",
"Insects", "MM"), row.names = c("X9", "X10", "X21", "X22", "X23",
"X31", "X50", "X51", "X62", "X63", "X64", "X65", "X82", "X101",
"X102", "X115", "X116", "X126", "X127", "X128", "X129", "X153",
"X154", "X169", "X170", "X178", "X179", "X180", "X199", "X200",
"X201", "X242", "X243", "X244", "X245", "X246", "X254", "X255",
"X256", "X257", "X295", "X296", "X303", "X305", "X306", "X307",
"X91", "X1021", "X213", "X224", "X235", "X316", "X518", "X629",
"X6310", "X6411", "X6512", "X8213", "X10114", "X10215", "X11516",
"X11617", "X12618", "X12719", "X12820", "X12921", "X15322", "X15423",
"X16925", "X17026", "X17827", "X17928", "X18029", "X19930", "X20031",
"X20132", "X24233", "X24334", "X24435", "X24536", "X24637", "X25438",
"X25539", "X25640", "X25741", "X29542", "X29643", "X30344", "X30545",
"X30747", "X948", "X1049", "X2150", "X2251", "X2352", "X3153",
"X5054", "X5155", "X6256", "X6357", "X6458", "X8260", "X10161",
"X10262", "X11563", "X11664", "X12665", "X12766", "X12867", "X12968",
"X15369", "X15470", "X15571", "X16972", "X17874", "X17975", "X18076",
"X19977", "X20078", "X20179", "X24280", "X24381", "X24482", "X24583",
"X24684", "X25485", "X25586", "X25687", "X25788", "X29589", "X29690",
"X30391", "X30592", "X30693", "X30794"), class = "data.frame")
Thanks
--
*************************************************************
___ Patrick Connolly
{~._.~} HortResearch Great minds discuss ideas;
_( Y )_ Mt Albert Average minds discuss events;
(:_~*~_:) Auckland Small minds discuss people.
(_)-(_) New Zealand .... Anon
Ph: +64-9 815 4200 x 7188
*************************************************************
______________________________________________________
The contents of this e-mail are privileged and/or confidential to the
named recipient and are not to be used by any other person and/or
organisation. If you have received this e-mail in error, please notify
the sender and delete all material pertaining to this e-mail.
______________________________________________________
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list