jointfit {mdhglm}

Fitting Multivariate Double Hierarchical Generalized Linear Models using h-likelihood Approach

Description

The jointfit is used to fit a multivariate double hierarchical generalized linear models (MDHGLMs) allowing different models for multivariate response variables where each response follow DHGLM. A variety of distributions and link functions for both response and the random effects are allowed. Fixed and random effects can also be fitted in both the mean and the dispersion components. To call the fitting function jointfit, models for the mean and dispersion must be specified by DHGLMMODELING object preferably created by calling the DHGLMMODELING function.

Usage

jointfit(RespDist="gaussian",BinomialDen=NULL, DataMain, MeanModel,DispersionModel=NULL,
PhiFix=NULL,LamFix=NULL,structure="correlated",mord=0,dord=1,convergence=1e-05,
Init_Corr=NULL, EstimateCorrelations=TRUE, ZZCorr=NULL, factor=NULL, REML=TRUE,order=1)

Arguments

RespDist

The distribution of the response is set by the option RespDist. The user can set it to: "gaussian" (default), "binomial", "poisson", or "gamma".

BinomialDen

When RespDist="binomial", one should use the option BinomialDen to specify the denominator for the binomial distribution. This should be "NULL" (default) or a numeric vector of length equal to the length of DataMain. When specified as BinomialDen=NULL and RespDist="binomial", the denominator is 1.

DataMain

The option DataMain determines the data frame to be used (non-optional).

MeanModel

For the mean model, this option requires DGHLMMODELING object which should specified by the option Model="mean".

DispersionModel

For the overdispersion model, this option requires DGHLMMODELING object which should be specified by the option Model="dispersion".

PhiFix

The option for overdispersion parameters (phi) to be estimated or maintaned constant. Specifying defaults such as PhiFix =NULL implies that phi is to be estimated. If not, phi is fixed at a value specified by PhiFix.

LamFix

The option for random-effect variance (lambda) to be estimated or maintaned constant. Specifying defaults such as LamFix =NULL implies that lambda is to be estimated. If not, lambda is fixed at a value specified by LamFix.

structure

The option structure determines structure of random effects. When structure="correlated" (or "shared"), correlated (or "shared") random-effects model is specified. Factor analysis can be fitted by specifying structure="factor". Furthremore, selection model for missing data can be also fitted by specifying structure="selection".

mord

The option mord specifies the order of Laplace approximation to the marginal likelihood for fitting the mean parameters. The choice is either 0 or 1 (default).

dord

The option dord specifies the order of adjusted profile likelihood for fitting the dispersion parameters. The choice is either 1 (default) or 2.

convergence

Setting this option determines the criterion for convergence, which is computed as the absolute difference between the values of all the estimated parameters in the previous and current iterations. The default criterion is 1e-06.

Init_Corr

Setting initial values of correlation (or shared parameters) between random effects

EstimateCorrelations

Correlation are estimated or fixed when EstimateCorrelations=TRUE (default) or EstimateCorrelations=FALSE

ZZCorr

List of model matrices for random effects

factor

factor structure when structure="factor"

REML

Giving REML estimates when REML=TRUE (default) or ML estimates when REML=FALSE

order

first order approximation when order=1 (default) or second order approximation when order=2 for factor analysis

Examples

data(eg)
eg1<-eg[1:100,] ## using sampled data to have faster results 

jm1<-DHGLMMODELING(Link="identity", LinPred=y1~dose+dose2+(1|litter),RandDist="gaussian")
jm2<-DHGLMMODELING(Link="logit", LinPred=y2~dose+dose2+(1|litter),RandDist="gaussian")

Init_Corr=list(c(0))
SSC=list(as.factor(c(eg1$litter,eg1$litter)),as.factor(c(eg1$litter,eg1$litter)))
EstimateOverDisp=c(TRUE,FALSE)
LaplaceFixed=c(TRUE,TRUE)
ZZ1<-model.matrix(~as.factor(eg1$litter)-1)
ZZCorr=list(ZZ1,ZZ1)

#### independent random-effects model ####
res_ind<-jointfit(RespDist=c("gaussian","binomial"),DataMain=list(eg1,eg1),
    MeanModel=list(jm1,jm2),structure="correlated",
    Init_Corr=Init_Corr,EstimateCorrelations=FALSE,convergence=1,ZZCorr=ZZCorr)

#### correlated random-effects model ####
res_corr<-jointfit(RespDist=c("gaussian","binomial"),DataMain=list(eg1,eg1),
    MeanModel=list(jm1,jm2),structure="correlated",
    Init_Corr=Init_Corr,convergence=1,ZZCorr=ZZCorr)

#### shared random-effects model ####
Init_Corr=c(1,-10)
res_saturated<-jointfit(RespDist=c("gaussian","binomial"),DataMain=list(eg1,eg1),
    MeanModel=list(jm1,jm2),structure="shared",
    Init_Corr=Init_Corr,convergence=1)