求助大神，如何自己写程序算出coxph中的loglik值

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我用coxph做了一个回归，算出来beta为（-0.383，-0.302），然后我写了一个cox的似然函数，将beta的值带入始终算不出来coxph$loglik里面beta=final的那个loglik值。我用beta=init带入我的程序可以得出beta=init的那个loglik值。求大神指点是否程序有问题。

1. set.seed(8943)
   
2. for (i in 1:100)
   
3. {x1<-c(runif(100,0,1))} #协变量x1
   
4. for (i in 1:100)
   
5. {x2<-c(rbinom(100,1,0.5))}#协变量X2
   
6. for (i in 1:100)
   
7. {s<-c(runif(100,0,1))}#生存时间随机数S
   
8. g<-c(sample(1:3,100,replace=T,prob=c(0.1,0.1,0.8)))
   
9. for (i in 1:100){
   
10. if (g[i]==3) g[i]<-0
    
11. else if (g[i]==2) g[i]<-5
    
12. else g[i]<- -5}
    
13. t<- ((-log(s))/exp(x1*1+x2*0.5+g))^2#生存时间t
    
14. lt<-c()
    
15. for (i in 1:100)
    
16. {lt<-c(runif(100,0,1)*2.990)}#删失时间lt
    
17. deta<-c()
    
18. for (i in 1:100){
    
19. if (t[i]<=lt[i]) deta[i]<-1 else deta[i]<-0}#deta删失指示,1dead，0alive
    
20. ht<-c()
    
21. for (i in 1:100)
    
22. {ht[i]<-min(t[i],lt[i])}
    
23. x<- cbind(x1,x2)#协变量
    
24. ga<- c(rep(0,100))
    
25. mn<- data.frame(x1,x2,ht,deta,ga,g)#建立数据
    
26. data<-mn[order(mn[,3]),]#将数据按ht的大小排序
    
27. library(survival)
    
28. xs<-coxph(Surv(ht, deta)~x, data)
    
29. beta<-xs$coef# beta0
    
30. xbl<-cbind(beta[1]*data[,1],beta[2]*data[,2],data[,5])
    
31.    co<-rowSums(xbl)
    
32.    fxj<- c(rep(0,100))
    
33.    for (d in 1:100){
    
34.      if (data[d,4]==1){
    
35.        fxj[d]<-log(sum(exp(co[d:100])))
    
36.      }
    
37.    }
    
38.    yy<- sum(data[,4]*(co-fxj))

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：求助大神 Log Cox xph final 程序 如何

另外，我还参照coxph源代码里的程序写了一个算loglik的程序

1. person<- seq(100,1,-1)
   
2. strata<- c(rep(0,100))
   
3. strata[100]<- 1
   
4. loglik<-0
   
5. covar<- cbind(data[,1],data[,2])
   
6. means<- cbind(mean(covar[,1]),mean(covar[,2]))
   
7. temp1<-cbind(covar[,1]-means[1],covar[,2]-means[2])
   
8. for (i in person){
   
9.   if (strata[i]==1){
   
10.     nrisk<- 0
    
11.     denom<- 0
    
12.   }
    
13.     ndead<- 0
    
14.     nrisk<- nrisk+1
    
15.   zbeta<- data[i,5]
    
16.     for (p in 1:2){
    
17.       zbeta<- zbeta+beta[p]*temp1[i,p]
    
18.     }
    
19.     risk<- exp(zbeta)
    
20.     denom<- denom+risk
    
21.   if (data[i,4]==1){
    
22.     ndead<- ndead+1
    
23.     loglik<- loglik+zbeta
    
24.   }
    
25.   if (ndead>0){
    
26.     loglik<-loglik-log(denom)
    
27.   }
    
28. }
    
29.   

复制代码

1. person<- seq(100,1,-1)
   
2. strata<- c(rep(0,100))
   
3. strata[100]<- 1
   
4. loglik<-0
   
5. for (i in person){
   
6.   if (strata[i]==1){
   
7.     nrisk<- 0
   
8.     denom<- 0
   
9.   }
   
10.     ndead<- 0
    
11.     nrisk<- nrisk+1
    
12.   zbeta<- data[i,5]
    
13.     for (p in 1:2){
    
14.       zbeta<- zbeta+beta[p]*data[i,p]
    
15.     }
    
16.     risk<- exp(zbeta)
    
17.     denom<- denom+risk
    
18.   if (data[i,4]==1){
    
19.     ndead<- ndead+1
    
20.     loglik<- loglik+zbeta
    
21.   }
    
22.   if (ndead>0){
    
23.     loglik<-loglik-log(denom)
    
24.   }
    
25. }
    
26.   
    
27.   

复制代码

上面两段程序跑出来的结果是一样，但是第一段程序zbeta用的是beta*x算的，而第二段程序zbeta用的是beta*（x-mean（x））算的，为什么这两段程序跑出来的结果会是一样的呢？