Nothing
#' @title Variance-covariance matrix of a CUBE model with covariates
#' @description Compute the variance-covariance matrix of parameter estimates of a CUBE model with covariates
#' for all the three parameters.
#' @aliases varcovcubecov
#' @usage varcovcubecov(m, ordinal, Y, W, Z, estbet, estgama, estalpha)
#' @param m Number of ordinal categories
#' @param ordinal Vector of ordinal responses
#' @param Y Matrix of covariates for explaining the uncertainty component
#' @param W Matrix of covariates for explaining the feeling component
#' @param Z Matrix of covariates for explaining the overdispersion component
#' @param estbet Vector of the estimated parameters for the uncertainty component, with length equal to
#' NCOL(Y)+1 to account for an intercept term (first entry)
#' @param estgama Vector of the estimated parameters for the feeling component, with length equal to
#' NCOL(W)+1 to account for an intercept term (first entry)
#' @param estalpha Vector of the estimated parameters for the overdispersion component, with length
#' equal to NCOL(Z)+1 to account for an intercept term (first entry)
#' @details The function checks if the variance-covariance matrix is positive-definite: if not,
#' it returns a warning message and produces a matrix with NA entries.
#' @keywords internal
#' @references
#' Piccolo, D. (2014), Inferential issues on CUBE models with covariates,
#' \emph{Communications in Statistics - Theory and Methods}, \bold{44}, DOI: 10.1080/03610926.2013.821487
varcovcubecov <-
function(m,ordinal,Y,W,Z,estbet,estgama,estalpha){
n<-length(ordinal)
Y<-as.matrix(Y); W<-as.matrix(W);Z<-as.matrix(Z);
if (ncol(W)==1){
W<-as.numeric(W)
}
if (ncol(Y)==1){
Y<-as.numeric(Y)
}
if (ncol(Z)==1){
Z<-as.numeric(Z)
}
p<-NCOL(Y);q<-NCOL(W);v<-NCOL(Z);
paivett<-logis(Y,estbet)
csivett<-logis(W,estgama)
phivett<-1/(-1+ 1/(logis(Z,estalpha)))
# Probability
probi<-paivett*(betabinomial(m,ordinal,csivett,phivett)-1/m)+1/m
uui<-1-1/(m*probi)
ubari<-uui+paivett*(1-uui)
### Matrix computations
mats1<-auxmat(m,csivett,phivett,1,-1,1,0,1)
mats2<-auxmat(m,csivett,phivett,0, 1,1,0,1)
mats3<-auxmat(m,csivett,phivett,1,-1,1,1,1)
mats4<-auxmat(m,csivett,phivett,0, 1,1,1,1)
mats5<-auxmat(m,csivett,phivett,1, 0,1,1,1)
### for D_i vectors
matd1<-auxmat(m,csivett,phivett,1,-1,2,0, 1)
matd2<-auxmat(m,csivett,phivett,0, 1,2,0,-1)
matd3<-auxmat(m,csivett,phivett,1,-1,2,1, 1)
matd4<-auxmat(m,csivett,phivett,0, 1,2,1,-1)
### for H_i vectors
math3<-auxmat(m,csivett,phivett,1,-1,2,2,-1)
math4<-auxmat(m,csivett,phivett,0, 1,2,2,-1)
math5<-auxmat(m,csivett,phivett,1, 0,2,2,-1)
###
S1<-S2<-S3<-S4<-D1<-D2<-D3<-D4<-H3<-H4<-rep(NA,m);
for(i in 1:n){
S1[i]<-sum(mats1[1:ordinal[i],i])-mats1[ordinal[i],i]
S2[i]<-sum(mats2[1:(m-ordinal[i]+1),i])-mats2[(m-ordinal[i]+1),i]
S3[i]<-sum(mats3[1:ordinal[i],i])-mats3[ordinal[i],i]
S4[i]<-sum(mats4[1:(m-ordinal[i]+1),i])-mats4[(m-ordinal[i]+1),i]
D1[i]<-sum(matd1[1:ordinal[i],i])-matd1[ordinal[i],i]
D2[i]<-sum(matd2[1:(m-ordinal[i]+1),i])-matd2[(m-ordinal[i]+1),i]
D3[i]<-sum(matd3[1:ordinal[i],i])-matd3[ordinal[i],i]
D4[i]<-sum(matd4[1:(m-ordinal[i]+1),i])-matd4[(m-ordinal[i]+1),i]
H3[i]<-sum(math3[1:ordinal[i],i])-math3[ordinal[i],i]
H4[i]<-sum(math4[1:(m-ordinal[i]+1),i])-math4[(m-ordinal[i]+1),i]
}
###### Attention !!!
S5<-colSums(mats5[1:(m-1),])
H5<-colSums(math5[1:(m-1),])
##########
CC<-S2-S1; EE<-S3+S4-S5;
DD<-D2-D1; FF<-D3+D4; GG<-H3+H4-H5;
### Computing vectors v and u
vibe<-uui*(1-paivett)
viga<-ubari*csivett*(1-csivett)*CC
vial<-ubari*phivett*EE #****
ubebe<-uui*(1-paivett)*(1-2*paivett)
ugabe<-ubari*csivett*(1-csivett)*(1-paivett)*CC
ualbe<-ubari*phivett*(1-paivett)*EE #****
ugaga<-ubari*csivett*(1-csivett)*((1-2*csivett)*CC+csivett*(1-csivett)*(CC^2+DD))
ualga<-ubari*phivett*csivett*(1-csivett)*(FF+CC*EE) #****
ualal<-ubari*phivett*(EE+phivett*(EE^2+GG)) #****
### Computing vectors g
gbebe<-Hadprod(vibe,vibe)-ubebe
ggabe<-Hadprod(viga,vibe)-ugabe
galbe<-Hadprod(vial,vibe)-ualbe
ggaga<-Hadprod(viga,viga)-ugaga
galga<-Hadprod(vial,viga)-ualga
galal<-Hadprod(vial,vial)-ualal
### Expanding matrices Y, W, Z
YY<-cbind(1,Y); WW<-cbind(1,W); ZZ<-cbind(1,Z);
### Elements of the Information matrix
infbebe<-t(YY)%*%(Hadprod(YY,gbebe))
infgabe<-t(WW)%*%(Hadprod(YY,ggabe))
infalbe<-t(ZZ)%*%(Hadprod(YY,galbe))
infgaga<-t(WW)%*%(Hadprod(WW,ggaga))
infalga<-t(ZZ)%*%(Hadprod(WW,galga))
infalal<-t(ZZ)%*%(Hadprod(ZZ,galal))
infbega<-t(infgabe)
infbeal<-t(infalbe)
infgaal<-t(infalga)
### Assembling Information matrix...
# # matinf<-rbind(
# cbind(infbebe,infbega,infbeal),
# cbind(infgabe,infgaga,infgaal),
# cbind(infalbe,infalga,infalal));
### Variance-covariance matrix
npai<-NCOL(Y)+1
ncsi<-NCOL(W)+1
nphi<-NCOL(Z)+1
nparam<-npai+ncsi+nphi
matinf<-matrix(NA,nrow=nparam,ncol=nparam)
for (i in 1:npai){
matinf[i,] <- matrix(c(infbebe[i,],infbega[i,],infbeal[i,]),nrow=1,byrow=TRUE)
}
for (i in (npai+1):(npai+ncsi)){
matinf[i,] <- matrix(c(infgabe[i-npai,],infgaga[i-npai,],infgaal[i-npai,]),nrow=1,byrow=TRUE)
}
for (i in (npai+ncsi+1):nparam){
matinf[i,] <- matrix(c(infalbe[i-npai-ncsi,],infalga[i-npai-ncsi,],infalal[i-npai-ncsi,]),nrow=1,byrow=TRUE)
}
if(any(is.na(matinf))==TRUE){
warning("ATTENTION: NAs produced")
varmat<-matrix(NA,nrow=nparam,ncol=nparam)
} else {
if(det(matinf)<=0){
warning("ATTENTION: Variance-covariance matrix NOT positive definite")
varmat<-matrix(NA,nrow=nparam,ncol=nparam)
} else {
varmat<-solve(matinf)
}
}
return(varmat)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.