The “R” codes | Ehsan Bahrami Samani

The “R” codes used for Analysis of mixed longitudinal ($k$,$l$)-Inflated power series, ordinal and continuous responses

F.Raziee, EhsanBahrami Samani, M. Ganjali ,(2017)

These Program present analysis of mixed longitudinal ($k$,$l$)-Inflated power series, ordinal and continuous responses with sensitivity analysis to non-ignorable missing mechanism. A simulation study is performed in which for count response ($k$,$l$)-inflated Poisson and ($k$,$l$)-inflated negative binomial distributions are considered

## Simulation Staudy for (0,1)-inflated Poisson and normal model

library(MASS)
library(nlme)

v=1000
vv=1000

b1=matrix(0,v,3)
b=array(0,dim=c(vv,3,v))

x11<-x12<-vector(“numeric”,v)
for(i in 1:v){

x11[i]=rnorm(1,0,1)
x12[i]=rnorm(1,0,1)

b1[i,]=rnorm(3,0,1)
}
for(i in 1:v){for(j in 1:vv){
b[j,,i]=rnorm(3,0,1)
}}
c1=c2=mc1=mc2=ppp1=ppp2=pp1=pp2=p1=p2=q1=q2=y1=y2=n1=nn1=n2=nn2=lla1=lla2=mm1=mm2=yo1=yo2=vector(“numeric”,v)

for(i in 1:v){

mm1[i]=x11[i]+b1[i,1]
mm2[i]=x12[i]+b1[i,1]

y1[i]=rnorm(1,mm1[i],1)
y2[i]=rnorm(1,mm2[i],1)

mc1[i]=x11[i]+b1[i,1]
mc2[i]=x12[i]+b1[i,1]

c1[i]=rnorm(1,mc1[i],1)
c2[i]=rnorm(1,mc2[i],1)

if(y1[i]<=-1){yo1[i]<-1}
if(-1<y1[i]&y1[i]<=1){yo1[i]<-2}
if(y1[i]>1){yo1[i]<-3}

if(y2[i]<=-1){yo2[i]<-1}
if(-1<y2[i]&y2[i]<=1){yo2[i]<-2}
if(y2[i]>1){yo2[i]<-3}

lla1[i]=exp(x11[i]+b1[i,1])
lla2[i]=exp(x12[i]+b1[i,1])

p1[i]<-exp(x11[i]+b1[i,2])/(1+exp(x11[i]+b1[i,2]))
p2[i]<-exp(x12[i]+b1[i,2])/(1+exp(x12[i]+b1[i,2]))
q1[i]<-exp(x11[i]+b1[i,3])/(1+exp(x11[i]+b1[i,3]))
q2[i]<-exp(x12[i]+b1[i,3])/(1+exp(x12[i]+b1[i,3]))

pp1[i]=(1-p1[i])*q1[i]
pp2[i]=(1-p2[i])*q2[i]

ppp1[i]=(1-p1[i])*(1-q1[i])
ppp2[i]=(1-p2[i])*(1-q2[i])

nn1[i]=rpois(1,lla1[i])
nn2[i]=rpois(1,lla2[i])

n1[i]<-sample(c(0,1,nn1[i]),1,prob=c(p1[i],pp1[i],ppp1[i]))
n2[i]<-sample(c(0,1,nn2[i]),1,prob=c(p2[i],pp2[i],ppp2[i]))
}
cbind(c1,c2,n1,n2,yo1,yo2)

################################################################################################
mcc1=mcc2=q11=q22=p11=pp11=ppp11=p22=pp22=ppp22=m1=m2=la1=la2=l1=l2=ll1=ll2=L1=lc1=lc2=matrix(0,v,vv)

L=vector(“numeric”,v)
L=L2=L4=L5=L6=vector(“numeric”,v)
aa=rep(0,3)
f<-function(p){

for(i in 1:v){for(j in 1:vv){
m1[i,j]=p[1]*x11[i]+p[11]*b[j,1,i]
m2[i,j]=p[2]*x12[i]+p[11]*b[j,1,i]

la1[i,j]=exp(p[3]*x11[i]+p[12]*b[j,1,i])
la2[i,j]=exp(p[4]*x12[i]+p[12]*b[j,1,i])

p11[i,j]=exp(p[5]*x11[i]+p[13]*b[j,2,i])/(1+exp(p[5]*x11[i]+p[13]*b[j,2,i]))
p22[i,j]=exp(p[6]*x12[i]+p[13]*b[j,2,i])/(1+exp(p[6]*x12[i]+p[13]*b[j,2,i]))
q11[i,j]=exp(p[7]*x11[i]+p[14]*b[j,3,i])/(1+exp(p[7]*x11[i]+p[14]*b[j,3,i]))
q22[i,j]=exp(p[8]*x12[i]+p[14]*b[j,3,i])/(1+exp(p[8]*x12[i]+p[14]*b[j,3,i]))

mcc1[i,j]=p[15]*x11[i]+p[17]*b[j,1,i]
mcc2[i,j]=p[16]*x12[i]+p[17]*b[j,1,i]
pp11[i,j]=q11[i,j]*(1-p11[i,j])
pp22[i,j]=q22[i,j]*(1-p22[i,j])

ppp11[i,j]=(1-p11[i,j])*(1-q11[i,j])
ppp22[i,j]=(1-p22[i,j])*(1-q22[i,j])

if(yo1[i]==1){l1[i,j]=pnorm((p[9]-m1[i,j]),0,1)}
if(yo1[i]==2){l1[i,j]=pnorm((p[10]-m1[i,j]),0,1)-pnorm((p[9]-m1[i,j]),0,1)}
if(yo1[i]==3){l1[i,j]=1-pnorm((p[10]-m1[i,j]),0,1)}
if(yo2[i]==1){ll1[i,j]=pnorm((p[9]-m2[i,j]),0,1)}
if(yo2[i]==2){ll1[i,j]=pnorm((p[10]-m2[i,j]),0,1)-pnorm((p[9]-m2[i,j]),0,1)}
if(yo2[i]==3){ll1[i,j]=1-pnorm((p[10]-m2[i,j]),0,1)}

if(n1[i]==0){l2[i,j]<-p11[i,j]+(ppp11[i,j])*dpois(0,la1[i,j])}
if(n1[i]==1){l2[i,j]<-pp11[i,j]+(ppp11[i,j])*dpois(1,la1[i,j])}
if((n1[i]>1)){l2[i,j]<-(ppp11[i,j])*dpois(n1[i],la1[i,j])}

if(n2[i]==0){ll2[i,j]<-p22[i,j]+(ppp22[i,j])*dpois(n2[i],la2[i,j])}
if(n2[i]==1){ll2[i,j]<-pp22[i,j]+(ppp22[i,j])*dpois(n2[i],la2[i,j])}
if((n2[i]>1)){ll2[i,j]<-(ppp22[i,j])*dpois(n2[i],la2[i,j])}

lc1[i,j]=dnorm(c1[i],mcc1[i,j],p[18])
lc2[i,j]=dnorm(c2[i],mcc2[i,j],p[18])

L1[i,j]=l1[i,j]*l2[i,j]*ll1[i,j]*ll2[i,j]*lc1[i,j]*lc2[i,j]

}
L2[i]=log(mean(L1[i,]))
}
L2=L2[L2!=-Inf]
Lfinal=sum(L2)
return(-Lfinal)

}

p<-c(1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,1,1)

f(p)
est<-nlminb(p,f,lower=c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0),upper=c(Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf))
est

h=fdHess(p,f)
h
ih=solve(h$Hessian)
ih
se<-sqrt(abs(diag(ih)))
se
upp=est$par+(1.96*se)
Low=est$par-(1.96*se)
est
##### Simulation Staudy for (0,1)-inflated Poisson and normal model with Missing data

library(MASS)
library(nlme)

v=1000
vv=1000
b1=matrix(0,v,3)
b=array(0,dim=c(vv,3,v))

x11<-x12<-vector(“numeric”,v)
for(i in 1:v){

x11[i]=rnorm(1,0,1)
x12[i]=rnorm(1,0,1)
}
for(i in 1:v){for(j in 1:vv){
b[j,,i]=rnorm(3,0,1)
}}
c1=c2=mc1=mc2=ppp1=ppp2=pp1=pp2=p1=p2=q1=q2=eyo1=eyo2=y1=y2=n1=nn1=n2=nn2=lla1=lla2=mm1=mm2=yo1=yo2=vector(“numeric”,v)

for(i in 1:v){

mm1[i]=x11[i]+b1[i,1]
mm2[i]=x12[i]+b1[i,1]

y1[i]=rnorm(1,mm1[i],1)
y2[i]=rnorm(1,mm2[i],1)

mc1[i]=x11[i]+b1[i,1]
mc2[i]=x12[i]+b1[i,1]

c1[i]=rnorm(1,mc1[i],1)
c2[i]=rnorm(1,mc2[i],1)

if(y1[i]<=-1){yo1[i]<-1}
if(-1<y1[i]&y1[i]<=1){yo1[i]<-2}
if(y1[i]>1){yo1[i]<-3}

if(y2[i]<=-1){yo2[i]<-1}
if(-1<y2[i]&y2[i]<=1){yo2[i]<-2}
if(y2[i]>1){yo2[i]<-3}

lla1[i]=exp(x11[i]+b1[i,1])
lla2[i]=exp(x12[i]+b1[i,1])

pp1[i]=(1-p1[i])*q1[i]
pp2[i]=(1-p2[i])*q2[i]

ppp1[i]=(1-p1[i])*(1-q1[i])
ppp2[i]=(1-p2[i])*(1-q2[i])

nn1[i]=rpois(1,lla1[i])
nn2[i]=rpois(1,lla2[i])

n1[i]<-sample(c(0,1,nn1[i]),1,prob=c(p1[i],pp1[i],ppp1[i]))
n2[i]<-sample(c(0,1,nn2[i]),1,prob=c(p2[i],pp2[i],ppp2[i]))
}
cbind(c1,c2,n1,n2,yo1,yo2)

P<-pnorm(.5*yo1)
R<-vector(“numeric”,v)
for(i in 1:v){
R[i]<-rbinom(1,1,P[i])
}

cbind(c1,c2,n1,n2,yo1,yo2,R)

################################################################################################
mcc1=mcc2=q11=q22=p11=pp11=ppp11=p22=pp22=ppp22=m1=m2=la1=la2=l1=l2=ll1=ll2=L1=lc1=lc2=matrix(0,v,vv)

L=vector(“numeric”,v)
L=L2=L4=L5=L6=vector(“numeric”,v)

f<-function(p){

for(i in 1:v){for(j in 1:vv){
m1[i,j]=p[1]*x11[i]+p[11]*b[j,1,i]
m2[i,j]=p[2]*x12[i]+p[11]*b[j,1,i]

la1[i,j]=exp(p[3]*x11[i]+p[12]*b[j,1,i])
la2[i,j]=exp(p[4]*x12[i]+p[12]*b[j,1,i])

mcc1[i,j]=p[15]*x11[i]+p[17]*b[j,1,i]
mcc2[i,j]=p[16]*x12[i]+p[17]*b[j,1,i]
pp11[i,j]=q11[i,j]*(1-p11[i,j])
pp22[i,j]=q22[i,j]*(1-p22[i,j])

ppp11[i,j]=(1-p11[i,j])*(1-q11[i,j])
ppp22[i,j]=(1-p22[i,j])*(1-q22[i,j])

if(n1[i]==0){l2[i,j]<-p11[i,j]+(ppp11[i,j])*dpois(0,la1[i,j])}
if(n1[i]==1){l2[i,j]<-pp11[i,j]+(ppp11[i,j])*dpois(1,la1[i,j])}
if((n1[i]>1)){l2[i,j]<-(ppp11[i,j])*dpois(n1[i],la1[i,j])}

if(n2[i]==0){ll2[i,j]<-p22[i,j]+(ppp22[i,j])*dpois(0,la2[i,j])}
if(n2[i]==1){ll2[i,j]<-pp22[i,j]+(ppp22[i,j])*dpois(1,la2[i,j])}
if((n2[i]>1)){ll2[i,j]<-(ppp22[i,j])*dpois(n2[i],la2[i,j])}

lc1[i,j]=dnorm(c1[i],mcc1[i,j],p[18])
lc2[i,j]=dnorm(c2[i],mcc2[i,j],p[18])

if(R[i]==1){L1[i,j]=pnorm((p[19]*yo1[i]),0,1)*l1[i,j]*l2[i,j]*ll1[i,j]*ll2[i,j]*lc1[i,j]*lc2[i,j]}
if(R[i]==0){aa1=(1-pnorm((p[19]*1),0,1))*(pnorm((p[9]-m1[i,j]),0,1))*l2[i,j]*ll1[i,j]*ll2[i,j]*lc1[i,j]*lc2[i,j]
aa2=(1-pnorm((p[19]*2),0,1))*(pnorm((p[10]-m1[i,j]),0,1)-pnorm((p[9]-m1[i,j]),0,1))*l2[i,j]*ll1[i,j]*ll2[i,j]*lc1[i,j]*lc2[i,j]
aa3=(1-pnorm((p[19]*3),0,1))*(1-pnorm((p[10]-m1[i,j]),0,1))*l2[i,j]*ll1[i,j]*ll2[i,j]*lc1[i,j]*lc2[i,j]
L1[i,j]=sum(aa1,aa2,aa3)
}
}
L2[i]=log(mean(L1[i,]))
}
L2=L2[L2!=-Inf]
Lfinal=sum(L2)
return(-Lfinal)

}

p<-c(1,1,1,1,1,1,1,1,-1,1,1,1,1,1,1,1,1,1,.5)

f(p)
est<-nlminb(p,f,lower=c(-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,0,-Inf),upper=c(Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf))
est

h=fdHess(p,f)
h
ih=solve(h$Hessian)
ih
se<-sqrt(abs(diag(ih)))
se
upp=est$par+(1.96*se)
Low=est$par-(1.96*se)
mse=((est$par-p)^2)+(se^2)
est
cbind(p,est$par,se,Low,upp,mse)
pa=est$par
f(pa)
Aici=2*(f(pa)+length(pa))
Aici