# R/mycor_jgee2.R In JGEE: Joint Generalized Estimating Equation Solver

#### Documented in mycor_jgee2

```mycor_jgee2 <-
function(N,nr,nt,y,X,K,family,beta_new,corstr1,Mv,corstr2,maxclsz,R1,R2,scale.fix,scale.value) {

aindex=nr*cumsum(nt)
index=c(0,aindex[-length(aindex)])

xnewmat=matrix(0,sum(nt)*nr,(nr*(K+1)))
for(i in 1:N){
xnewmat[(1+index[i]):(nt[i]*nr+index[i]),]=kronecker(diag(nr),X[(1+(i-1)*nt[i]):(i*nt[i]),])
#cat((1+index[i]):(nt[i]*nr+index[i]),"\n")
}

eta=xnewmat%*%beta_new

if (is.character(family)) family <- get(family)
if (is.function(family))  family <- family()

mu=sapply(1:nr, function(r) family[[r]]\$linkinv(eta[as.vector(sapply(1:N, function(i) sapply(1:nt[i], function(j) j+index[i]+nt[i]*(r-1))))]))

muy=matrix(0,sum(nt)*nr,1)
yy=matrix(0,sum(nt)*nr,1)

for(i in 1:N) {
for(j in 1:nt[i]){
for(r in 1:nr){
yy[j+index[i]+nt[i]*(r-1)]=y[j+(i-1)*nt[i],r]
muy[j+index[i]+nt[i]*(r-1)]=mu[j+(i-1)*nt[i],r]
}
}
}

mu=muy

sd=sapply(1:nr, function(r) sqrt(family[[r]]\$variance(mu[as.vector(sapply(1:N, function(i) sapply(1:nt[i], function(j) j+index[i]+nt[i]*(r-1))))])))
sdy=matrix(0,sum(nt)*nr,1)
for(i in 1:N) {
for(j in 1:nt[i]){
for(r in 1:nr){
sdy[j+index[i]+nt[i]*(r-1)]=sd[j+(i-1)*nt[i],r]
}
}
}

sd=sdy

res<-(as.vector(yy)-mu)/sd

#cat("eta",eta,"yy",yy,"mu",mu,"sd",sd,"res",res,"\n")

if(scale.fix==0) {
fi<-sum(res^2)/(nr*sum(nt))
} else
if(scale.fix==1) {
fi<-scale.value
}

if (corstr1=="independence")
{alfa_hat<-0} else
if (corstr1=="exchangeable")
{
sum1<-0
sum3<-0
for ( i in  1:N)         {
for ( r in 1:nr)         {
for ( j in 1:nt[i])      {
for ( jj in 1:nt[i])     {
if  ( j!=jj)             {
#cat("i",i,"r",r,"j",j,"jj",jj,"\n")
sum2<-res[j+index[i]+nt[i]*(r-1)]*res[jj+index[i]+nt[i]*(r-1)]
sum1<-sum1+sum2
#cat("i",i,"r",r,"j",j,"jj",jj,"sum2",sum2,"sum1",sum1,"\n")
}
}
}
}
sum4<-nt[i]*(nt[i]-1)
sum3<-sum3+sum4
} #i
alfa_hat<-sum1/(nr*sum3*fi)
} else
if (corstr1=="AR-1")
{
sum5<-0
sum6<-0
for ( i in  1:N)           {
for ( r in  1:nr)          {
for ( j in  1:nt[i])       {
for ( jj in 1:nt[i])       {
if( j>jj && abs(j-jj)==1)  {
#cat("i",i,"r",r,"j",j,"jj",jj,"\n")
sum7<-res[j+index[i]+nt[i]*(r-1)]*res[jj+index[i]+nt[i]*(r-1)]
sum5<-sum5+sum7
#cat("i",i,"r",r,"j",j,"jj",jj,"sum7",sum7,"sum5", sum5, "\n")
}
}
}
}
sum8<-nt[i]-1
sum6<-sum6+sum8
}
alfa_hat<-sum5/(nr*sum6*fi)
}  else
if (corstr1=="stat_M_dep")
{
alfa_hat=matrix(0,Mv,1)
for(m in 1:Mv) {
sum12<-0
sum14<-0
for ( i in  1:N)           {
for ( r in  1:nr)          {
for ( j in  1:nt[i])       {
for ( jj in 1:nt[i])       {
if( j>jj && abs(j-jj)==m)  {
#cat("m",m,"i",i,"r",r,"j",j,"jj",jj,"\n")
sum11<-res[j+index[i]+nt[i]*(r-1)]*res[jj+index[i]+nt[i]*(r-1)]
sum12<-sum12+sum11
#cat("m",m,"i",i,"r",r,"j",j,"jj",jj,"sum11",sum11,"sum12", sum12, "\n")
} #if
}
}
}
sum13<-nt[i]-1
sum14<-sum14+sum13
}
alfa_hat[m]<-sum12/(nr*sum14*fi)
} #m
} else
if (corstr1=="non_stat_M_dep")
{
alfa_hat<-matrix(0,nt[1],nt[1]) #not allowed for unequal number of cluster sizes.
for( m in 1:Mv)            {
for ( j in  1:nt[1])       {
for ( jj in 1:nt[1])       {
if( j>jj && abs(j-jj)==m)  {
sum16<-0
for ( i  in 1:N)     {
for ( r in 1:nr)     {
#cat("m",m,"j",j,"jj",jj,"i",i,"r",r,"\n")
sum15<-res[j+index[1]+nt[i]*(r-1)]*res[jj+index[1]+nt[i]*(r-1)]
sum16<-sum15+sum16
#cat("j",j,"jj",jj,"i",i,"r",r,"sum15",sum15,"sum16",sum16,"\n")
}
}
#cat("j",j,"jj",jj,"sum16",sum16,"\n")
alfa_hat[j,jj]<-sum16/(nr*N*fi)
}
}
}
}
} else
if (corstr1=="unstructured")
{
alfa_hat<-matrix(0,nt[1],nt[1]) #not allowed for unequal number of cluster sizes.
for ( j in 1:nt[1])  {
for ( jj in 1:nt[1]) {
sum20<-0
if (j > jj)          {
for ( i  in 1:N )    {
for ( r in 1:nr)     {
#cat("i",i,"r",r,"j",j,"jj",jj,"\n")
sum21<-res[j+index[i]+nt[i]*(r-1)]*res[jj+index[i]+nt[i]*(r-1)]
sum20<-sum21+sum20
}
}
#cat("j",j,"jj",jj,"sum20",sum20,"\n")
alfa_hat[j,jj]<-sum20/(nr*N*fi)
}
}
}
} #end

if (corstr2=="independence")
{gamma_hat<-0} else
if (corstr2=="exchangeable")
{
sum30<-0
for ( i in 1:N)    {
for ( j in 1:nt[i]){
for ( r in 1:nr)   {
for ( rr in 1:nr)  {
if  ( r != rr)     {
#cat("i",i,"j",j,"r",r,"rr",rr,"\n")
sum31<-res[j+index[i]+nt[i]*(r-1)]*res[j+index[i]+nt[i]*(rr-1)]
sum30<-sum30+sum31
#cat("i",i,"j",j,"r",r,"rr",rr,"sum31",sum31,"sum32",sum32,"\n")
}
}
}
}
}
gamma_hat<-sum30/(nr*(nr-1)*sum(nt)*fi)
} else
if (corstr2=="unstructured")
{
gamma_hat<-matrix(0,nr,nr)
for ( r in 1:nr)   {
for ( rr in 1:nr)  {
if  ( r > rr)      {
sum40<-0
for ( i in 1:N)       {
for ( j in 1:nt[1])   {
sum41<-res[j+index[i]+nt[i]*(r-1)]*res[j+index[i]+nt[i]*(rr-1)]
#cat( "i",i,"j",j,"r",r,"rr",rr,"\n")
sum40<-sum40+sum41
#cat("i",i,"j",j,"r",r,"rr",rr,"sum41",sum41,"sum40",sum40,"\n")
}
}
gamma_hat[r,rr]<-sum40/(sum(nt)*fi)
}
}
}
} #

cor2<-matrix(0,nr,nr)
if (corstr2=="independence")
{cor2<-diag(nr) } else
if (corstr2=="exchangeable")
{ for (r1 in 1:nr) {
for (r2 in 1:nr) {
if (r1!=r2)
{cor2[r1,r2]<-gamma_hat} else
{cor2[r1,r2]<-1}
}
}
} else
if (corstr2=="unstructured")
{ cor2<-gamma_hat+t(gamma_hat)
diag(cor2)=1
}else
if (corstr2=="fixed")
{cor2=R2
}

#cat("alfa_hat",alfa_hat,"gamma_hat",gamma_hat,"\n")

C1hat<-array(0,c(maxclsz,maxclsz,N))
Ehat<-array(0,c(maxclsz*nr,maxclsz*nr,N))

for(i in 1:N){

cor1<-matrix(0,nt[i],nt[i])

if (corstr1=="independence")
{cor1<-diag(nt[i])} else
if (corstr1=="exchangeable")
{ for (t1 in 1:nt[i]) {
for (t2 in 1:nt[i]) {
if (t1!=t2)
{cor1[t1,t2]<-alfa_hat} else
{cor1[t1,t2]<-1}
}
}
} else
if (corstr1=="AR-1")
{ for (t1 in 1:nt[i]) {
for (t2 in 1:nt[i]) {
cor1[t1,t2]<-alfa_hat^abs(t1-t2)
}
}
}  else
if (corstr1=="stat_M_dep")
{ for (t1 in 1:nt[1]) {
for (t2 in 1:nt[1]) {
if (abs(t1-t2)==0)
{cor1[t1,t2]<-1} else
for(m in 1:Mv) {
if (abs(t1-t2)==m)
{cor1[t1,t2]<-alfa_hat[m]}
}
}
}
} else
if (corstr1=="non_stat_M_dep")
{
cor1=alfa_hat+t(alfa_hat)
diag(cor1)=1
} else
if ( corstr1=="unstructured")
{cor1=alfa_hat+t(alfa_hat)
diag(cor1)=1
} else
if (corstr1=="fixed")
{cor1=R1
}

Ehat[1:(nr*nt[i]),1:(nr*nt[i]),i]<-kronecker(cor2,cor1)
C1hat[1:nt[i],1:nt[i],i]<-cor1
}

return(list("Ehat"=Ehat,"fi"=fi,"cor1"=C1hat,"cor2"=cor2,"yy"=yy,"et"=eta,"mm"=mu))

}
```

## Try the JGEE package in your browser

Any scripts or data that you put into this service are public.

JGEE documentation built on May 2, 2019, 11:51 a.m.