[R] betareg help

[Ben Bolker]

   The problem seems to be that the algorithm for coming up with a
starting guess for the phi (dispersion) parameter is getting a negative
number.  It's not all that easy to figure this out ... The data set is a
little bit nasty (lots of points stacked on the equivalent of (0,0)),
but not pathological.
   I'm cc'ing the maintainer of the package --

adding the lines

  if (sum(phi_y)<0) {
       stop("bad estimated start value for phi: consider setting start
values manually\n",
            "(see ?betareg.control)\n",
            "Estimated starting values for
mean:\n",paste(beta,collapse=","))
    }

at line 162 of betareg.R in the current CRAN version provides
a more informative error message in this case.

  See below for solutions.
=============


results <- read.csv2("betareg_tmp.csv")
results$drugcat <- cut(results$drug,c(0,0.005,0.06,0.17))
table(results$drug)
table(results$alcoh)
table(results$cond)
## shows that fairly large fractions of the data are
##  in the lowest category:
## 165/209=0.001 'drug'
## 54/209=0.001 'alcoh'
## 38/209=0.001 'cond'
## so this will be a fairly challenging problem in any case

library(ggplot2)

ggplot(results,aes(x=alcoh,y=cond))+stat_sum(aes(size=..n..),alpha=0.7)+
  facet_wrap(~drugcat)+theme_bw()


library(betareg)
## set phi link to logarithmic
## basic problem (digging through betareg.fit etc.) is
##   that initial estimate of phi, based on
##   linear model of logit(cond) ~ alcoh + cond, is NEGATIVE ...
## doesn't seem to be any way to override this starting value
## brute force

try(gyl<-betareg(cond ~ alcoh + drug, data=results,
             link.phi="log"))

## pick through, debugging ... find starting values used

svec <- c(-1.6299469,0.8048446,1.7071124,0)
gyl<-betareg(cond ~ alcoh + drug, data=results,
             link.phi="log",
             control=betareg.control(start=svec))

## would work fine with more generic starting values

svec2 <- c(qlogis(mean(results$cond)),0,0,0)
gyl2<-betareg(cond ~ alcoh + drug, data=results,
             link.phi="log",
             control=betareg.control(start=svec2))

## before I got that to work, I also tried this (which
##  will be slower and less efficient but is a useful
## alternative

library(bbmle)

## define a variant parameterization of the beta distribution with
##  m=a/(a+b), phi=(a+b)
dbeta2 <- function(x,m,phi,log=FALSE) {

  a <- m*phi
  b <- phi*(1-m)
  dbeta(x,shape1=a,shape2=b,log=log)
}

m1 <- mle2(cond~dbeta2(m=plogis(mu),phi=exp(logphi)),
     parameters=list(mu~alcoh+drug),
     data=results,
     start=list(mu=qlogis(mean(results$cond)),logphi=0))
summary(m1)
p1 <- profile(m1)
plot(p1,show.points=TRUE)
confint(p1)
confint(m1,method="quad") ## not much difference

coef(m1)
coef(gyl)

On 11-03-13 12:59 PM, Vlatka Matkovic Puljic wrote:
> Sorry, here is my data (attached).
> 
> 2011/3/12 Ben Bolker <bbolker at gmail.com <mailto:bbolker at gmail.com>>
> 
>     Vlatka Matkovic Puljic <v.matkovic.puljic <at> gmail.com
>     <http://gmail.com>> writes:
> 
>     >
>     > That was also my first  thought.
>     > But I guess  it has something to do with W and phihat
>     > (which I'm struggling to check
> 
>      Again, it would help to post a reproducible example ...
>     hard to debug/diagnose by remote control.  If you can't
>     possibly post the data to the list, or put them on a web
>     site somewhere, or randomize them a bit so you're not
>     giving anything away, or find a simulated example that
>     shows the same problem, you could as a last resort send them
>     to me.
> 
>      Ben Bolker
> 
>     ______________________________________________
>     R-help at r-project.org <mailto:R-help at 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.
> 
> 
> 
> 
> -- 
> **************************
> Vlatka Matkovic Puljic
> +32/ 485/ 453340
>