forked from moreno-betancur/survMCOD
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlikelihoods.r
151 lines (119 loc) · 5.11 KB
/
likelihoods.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
########################## Likelihoods functions ##########################
# All of these functions assume that datamat has in its first columns Z1,Z2,..,ZK X1,X2,...,XJ,
# in that order and that data is sorted in ascending failure time
#' Loglikelihood L as a function of (xi,rho,phi) to get starting values
#'
#' Internal function used in estimation procedure
#' @export
#' @keywords internal
loglik=function(par,nxi,nX,nZ,pset,datamat,psetTT,txi,tt,datInd) {###par contains (xi,rho,phi)
xi=par[1:nxi]
if(nX==0){rho=as.matrix(0)}else{rho=as.matrix(c(rep(0,times=nZ),par[(nxi+1):(nxi+nX)]))}
if(nZ==0){phi=as.matrix(0)}else{phi=as.matrix(c(par[(nX+nxi+1):(nZ+nX+nxi)],rep(0,nX)))}
ll<-0
for(pp in pset)
{
q<-which(psetTT==pp)
datXi<-NULL
for(ff in 1:nxi)
datXi<-cbind(datXi,pp*exp(datamat[,1:(nZ+nX)]%*%rho)+(1-pp)*exp(-xi[ff]+datamat[,1:(nZ+nX)]%*%phi))
# get denominator: sum of contributions of individuals at risk for each uniqute type-q event time
ttq<-unique(datamat[datamat[,"status"]==q,"AgeExit"])
cum<-vector()
for(ff in 1:nxi)
{tind<-(tt%in%ttq)&(txi[ff]<=tt)&(tt<ifelse(is.na(txi[ff+1]),Inf,txi[ff+1]))
datq<-datInd[,tind]
cum<-c(cum,as.matrix(t(datXi[,ff])%*%datq))}
# get power of denominator: number of ties for each uniqute type-q event time
ttq<-datamat[datamat[,"status"]==q,"AgeExit"]
tt2<-table(ttq)
tt2<-data.frame(ID=1:nrow(tt2),ttq=as.numeric(dimnames(tt2)$ttq),tt2,row.names=NULL)[,c(1,2,4)]
ll<-ll-sum(log((cum)^tt2$Freq))
for(ff in 1:nxi)
{
dp<-(datamat[,"status"]==q&txi[ff]<=datamat[,"AgeExit"]&datamat[,"AgeExit"]<ifelse(is.na(txi[ff+1]),Inf,txi[ff+1]))
ll<-ll+sum(log(datXi[dp,ff]))
}
}
return(-ll) #ll is negative because optim performs minimization
}
#' Loglikelihood L as a function of (rho,phi) only, with xi known from previous iteration
#' Exactly the same as loglik except now xi is not an argument to the function so the first three
#' lines change
#'
#' Internal function used in estimation procedure
#' @export
#' @keywords internal
loglikopt=function(par,nxi,nX,nZ,pset,datamat,psetTT,txi,tt,datInd,opt,j){ ###par contains (rho,phi)
xi=opt[j,1:nxi]
if(nX==0){rho=as.matrix(0)}else{rho=as.matrix(c(rep(0,nZ),par[1:nX]))}
if(nZ==0){phi=as.matrix(0)}else{phi=as.matrix(c(par[(nX+1):(nZ+nX)],rep(0,nX)))}
ll<-0
for(pp in pset)
{
q<-which(psetTT==pp)
datXi<-NULL
for(ff in 1:nxi)
datXi<-cbind(datXi,pp*exp(datamat[,1:(nZ+nX)]%*%rho)+(1-pp)*exp(-xi[ff]+datamat[,1:(nZ+nX)]%*%phi))
# get denominator: sum of contributions of individuals at risk for each uniqute type-q event time
ttq<-unique(datamat[datamat[,"status"]==q,"AgeExit"])
cum<-vector()
for(ff in 1:nxi)
{tind<-(tt%in%ttq)&(txi[ff]<=tt)&(tt<ifelse(is.na(txi[ff+1]),Inf,txi[ff+1]))
datq<-datInd[,tind]
cum<-c(cum,as.matrix(t(datXi[,ff])%*%datq))}
# get power of denominator: number of ties for each uniqute type-q event time
ttq<-datamat[datamat[,"status"]==q,"AgeExit"]
tt2<-table(ttq)
tt2<-data.frame(ID=1:nrow(tt2),ttq=as.numeric(dimnames(tt2)$ttq),tt2,row.names=NULL)[,c(1,2,4)]
ll<-ll-sum(log((cum)^tt2$Freq))
for(ff in 1:nxi)
{
dp<-(datamat[,"status"]==q&txi[ff]<=datamat[,"AgeExit"]&datamat[,"AgeExit"]<ifelse(is.na(txi[ff+1]),Inf,txi[ff+1]))
ll<-ll+sum(log(datXi[dp,ff]))
}
}
return(-ll) #ll is negative because optim performs minimization
}
#' Loglikelihood L* as a function of xi only, with rho and phi known from previous iteration
#' Only difference with L is denominator
#'
#' Internal function used in estimation procedure
#' @export
#' @keywords internal
#'
logliksteropt=function(par,nxi,nX,nZ,pset,datamat,psetTT,txi,tt,datInd,opt,i){ ###par contains (xi)
xi=par
if(nX==0){rho=as.matrix(0)}else{rho=as.matrix(c(rep(0,nZ),opt[i,(nxi+1):(nxi+nX)]))}
if(nZ==0){phi=as.matrix(0)}else{phi=as.matrix(c(opt[i,(nX+nxi+1):(nZ+nX+nxi)],rep(0,nX)))}
ll<-0
Q<-length(pset)
datXiT<-NULL
for(ff in 1:nxi)
datXiT<-cbind(datXiT,sum(pset)*exp(datamat[,1:(nZ+nX)]%*%rho)+(Q-sum(pset))*exp(-xi[ff]+datamat[,1:(nZ+nX)]%*%phi)) #aT
for(pp in pset)
{
q<-which(psetTT==pp)
datXi<-NULL
for(ff in 1:nxi)
datXi<-cbind(datXi,pp*exp(datamat[,1:(nZ+nX)]%*%rho)+(1-pp)*exp(-xi[ff]+datamat[,1:(nZ+nX)]%*%phi))
# get denominator: sum of contributions of individuals at risk for each uniqute type-q event time
ttq<-unique(datamat[datamat[,"status"]==q,"AgeExit"])
cum<-vector()
for(ff in 1:nxi)
{tind<-(tt%in%ttq)&(txi[ff]<=tt)&(tt<ifelse(is.na(txi[ff+1]),Inf,txi[ff+1]))
datq<-datInd[,tind]
cum<-c(cum,as.matrix(t(datXiT[,ff])%*%datq))}
# get power of denominator: number of ties for each uniqute type-q event time
ttq<-datamat[datamat[,"status"]==q,"AgeExit"]
tt2<-table(ttq)
tt2<-data.frame(ID=1:nrow(tt2),ttq=as.numeric(dimnames(tt2)$ttq),tt2,row.names=NULL)[,c(1,2,4)]
ll<-ll-sum(log((cum)^tt2$Freq))
for(ff in 1:nxi)
{
dp<-(datamat[,"status"]==q&txi[ff]<=datamat[,"AgeExit"]&datamat[,"AgeExit"]<ifelse(is.na(txi[ff+1]),Inf,txi[ff+1]))
ll<-ll+sum(log(datXi[dp,ff]))
}
}
return(-ll)
}