At 19:34 29/03/2012, JimeBoHe wrote:
I am a new user in R, so I am sorry if this is a basic question but I am kind
of lost here and I really would appreciatte your help.
Dear Jimena
Comments in-line
I have behavioral information of sea lions and I've done a binomial
Generalized Linear model using Rcmdr, to understand what variables are
affecting the place where Male's fights are being made (water or land). One
of my possible models includes Agression Type, number of females around the
fight and temperature in the moment.
*******************
> GLM.2 <- glm(Place ~ AgType + Females + Temp, family=binomial(cloglog),
+ data=June)
> summary(GLM.2)
Call:
glm(formula = Place ~ AgType + Females + Temp, family = binomial(cloglog),
data = June)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4448 -0.7248 0.4220 0.6741 1.8798
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.43071 0.82164 4.175 2.97e-05 ***
AgType -1.23463 0.37166 -3.322 0.000894 ***
Females 0.05812 0.01320 4.404 1.06e-05 ***
Temp -0.08614 0.02453 -3.512 0.000444 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 175.79 on 141 degrees of freedom
Residual deviance: 125.33 on 138 degrees of freedom
AIC: 133.33
Number of Fisher Scoring iterations: 8
****************************************
When I ran hosmerlem test it seems that I have a good fit to the data.
************************
> hosmerlem(Place, fitted(GLM.2))
X^2 Df P(>Chi)
7.5123428 8.0000000 0.4824925
*****************************
However if I run only AgType, or Agtype + Temp, when running Hosmerlem test;
it gives me an error that says: * ERROR: 'breaks' are not unique*.
I've been having the same error when I try to run other models or when I use
different link functions like logit.
So... right now my questions are:
1. What is this error means? as far as I have seen in the Hosmerlem test
code that I am using, it seems to depend on the yhat probs, but to be honest
I don't really understand what is this "yhat". The code I am using is:
yhat is the predicted value of y (the outcome variable). You do not
provide us with a reproducible example so this is only a guess but I
suspect if you look at the values of yhat for the failing models you
will find you have many replicated values and so at least two values
of the breakpoints are identical. But, as I say, this is at best a guess.
hosmerlem <-
function (y, yhat, g = 10)
{
cutyhat <- cut(yhat, breaks = quantile(yhat, probs = seq(0,
1, 1/g)), include.lowest = T)
obs <- xtabs(cbind(1 - y, y) ~ cutyhat)
expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
chisq <- sum((obs - expect)^2/expect)
P <- 1 - pchisq(chisq, g - 2)
c("X^2" = chisq, Df = g - 2, "P(>Chi)" = P)
}
2. I've seen in older posts that Cessie - van Houwelingen - Copas - Hosmer
unweighted sum of squares test for global goodness of fit test; seems to be
a better alternative. However I am probably doing something wrong because I
can't get it work.
I installed and loaded rms package and I ran the following command:
> resid(GLM.2, 'gof')
As far as I understand 'gof' is another package with this description:
Implementation of model-checking technique based on cumulative residuals.
So, for the command I also installed and loaded 'gof' package
You seem to be so far away from success here that I think remote
advice is not going to help.
However I get this: *ERROR: 'arg' should be one of "deviance", "pearson",
"working", "response", "partial"*
If I use any of the arguments that it suggests, it gives me a list of
residuals, but I doesn't show the Z or the P.
That is because it is doing what it is documented to. I am not
familiar with the gof function nor do I see it in the index for rms
(as you implied) so it is a bit hard to know what you are trying to do here.
So, if someone could help me understand what I am doing wrong I really would
appreciate it.
--
View this message in context:
http://r.789695.n4.nabble.com/Hosmer-Lemeshow-goodness-of-fit-tp3508127p4516441.html
Sent from the R help mailing list archive at Nabble.com.
Michael Dewey
i...@aghmed.fsnet.co.uk
http://www.aghmed.fsnet.co.uk/home.html
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
