[R-sig-ME] bivariate logistic/probit with hierarchical structure
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sun May 15 18:39:46 CEST 2011
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.
More information about the R-sig-mixed-models
mailing list