Description Usage Arguments Details Examples
If clusters contain more than two times, the algoritm uses a compososite likelihood based on the pairwise bivariate models.
1 2 3 4 5 6 7 | easy.binomial.twostage(margbin = NULL, data = sys.parent(),
score.method = "nlminb", response = "response", id = "id", Nit = 60,
detail = 0, silent = 1, weights = NULL, control = list(),
theta = NULL, theta.formula = NULL, desnames = NULL, deshelp = 0,
var.link = 1, iid = 1, step = 0.5, model = "plackett",
marginal.p = NULL, strata = NULL, max.clust = NULL,
se.clusters = NULL)
|
margbin |
Marginal binomial model |
data |
data frame |
response |
name of response variable in data frame |
id |
name of cluster variable in data frame |
score.method |
Scoring method |
Nit |
Number of iterations |
detail |
Detail for more output for iterations |
silent |
Debug information |
weights |
Weights for log-likelihood, can be used for each type of outcome in 2x2 tables. |
control |
Optimization arguments |
theta |
Starting values for variance components |
theta.formula |
design for depedence, either formula or design function |
desnames |
names for dependence parameters |
deshelp |
if 1 then prints out some data sets that are used, on on which the design function operates |
var.link |
Link function for variance |
iid |
Calculate i.i.d. decomposition |
step |
Step size |
model |
model |
marginal.p |
vector of marginal probabilities |
strata |
strata for fitting |
max.clust |
max clusters |
se.clusters |
clusters for iid decomposition for roubst standard errors |
The reported standard errors are based on the estimated information from the likelihood assuming that the marginals are known. This gives correct standard errors in the case of the plackett distribution (OR model for dependence), but incorrect for the clayton-oakes types model. The OR model is often known as the ALR model, but our fitting procedures gives correct standard errors and is quite a bit quicker.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 | data(twinstut)
twinstut0 <- subset(twinstut, tvparnr<2300000)
twinstut <- twinstut0
theta.des <- model.matrix( ~-1+factor(zyg),data=twinstut0)
margbin <- glm(stutter~factor(sex)+age,data=twinstut0,family=binomial())
bin <- binomial.twostage(margbin,data=twinstut0,
clusters=twinstut0$tvparnr,theta.des=theta.des,detail=0,
score.method="fisher.scoring")
summary(bin)
twinstut0$cage <- scale(twinstut0$age)
theta.des <- model.matrix( ~-1+factor(zyg)+cage,data=twinstut0)
bina <- binomial.twostage(margbin,data=twinstut0,
clusters=twinstut0$tvparnr,theta.des=theta.des,detail=0,
score.method="fisher.scoring")
summary(bina)
theta.des <- model.matrix( ~-1+factor(zyg)+factor(zyg)*cage,data=twinstut0)
bina <- binomial.twostage(margbin,data=twinstut0,
clusters=twinstut0$tvparnr,theta.des=theta.des,detail=0,
score.method="fisher.scoring")
summary(bina)
twinstut0$binstut <- (twinstut0$stutter=="yes")*1
out <- easy.binomial.twostage(stutter~factor(sex)+age,data=twinstut0,
response="binstut",id="tvparnr",
theta.formula=~-1+factor(zyg1),
score.method="fisher.scoring")
summary(out)
## refers to zygosity of first subject in eash pair : zyg1
## could also use zyg2 (since zyg2=zyg1 within twinpair's))
desfs <- function(x,num1="zyg1",namesdes=c("mz","dz","os"))
c(x[num1]=="mz",x[num1]=="dz",x[num1]=="os")*1
out3 <- easy.binomial.twostage(binstut~factor(sex)+age,
data=twinstut0, response="binstut",id="tvparnr",
score.method="fisher.scoring",
theta.formula=desfs,desnames=c("mz","dz","os"))
summary(out3)
n <- 10000
set.seed(100)
dd <- simBinFam(n,beta=0.3)
binfam <- fast.reshape(dd,varying=c("age","x","y"))
## mother, father, children (ordered)
head(binfam)
########### ########### ########### ########### ########### ###########
#### simple analyses of binomial family data
########### ########### ########### ########### ########### ###########
desfs <- function(x,num1="num1",num2="num2")
{
pp <- 1*(((x[num1]=="m")*(x[num2]=="f"))|(x[num1]=="f")*(x[num2]=="m"))
pc <- (x[num1]=="m" | x[num1]=="f")*(x[num2]=="b1" | x[num2]=="b2")*1
cc <- (x[num1]=="b1")*(x[num2]=="b1" | x[num2]=="b2")*1
c(pp,pc,cc)
}
ud <- easy.binomial.twostage(y~+1,data=binfam,
response="y",id="id",
score.method="fisher.scoring",deshelp=0,
theta.formula=desfs,desnames=c("pp","pc","cc"))
summary(ud)
udx <- easy.binomial.twostage(y~+x,data=binfam,
response="y",id="id",
score.method="fisher.scoring",
theta.formula=desfs,desnames=c("pp","pc","cc"))
summary(udx)
########### ########### ########### ########### ########### ###########
#### now allowing parent child POR to be different for mother and father
########### ########### ########### ########### ########### ###########
desfsi <- function(x,num1="num1",num2="num2")
{
pp <- (x[num1]=="m")*(x[num2]=="f")*1
mc <- (x[num1]=="m")*(x[num2]=="b1" | x[num2]=="b2")*1
fc <- (x[num1]=="f")*(x[num2]=="b1" | x[num2]=="b2")*1
cc <- (x[num1]=="b1")*(x[num2]=="b1" | x[num2]=="b2")*1
c(pp,mc,fc,cc)
}
udi <- easy.binomial.twostage(y~+1,data=binfam,
response="y",id="id",
score.method="fisher.scoring",
theta.formula=desfsi,desnames=c("pp","mother-child","father-child","cc"))
summary(udi)
##now looking to see if interactions with age or age influences marginal models
##converting factors to numeric to make all involved covariates numeric
##to use desfai2 rather then desfai that works on binfam
nbinfam <- binfam
nbinfam$num <- as.numeric(binfam$num)
head(nbinfam)
desfsai <- function(x,num1="num1",num2="num2")
{
pp <- (x[num1]=="m")*(x[num2]=="f")*1
### av age for pp=1 i.e parent pairs
agepp <- ((as.numeric(x["age1"])+as.numeric(x["age2"]))/2-30)*pp
mc <- (x[num1]=="m")*(x[num2]=="b1" | x[num2]=="b2")*1
fc <- (x[num1]=="f")*(x[num2]=="b1" | x[num2]=="b2")*1
cc <- (x[num1]=="b1")*(x[num2]=="b1" | x[num2]=="b2")*1
agecc <- ((as.numeric(x["age1"])+as.numeric(x["age2"]))/2-12)*cc
c(pp,agepp,mc,fc,cc,agecc)
}
desfsai2 <- function(x,num1="num1",num2="num2")
{
pp <- (x[num1]==1)*(x[num2]==2)*1
agepp <- (((x["age1"]+x["age2"]))/2-30)*pp ### av age for pp=1 i.e parent pairs
mc <- (x[num1]==1)*(x[num2]==3 | x[num2]==4)*1
fc <- (x[num1]==2)*(x[num2]==3 | x[num2]==4)*1
cc <- (x[num1]==3)*(x[num2]==3 | x[num2]==4)*1
agecc <- ((x["age1"]+x["age2"])/2-12)*cc ### av age for children
c(pp,agepp,mc,fc,cc,agecc)
}
udxai2 <- easy.binomial.twostage(y~+x+age,data=binfam,
response="y",id="id",
score.method="fisher.scoring",deshelp=0,detail=0,
theta.formula=desfsai,
desnames=c("pp","pp-age","mother-child","father-child","cc","cc-age"))
summary(udxai2)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.