On Oct 18, 2011, at 20:12 , aazaff wrote:
>
> # My understanding of how this is supposed to work is that the final line of
> the HOF function:
> "sum(-dbinom(y,M,fv,log=TRUE))" is supposed to give the negative log-odds of
> that function being a good fit.
>
>> HOF(p,x=Sample,y=Aequipecten,model=4)
> [1] 63.38763
>
> # Notice, however, that it returns a positive value, when it is supposed to
> be negative. What does this mean? What am I doing wrong?
Not reading up on the theory.
a: That looks like a negative log likelihood. "Log odds of being a good fit"
makes no sense (to me at least).
b: Likelihoods for discrete data are probabilities, hence the log likelihood is
negative and the negative log likelihood is positive.
c: You usually want to maximize the likelihood, hence minimize the negative log
likelihood.
d: Your "gaussian logistic regression" below, isn't. It's a binary logistic
regression, basically doing what HOF.start does, but calculating a different
set of transformed parameters.
e: The names on the return value from HOF.start doesn't seem to match the
function definition???
> #Furthermore, if I attempt a non-linear maximization on this function, I
> continue to get a positive negative log-odds and what I think are
> unreasonable looking parameter estimates. What does this mean?
>
>> sol<-nlm(HOF,p,x=Sample,y=Aequipecten,model=4)
>> sol
> $minimum
> [1] 32.72701
>
> $estimate
> [1] -8.751336 12.687632 8.281556
>
> $gradient
> [1] -4.640475e-07 -5.290583e-07 5.618925e-08
>
> $code
> [1] 1
>
> $iterations
> [1] 16
>
>
> # If I compare these results to simply doing a straight gaussian logistic
> regression on the vectors I get quite a different result. I'm very confident
> that this equation is working properly. I would expect the outputs of this
> function to be roughly equivalent to the results from object "p" as seen
> above, but they don't look similar at all. I'm not sure what this means.
>
> tBLparamsForOneSpecies <- function(theData) {
> x <- theData[ ,1]
> y <- theData[ ,2]
> logitReg <- glm(y ~ x + I(x^2), family=binomial)
> b0 <- logitReg$coefficients[1]
> b1 <- logitReg$coefficients[2]
> b2 <- logitReg$coefficients[3]
> opt <- (-b1)/(2*b2)
> tol <- 1 / sqrt(-2*b2)
> pmax<-1/(1+exp(b1^2/(4*b2)-b0))
> theParams <- cbind(opt, tol, pmax)
> theParams
> }
>
>> theData<-cbind(Sample,Aequipecten)
>> tBLparamsForOneSpecies(theData)
>> tBLparamsForOneSpecies(theData)
> opt tol pmax
> x 0.6337473 0.1468823 0.2433407
>
> #Nor do to the coefficients look similar to p either.
>
>> glm(Aequipecten~Sample+I(Sample^2),family=binomial)
>
> Call: glm(formula = Aequipecten ~ Sample + I(Sample^2), family = binomial)
>
> Coefficients:
> (Intercept) Sample I(Sample^2)
> -10.44 29.37 -23.18
>
> Degrees of Freedom: 86 Total (i.e. Null); 84 Residual
> Null Deviance: 73.38
> Residual Deviance: 68.07 AIC: 74.07
>
>
> # So anyway, to sum up my problem, I'm concerned that the results of my
> function HOF and of my non-linear maximization of HOF keep giving me a
> positive log-odds. I'm not sure if this is a result of a bug in the code
> somewhere, or if I actually just don't understand what I'm doing in terms of
> the statistics. I've tried to provide as much information as possible, but I
> really don't know where or what the problem is.
>
> Thank you for any help,
>
> Andrew
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/Non-linear-maximization-function-in-R-tp3916240p3916240.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd....@cbs.dk Priv: pda...@gmail.com
______________________________________________
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.
