Nothing
# Density of the univariate marginal distribution
# input:
# y the vector of (non-negative integer) quantiles.
# mu the mean parameter of the univariate distribution.
# gam the parameter gamma of the negative binomial distribution.
# output:
# the density of the univariate marginal distribution
dmargmodel.ord<-function(y,mu,gam,link)
{ cuts<-c(-10,gam,10)
lb=cuts[y]+mu
ub=cuts[y+1]+mu
if(link=="probit") res<-pnorm(ub)-pnorm(lb) else res<-plogis(ub)-plogis(lb)
res[y<1]<-0
res
}
# CDF of the univariate marginal distribution
# input:
# y the vector of (non-negative integer) quantiles.
# mu the mean parameter of the univariate distribution.
# gam the parameter gamma of the negative binomial distribution.
# output:
# the cdf of the univariate marginal distribution
pmargmodel.ord<-function(y,mu,gam,link)
{ cuts<-c(-10,gam,10)
ub=cuts[y+1]+mu # for mprobit
#ub=cuts[y+1]-mu # for polr
if(link=="probit") res<-pnorm(ub) else res<-plogis(ub)
res[y<1]<-0
res
}
# Bivariate composite likelihood for multivariate normal copula with ordinal regression.
# input:
# r the vector of normal copula parameters
# b the regression coefficients
# gam the gamma parameter
# xdat the matrix of covariates (use the constant 1 for the first covariate)
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# link is the link function. Choices are ?log? for the log link function,
# ?logit? for the logit link function, and ?probit? for the probit link function.
# corstr indicates the latent correlation structure of normal copula.
# Choices are exch, ar, and unstr for exchangeable, ar(1) and
# unstrucutred correlation structure, respectively.
# output:
# negative bivariate composite likelihood for multivariate normal copula
# with ordinal regression.
bcl.ord<-function(r,b,gam,xdat,ydat,id,tvec,corstr,link)
{ s<-0
uid<-unique(id)
d<-id.size(id)
maxd<-max(d)
tvec<-id.time(tvec,d)
n<-1:length(ydat) #### change here
pairs<-maxpairs(d)
rmat<-cormat(maxd,r,pairs,corstr)
c1<-(r[1] > 1 | r[1] < (-1/(maxd-1)))
c2<-(r[1] > 1 | r[1] < -1)
c3<-(det(rmat)<0 | min(rmat)< -1 | max(rmat)>1)
if((corstr=="exch" & c1) | (corstr=="ar" & c2 ) | (corstr=="unstr" & c3)) {s<-1e10}
else {
for(i in uid)
{ cases<-id==i
irow=n[cases]
yi<-ydat[irow]
ti<-tvec[irow]
newyi<-rep(NA,maxd)
newmui<-rep(NA,maxd)
newyi[ti]<-yi
if(is.vector(xdat)) x<-xdat[cases] else x<-xdat[cases,] ###change here
mui<-ordreg.mu(x,b)
newmui[ti]<-mui
vlow<-pmargmodel.ord(newyi-1,newmui,gam,link)
tem<-dmargmodel.ord(newyi,newmui,gam,link)
vupp<-vlow+tem
zlow=qnorm(vlow)
zupp=qnorm(vupp)
for(j in 1:dim(pairs)[1])
{ k1<-pairs[j,][1]
k2<-pairs[j,][2]
if((sum(k1==ti)==1) & (sum(k2==ti)==1))
{ if(corstr=="exch")
{ s<-s +bivlik(zlow[c(k1,k2)],zupp[c(k1,k2)],r) }
else
{ if(corstr=="ar")
{ s<-s +bivlik(zlow[c(k1,k2)],zupp[c(k1,k2)],r^(k2-k1)) }
else
{ # corstr=="unstr"
s<-s +bivlik(zlow[c(k1,k2)],zupp[c(k1,k2)],r[j])
}}}
else { s<-s }
}
}
-s
}
}
# optimization routine for composite likelihood for MVN copula
# b the regression coefficients
# gam the gamma parameter
# xdat the matrix of covariates (use the constant 1 for the first covariate)
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# link is the link function. Choices are
# ?logit? for the logit link function, and ?probit? for the probit link function.
# corstr indicates the latent correlation structure of normal copula.
# Choices are ?exch?, ?ar?, and ?unstr? for exchangeable, ar(1) and
# unstrucutred correlation structure, respectively.
# output: A list containing the following components:
# minimum the value of the estimated minimum of CL1 for MVN copula
# estimate the CL1 estimates
# gradient the gradient at the estimated minimum of CL1
# code an integer indicating why the optimization process terminated, see nlm.
cl1.ord<-function(b,gam,xdat,ydat,id,tvec,corstr,link)
{ d<-id.size(id)
pairs<-maxpairs(d)
if(corstr=="unstr")
{ nom<-dim(pairs)[1]
nlm(bcl.ord,rep(0.1,nom),b,gam,xdat,ydat,id,tvec,corstr,link)
}
else
{ nlm(bcl.ord,0.1,b,gam,xdat,ydat,id,tvec,corstr,link) }
}
# derivative of the ordinal loglikelihood with respect to gamma
# input:
# gam the gamma parameter
# mu the mean parameter
# output:
# the vector with the derivatives of the NB loglikelihood with respect to gamma
derlik.gam.ord<-function(mu,gam,u,link)
{ K<-length(gam)+1
k<-1:K
cuts<-c(-10,gam,10)
lb=cuts[k]+mu
ub=cuts[k+1]+mu
if(link=="probit") { dlatent=dnorm; platent=pnorm } else { dlatent=dlogis; platent=plogis }
dlatentub<-dlatent(ub)
dlatentlb<-dlatent(lb)
den<-platent(ub)-platent(lb)
res<-rep(NA,K)
for(i in 1:K)
{ if(u==i)
{ res[i]=dlatentub[i]/den[i] }
else
{ if(u==i-1)
{ res[i]=-dlatentlb[i]/den[i] }
else {res[i]=0}}
}
res
}
# derivative of the NB loglikelihood with respect to gamma
# input:
# gam the gamma parameter
# mu the mean parameter
# y the value of a non-negative integer quantile
# output:
# the derivative of the NB loglikelihood with respect to gamma
iderlik.gam.ord<-function(mu,gam,y,u,link)
{ cuts<-c(-10,gam,10)
lb=cuts[y]+mu
ub=cuts[y+1]+mu
if(link=="probit") { dlatent=dnorm; platent=pnorm } else { dlatent=dlogis; platent=plogis }
den<-platent(ub)-platent(lb)
if(u==y) dlatent(ub)/den
else if(u==y-1) -dlatent(lb)/den else 0
}
der.dnorm<-function(x)
{ -x*dnorm(x) }
der.dlogis<-function(x)
{ expx=exp(x)
expx*(1-expx)/(1+expx)^3
}
# minus expectation of the second derivative of the marginal ordinal loglikelihood
# with resect to gamma
# input:
# mu the mean parameter
# gam the gamma parameter
# u the univariate cdfs
# output:
# the vector with the minus expectations of the margmodel loglikelihood
# with respect to gamma
fisher.gam.ord<-function(mu,gam,u,v,link)
{ cuts<-c(-10,gam,10)
K<-length(gam)+1
k<-1:K
lb=cuts[k]+mu
ub=cuts[k+1]+mu
if(link=="probit") { dlatent=dnorm; platent=pnorm; der.dlatent=der.dnorm } else {
dlatent=dlogis; platent=plogis; der.dlatent=der.dlogis }
den<-platent(ub)-platent(lb)
dlatentub<-dlatent(ub)
dlatentlb<-dlatent(lb)
der.dlatentub<-der.dlatent(ub)
der.dlatentlb<-der.dlatent(lb)
h<-rep(NA,K)
for(k in 1:K)
{ if(u==k & v==k)
{ num1<-der.dlatentub[k]
num2<-dlatentub[k]
tem<-num2/den[k]
h[k]<--num1/den[k]+tem*tem
}
else
{ if((u==k & v==k-1) | (u==k-1 & v==k))
{ h[k]<--dlatentub[k]*dlatentlb[k]/den[k]/den[k]
}
else
{ if(u==k-1 & v==k-1)
{ num1<-der.dlatentlb[k]
num2<-dlatentlb[k]
tem<-num2/den[k]
h[k]<-num1/den[k]+tem*tem
} else h[k]<-0}}
}
sum(h*den)
}
# the mean values of the univariate marginal distribution
# corresonding to the used link function
# input:
# x the matix of the covariates
# b the vector with the regression coefficients
# output:
# the mean values of the univariate marginal distribution
ordreg.mu<-function(x,b)
{ if(length(b)!=1) mu<-x %*% b else mu<-x*b }
# select the present column and lines in Omega, Delta and X matrices
# for unbalanced data
# input:
# tvec a vector of the time for an individual
subselect.ord<-function(tvec,q)
{ sel<-NULL
for(i in 1:length(tvec))
{ tm<-tvec[i]
k<-((tm-1)*q+1):((tm-1)*q+q)
sel<-c(sel,k)
}
sel
}
# Calculating the number of categories
# input:
# gam the cutpoints
# output:
# the number of categories
noCategories<-function(gam)
{ length(gam)+1 }
# covariance matrix of the scores Omega_i
# input:
# scgam the array of the score functions with respect to gam
# index the bivariate pair
# pmf the matrix of rectangle probabilities
scoreCov.ord<-function(scgam,pmf,index)
{ j1<-index[1]
j2<-index[2]
q<-dim(scgam)[2]
cov22<-matrix(NA,q,q)
for(i1 in 1:q)
{ for(i2 in 1:q)
{ cov22[i1,i2]<-t(scgam[,i1,j1])%*%pmf%*%scgam[,i2,j2] }
}
cov22
}
# weight matrix fixed at values from the CL1 estimator
# input:
# b the regression coefficients
# gam the gamma parameter
# xdat the matrix of covariates (use the constant 1 for the first covariate)
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# rh the vector with CL1 estimates
# link is the link function. Choices are
# ?logit? for the logit link function, and ?probit? for the probit link function.
# corstr indicates the latent correlation structure of normal copula.
# Choices are ?exch?, ?ar?, and ?unstr? for exchangeable, ar(1) and
# unstrucutred correlation structure, respectively.
# output: A list containing the following components:
# omega the array with the Omega matrices
# delta the array with the Delta matrices
# X the array with the X matrices
weightMat.ord<-function(b,gam,rh,xdat,ydat,id,tvec,corstr,link)
{ uid<-unique(id)
d<-id.size(id)
maxd<-max(d)
q<-length(gam)
lid<-length(uid)
if(is.matrix(xdat))
{ dim<-dim(xdat)
n<-1:dim[1]
p<-dim[2]
} else {n<-1:length(xdat); p<-1}
pairs<-maxpairs(d)
tvec<-id.time(tvec,d)
omega<-array(NA,c(maxd*q,maxd*q,lid))
X<-array(NA,c(p+q,maxd*q,lid))
delta<-array(NA,c(maxd*q,maxd*q,lid))
dom<-maxd*q
if(q>1) pos<-seq(1,dom-1,by=q) else pos<-seq(1,dom) #for binary
m<-0
for(i in uid)
{ m<-m+1
cases<-id==i
irow=n[cases]
newyi<-rep(NA,maxd)
newmui<-rep(NA,maxd)
ti<-tvec[irow]
if(is.matrix(xdat))
{ x<-xdat[cases,]
newx<-matrix(NA,maxd,p)
newx[ti,]<-x} else {
x<-xdat[cases]
newx<-rep(NA,maxd)
newx[ti]<-x }
mui<-ordreg.mu(x,b)
newmui[ti]<-mui
ub<-noCategories(gam)
du<-matrix(NA,ub,maxd)
scgam<-array(NA,c(ub,ub-1,maxd))
for(j in ti)
{ du[,j]<-dmargmodel.ord(1:ub,newmui[j],gam,link)
for(k in 1:(ub-1))
{ scgam[,k,j]<-derlik.gam.ord(newmui[j],gam,k,link) }
}
u<-apply(du,2,cumsum)
z<-qnorm(u)
z[is.nan(z)]<-7
z[z>10]<-7
z<-rbind(-7,z)
pz<-pnorm(z)
dz<-dnorm(z)
zs<-z*z
zc<-z*zs
zf<-zs*zs
xi<-NULL
diagonali<-array(NA,c(q,q,maxd))
for(j in 1:maxd)
{ if(is.matrix(newx))
{ temp1<-matrix(rep(newx[j,],each=q),q) } else {
temp1<-matrix(rep(newx[j],each=q),q) }
xi<-cbind(xi,t(cbind(temp1,diag(q))))
fisher<-matrix(NA,ub-1,ub-1)
for(k1 in 1:(ub-1))
{ for(k2 in 1:(ub-1))
{ fisher[k1,k2]<-fisher.gam.ord(newmui[j],gam,k1,k2,link)
}
}
diagonali[,,j]<-fisher
}
deltai<-matrix(0,dom,dom)
minus<-0
for(j in pos)
{ deltai[j:(j+q-1),j:(j+q-1)]<-diagonali[,,j-minus]
if(q>1) minus<-minus+q-1 else minus<-0
}
offi<-array(NA,c(q,q,dim(pairs)[1]))
for(k in 1:dim(pairs)[1])
{ k1<-pairs[k,][1]
k2<-pairs[k,][2]
if((sum(k1==ti)==1) & (sum(k2==ti)==1))
{ x1=z[,k1]; x2=z[,k2]
x1s=zs[,k1]; x2s=zs[,k2]
x1c=zc[,k1]; x1f=zf[,k1]
x2c=zc[,k2]; x2f=zf[,k2]
t1=outer(pz[,k1],pz[,k2])
t2=outer(dz[,k1],dz[,k2])
if(corstr=="exch")
{ cdf<-approxbvncdf(rh,x1,x2,x1s,x2s,x1c,x2c,x1f,x2f,t1,t2) } else {
if(corstr=="ar")
{ cdf<-approxbvncdf(rh^(k2-k1),x1,x2,x1s,x2s,x1c,x2c,x1f,x2f,t1,t2) } else { # corstr=="unstr"
cdf<-approxbvncdf(rh[k],x1,x2,x1s,x2s,x1c,x2c,x1f,x2f,t1,t2)
}
}
cdf1=apply(cdf,2,diff)
pmf=apply(t(cdf1),2,diff)
pmf=t(pmf)
offi[,,k]<-scoreCov.ord(scgam,pmf,pairs[k,])
}
}
omegai<-deltai
ch1<-0
ch2<-0
ch3<-0
if(d[m]>1)
{
for(j in 1:(d[m]-1))
{ for(r in pos[-(1:j)])
{ #print(c(j,r))
omegai[(1+(j-1)*q):(j*q),r:(r+q-1)]<-offi[,,(j+ch2-ch1)]
omegai[r:(r+q-1),(1+(j-1)*q):(j*q)]<-t(offi[,,(j+ch2-ch1)])
ch2<-ch2+1
}
ch1<-ch1+1
}
}
X[,,m]<-xi
delta[,,m]<-deltai
omega[,,m]<-omegai
}
list(omega=omega,X=X,delta=delta)
}
# the weigted scores equations
# input:
# param the vector of regression and not regression parameters
# xdat the matrix of covariates
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# link is the link function. Choices are
# ?logit? for the logit link function, and ?probit? for the probit link function.
# corstr indicates the latent correlation structure of normal copula.
# Choices are ?exch?, ?ar?, and ?unstr? for exchangeable, ar(1) and
# unstrucutred correlation structure, respectively.
# WtScMat is a list containing the following components:
# omega the array with the Omega matrices
# delta the array with the Delta matrices
# X the array with the X matrices
# output
# the weigted scores equations
wtsc.ord<-function(param,WtScMat,xdat,ydat,id,tvec,link)
{ uid<-unique(id)
d<-id.size(id)
maxd<-max(d)
if(is.matrix(xdat))
{ dim<-dim(xdat)
n<-1:dim[1]
p<-dim[2]
} else {n<-1:length(xdat); p<-1}
tvec<-id.time(tvec,d)
b<-param[1:p]
q<-length(unique(ydat))-1
gam<-param[(p+1):(p+q)]
g<-0
m<-0
for(i in uid)
{ #print(i)
cases<-id==i
irow=n[cases]
m<-m+1
if(is.matrix(xdat)) x<-xdat[cases,] else x<-xdat[cases]
mu<-ordreg.mu(x,b)
ub<-noCategories(gam)
scgam<-array(NA,c(ub,ub-1,d[m]))
for(j in 1:d[m])
{ for(k in 1:(ub-1))
{ scgam[,k,j]<-derlik.gam.ord(mu[j],gam,k,link) }
}
y<-ydat[irow]
sci<-NULL
for(j in 1:d[m])
{ scgami<-NULL
for(k in 1:(ub-1))
{ scgami<-c(scgami,scgam[y[j],k,j])}
sci<-c(sci,scgami)
}
ti<-tvec[irow]
seli<-subselect.ord(ti,q)
Xi<-WtScMat$X[,,m]
Xi<-Xi[,seli]
deltai<-WtScMat$delta[,,m]
deltai<-deltai[seli,seli]
omegai<-WtScMat$omega[,,m]
omegai<-omegai[seli,seli]
gi<-Xi%*%t(deltai)%*%solve(omegai,sci)
g<-g+gi
}
g
}
# solving the weigted scores equations
# input:
# start the starting values (IEE estimates) for the vector of
# regression and not regression parameters
# xdat the matrix of covariates (use the constant 1 for the first covariate)
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# link is the link function. Choices are
# ?logit? for the logit link function, and ?probit? for the probit link function.
# corstr indicates the latent correlation structure of normal copula.
# Choices are ?exch?, ?ar?, and ?unstr? for exchangeable, ar(1) and
# unstrucutred correlation structure, respectively.
# WtScMat is a list containing the following components:
# omega the array with the Omega matrices
# delta the array with the Delta matrices
# X the array with the X matrices
# output:
# the weighted scores estimates
solvewtsc.ord<-function(start,WtScMat,xdat,ydat,id,tvec,link)
{ multiroot(f=wtsc.ord,start,atol=1e-4,rtol=1e-4,ctol=1e-4,
WtScMat=WtScMat,xdat=xdat,ydat=ydat,id=id,tvec=tvec,link=link)
}
# inverse Godambe matrix with Delta and Omega evaluated at IEE estimator
# input:
# b the regression coefficients
# gam the gamma parameter
# xdat the matrix of covariates (use the constant 1 for the first covariate)
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# rh the vector with CL1 estimates
# WtScMat a list containing the following components:
# omega the array with the Omega matrices
# delta the array with the Delta matrices
# X the array with the X matrices
# link is the link function. Choices are
# ?logit? for the logit link function, and ?probit? for the probit link function.
# output:
# the inverse Godambe matrix
godambe.ord<-function(param,WtScMat,xdat,ydat,id,tvec,link)
{ uid<-unique(id)
d<-id.size(id)
maxd<-max(d)
if(is.matrix(xdat))
{ dim<-dim(xdat)
n<-1:dim[1]
p<-dim[2]
} else {n<-1:length(xdat); p<-1}
tvec<-id.time(tvec,d)
b<-param[1:p]
q<-length(unique(ydat))-1
gam<-param[(p+1):(p+q)]
v<-v1<-v2<-v3<-0
fv<-fv1<-fv2<-fv3<-0
m<-0
for(i in uid)
{ m<-m+1
cases<-id==i
irow=n[cases]
ti<-tvec[irow]
if(is.matrix(xdat))
{ x<-xdat[cases,]
newx<-matrix(NA,maxd,p)
newx[ti,]<-x} else {
x<-xdat[cases]
newx<-rep(NA,maxd)
newx[ti]<-x }
newy<-rep(NA,maxd)
newmui<-rep(NA,maxd)
mui<-ordreg.mu(x,b)
newmui[ti]<-mui
ub<-noCategories(gam)
scgam<-array(NA,c(ub,ub-1,maxd))
for(j in ti)
{ for(k in 1:(ub-1))
{ scgam[,k,j]<-derlik.gam.ord(newmui[j],gam,k,link) }
}
y<-ydat[irow]
newy[ti]<-y
sci<-NULL
for(j in 1:maxd)
{ scgami<-NULL
for(k in 1:(ub-1))
{ scgami<-c(scgami,scgam[newy[j],k,j])}
sci<-c(sci,scgami)
}
sci<-sci[!is.na(sci)]
seli<-subselect.ord(ti,q)
Xi<-WtScMat$X[,,m]
Xi<-Xi[,seli]
fdeltai<-WtScMat$delta[,,m]
fdeltai<-fdeltai[seli,seli]
fomegai<-WtScMat$omega[,,m]
fomegai<-fomegai[seli,seli]
fci<-fdeltai%*%t(Xi)
nrhs1=ncol(fdeltai)
nrhs2=ncol(fci)
tem.lin=solve(fomegai,cbind(fdeltai,fci))
finvai=t(tem.lin[,1:nrhs1])
xi.invai=Xi %*% finvai
fai<-solve(finvai)
tem=solve(t(fai),t(Xi))
fv1i<-xi.invai %*% fci
fv2i<-xi.invai %*% (sci%*% t(sci)) %*% tem
fv3i<-t(fci) %*% tem
fv1<-fv1+fv1i
fv2<-fv2+fv2i
fv3<-fv3+fv3i
}
solve(fv1,fv2)%*%solve(fv3)
}
# the weighted scores wrapper function: handles all the steps in the weighted scores
# xdat the matrix of covariates (use the constant 1 for the first covariate)
# ydat the vector with the response
# id the the vector with the id
# tvec the time related vector
# rh the vector with CL1 estimates
# WtScMat a list containing the following components:
# omega the array with the Omega matrices
# delta the array with the Delta matrices
# X the array with the X matrices
# corstr indicates the latent correlation structure of normal copula.
# Choices are ?exch?, ?ar?, and ?unstr? for exchangeable, ar(1) and
# unstrucutred correlation structure, respectively.
# link is the link function. Choices are
# ?logit? for the logit link function, and ?probit? for the probit link function.
# iprint indicates printing of some intermediate results, default FALSE
wtsc.ord.wrapper<-function(xdat,ydat,id,tvec,corstr,link,iprint=FALSE)
{
i.est<-iee.ord(xdat,ydat,link)
if(iprint)
{ cat("\niest: IEE estimates\n")
print(c(i.est$reg,i.est$gam))
}
est.rho<-cl1.ord(b=i.est$reg,gam=i.est$gam,xdat,ydat,id,tvec,corstr,link)
if(iprint)
{ cat("\nest.rho: CL1 estimates\n")
print(est.rho$e)
cat("\nest.rho: CL1 likelihood\n")
print(-est.rho$m)
}
WtScMat<-weightMat.ord(b=i.est$reg,gam=i.est$gam,rh=est.rho$e,
xdat,ydat,id,tvec,corstr,link)
ws<-solvewtsc.ord(start=c(i.est$reg,i.est$gam),WtScMat,xdat,ydat,id,
tvec,link)
if(iprint)
{ cat("ws=parameter estimates\n")
print(ws$r)
}
acov<-godambe.ord(ws$r,WtScMat,xdat,ydat,id,tvec,link)
se<-sqrt(diag(acov))
if(iprint)
{ cat("\nacov: inverse Godambe matrix with W based on first-stage wt matrices\n")
print(acov)
cat("\nse: robust standard errors\n")
print(se)
res<-round(cbind(ws$r,se),3)
cat("\nres: Weighted scores estimates and standard errors\n")
print(res)
}
list(IEEest=c(i.est$reg,i.est$gam),CL1est=est.rho$e,CL1lik=-est.rho$m,WSest=ws$r, asympcov=acov)
}
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.