bivariate logistic/probit with hierarchical structure

Sun, May 15, 2011 9:39 AM

Hi David,

Yes, I see, it can be done that way, but care has to be taken. In your  
example the observation-level correlation is actually negative which  
means the group variance gets set to zero. You could reverse the  
scoring for one of the responses though:

sanction2$response2<-sanction2$response

sanction2$response2[which(sanction2$type=="export")]<-as.numeric(sanction2$response2[which(sanction2$type=="export")]==0)

m2<-lmer(response2 ~ type + type:coop + type:cost + type:target +  
(1|group), data=sanction2, family=binomial())

which gives a positive estimate of the group variance (i.e. a negative  
observation-level correlation).

Cheers,

Jarrod





Quoting David Duffy <davidD at qimr.edu.au> on Thu, 12 May 2011 08:48:52  
+1000 (EST):

On Wed, 11 May 2011, Jarrod Hadfield wrote:

Quoting David Duffy <davidD at qimr.edu.au> on Wed, 11 May 2011 16:10:04

On Fri, 6 May 2011, Arthur Charpentier wrote:

I would like to fit a bivariate probit model, i.e. I observe Y=(Y1,Y2),
I have some covariates X1,...Xk, and I assume that there is a hierarchical
structure (and random effects), i.e. Y_i is actually a Y_{i,j}  
where j is a
subgroup index.
So this yield two correlations levels,
- within subgroup correlation between Y1_i's (if they belong to the same
subgroup)

Since you are doing "only" a bivariate, you might consider  
rewriting it as a univariate problem with appropriate random  
effects (that is Y1 and Y2 are in a lower level with variable  
specific random intercepts etc).

The `residual' variances would not be identifiable even though the  
residual correlation is, which I think would cause problems unless  
the variances are fixed at some value (may be possible with the  
development version of lmer? )

I'm a bit slow, so I don't see where that would be a problem fitting  
in lmer per se.  Anyway, here is an example

library(Zelig)
require(VGAM)
data(sanction)
s1 <- sanction[,-5]
s2 <- sanction[,-4]
s1$group <- rownames(s1)
s2$group <- rownames(s2)
s1$type   <- "import"
s2$type   <- "export"
names(s2)[4] <- names(s1)[4] <- "response"
sanction2 <- rbind(s1,s2)
sanction2 <- sanction2[order(sanction2$group, sanction2$type),]
sanction2$type <- factor(sanction2$type)

lmer(response ~ coop + cost + target + type + (1|group),
     data=sanction2, family=binomial())
summary(vglm(cbind(import, export) ~ coop + cost + target,
             binom2.rho(exchangeable = TRUE), data = sanction))

lmer(response ~ type + type:coop + type:cost + type:target +  
(1|group), data=sanction2, family=binomial())
summary(vglm(cbind(import, export) ~ coop + cost + target,
             binom2.rho(exchangeable = FALSE), data = sanction))



-- 
| David Duffy (MBBS PhD)                                         ,-_|\
| email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

bivariate logistic/probit with hierarchical structure

Thread (6 messages)