mets (version 1.2.3.1)

biprobit: Bivariate Probit model

Description

Bivariate Probit model

Usage

biprobit(x, data, id, rho = ~1, num = NULL, strata = NULL,
  eqmarg = TRUE, indep = FALSE, weights = NULL, biweight,
  samecens = TRUE, randomeffect = FALSE, vcov = "robust",
  pairs.only = FALSE, allmarg = samecens & !is.null(weights),
  control = list(trace = 0), messages = 1, constrain = NULL,
  table = pairs.only, p = NULL, ...)

Arguments

x

formula (or vector)

data

data.frame

id

The name of the column in the dataset containing the cluster id-variable.

rho

Formula specifying the regression model for the dependence parameter

num

Optional name of order variable

strata

Strata

eqmarg

If TRUE same marginals are assumed (exchangeable)

indep

Independence

weights

Weights

biweight

Function defining the bivariate weight in each cluster

samecens

Same censoring

randomeffect

If TRUE a random effect model is used (otherwise correlation parameter is estimated allowing for both negative and positive dependence)

vcov

Type of standard errors to be calculated

pairs.only

Include complete pairs only?

allmarg

Should all marginal terms be included

control

Control argument parsed on to the optimization routine. Starting values may be parsed as 'start'.

messages

Control amount of messages shown

constrain

Vector of parameter constraints (NA where free). Use this to set an offset.

table

Type of estimation procedure

p

Parameter vector p in which to evaluate log-Likelihood and score function

...

Optional arguments

Examples

# NOT RUN {
data(prt)
prt0 <- subset(prt,country=="Denmark")
a <- biprobit(cancer~1+zyg, ~1+zyg, data=prt0, id="id")
b <- biprobit(cancer~1+zyg, ~1+zyg, data=prt0, id="id",pairs.only=TRUE)
predict(b,newdata=lava::Expand(prt,zyg=c("MZ")))
predict(b,newdata=lava::Expand(prt,zyg=c("MZ","DZ")))

# }
# NOT RUN {
 ## Reduce Ex.Timings
library(lava)
m <- lvm(c(y1,y2)~x)
covariance(m,y1~y2) <- "r"
constrain(m,r~x+a+b) <- function(x) tanh(x[2]+x[3]*x[1])
distribution(m,~x) <- uniform.lvm(a=-1,b=1)
ordinal(m) <- ~y1+y2
d <- sim(m,1000,p=c(a=0,b=-1)); d <- d[order(d$x),]
dd <- fast.reshape(d)

a <- biprobit(y~1+x,rho=~1+x,data=dd,id="id")
summary(a, mean.contrast=c(1,.5), cor.contrast=c(1,.5))
with(predict(a,data.frame(x=seq(-1,1,by=.1))), plot(p00~x,type="l"))

pp <- predict(a,data.frame(x=seq(-1,1,by=.1)),which=c(1))
plot(pp[,1]~pp$x, type="l", xlab="x", ylab="Concordance", lwd=2, xaxs="i")
confband(pp$x,pp[,2],pp[,3],polygon=TRUE,lty=0,col=Col(1))

pp <- predict(a,data.frame(x=seq(-1,1,by=.1)),which=c(9)) ## rho
plot(pp[,1]~pp$x, type="l", xlab="x", ylab="Correlation", lwd=2, xaxs="i")
confband(pp$x,pp[,2],pp[,3],polygon=TRUE,lty=0,col=Col(1))
with(pp, lines(x,tanh(-x),lwd=2,lty=2))

xp <- seq(-1,1,length.out=6); delta <- mean(diff(xp))
a2 <- biprobit(y~1+x,rho=~1+I(cut(x,breaks=xp)),data=dd,id="id")
pp2 <- predict(a2,data.frame(x=xp[-1]-delta/2),which=c(9)) ## rho
confband(pp2$x,pp2[,2],pp2[,3],center=pp2[,1])


# }
# NOT RUN {
## Time
# }
# NOT RUN {
    a <- biprobit.time(cancer~1, rho=~1+zyg, id="id", data=prt, eqmarg=TRUE,
                       cens.formula=Surv(time,status==0)~1,
                       breaks=seq(75,100,by=3),fix.censweights=TRUE)

    a <- biprobit.time2(cancer~1+zyg, rho=~1+zyg, id="id", data=prt0, eqmarg=TRUE,
                       cens.formula=Surv(time,status==0)~zyg,
                       breaks=100)

    a1 <- biprobit.time2(cancer~1, rho=~1, id="id", data=subset(prt0,zyg=="MZ"), eqmarg=TRUE,
                       cens.formula=Surv(time,status==0)~1,
                       breaks=100,pairs.only=TRUE)

    a2 <- biprobit.time2(cancer~1, rho=~1, id="id", data=subset(prt0,zyg=="DZ"), eqmarg=TRUE,
                        cens.formula=Surv(time,status==0)~1,
                        breaks=100,pairs.only=TRUE)

    prt0$trunc <- prt0$time*runif(nrow(prt0))*rbinom(nrow(prt0),1,0.5)
    a3 <- biprobit.time(cancer~1, rho=~1, id="id", data=subset(prt0,zyg=="DZ"), eqmarg=TRUE,
                        cens.formula=Surv(trunc,time,status==0)~1,
                        breaks=100,pairs.only=TRUE)

    plot(a,which=3,ylim=c(0,0.1))
# }