#Poisson regression
no=c(0,1,2,3,1,5,10,17,23,31,20,25,37,45)
time=c(1:14)
poi=glm(no~time, family=poisson(link=log))
summary(poi)
par=poi$coeff
names(par)=NULL
par
beta0=par[1]
beta1=par[2]
#par(mfrow=c(2,2))
c1=function(time) exp(beta0+beta1*time)
plot(time,no,xlab="3 month periods: Jan-March 1983 to April-June 1986",ylab="Number of AIDS deaths", cex=2,pch=20,col='black' )
curve(c1,1,14,add=TRUE,col='blue',lwd=2)
abline(lsfit(time,no),col='red',lwd=2)
legend(0.1*max(time),0.9*max(no),text.width=1.2,y.intersp=1.2,bty="n",c("linear regression","Poisson regression"),col = c("red","blue"),lty = c(1,1),lwd=c(2,2),cex=1)  
#title("Poisson regression")

#Logistic regression
y=c(0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1)
dose=c(0:25)/10
log=glm(y~dose, family=binomial(link=logit))
summary(log)
par=log$coeff
names(par)=NULL
par
beta0=par[1]
beta1=par[2]
c1=function(dose) exp(beta0+beta1*dose)/(1+exp(beta0+beta1*dose))
plot(dose,y,ylab="Response (1) or not (0)", pch=20,col='black')
curve(c1,0,2.5,col='blue',add=TRUE,lwd=2)
abline(lsfit(dose,y),col='red',lwd=2)
legend(0.1,0.9,text.width=1.2,y.intersp=1.2,bty="n",c("linear regression","logistic regression"),col = c("red","blue"),lty = c(1,1),lwd=c(2,2),cex=1)  
#title("Logistic regression")


#maximization of using optim, it works with minization
f=function(x) -(20*(3*log(1-x)+log(x)))
optim(c(0.1),f,method = "BFGS")


#maximization of logestic regression
y=c(0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1)
dose=c(0:25)/10

gb<-function(pars){
a<-pars[1]
b<-pars[2]
logl = sum((a+b*dose)*y-log(1+exp(a+b*dose)))
return(-logl)
}

optim(c(0.1,0.1),gb,method="BFGS")


#Fisher Score
n=20
ym=3
pi=0.1  #starting value
result=matrix(0,10,7)

for (i in 1:10) {
dl=n*(1/pi-ym/(1-pi))
edl2=-n/(pi^2*(1-pi))
pi=pi-dl/edl2
#pi=pi+(1-pi-pi*ym)*pi
se=sqrt(-1/edl2)
l=n*(ym*log(1-pi)+log(pi))
step=(1-pi-pi*ym)*pi
result[i,]=c(i,pi,se,l,dl,edl2,step)
}
colnames(result)=c("Iter","pi","se","l","dl","edl2","step")
result


#Newton-Raphson
n=20
ym=3
pi=0.1  #starting value 
result=matrix(0,10,7)

for (i in 1:10) {
dl=n*(1/pi-ym/(1-pi))
dl2=-n*(1/pi^2+ym/(1-pi)^2)
pi=pi-dl/dl2
#pi=pi+(1-pi)*(1-pi-pi*ym)*pi/(1-2*pi+(1+ym)*pi^2)
se=sqrt(-1/dl2)
l=n*(ym*log(1-pi)+log(pi))
step=(1-pi)*(1-pi-pi*ym)*pi/(1-2*pi+(1+ym)*pi^2)
result[i,]=c(i,pi,se,l,dl,dl2,step)
}
colnames(result)=c("Iter","pi","se","l","dl","dl2","step")
result


#likelihood, score and expected information fcn
n=20
ym=3
pi=c(1:100)/100
#out=cbind(pi,logl)
#out
logl=function(pi) n*(ym*log(1-pi)+log(pi))
score=function(pi) n*(1/pi-ym/(1-pi))
Ie=function(pi) n/(pi^2*(1-pi))
Io=function(pi) n*(1/pi^2+ym/(1-pi)^2)

logl1=n*(ym*log(1-pi)+log(pi))
score1=n*(1/pi-ym/(1-pi))
Ie1=n/(pi^2*(1-pi))
Io1=n*(1/pi^2+ym/(1-pi)^2)
c(pi[logl1==max(logl1)],pi[score1==0],max(logl1))
c(Io[pi==0.15],Ie1[pi==0.15],Io1[pi==0.25],Ie1[pi==0.25])

n=20
ym=3
pi=c(1:100)/100
#out=cbind(pi,logl)
#out
logl=function(pi) n*(ym*log(1-pi)+log(pi))
score=function(pi) n*(1/pi-ym/(1-pi))
Ie=function(pi) n/(pi^2*(1-pi))
Io=function(pi) n*(1/pi^2+ym/(1-pi)^2)

logl1=n*(ym*log(1-pi)+log(pi))
score1=n*(1/pi-ym/(1-pi))
Ie1=n/(pi^2*(1-pi))
Io1=n*(1/pi^2+ym/(1-pi)^2)
c(pi[logl1==max(logl1)],pi[score1==0],max(logl1)) #MLE over grid
c(Io1[pi==0.15],Ie1[pi==0.15],Io1[pi==0.25],Ie1[pi==0.25])

par(mfrow=c(2,2))
plot(logl, col='red',xlab="pi",ylab="logl(pi)")
points(pi[score1==0],logl1[pi==pi[score1==0]],pch=2,col="red",cex=0.6)
title("log-likelihood function")
plot(score, col='red',xlab="pi",ylab="score(pi)")
abline(h = 0)
points(pi[score1==0],0,pch=2,col="red",cex=0.6)
title("score function")
plot(Ie, col='red',xlab="pi",ylab="Ie(pi)")
title("Expected information function")
plot(Io, col='red',xlab="pi",ylab="Io(pi)")
title("Observed information function")

#ML using maximizer
x=c(3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,4,4,4,4,4)
logl=function(pi,x) -1*sum(x*log(1 - pi)+log(pi))
pi.start=0.1
out=nlm(logl, pi.start,x=x)
pi.hat = out$estimate
pi.hat


#ML using maximizer
logl = function(pi) -20*(3*log(1-pi)+log(pi))
pi.hat = optimize(logl, c(0, 1), tol = 0.0001)
pi.hat


#EM mixture 
y=c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75)
n=length(y)
y1=y[3:n]
y2=y[1:2]
p=3  #no. of par
iterE=5
iterM=10
dim1=iterE*iterM
dim2=2*p+3
dl=c(rep(0,p))
dl2=matrix(0,p,p)
result=matrix(0,dim1,dim2)
theta=c(mean(y1),mean(y2),var(y1))  #starting values

  for (k in 1:iterE) {   # E-step
  ew1=0.9*exp(-0.5*(y-theta[1])^2/theta[3])
  ew2=0.1*exp(-0.5*(y-theta[2])^2/theta[3])
  w1=ew1/(ew1+ew2)
  w1m=mean(w1)
  w2=1-w1
  sw1=sum(w1)
  sw2=sum(w2) 

    for (i in 1:iterM) {   # M-step
    r1=y-theta[1]
    r2=y-theta[2]
    s1=r1^2-theta[3]
    s2=r2^2-theta[3]
    dl[1]=sum(w1*r1)/theta[3]
    dl[2]=sum(w2*r2)/theta[3]
    dl[3]=0.5*(sum(w1*s1)+sum(w2*s2))/theta[3]^2
    dl2[1,1]=-sw1/theta[3]
    dl2[2,2]=-sw2/theta[3]
    dl2[3,3]=-(sum(w1*s1)+sum(w2*s2))/theta[3]^3-0.5*n/theta[3]^2
    dl2[3,1]=dl2[1,3]=-sum(w1*r1)/theta[3]^2
    dl2[3,2]=dl2[2,3]=-sum(w2*r2)/theta[3]^2
    dl2i=solve(dl2)
    theta=theta-dl2i%*%dl
    se=sqrt(diag(-dl2i))
    l=log(0.9)*sw1+log(0.1)*sw2-n*log(2*pi*theta[3])/2-(sum(w1*r1^2)+sum(w2*r2^2))/(2*theta[3]) 
    row=(k-1)*10+i 
    result[row,]=c(k,i,theta[1],se[1],theta[2],se[2],theta[3],se[3],l)  
    }
  }
colnames(result)=c("iE","iM","mu1","se","mu2","se","sigma2","se","logL") 
result
w=cbind(w1,w2)
w
mean(w1)

library(flexmix) 
mm = flexmix(y ~ 1, k = 2, model = FLXglm(family = "gaussian"))
parameters(mm)
mm.refit=refit(mm)
summary(mm.refit)
 

#EM plots
x=rep(-0.001,n)
x1=seq(-120,100,0.1)
#fx1=exp(-0.5*(x1-theta[1])^2/theta[3])/(sqrt(2*pi*theta[3]))
#fx2=exp(-0.5*(x1-theta[2])^2/theta[3])/(sqrt(2*pi*theta[3]))
fx1=dnorm(x1,theta[1],sqrt(theta[3])) 
fx2=dnorm(x1,theta[2],sqrt(theta[3])) 
fx=0.9*fx1+0.1*fx2
plot(x1, fx1, xlab="x", ylab="f(x)", ylim=c(-0.001,0.025), xlim=c(-120,100), pch=20, col="red",cex=0.7)
points(x1,fx2,pch=20,col='blue',cex=0.5)
points(x1,fx,pch=20,cex=0.5)
points(y,x,pch=20,cex=0.5)
title("Mixture of normals for Darwin data")


#EM censored data
library("msm")
cy=c(6,8,14,16,23,24,28,29,41,49,56,60,75) #first 4 obs are censor time
w=c(0,0,0,0,1,1,1,1,1,1,1,1,1)
n=length(cy)
S=10000  #sim 10000 z for z hat 
m=4   #first 4 censored
p=2   #2 parmeters
cen=cy[1:m]  #censored obs
y=cy[(m+1):n]  #uncensored obs
iterE=10
dim=p+m+1
result=matrix(0,iterE,dim)
z=rep(0,m)
simz=matrix(0,m,S)
theta=c(mean(cy),var(cy))   #mu & var

  for (k in 1:iterE) {   #E-step

    for (j in 1:m) {
    simz[j,]=rtnorm(S,mean=theta[1],sd=sqrt(theta[2]),lower=cen[j],upper=Inf) 
    z[j]=mean(simz[j,])
    
#    cz=(cen[j]-theta[1])/sqrt(theta[2]) 
#    z[j]=theta[1]+dnorm(cz)*sqrt(theta[2])/(1-pnorm(cz))

                   } 
    yr=c(z,y)
    theta[1]=mean(yr)
    theta[2]=(sum(yr^2)-sum(yr)^2/n)/n
    result[k,]=c(k,theta[1],theta[2],z[1],z[2],z[3],z[4])  

      }
colnames(result)=c("iE","mu","sigma2","ez1","ez2","ez3","ez4") 
print(result,digit=5)
print(htheta,digit=5)


#MCEM censored data
library("msm")
cy=c(6,8,14,16,23,24,28,29,41,49,56,60,75) #first 4 obs are censor time
w=c(0,0,0,0,1,1,1,1,1,1,1,1,1)
mean(cy)
n=length(cy)
T=10000  #sim 10000 z for z hat 
m=4   #first 4 censored
p=2   #2 pars
cen=cy[1:m]  #censor obs
y=cy[(m+1):n]  #uncensored obs
iterE=5
iterM=10
dim1=iterE*iterM
dim2=2*p+7
dl=c(rep(0,p))
dl2=matrix(0,p,p)
result=matrix(0,dim1,dim2)
simz=matrix(0,m,T)
z=matrix(0,m,1)
rz=rep(0,m)
r2z=rep(0,m)
theta=c(mean(y),var(y))   #starting values for mu & var

  for (k in 1:iterE) {   #E-step

    for (j in 1:m) {
    simz[j,]=rtnorm(T,mean=theta[1],sd=sqrt(theta[2]),lower=cen[j],upper=Inf) 
    z[j]=mean(simz[j,])  #mean over sim
                   } 

    for (i in 1:iterM) {   #M-step
    rz=z-theta[1] 
    ry=y-theta[1]
    r=c(rz,ry)
    r2z=apply((simz-theta[1])^2,1,mean) #sq 1st; then mean over sim 
    r2=c(r2z,ry^2)
    s2=r2-theta[2]
    l=-n*log(2*pi*theta[2])/2-sum(s2)/(2*theta[2]) 
    dl[1]=sum(r)/theta[2]
    dl[2]=0.5*sum(s2)/theta[2]^2

    dl2[1,1]=-n/theta[2]
    dl2[2,2]=-sum(s2)/theta[2]^3-0.5*n/theta[2]^2
    dl2[2,1]=dl2[1,2]=-sum(r)/theta[2]^2
    dl2i=solve(dl2) 
    theta=theta-dl2i%*%dl
    se=sqrt(diag(-dl2i))
    l=-n*log(2*pi*theta[2])/2-sum(r^2)/(2*theta[2]) 
    row=(k-1)*10+i 
    result[row,]=c(k,i,theta[1],se[1],theta[2],se[2],l,z[1],z[2],z[3],z[4])  
    }
  }

colnames(result)=c("iE","iM","mu","se","sigma2","se","logL","ez1","ez2","ez3","ez4") 
print(result,digit=5)


#ECM Darwin data
y=c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75)
n=length(y)
p=3  #no. of par
iterE = 1000
iterM = 20
dim1=iterE
dim2=p+3
result=matrix(0,dim1,dim2)
tol=1E-7

ECM.method="ECME"	
mu=mean(y)  #starting values
sig2=var(y)
df=5

#### start ECM algorithm #####

for (k in 1:iterE) { 
  
### E-step ###
z2={y-mu}^2/sig2
elam={df+1}/{z2+df}

slamy <- sum(y*elam)
slam <- sum(elam)
	
### CM-step 1 ###
mu=slamy/slam
sig2=sum(elam*{y-mu}^2)/n
	
### E-step ### 
z2={y-mu}^2/sig2 
elam={df+1}/{z2+df}
eloglam=digamma({df+1}/2) - log({df+z2}/2)
	
slam=sum(elam)
sloglam=sum(eloglam)
	
### CM-step2 ###	using Newton Rapson method
  if (ECM.method == "MCECM"){

  for (i in 1:iterM) {		
  dl=n/2*{log(df/2) + 1 - digamma(df/2)} + 0.5*{sloglam - slam}
  dl2=n/2*{1/df - 0.5*trigamma(df/2)}
  df=df - dl/dl2 

  }
  }

  if (ECM.method == "ECME"){

  for (i in 1:iterM) {		
  dl =n/2*{digamma({1+df}/2) - digamma(df/2) - 1/df} - 0.5*sum(log(1+z2/df)) + {1+df}/2/df^2*sum(z2/{1+z2/df})
  dl2=n/4*{2/df^2 + trigamma({1+df}/2) - trigamma(df/2)} + {df+1}/2/df^4*sum(z2^2/{1+z2/df}^2) - 1/df^3*sum(z2/{1+z2/df})
  df=max(df - dl/dl2 , 0.5)
  df=min(df, 20)
  }
  }
	
l=n*{lgamma(df/2+0.5) - lgamma(df/2) - log(pi*df*sig2)/2} - {df/2+0.5}*sum(log(1+z2/df)) 
result[k,]=c(k,i,mu,sig2,df,l)

if (k > 1){if(max(abs(result[k,3:{2+p}]-result[k-1,3:{2+p}])) < tol) break}

}

colnames(result) <- c("iE","iM","mu","sigma2","df","logL"); print(result[1:k,])

par(mfrow=c(2,3))
plot(result[1:k,colnames(result)=="mu"], ylab="", main="mu",pch=20,cex=1.2)
plot(result[1:k,colnames(result)=="sigma2"], ylab="", main="sigma2",pch=20,cex=1.2)
plot(result[1:k,colnames(result)=="df"], ylab="", main="df",pch=20,cex=1.2)
plot(result[1:k,colnames(result)=="logL"], ylab="", main="logL",pch=20,cex=1.2)
dt <- function(x,mu,sigma,df) 1/beta(df/2,0.5)/sqrt(df)/sigma * {1 + {x-mu}^2/df/sigma^2}^{-{df+1}/2}
plot(density(y), lwd = 2,ylim=c(0,0.017),main="obs & fitted pdf",)
curve(dt(x, mu, sqrt(sig2), df1), add=TRUE, col="red", lwd = 2)
points(y,rep(0,n),pch=20,cex=1.2)


#check var
x1=2
x2=3
x3=5
s1=0.1
s2=0.4
s3=0.3
x=matrix(c(1,1,1,x1,x2,x3),3,2)
s=matrix(c(s1,0,0,0,s2,0,0,0,s3),3,3)
s1=solve(s)
xx=t(x)%*%x
xx1=solve(xx)
xsx=t(x)%*%s%*%x
xs1x=t(x)%*%s1%*%x
m1=xx1%*%xsx%*%xx1
m2=solve(xs1x)
m1
m2



tol = 0.0001

n=length(y)
p=3  #no. of par
iterE=5000
iterM=15
dim1=iterE
dim2=p+3
result=matrix(0,dim1,dim2)

df=5   
mu=mean(y)  #starting values
sig2=var(y)*(n-1)/n

for (k in 1:iterE) {   

### E-step ###
z2=(y-mu)^2/sig2 
elam=(df+1)/(z2+df)  #2nd expectation
elnlam=digamma({df+1}/2) - log({df+z2}/2)

slamy <- sum(y*elam)
slam <- sum(elam)
slamr2 <- sum(elam*(y-mu1)^2)

mu=slamy/slam
sig2=slamr2/n
z2=(y-mu)^2/sig2 

  ### M-Step ###    using NR method
  for (i in 1:iterM) {        
  dl=log(df/(df+1))+1-digamma(df/2)+digamma((df+1)/2)+slamln/n-slam/n
  dl2=n/2*{1/df - 1/2*trigamma(df/2)}
  df=df - dl/dl2
  }

l <- n*(log(gamma(df/2+0.5))-log(gamma(df/2))-log(pi*df*sig21)/2)-(df/2+0.5)*sum(log(1+z2/df)) 

result[k,]=c(k,i,mu1,sqrt(sig21),df,l1)
if (k > 1) {if(max(abs(result[k,3:5]-result[k-1,3:5])) < tol) break }
}

colnames(result) <- c("iE","iM","mu","sigma","df","logL"); print(result[-{k+1:dim1},])


#F test
y=c(1,4,8,9,3,8,9)
x1=c(-1,1,-1,1,0,0,0)
x2=c(-1,-1,1,1,0,1,2)
x3=x1*x1
I=matrix(0,7,7)
for (i in 1:7){  I[i,i]=1 }
X=matrix(0,7,4)
X[,1]=c(rep(1,7))
X[,2]=x1
X[,3]=x2
X[,4]=x3
X
XX=t(X)%*%X
XX
XY=t(X)%*%y
XY
XXI=solve(XX)
beta=c(0,0,0,0)
beta=XXI%*%XY
beta

SSB=t(beta)%*%XY
SST=t(y)%*%y
SSR=SST-SSB
c(SSB,SST,SSR)
Z=matrix(0,7,2)
Z[,1]=c(rep(1,7))
Z[,2]=x1+x2
ZZ=t(Z)%*%Z
ZZ
ZY=t(Z)%*%y
ZY
ZZI=solve(ZZ)
beta1=c(0,0)
beta1=ZZI%*%ZY
beta1
SSB1=t(beta1)%*%ZY
SSR1=SST-SSB1
DIFF=SSB-SSB1
c(SSB1,SSR1,DIFF)
F0=(DIFF/2)/(SSR/3)
F0

#Cook's distance
fit= lm(y ~ x1 + x2 + x3)
cooks.distance(fit)
n=dim(X)[1]
p=dim(X)[2]
H=X%*%XXI%*%t(X)
H
R=y-X%*%beta
R
SSR=t(R)%*%R
SSR
sigma2=SSR/3
sigma2
D=c(rep(0,7))
for (i in 1:7){ D[i]=R[i]^2*H[i,i]/((1-H[i,i])^2*4*sigma2) }
D
h=c(rep(0,7))
for (i in 1:7){ h[i]=H[i,i] }
h
hlimit=2*p/n
Dlimit=qf(0.5,p,n-p)
c(hlimit,Dlimit)

#ridge regression
ridge = function (x,y,k) {
n = dim(x)[1]
#intercept = rep(1, n) # an intercept is assumed
#x = cbind(intercept, x)
p = dim(x)[2]
xxk = crossprod(x,x)+diag(rep(k,p))
xy = crossprod(x,y)
betahat = qr.solve(xxk, xy)
resid = y - crossprod(t(x), betahat)
SSE = crossprod(resid, resid)
Syy = var(as.vector(y))*(length(as.vector(y))-1)
R2 = (Syy-SSE)/Syy
list(coefficients=betahat, SSE=as.vector(SSE), R2=R2)
}


p=dim(X)[2]
cor(X)
out=matrix(0,51,5)
for (i in 1:51){
k=i-1
XXk=t(X)%*%X+diag(rep(k,p))
betah=solve(XXk)%*%XY
out[i,]=c(k,betah[1],betah[2],betah[3],betah[4])
}
khat=p*sigma2/(t(beta)%*%beta)
khat
out
par(mfrow=c(2,2))
plot(out[,1],out[,2],xlab="k",ylab="beta*",ylim=c(0,4),xlim=c(0,51),pch=20,col="green")
points(out[,1],out[,3],pch=20,col="blue")
points(out[,1],out[,4],pch=20,col="red")
title("Ridge Traces")
khat=p*sigma2/(t(beta)%*%beta)
khat
ridge1 <- ridge(X,y,k=khat)
ridge1$coefficients


#WLS in GLM
y=c(0,1,2,3,1,5,10,17,23,31,20,25,37,45)
x=c(1:14)
n=length(x)
one=c(rep(1,n))
X=cbind(one,x)
beta=matrix(0.3,2,1)
for (i in 1:10){
eta=X%*%beta
mu=exp(eta)
va=mu
W=matrix(0,n,n)
+ for (j in 1:n)  {W[j,j]=va[j]}
z=eta+(y-mu)/va
XWX=t(X)%*%W%*%X
XWXI=solve(XWX)
XWZ=t(X)%*%W%*%z
beta=XWXI%*%XWZ
VA=diag(XWXI)
se=sqrt(VA)
names(se) = NULL
result=c(i,beta[1],se[1],beta[2],se[2])
print(result)
}



#GEE #binary data #AR(1)
x=c(0:25)/10
y=c(rep(0,5),1,rep(0,4),1,0,1,0,0,rep(1,11))
n=length(x)
p=2
iter=10
one=c(rep(1,n))
X=cbind(one,x)
results=matrix(0,iter,6)
beta=matrix(0.1,2,1)  #starting values

for (k in 1:iter){
eta=X%*%beta
mu=exp(eta)/(1+exp(eta))
se=sqrt(mu*(1-mu))
S=matrix(0,n,n)
D=matrix(0,n,p)

C=diag(one)
rhov=c(rep(0,n-1))

for (i in 1:n) {S[i,i]=se[i]}

for (i in 1:n) 
{ for (j in 1:p) { D[i,j]=se[i]^2*X[i,j] } }  #D matrix

for (i in 1:n-1) 
{ rhov[i]=(y[i]-mu[i])*(y[i+1]-mu[i+1])/(se[i]*se[i+1]) }  #rho hat
rho=sum(rhov)/(n-1)  #rho hat

for (i in 1:n) 
{ for (j in 1:n) { C[i,j]=rho^(abs(i-j)) } }  #working corr. matrix C (AR1)

V=S%*%C%*%S  #working cov matrix
VI=solve(V)
z=y-mu
DVID=t(D)%*%VI%*%D
DVIDI=solve(DVID)
DVIZ=t(D)%*%VI%*%z
beta=beta+DVIDI%*%DVIZ
VA=diag(DVIDI)
se=sqrt(VA)
names(se) = NULL
results[k,]=c(k,beta[1],se[1],beta[2],se[2],rho)
}
results



#GLMM
dat=read.csv("data/glm/random.csv")
dat
y=dat[,1]
n=length(y)
q=10
one=c(rep(1,q))
X=c(rep(1,n))
Z=as.matrix(dat[,-1])
S=diag(one)
SI=solve(S)
I=diag(X)
IZ=I+Z%*%S%*%t(Z)
IZ
IZ=solve(IZ)
XIZY=t(X)%*%IZ%*%y
XIZXI=solve(t(X)%*%IZ%*%X)
beta=XIZXI%*%XIZY
beta
ZR=t(Z)%*%(y-X*beta)
ZZSI=t(Z)%*%Z+solve(S)
ZZSI
U=solve(ZZSI)%*%ZR
U
sum(U)
Res=y-X%*%beta-Z%*%U
sigma2=(t(Res)%*%Res+t(U)%*%SI%*%U)/(n+q)
sigma2


#Poisson regression
data=read.table("CEBMeanVar.txt")
mug=data[,1]
var=data[,2]
plot(log(mean),log(var),pch=16,cex=0.7)


rm(list=ls())
data=read.table("CEB.txt",skip=1)
mu=data[,1]
ni=data[,3]
ln=log(ni)
dur=as.factor(data[,4])
resid=as.factor(data[,5])
edu=as.factor(data[,6])
y=round(mu*ni,0)
d0=glm(y~offset(ln),family=poisson)$dev
d1=glm(y~offset(ln)+dur,family=poisson)$dev
d2=glm(y~offset(ln)+resid,family=poisson)$dev
d3=glm(y~offset(ln)+edu,family=poisson)$dev
d4=glm(y~offset(ln)+dur+resid,family=poisson)$dev
d5=glm(y~offset(ln)+dur+edu,family=poisson)$dev
d6=glm(y~offset(ln)+dur*resid,family=poisson)$dev
d7=glm(y~offset(ln)+dur*edu,family=poisson)$dev
d8=glm(y~offset(ln)+dur+edu+resid,family=poisson)$dev
d9=glm(y~offset(ln)+dur+edu*resid,family=poisson)$dev
d10=glm(y~offset(ln)+dur*resid+edu,family=poisson)$dev
d11=glm(y~offset(ln)+dur*edu+resid,family=poisson)$dev
d12=glm(y~offset(ln)+dur*resid+edu*resid,family=poisson)$dev
d13=glm(y~offset(ln)+dur*edu+resid*edu,family=poisson)$dev
d14=glm(y~offset(ln)+dur*resid+dur*edu,family=poisson)$dev
d15=glm(y~offset(ln)+dur*resid+dur*edu+resid*edu,family=poisson)$dev
c(d0,d1,d2,d3,d4,d5,d6)
c(d7,d8,d9,d10,d11,d12,d13)
summary(glm(y~offset(ln)+dur+resid+edu,family=poisson))
library(MASS)
glm.nb(y ~ offset(ln)+dur+resid+edu)
mean(y)
var(y)

#AIDS deaths
no=c(0,1,2,3,1,5,10,17,23,31,20,25,37,45)
time=c(1:14)

#NB mean regression
logl1 <- function(pars) {
b0<-pars[1]
b1<-pars[2]
r<-pars[3]
p=exp(b0+b1*time)/(r+exp(b0+b1*time))
l1=sum(lgamma(no+r)-lgamma(no+1)-lgamma(r)+r*log(1-p)+no*log(p))
return(-l1)
}
mreg=optim(c(0.37571,0.25365,10),logl1,method = "BFGS",hessian=T)
OI1<-solve(mreg$hessian)

se1<-sqrt(diag(OI1))
par1=mreg$par
rbind(par1,se1)
var1=exp(par1[1]+par1[2]*time)*(par1[3]+exp(par1[1]+par1[2]*time))/par1[3]
AIC1=mreg$value*2+length(par1)*2
AIC1

summary(glm.nb(no~time, link = log))

#NB logistic regression same as mean reg, mean=exp(lnr + xbeta); checking   
logl2 <- function(pars) {
b0<-pars[1]
b1<-pars[2]
r<-pars[3]
p=exp(b0+b1*time)/(1+exp(b0+b1*time))
l2=sum(lgamma(no+r)-lgamma(no+1)-lgamma(r)+r*log(1-p)+no*log(p))
return(-l2)
}
lreg=optim(c(0.3,0.25,10),logl2,method = "BFGS",hessian=T)
OI2<-solve(lreg$hessian)

se2<-sqrt(diag(OI2))
par2=lreg$par
rbind(par2,se2)
AIC2=lreg$value*2+length(par2)*2
AIC2
c(par1[1],log(par2[3])+par2[1]) #2 pars are the same

#GP regression
logl3 <- function(pars) {
b0<-pars[1]
b1<-pars[2]
r<-pars[3]
lam=(1-r)*exp(b0+b1*time)
l3=sum(log(lam)+(no-1)*log(lam+r*no)-(lam+r*no)-lgamma(no+1))
return(-l3)
}
gpreg=optim(c(0.3,0.25,0.01),logl3,method = "BFGS",hessian=T)
OI3<-solve(gpreg$hessian)

se3<-sqrt(diag(OI3))
par3=gpreg$par
rbind(par3,se3)
var3=exp(par3[1]+par3[2]*time)/(1-par3[3])^2
AIC3=gpreg$value*2+length(par3)*2
AIC3

preg=glm(no~time, family=poisson(link=log))
par0=preg$coeff

c0=function(time) exp(par0[1]+par0[2]*time)
c1=function(time) exp(par1[1]+par1[2]*time)
c1v=function(time) exp(par1[1]+par1[2]*time)*(par1[3]+exp(par1[1]+par1[2]*time))/par1[3]
c3=function(time) exp(par3[1]+par3[2]*time)
c3v=function(time) exp(par3[1]+par3[2]*time)/(1-par3[3])^2
plot(time,no,xlab="3 month periods: Jan-March 1983 to April-June 1986",ylab="Number of AIDS deaths", cex=2,pch=20,col='black' )
curve(c0,1,14,add=TRUE,col='blue',lwd=2)
curve(c1,1,14,add=TRUE,col='gray50',lwd=2)
curve(c1v,1,14,add=TRUE,col='gray50',lwd=2,lty=2)
curve(c3,1,14,add=TRUE,col='magenta',lwd=2)
curve(c3v,1,14,add=TRUE,col='magenta',lwd=2,lty=2)
abline(lsfit(time,no),col='red',lwd=2)
legend(0.1*max(time),0.9*max(no),text.width=1.2,y.intersp=1.2,bty="n",c("linear mean","Poisson mean","NB mean","NB var","GP mean","GP var"),col = c("red","blue","gray50","gray50","magenta","magenta"),lty = c(1,1,1,2,1,2),lwd=c(2,2,2,2,2,2),cex=1)  
#title("Poisson regression")

#plot GP pmf
#across r
r=c(-0.1,0.3,0.6,0.9)
N=length(r)
y=0:50
y1=y+1
n1=length(y)
n=n1-1
lam=5
mean=lam/(1-r)
mean
var=lam/(1-r)^3
var

f=matrix(0,N,n1)

for (i in 1:N) {
f[i,]=lam*(lam+r[i]*y)^(y-1)*exp(-1*(lam+r[i]*y))/factorial(y)
  }

rownames(f)=c("r=-0.1","r=0.3","r=0.6","r=0.9")
colnames(f)=0:n
barplot(f,beside=T,col=c("black","red","blue","gray50"),border=NA)
legend(100,0.2,text.width=0.5,y.intersp = 1,bty="n",c("r=-0.1, m=4.5, v=3.8","r=0.3, m=7.1, v=14.6","r=0.6, m=12.5, v=78.1","r=0.9, m=50, v=5000"),col=c("black","red","blue","gray50"),lwd=c(1,1,1,1),lty=c(1,1,1,1),cex=1)
title(main ="pmf of GP across r when lam=5", font.main = 1.5)

#across lam
lam=c(1,3,7,10)
N=length(lam)
y=0:50
y1=y+1
n1=length(y)
n=n1-1
r=0.4
mean=lam/(1-r)
mean
var=lam/(1-r)^3
var

f=matrix(0,N,n1)

for (i in 1:N) {
f[i,]=lam[i]*(lam[i]+r*y)^(y-1)*exp(-1*(lam[i]+r*y))/factorial(y)
  }

rownames(f)=c("lam=1","lam=2","lam=5","lam=7")
colnames(f)=0:n
barplot(f,beside=T,col=c("black","red","blue","gray50"),border=NA)
legend(100,0.35,text.width=0.5,y.intersp = 1,bty="n",c("lam=1, m=1.7, v=4.6","lam=2, m=5, v=13.9","lam=5, m=11.7, v=32.4","lam=7, m=16.7, v=46.3"),col=c("black","red","blue","gray50"),lwd=c(1,1,1,1),lty=c(1,1,1,1),cex=1)
title(main ="pmf of GP across lam when r=0.4", font.main = 2)


#Cannabis data; Poisson mixture model
dat=read.table("CannabisData.txt")
y=dat[,2]
x=dat[,1]
ym=read.table("CannabisData2.txt")
ym=ym[-4,] #delete outlier
m=nrow(ym)
time=1:8
l1=rep(0,m)
l2=rep(0,m)
l=rep(0,m)
ind=rep(0,m)

logl <- function(pars) {
b01=pars[1]; b11=pars[2]; b02=pars[3]; b12=pars[4]; pi=pars[5]
mu1=exp(b01+b11*time)
mu2=exp(b02+b12*time)
for (i in 1:m) {
l1[i]=prod(exp(-mu1)*mu1^ym[i,]/factorial(ym[i,]))
l2[i]=prod(exp(-mu2)*mu2^ym[i,]/factorial(ym[i,]))
l[i]=pi*l1[i]+(1-pi)*l2[i]
}
ll=sum(log(l))
return(-ll)
}

mixreg=optim(c(1.87,-0.25,1.00,-0.05,0.7),logl,method = "BFGS",hessian=T)
OIm<-solve(mixreg$hessian)

sem<-sqrt(diag(OIm))
parm=mixreg$par
rbind(parm,sem)
AICm=mixreg$value*2+length(parm)*2
AICm

b01=parm[1]; b11=parm[2]; b02=parm[3]; b12=parm[4]; pi=parm[5]
mu1=exp(b01+b11*time)
mu2=exp(b02+b12*time)
for (i in 1:m) {
l1[i]=prod(exp(-mu1)*mu1^ym[i,]/factorial(ym[i,]))
l2[i]=prod(exp(-mu2)*mu2^ym[i,]/factorial(ym[i,]))
ind[i]=pi*l1[i]/(pi*l1[i]+(1-pi)*l2[i])
}
ind=round(ind,2)
wt=sum(ind)/m
wt
mi=1:m
gp1=mi[ind==1]
gp2=mi[ind==0]
n1=length(gp1)
n2=length(gp2)
y1=ym[gp1,]
y2=ym[gp2,]
ym1=apply(y1,2,mean)
ym2=apply(y2,2,mean)
yv1=apply(y1,2,var)
yv2=apply(y2,2,var)

c1=function(time) exp(parm[1]+parm[2]*time)
c2=function(time) exp(parm[3]+parm[4]*time)
plot(time,ym1,ylim=c(0,max(y1)),xlab="time",ylab="Number of cannabis offences",cex=2,pch=20,col='black' )
points(time,ym2,type="p",cex=2,pch=20,col='gray50')
curve(c1,1,8,add=TRUE,col='black',lwd=2)
curve(c2,1,8,add=TRUE,col='gray50',lwd=2)
for (i in 1:n1) {points(time,y1[i,],type="p",cex=0.5,pch=16,col='black')}
for (i in 1:n2) {points(time,y2[i,],type="p",cex=0.5,pch=16,col='gray50')}
points(time,yv1,type="p",cex=1.2,pch=17,col='black')
points(time,yv2,type="p",cex=1.2,pch=17,col='gray50')
#par(xpd=NA,oma=c(0,0,0,0))
#par(xpd=NA,mar= c(3, 3, 2, 8))
#legend("right",text.width=1.2,y.intersp=1.2,bty="n",c("high obs mean","high fitted mean","low obs mean","low fitted mean"),col = c("black","black","gray50","gray50"),lty = c(1,NA,1,NA),lwd=c(2,NA,2,NA),pch=c(NA,16,NA,16),cex=1)  
#legend("bottom", xpd =NA, horiz = TRUE,inset=c(0,-0.2),text.width=1.2,y.intersp=1.2,bty="n",c("high obs mean","high fitted mean","low obs mean","low fitted mean"),col = c("black","black","gray50","gray50"),lty = c(1,NA,1,NA),lwd=c(2,NA,2,NA),pch=c(NA,16,NA,16),cex=1)  
legend(0.47*max(time),1.05*max(y1),text.width=0.8,y.intersp=0.8,bty="n",c("high fitted mean","high obs mean","high obs var","low fitted mean","low obs mean","low obs var"),col = c("black","black","black","gray50","gray50","gray50"),lty = c(1,NA,NA,1,NA,NA),lwd=c(2,NA,NA,2,NA,NA),pch=c(NA,20,17,NA,20,17),pt.cex=c(NA,1.2,1.2,NA,1.2,1.2),cex=0.9)  


#library(flexmix) #not work for cluster data
#mm = flexmix(y ~ x, k = 2, model = FLXglm(family = "poisson"))

#plot pmf
ey1=exp(parm[1]+parm[2]*time)
ey2=exp(parm[3]+parm[4]*time)
rbind(ey1,ey2)
T=length(ey1)
n=30
yy=0:n
nyy=length(yy)
f1=matrix(0,T,nyy)
f2=matrix(0,T,nyy)

for (i in 1:T) {
f1[i,]=dpois(yy,ey1[i])
f2[i,]=dpois(yy,ey2[i])
  }
fw=wt*f1+(1-wt)*f2
fr=rbind(f1[2,],f1[4,],f1[6,],f1[8,],f2[2,],f2[4,],f2[6,],f2[8,])
frr=rbind(f1[1,],f2[1,],fw[1,])

#plot pmf for 2 gp at t=2,4,6,8
colnames(fr)=0:n
barplot(fr,beside=T,col=c("black","red","blue","gray50","black","red","blue","gray50"),border=NA)
legend(200,0.25,text.width=0.5,y.intersp = 1,bty="n",c("t=2","t=4","t=6","t=8"),col=c("black","red","blue","gray50"),lwd=c(1,1,1,1),lty=c(1,1,1,1),cex=1)
title(main ="pmf of the 2 groups for time=2,4,6,8", font.main = 1.5)

#plot pmf for 2 gp at t=2,4,6,8
colnames(frr)=0:n
barplot(frr,beside=T,ylim=c(0,max(frr)),col=c("red","blue","black"),border=NA)
legend(80,0.95*max(frr),text.width=0.5,y.intersp = 1,bty="n",c("high level group","low level group","overall"),col=c("red","blue","black"),lwd=c(1,1,1),lty=c(1,1,1),cex=1)
title(main ="pmf of the 2 groups and overall for time=1", font.main = 2)


#GLM #Disturbed Dream data
y=c(7,3,4,7,13,11,15,10,7,11,9,23,10,12,9,28,3,4,5,32)
age=factor(c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5))
rate=factor(c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4))
glm1=glm(y~age+rate,family=poisson)
summary(glm1)
ey=glm1$fitted.values
ey
n=length(y)
mu=mean(y)
mu
D.null=2*sum(y*log(y/mu)-(y-mu))
D.null
D.model=2*sum(y*log(y/ey)-(y-ey))
D.model
AIC=-2*sum(y*log(ey)-ey-log(factorial(y)))+2*(1+3+4)
AIC
r1=factor(c(1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2,1,1,1,2))
glm2=glm(y~age+r1,family=poisson)
summary(glm2)
glm3=glm(y~age*r1,family=poisson)
summary(glm3)
r=residuals.glm(glm3,type="pearson")
par=glm3$coeff
r1=(y[1]-exp(par[1]))/sqrt(exp(par[1]))
r1
glm=glm(y~age*rate,family=poisson)
glm


#logistic regression using different links
#Contraceptive use
y=c(rep(1,507),rep(0,1607-507))
prop.sam=507/1607
prop.sam
glm1=glm(y~1, family=binomial(link=logit))
summary(glm1)
eta1=glm1$coeff
eta1s=eta1/(pi/sqrt(3)) #standardize
pi1=exp(eta1)/(1+exp(eta1))
glm2=glm(y~1, family=binomial(link=probit))
summary(glm2)
eta2=glm2$coeff
pi2=pnorm(eta2,0,1)
glm3=glm(y~1, family=binomial(link=cloglog))
summary(glm3)
eta3=glm3$coeff
eta3s=eta3/(pi/sqrt(6)) #standardize
pi3=1-exp(-exp(eta3))
c(pi1,pi2,pi3)
c(eta1,eta3,eta1s,eta2,eta3s) #eta1,3s comparable to eta2


#logistic regression for 2x2 table 
#Smoker data
smoke=c(32,60)
nonsmoke=c(11,3)
total=smoke+nonsmoke
y=cbind(smoke,nonsmoke)
y
cancer=factor(c(1,2))
glm1=glm(y~cancer, family=binomial)
summary(glm1)
glm1$fitted
smoke/total
beta=glm1$coeff
beta
phat=c(exp(beta[1])/(1+exp(beta[1])),exp(beta[1]+beta[2])/(1+exp(beta[1]+beta[2])))
phat
AIC=-2*sum(smoke*log(phat/(1-phat))+total*log(1-phat)+log(choose(total,smoke)))+2*2
AIC
Dev=2*sum(smoke*log(smoke/(total*phat))+nonsmoke*log(nonsmoke/(total-total*phat)))
Dev
glm0=glm(y~1, family=binomial)
summary(glm0)
beta0=glm0$coeff
phat0=c(exp(beta0)/(1+exp(beta0)),exp(beta0)/(1+exp(beta0)))
phat0
Dev0=2*sum(smoke*log(smoke/(total*phat0))+nonsmoke*log(nonsmoke/(total-total*phat0)))
Dev0
Dev0i=sqrt(2*(smoke*log(smoke/(total*phat0))+nonsmoke*log(nonsmoke/(total-total*phat0))))
Dev0i
#options(contrasts=c("contr.treatment"))
cancer=c(3,60)
nocancer=c(11,32)
total=cancer+nocancer
y=cbind(cancer,nocancer)
smoke=factor(c(1,2))
glm2=glm(y~smoke,family=binomial)
summary(glm2)
glm2$fitted
cancer/total
beta=glm2$coeff
beta



smoke=c(32,60)
nonsmoke=c(11,3)
total=smoke+nonsmoke
x=c(0,1)
y=smoke
n=2
one=c(rep(1,n))
X=cbind(one,x)
beta=matrix(1,2,1)
for (i in 1:10){
eta=X%*%beta
pi=exp(eta)/(1+exp(eta))
mu=total*pi
va=total*pi*(1-pi)
W=matrix(0,n,n)
+ for (j in 1:n)  {W[j,j]=va[j]}
z=eta+(y-mu)/va
XWX=t(X)%*%W%*%X
XWXI=solve(XWX)
XWZ=t(X)%*%W%*%z
beta=XWXI%*%XWZ
VA=diag(XWXI)
se=sqrt(VA)
names(se) = NULL
result=c(i,beta[1],se[1],beta[2],se[2])
print(result)
}
XWXI  #Covariance matrix





#Detergent data
X=c(68,42,37,24,66,33,47,23,63,29,57,19)
T=c(110,72,89,67,116,56,102,70,116,56,106,48)
M=T-X
y=cbind(X,M)
soft=factor(c(1,1,1,1,2,2,2,2,3,3,3,3))
user=factor(c(1,1,2,2,1,1,2,2,1,1,2,2))
temp=factor(c(1,2,1,2,1,2,1,2,1,2,1,2))
dat=data.frame(X,T,soft,user,temp)
glm2=glm(y~soft*user*temp,family=binomial)
summary(glm2)
anova(glm2,test="Chisq")
glm3=glm(y~user+temp,family=binomial)
summary(glm3)
anova(glm3,test="Chisq")
beta3=glm3$coeff
beta3
p1=exp(beta3[1])/(1+exp(beta3[1]))
p2=exp(beta3[1]+beta3[2])/(1+exp(beta3[1]+beta3[2]))
c(p1,p2)


#Matched pair 
y=c(2,125,121,8254)
u=c(0,1,-1,0)
f=c(0,0,0,0)
vet=cbind(y,f)
glm1=glm(vet~u-1,family=binomial)
summary(glm1)


#logistic regression 
#Contraceptive use
#dat=read.csv("data/Contra1.csv")
dat=read.csv("data/Contra.csv")
attach(dat)
yesc=c(rep(0,8))
noc=c(rep(0,8))
for (i in 1:4) {
yesc[2*i-1]=yes[4*(i-1)+1]+yes[4*(i-1)+3]
yesc[2*i]=yes[4*(i-1)+2]+yes[4*(i-1)+4]
noc[2*i-1]=no[4*(i-1)+1]+no[4*(i-1)+3]
noc[2*i]=no[4*(i-1)+2]+no[4*(i-1)+4]
}
agec=c(rep('<25',2),rep('25-29',2),rep('30-39',2),rep('40-49',2))
desc=c(rep(c('Yes','No'),4))
yesc=as.numeric(yesc)
noc=as.numeric(noc)
agec=as.factor(agec)
desc=as.factor(desc)
datc=cbind(agec,desc,yesc,noc)
datc
yc=cbind(yesc,noc)
d0=glm(yc~1,family=binomial(link=logit))$dev
d1=glm(yc~agec,family=binomial(link=logit))$dev
d2=glm(yc~desc,family=binomial(link=logit))$dev
d3=glm(yc~agec+desc,family=binomial(link=logit))$dev
d4=glm(yc~agec*desc,family=binomial(link=logit))$dev
c(d0,d1,d2,d3,d4)
summary(glm(yc~agec+desc,family=binomial(link=logit)))
glm(yc~agec*desc,family=binomial(link=logit))$coeff
ager=c(rep(12.5,2),rep(27.5,2),rep(35,2),rep(45,2))
summary(glm(yc~ager,family=binomial(link=logit)))
summary(glm(yc~ager+desc,family=binomial(link=logit)))
summary(glm(yc~ager*desc,family=binomial(link=logit)))
ager2=ager*ager

glm(yc~ager*desc+ager2*desc,family=binomial(link=logit))$dev
glm(yc~ager*desc+ager2*desc,family=binomial(link=logit))$coeff


dat=read.csv("data/Contra.csv")
attach(dat)
y=cbind(yes,no)
age=as.factor(age)
d0=glm(y~1,family=binomial(link=logit))$dev
d1=glm(y~age,family=binomial(link=logit))$dev
d2=glm(y~edu,family=binomial(link=logit))$dev
d3=glm(y~des,family=binomial(link=logit))$dev
d4=glm(y~age+edu,family=binomial(link=logit))$dev
d5=glm(y~age+des,family=binomial(link=logit))$dev
d6=glm(y~edu+des,family=binomial(link=logit))$dev
d7=glm(y~age*edu,family=binomial(link=logit))$dev
d8=glm(y~age*des,family=binomial(link=logit))$dev
d9=glm(y~edu*des,family=binomial(link=logit))$dev
c(d0,d1,d2,d3,d4,d5,d6,d7,d8,d9)
d10=glm(y~age+edu+des,family=binomial(link=logit))$dev
d11=glm(y~age*edu+des,family=binomial(link=logit))$dev
d12=glm(y~age*des+edu,family=binomial(link=logit))$dev
d13=glm(y~age+des*edu,family=binomial(link=logit))$dev
d14=glm(y~age*edu+age*des,family=binomial(link=logit))$dev
d15=glm(y~age*edu+edu*des,family=binomial(link=logit))$dev
d16=glm(y~age*des+edu*des,family=binomial(link=logit))$dev
d17=glm(y~age*edu+age*des+edu*des,family=binomial(link=logit))$dev
d18=glm(y~age*des*edu,family=binomial(link=logit))$dev
c(d10,d11,d12,d13,d14,d15,d16,d17,d18)

1-pchisq(6.97,1)
1-pchisq(6.77,3)
1-pchisq(17.29,3)
1-pchisq(6.9,1)
1-pchisq(6.83,3)
1-pchisq(1.81,1)
1-pchisq(2.44,3)
1-pchisq(3.36,1)
1-pchisq(11.32,3)
1-pchisq(8.38,3)
1-pchisq(17.35,3)

summary(glm(y~age*des+edu,family=binomial(link=logit)))


#survival analysis
#Leukaemia Example #KM estimator
#load package survival
t=c(6,6,6,6,7,9,10,10,11,13,16,17,19,20,22,23,25,32,32,34,35,
1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)
w=c(0,1,1,1,1,0,0,1,0,1,1,0,0,0,1,1,0,0,0,0,0,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
c=c(rep(0,21),rep(1,21))
l=log(t)
data=cbind(t,w,c,l)
data
surv1=survfit(Surv(t,w)~c)
surv0=survfit(Surv(t,w)~1)
summary(surv1)
summary(surv0)
plot(surv0,xlab="time of remission(weeks)", ylab="survival",main="Leukaemia Data")
plot(surv1,xlab="time of remission(weeks)", ylab="survival",main="Leukaemia Data")
survdiff(Surv(t,w)~c)
#options(contrasts=c("contr.treatment","contr.poly"))
#options(digits=8)


#Leukaemia Example #Proportional hazard model with exp dist.
summary(glm(w~offset(l),family=poisson)) #null model
summary(glm(w~offset(l)+c,family=poisson)) #exp with gp eff.
summary(survreg(Surv(t,w)~c,dist="exp"))  #exp with gp eff.

par=glm(w~offset(l)+gp,family=poisson)$coeff  #plots
loglam=par[1]
lam=exp(par[1])
beta=par[2]
haz=c(lam,lam*exp(beta))
names(haz) = NULL
haz
time=seq(0,35,0.1)
surv1=exp(-lam*time)
surv2=exp(-lam*time)^exp(beta)

par(mfrow=c(2,2))
plot(100,100,xlab="time of remission(weeks)", ylab="hazard",xlim=c(0,35),ylim=c(0,0.5),main="Hazard functions")
abline(haz[1],0,col="blue")
abline(haz[2],0,col="red")
plot(time,surv1,xlab="time of remission(weeks)", ylab="survival",xlim=c(0,35),ylim=c(0,1),col="blue",cex=0.2,main="Survival functions")
points(time,surv2,pch=20,col="red",cex=0.2)


#Leukaemia Example #Proportional hazard model with Weibull dist.
summary(survreg(Surv(t,w)~c,dist="weibull"))
summary(survreg(Surv(l,w)~c,dist="extreme"))

par=survreg(Surv(t,w)~c,dist="weibull")$coeff  #plots
scale=survreg(Surv(t,w)~c,dist="weibull")$scale
alpha=1/scale
beta0=-par[1]*alpha
beta1=-par[2]*alpha
time=seq(0,35,0.1)
haz1=alpha*time^(alpha-1)*exp(beta0)
haz2=alpha*time^(alpha-1)*exp(beta0+beta1)
surv1=exp(-time^alpha*exp(beta0))
surv2=exp(-time^alpha*exp(beta0+beta1))

par(mfrow=c(2,2))
plot(time,haz1,xlab="time of remission(weeks)", ylab="hazard",xlim=c(0,35),ylim=c(0,0.5),col="blue",cex=0.2,main="Hazard functions")
points(time,haz2,pch=20,col='red',cex=0.2)
plot(time,surv1,xlab="time of remission(weeks)", ylab="survival",xlim=c(0,35),ylim=c(0,1),col="blue",cex=0.2,main="Survival functions")
points(time,surv2,pch=20,col='red',cex=0.2)


#estimate alpha & beta alternatively with Weibull dist.

omega=sum(w)
alpha=1
iter=15
result=matrix(0,iter,3)
for (i in 1:iter) {
par=glm(w~offset(alpha*l)+c,family=poisson)$coeff
beta0=par[1]
beta1=par[2]
result[i,]=c(alpha,par[1],par[2])
mu=t^alpha*exp(beta0+c*beta1) 
sum1=sum((mu-w)*(log(t)+beta0/alpha))
alpha=0.5*(alpha+omega/sum1)
}
colnames(result)=c("alpha","beta0","beta1")
result

lam=exp(beta0/alpha)
lam


#Leukaemia Example #Cox proportional hazard model
d2=c(2,2,1,2,2,0,0,4,0,2,2,0,1,0,1,1,1)
d=c(2,2,1,2,2,3,1,4,1,2,2,1,1,1,1,2,2)
r2=c(21,19,17,16,14,12,12,12,8,8,6,4,4,3,3,2,1)
r1=c(21,21,21,21,21,21,17,16,15,13,13,12,11,11,10,7,6)
l=log(r2/r1)
d1=d-d2
dd=cbind(d2,d1)
summary(glm(dd~offset(l),family=binomial))
coxph( Surv(t,w)~c)   #alternative module


#alternate way to read the data
a=c(2,2,21,21,2,2,19,21,1,1,17,21,2,2,16,21,2,2,14,21,0,3,12,21,0,1,12,17,4,4,12,
16,0,1,8,15,2,2,8,13,2,2,6,12,0,1,4,12,1,1,4,11,0,1,3,11,1,1,3,10,1,2,2,7,1,2,1,6)
a1=matrix(a,ncol=4,byrow=T)
d2=a1[,1]
d=a1[,2]
r2=a1[,3]
r1=a1[,4]
l=log(r2/r1)
cbind(d2,d,r2,r1,l)


#survival Infant and Child Mortality data
w=c(168,48,63,89,102,81,40,197,48,62,81,97,103,39,195,55,58,85,87,70,10)
t=c(278.4,538.8,794.4,1550.8,3006,8743.5,14270,403.2,786,1165.3,2294.8,
4500.5,13201.5,19525,495.3,956.7,1381.4,2604.5,4618.5,9814.5,5802.5)
age=as.factor(rep(c('A:0-1m','B:1-3m','C:3-6m','D:6-12m','E:1-2y','F:2-5y','G:5-10y'),3))
cohort=as.factor(rep(c('1941-59','1960-67','1968-76'),c(7,7,7)))
l=log(t)
d0=glm(w~offset(l),family=poisson)$dev
d1=glm(w~offset(l)+age,family=poisson)$dev
d2=glm(w~offset(l)+cohort,family=poisson)$dev
d3=glm(w~offset(l)+age+cohort,family=poisson)$dev
d4=glm(w~offset(l)+age*cohort,family=poisson)$dev
round(c(d0,d1,d2,d3,d4),3)
summary(glm(w~offset(l)+age+cohort,family=poisson))
par=glm(w~offset(l)+age+cohort,family=poisson)$coeff
width=c(1/12,2/12,3/12,1/2,1,3,5)
npi=101
widthp=width/npi  #for plot
#width=c(30,60,91,182,365,1095,1825)
month=c(1,3,6,12,24,60)
loghaz=c(0,0,0,0,0,0,0)
loghaz[1]=par[1]
loghaz[2]=par[1]+par[2]
loghaz[3]=par[1]+par[3]
loghaz[4]=par[1]+par[4]
loghaz[5]=par[1]+par[5]
loghaz[6]=par[1]+par[6]
loghaz[7]=par[1]+par[7]
haz=exp(loghaz)
hazless60=haz
haz60_67=haz*exp(par[8])
haz68_76=haz*exp(par[9])
width_haz=width*haz
cumhaz=c(width_haz[1],0,0,0,0,0,0)
cumwid=c(width[1],0,0,0,0,0,0)
for (i in 2:7) {cumhaz[i]=cumhaz[i-1]+width_haz[i]}
for (i in 2:7) {cumwid[i]=cumwid[i-1]+width[i]}
surless60=exp(-cumhaz)
sur60_67=exp(-cumhaz)^(exp(par[8]))
sur68_76=exp(-cumhaz)^(exp(par[9]))
surless60=c(1,surless60)
sur60_67=c(1,sur60_67)
sur68_76=c(1,sur68_76)
cumwid0=c(0,cumwid)
cumhaz0=c(0,cumhaz)

result1=cbind(width,loghaz,haz,hazless60,haz60_67,haz68_76)
rownames(result1)=c('0-1m','1-3m','3-6m','6-12m','1-2y','2-5y','5-10y') 
round(result1,digits=4) 
result2=cbind(width_haz,cumhaz,surless60,sur60_67,sur68_76)
rownames(result2)=c('0-1m','1-3m','3-6m','6-12m','1-2y','2-5y','5-10y') 
round(result2,digits=4)
 
hazfun1=stepfun(month,hazless60,f=0)
hazfun2=stepfun(month,haz60_67,f=0)
hazfun3=stepfun(month,haz68_76,f=0)
plot(hazfun1,month,xlab="age(months)", ylab="hazard function",xlim=c(0,120),ylim=c(0,1),main="Piece-wise hazard function")
plot(hazfun2,month,add=TRUE,xlim=c(0,120))
plot(hazfun3,month,add=TRUE,xlim=c(0,120))

plot(cumwid,surless60,xlab="age(year)", ylab="survival function",xlim=c(0,10),ylim=c(0.8,1),type="o",lty=1,lwd=2,col="blue",main="Piece-wise survival function")
points(cumwid,sur60_67,type="o",lty=1,lwd=2,col="darkgreen",xlim=c(0,120))
points(cumwid,sur68_76,type="o",lty=1,lwd=2,col="red",xlim=c(0,120))

#add more points but the curve still not smooth
for (i in 1:7) {
    if (i==1) {
      cumwidp=seq(cumwid0[i],cumwid0[i+1],length=npi)
      cumwidp=cumwidp[1:(length(cumwidp)-1)] 
      cumhazp=seq(cumhaz0[i],cumhaz0[i+1],length=npi)
      cumhazp=cumhazp[1:(length(cumhazp)-1)] }

    if (i>1) { 
      cumwidp=c(cumwidp,seq(cumwid0[i],cumwid0[i+1],length=npi))
      cumwidp=cumwidp[1:(length(cumwidp)-1)] 
      cumhazp=c(cumhazp,seq(cumhaz0[i],cumhaz0[i+1],length=npi))
      cumhazp=cumhazp[1:(length(cumhazp)-1)]       }
}

cumwidp
cumhazp
cumwid0
cumhaz0
surless60p=exp(-cumhazp)
sur60_67p=exp(-cumhazp)^(exp(par[8]))
sur68_76p=exp(-cumhazp)^(exp(par[9]))


plot(cumwidp,surless60p,xlab="age(year)", ylab="survival function",xlim=c(0,10),ylim=c(0.8,1),type="l",lty=1,lwd=2,col="blue",main="Piece-wise survival function")
points(cumwidp,sur60_67p,type="l",lty=1,lwd=2,col="darkgreen")
points(cumwidp,sur68_76p,type="l",lty=1,lwd=2,col="red")
points(cumwid0,surless60,type="p",lty=1,lwd=2,col="blue")
points(cumwid0,sur60_67,type="p",lty=1,lwd=2,col="darkgreen")
points(cumwid0,sur68_76,type="p",lty=1,lwd=2,col="red")

agec=month/12
agec=c(agec,10)
agec=c(agec,agec,agec)
summary(glm(y~offset(l)+agec+cohort,family=poisson))
par=glm(y~offset(l)+agec+cohort,family=poisson)$coeff
par


#log-linear model 2x2 table
smoke=c(32,60)
nonsmoke=c(11,3)
total=smoke+nonsmoke
y=cbind(smoke,nonsmoke)
yr=c(nonsmoke,smoke)
smokef=factor(c(0,0,1,1))
cancerf=factor(c(0,1,0,1))
cancer=factor(c(1,2))
glm(y~1, family=binomial)$dev
glm(yr~smokef+cancerf,family=poisson)$dev


#Incomplete Contingency Table Example 1
y1=c(9,6,5,0,8,5,7,16)
a=factor(c(1,1,2,2,1,1,2,2))
b=factor(c(1,2,1,2,1,2,1,2))
c=factor(c(1,1,1,1,2,2,2,2))
glm1=glm(y1~a*b+a*c+b*c,family=poisson)
summary(glm1)
glm1$fitted
#Example with marginal 0 total
y2=c(9,6,5,0,8,5,7,0)
glm2=glm(y2~a*b+a*c+b*c,family=poisson)
summary(glm2)
glm2$fitted


#Teenage Health Concern Data
y=c(4,2,9,7,42,7,19,10,57,20,71,31,4,8)
H=factor(c(1,1,1,1,2,2,2,2,3,3,3,3,4,4))
A=factor(c(1,2,1,2,1,2,1,2,1,2,1,2,1,2))
S=factor(c(1,1,2,2,1,1,2,2,1,1,2,2,2,2))
d0=glm(y~H*A*S,family=poisson)$dev
d1=glm(y~H*A+H*S+A*S,family=poisson)$dev
d2=glm(y~H*A+H*S,family=poisson)$dev
d3=glm(y~H*A+A*S,family=poisson)$dev
d4=glm(y~H*S+A*S,family=poisson)$dev
d5=glm(y~A*S+H,family=poisson)$dev
d6=glm(y~H*S+A,family=poisson)$dev
d7=glm(y~H*A+S,family=poisson)$dev
d8=glm(y~H+S+A,family=poisson)$dev
c(d0,d1,d2,d3,d4,d5,d6,d7,d8)
summary(glm(y~H*A*S,family=poisson))
summary(glm(y~H*A+H*S+A*S,family=poisson))
summary(glm(y~H*A+H*S,family=poisson))
summary(glm(y~H*A+A*S,family=poisson))
summary(glm(y~H*S+A*S,family=poisson))
summary(glm(y~A*S+H,family=poisson))
summary(glm(y~H*S+A,family=poisson))
summary(glm(y~H*A+S,family=poisson))
summary(glm(y~H+S+A,family=poisson))


#British Election Data
y=c(157,4,17,9,16,159,13,9,11,9,51,1,18,12,11,15)
a=factor(c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4))
b=factor(c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4))
summary(glm(y~a+b,family=poisson))
yh=glm(y~a+b,family=poisson)$fitted
chi2=sum((y-yh)^2/yh)
chi2

yr=c(4,17,9,16,13,9,11,9,1,18,12,11)
ar=factor(c(1,1,1,2,2,2,3,3,3,4,4,4))
a1r=factor(c(1,1,1,1,1,1,2,2,2,1,1,1))
br=factor(c(2,3,4,1,3,4,1,2,4,1,2,3))
summary(glm(yr~ar+br,family=poisson))
yh=glm(yr~ar+br,family=poisson)$fitted
chi2=sum((yr-yh)^2/yh)
chi2


#Bayesian
model=c(0:10)/10
prior=c(0.05,0.10,0.15,0.19,0.18,0.15,0.09,0.05,0.02,0.01,0.01)
n=10
y=6
like=pbinom(y,n,model)-pbinom(y-1,n,model)
priorxlike=prior*like
posterior=priorxlike/sum(priorxlike)
result=cbind(model,prior,like,priorxlike,posterior)
print(formatC(result,dig=5,format="f")) 
smean=y/n
priormean=sum(model*prior)
postmean=sum(model*posterior)
c(smean,priormean,postmean)
par(mfrow=c(2,2))
plot(model,prior,xlab="model (value of p)",ylab="probability*",ylim=c(0,0.4),xlim=c(0,1),pch=20,col="red")
points(model,posterior,pch=20,col="blue")
title("Distrobution of p")


#Bayesian WinBUGS Darwin data #2 var, 1 mean
#Model: y_i ~ p*N(mu,sigma1^2)+(1-p)*N(mu,sigma2^2)
model
{ for (i in 1:N) {
x[i] ~ dnorm(mu,invar[i])
invar[i] <- ind[i]*insigma12+(1-ind[i])*insigma22
ind[i] ~ dbern(p)
}
p <- 0.9
#p ~ dunif(0,1)
mu ~ dnorm(0,0.0001)
#k <- 10
insigma12 ~ dgamma(0.001,0.001)
insigma22 ~ dgamma(0.001,0.001)
sigma12 <- 1/sqrt(insigma12)
sigma22 <- 1/sqrt(insigma22)
}


#Bayesian WinBUGS Darwin data #2 mean, 1 var 
#Model: y_i ~ p*N(mu1,sigma^2)+(1-p)*N(mu2,sigma^2)
model
{ for (i in 1:N) {
x[i] ~ dnorm(mean[i],insigma2)
mean[i]<- w1[i]*mu1+(1-w1[i])*mu2
w1[i] ~ dbern(p)
}
p <- 0.9
#p ~ dunif(0,1)
mu1 ~ dnorm(0,0.0001)
mu2 ~ dnorm(0,0.0001)
insigma2 ~ dgamma(0.001,0.001)
sigma2 <- 1/insigma2
}

#Bayesian WinBUGS Darwin data #2 mean, 1 var, t as SMN
#Model: y_i ~ p*t(mu1,sigma2,df)+(1-p)*t(mu2,sigma2,df)
model
{ for (i in 1:N) {
x[i] ~ dnorm(mean[i],tau2[i])
tau2[i]<- insigma2*lam[i]
mean[i]<- w1[i]*mu1+(1-w1[i])*mu2
w1[i] ~ dbern(p)
lam[i] ~ dgamma(df2,df2)
}
p <- 0.9
#p ~ dunif(0,1)
df <- 5
#df ~ dgamma(0.001,0.001)
df2 <- df/2
mu1 ~ dnorm(0,0.0001)
mu2 ~ dnorm(0,0.0001)
insigma2 ~ dgamma(0.001,0.001)
sigma2 <- 1/insigma2
}


#Bayesian WinBUGS Darwin data #1 mean, 1 var, normal
#Model: y_i ~ N(mu,sigma2)
model
{ for (i in 1:N) {
x[i] ~ dnorm(mu,insigma2)
}
mu ~ dnorm(0,0.0001)
insigma2 ~ dgamma(0.001,0.001)
sigma2 <- 1/insigma2
}


#Bayesian WinBUGS Darwin data #1 mean, 1 var, t
#Model: y_i ~ t(mu,sigma2,nu) as SMN
model
{ for (i in 1:N) {
x[i] ~ dnorm(mu,tau2[i])
tau2[i]<- insigma2*lam[i]
lam[i] ~ dgamma(df2,df2)
}
df <- 5
df2 <- df/2
mu ~ dnorm(0,0.0001)
insigma2 ~ dgamma(0.001,0.001)
sigma2 <- 1/insigma2
}

list(x=c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75),N=15) #data

list(mu=25,insigma12=0.001,insigma12=0.001) #initials 2 var, 1 mean 
list(mu1=25,mu2=-25,insigma2=0.001) #initials 2 mean, 1 var 
list(mu=25,insigma2=0.001)  #initials 1 var, 1 mean



#Bayesian WinBUGS AIDS death data
#Model: y_i ~ Poi(b0+b1*i)
model
{ for (i in 1:N) {
y[i] ~ dpois(mean[i])
mean[i]<- exp(b0+b1*i)
}
b0 ~ dnorm(0,0.0001)
b1 ~ dnorm(0,0.0001)
}

list(y=c(0,1,2,3,1,5,10,17,23,31,20,25,37,45),N=14) #data
list(b0=1,b1=1) #initials
list(b0=1)


#Bayesian WinBUGS random effect data
#Model: y_ij ~ N(b0+ui,sigma2), ui ~ N(0,sigma2)
model
{ for (ij in 1:N) {
y[ij] ~ dnorm(mu[ij],insigma2)
mu[ij]<- b0+u[g[ij]]
}

for (i in 1:G) {
u[i] ~ dnorm(0,insigma2)
}

b0 ~ dnorm(0,0.0001)
insigma2 ~ dgamma(0.001,0.001)
sigma2 <- 1/insigma2
}

list(b0=1,insigma2=0.1,u=c(0,0,0,0,0,0,0,0,0,0)) #initials
list(N=60,G=10)  #data
y[] g[]
5.058	1
8.611	1
3.739	1
5.201	1
6.944	1
6.387	1
3.901	2
6.735	2
3.03	2
6.539	2
5.508	2
3.345	2
7.623	3
6.876	3
3.451	3
5.468	3
4.942	3
2.405	4
5.007	4
3.153	4
2.856	4
6.309	4
7.839	5
5.476	5
2.892	5
6.143	5
5.605	5
3.957	6
1.025	6
5.463	6
3.399	6
6.285	6
5.28	6
4.602	6
9.057	7
6.325	7
6.394	7
3.464	7
7.541	7
6.603	7
5.909	8
7.404	8
7.937	8
4.275	8
6.391	8
4.023	8
8.56	8
17.779	9
13.841	9
14.033	9
14.526	9
15.599	9
15.848	9
11.967	9
17.839	10
15.476	10
12.894	10
16.143	10
15.604	10
14.678	10
END

#Bayesian plot mixtures of normal and mixtures of t 
y=c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75)
n=length(y)
df=5
x=rep(-0.001,n)
theta=c(0,0,0)
thetat=c(0,0,0)
x1=seq(-120,100,0.1)
theta[1]=31.23 #normal
theta[2]=-52.43
theta[3]=572.7
fx1=dnorm(x1,theta[1],sqrt(theta[3])) 
fx2=dnorm(x1,theta[2],sqrt(theta[3])) 
fx=0.9*fx1+0.1*fx2 #pdf for normal

thetat[1]=29.17 #t 
thetat[2]=-49.37
thetat[3]=521.6
zt1=(x1-thetat[1])/sqrt(thetat[3]) 
zt2=(x1-thetat[2])/sqrt(thetat[3])
ftx1=dt(zt1,df)/sqrt(thetat[3]) 
ftx2=dt(zt2,df)/sqrt(thetat[3]) 
ftx=0.9*ftx1+0.1*ftx2 #pdf for t

plot(x1, fx, xlab="x", ylab="f(x)", ylim=c(-0.001,0.018), 
xlim=c(-120,100), pch=20, col="red",cex=0.5) #plot
points(x1,ftx,pch=20,col='blue',cex=0.5)
points(y,x,pch=20,cex=0.8)
title("Mixtures of normal and t for Darwin data")

fy1=dnorm(y,theta[1],sqrt(theta[3])) #AIC for normal
fy2=dnorm(y,theta[2],sqrt(theta[3]))
aic=-2*sum(log(0.9*fy1+0.1*fy2))+6
zty1=(y-thetat[1])/sqrt(thetat[3]) #AIC for t
zty2=(y-thetat[2])/sqrt(thetat[3])
fty1=dt(zty1,df)/sqrt(thetat[3]) 
fty2=dt(zty2,df)/sqrt(thetat[3])
aict=-2*sum(log(0.9*fty1+0.1*fty2))+6
c(aic,aict)


#Bayesian without the use of WinBUGS 
#Darwin data, Y ~ N(mu,sigma2) 
y=c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75)
ym=mean(y)
n=length(y)
samsize=2000  #posterior sample size
burn=2000
skip=10
iters=samsize*skip+burn  #total iterations

tau2=1000  #hyperparameters
mu0=0
alpha=0.001
beta=0.001

sigma2=1  #starting value
mu=0     #starting value

mu.sam=matrix(0,samsize)
sigma2.sam=matrix(0,samsize)
dic.sam=matrix(0,samsize)
t=matrix(0,samsize)
result=matrix(0,3,5)
count=0

for (i in 1:iters) {
vari=1/(sigma2/n)+1/tau2
mean=(ym/(sigma2/n)+mu0/tau2)/vari
sd=1/sqrt(vari)
mu=rnorm(1,mean,sd)
shape=alpha+n/2
scalee=beta+0.5*sum((y-mu)^2)
scalei=1/scalee
sigma2i=rgamma(1,shape,rate=scalee,scale=scalei)
sigma2=1/sigma2i
if (i>burn) {
if ((i-burn)%%skip==0) {
count=count+1
t[count]=i
mu.sam[count]=mu
sigma2.sam[count]=sigma2
dic.sam[count]=-2*sum(log(dnorm(y,mu,sqrt(sigma2))))
}
}
}

par(mfrow=c(2,1))  #history plot
plot(t,mu.sam,type="l",col="red",main="Posterior sample of mu")
plot(t,sigma2.sam,type="l",col="red",main="Posterior sample of sigma2")
emu=mean(mu.sam)
esigma2=mean(sigma2.sam)
Dbar=mean(dic.sam)
Dhat=-2*sum(log(dnorm(y,emu,sqrt(esigma2))))
pD=Dbar-Dhat
rownames(result)=c('mu','sigma2','dic')  
colnames(result)=c('est','sd','median','2.5%','97.5%') 
result[1,]=c(mean(mu.sam),sd(mu.sam),median(mu.sam),quantile(mu.sam,0.025),quantile(mu.sam,0.975))
result[2,]=c(mean(sigma2.sam),sd(sigma2.sam),median(sigma2.sam),quantile(sigma2.sam,0.025),quantile(sigma2.sam,0.975))
result[3,]=c(Dbar,Dhat,pD,Dbar+pD,NA)
round(result,digits=4)



x=seq(-4,4,0.01)
c1=pi/sqrt(3)
c3=pi/sqrt(6)
l1=exp(x*c1)/(1+exp(x*c1))
l2=pnorm(x)
l3=1-exp(-exp(x*c3))

plot(x,l1,xlab="link", ylab="probability",xlim=c(-4,4),ylim=c(0,1),col="blue",lwd=2,cex=0.3,cex.lab=1.5,cex.axis=1.5,cex.main=1.5,main="Link functions")
points(x,l2,pch=20,col="red",lwd=2,cex=0.3)
points(x,l3,pch=20,col="darkgreen",lwd=2,cex=0.3)
legend(0,0.3,text.width=2,bty="n",c("logit","probit","com-loglog"),col = c("blue","red","darkgreen"),lty = c(1,1,1),lwd=c(2,2,2),cex=1.5)  #end skip