-
Notifications
You must be signed in to change notification settings - Fork 0
/
multilikelihood.R
93 lines (66 loc) · 5.19 KB
/
multilikelihood.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
multiloglikelihood = function(c,Y1,Y2,D1,D2,That,K, beta, ro, r, si,sig2.dat,gmmx1, gmmx2, regy1, regy2 ) {
GMMlikelihood <- c(0)
AFTlikelihood <- c(0)
loglikelihood <- c(0)
gmmx1 <- gmmx1
gmmx2 <- gmmx2
regy1 <- regy1
regy2 <- regy2
disclass <- table(factor(c, levels = 1:K))
activeclass <- which(disclass!=0)
Ytemp1 <- Y1
Ytemp2 <- Y2
N = nrow(Y1)
Ytemp1.scaled <- matrix(NA, nrow = N, ncol = D1)
Ytemp2.scaled <- matrix(NA, nrow = N, ncol = D2)
for ( i in 1:length(activeclass)) {
clust <- which(c == activeclass[i])
if (length(clust)==1){
Ytemp1.scaled[clust,1:D1] <- matrix(0, nrow =1, ncol =D1)
Ytemp2.scaled[clust,1:D2] <- matrix(0, nrow =1, ncol =D2)
} else {
Ytemp1.scaled[clust,1:D1] <- scale(Ytemp1[clust,1:D1], center = TRUE, scale = TRUE)
Ytemp2.scaled[clust,1:D2] <- scale(Ytemp2[clust,1:D2], center = TRUE, scale = TRUE)
}
}
loglikelihood <-c(0)
for (j in 1:length(activeclass)) {
loglikelihood[j] <- 0
# clust <- which(c==activeclass[j])
# if (length(clust)==1){
# loglikelihood[j] <- 0
# } else {
# for ( l in 1:length(clust)) {
# if (Time[clust[l],2]==1){
# loglikelihood[j] <- loglikelihood[j] + dMVN(x = as.vector(t(Ytemp1[clust[l],])), mean = gmmx1$mu[c[clust[l]],1:D1], Q = gmmx1$S[c[clust[l]],1:D1,1:D1], log =TRUE) + dMVN(x = as.vector(t(Ytemp2[clust[l],])), mean = gmmx2$mu[c[clust[l]],1:D2], Q = gmmx2$S[c[clust[l]],1:D2,1:D2], log = TRUE) + dnorm(x = Time[clust[l],1], mean = regy1$beta0[c[clust[l]]] + regy1$betahat[c[clust[l]],] %*% as.vector(t(Ytemp1.scaled[clust[l],])), sd = sqrt(regy1$sigma2[c[clust[l]]]), log =TRUE) + dnorm(x = Time[clust[l],1], mean = regy2$beta0[c[clust[l]]] + regy2$betahat[c[clust[l]],] %*% as.vector(t(Ytemp2.scaled[clust[l],])), sd = sqrt(regy2$sigma2[c[clust[l]]]), log = TRUE)
# } else{
# loglikelihood[j] <- loglikelihood[j] + dMVN(x = as.vector(t(Ytemp1[clust[l],])), mean = gmmx1$mu[c[clust[l]],1:D1], Q = gmmx1$S[c[clust[l]],1:D1,1:D1], log =TRUE) + dMVN(x = as.vector(t(Ytemp2[clust[l],])), mean = gmmx2$mu[c[clust[l]],1:D2], Q = gmmx2$S[c[clust[l]],1:D2,1:D2], log =TRUE) + log(dtruncnorm(x = Time[clust[l],1], a = Time[clust[l],1], b = Inf , mean = regy1$beta0[c[clust[l]]] + regy1$betahat[c[clust[l]],] %*% as.vector(t(Ytemp1.scaled[clust[l],])), sd = sqrt(regy1$sigma2[c[clust[l]]]))) + log(dtruncnorm(x = Time[clust[l],1], a = Time[clust[l],1], b = Inf, mean = regy2$beta0[c[clust[l]]] + regy2$betahat[c[clust[l]],] %*% as.vector(t(Ytemp2.scaled[clust[l],])), sd = sqrt(regy2$sigma2[c[clust[l]]])))
# }
# }
# }
clust <- which(c == activeclass[j])
luk1 <- c(0)
luk2 <- c(0)
for ( l in 1:length(clust)) {
luk1[l] <- dMVN(x = as.vector(t(Ytemp1[clust[l],])), mean = gmmx1$mu[c[clust[l]],1:D1], Q = gmmx1$S[c[clust[l]],1:D1,1:D1], log =TRUE) + dMVN(x = as.vector(t(Ytemp2[clust[l],])), mean = gmmx2$mu[c[clust[l]],1:D2], Q = gmmx2$S[c[clust[l]],1:D2,1:D2], log = TRUE)
luk2[l] <- dnorm(x = That[clust[l]], mean = regy1$beta0[c[clust[l]]] + regy1$betahat[c[clust[l]],] %*% as.vector(t(Ytemp1.scaled[clust[l],])), sd = sqrt(regy1$sigma2[c[clust[l]]]), log =TRUE) + dnorm(x = That[clust[l]], mean = regy2$beta0[c[clust[l]]] + regy2$betahat[c[clust[l]],] %*% as.vector(t(Ytemp2.scaled[clust[l],])), sd = sqrt(regy2$sigma2[c[clust[l]]]), log = TRUE)
}
GMMlikelihood[j] <- sum(luk1)
AFTlikelihood[j] <- sum(luk2)
loglikelihood[j] <- GMMlikelihood[j] + AFTlikelihood[j]
}
list('loglikelihood' = sum(loglikelihood),'GMMlikelihood' = sum(GMMlikelihood), 'AFTlikelihood' = sum(AFTlikelihood))
}
# clust <- which(c == activeclass[j])
# luk1 <- c(0)
# luk2 <- c(0)
# for ( l in 1:length(clust)) {
# if (Time[clust[l],2]==1){
# luk1[l] <- dMVN(x = as.vector(t(Ytemp1[clust[l],])), mean = gmmx1$mu[c[clust[l]],1:D1], Q = gmmx1$S[c[clust[l]],1:D1,1:D1], log =TRUE) + dMVN(x = as.vector(t(Ytemp2[clust[l],])), mean = gmmx2$mu[c[clust[l]],1:D2], Q = gmmx2$S[c[clust[l]],1:D2,1:D2], log = TRUE)
# luk2[l] <- dnorm(x = That[clust[l]], mean = regy1$beta0[c[clust[l]]] + regy1$betahat[c[clust[l]],] %*% as.vector(t(Ytemp1.scaled[clust[l],])), sd = sqrt(regy1$sigma2[c[clust[l]]]), log =TRUE) + dnorm(x = That[clust[l]], mean = regy2$beta0[c[clust[l]]] + regy2$betahat[c[clust[l]],] %*% as.vector(t(Ytemp2.scaled[clust[l],])), sd = sqrt(regy2$sigma2[c[clust[l]]]), log = TRUE)
# } else{
# luk1[l] <- dMVN(x = as.vector(t(Ytemp1[clust[l],])), mean = gmmx1$mu[c[clust[l]],1:D1], Q = gmmx1$S[c[clust[l]],1:D1,1:D1], log =TRUE) + dMVN(x = as.vector(t(Ytemp2[clust[l],])), mean = gmmx2$mu[c[clust[l]],1:D2], Q = gmmx2$S[c[clust[l]],1:D2,1:D2], log =TRUE)
# luk2[l] <- log(dtruncnorm(x = Time[clust[l],1], a = Time[clust[l],1], b = Inf, mean = regy1$beta0[c[clust[l]]] + regy1$betahat[c[clust[l]],] %*% as.vector(t(Ytemp1.scaled[clust[l],])), sd = sqrt(regy1$sigma2[c[clust[l]]]))) + log(dtruncnorm(x = Time[clust[l],1], a = Time[clust[l],1], b = Inf, mean = regy2$beta0[c[clust[l]]] + regy2$betahat[c[clust[l]],] %*% as.vector(t(Ytemp2.scaled[clust[l],])), sd = sqrt(regy2$sigma2[c[clust[l]]])))
# }
# }
# sum(luk1) + sum(luk2)