Dear R-list
I have been trying to make survreg fit a normal regression model with left
truncated data, but unfortunately I am not able to figure out how to do it.
The following survreg-call seems to work just fine when the observations are
right censored:
library(survival)
n<-100000
#censored observations
x<-rnorm(n)
y<-rnorm(n,mean=x)
d<-data.frame(x,y)
d$ym<-pmin(y,0.5)
d$c<-ifelse(d$ym==0.5,0,1)
s<-survreg(Surv(ym,c)~x, dist="gaussian",data=d)
The truncated observations can be fitted with the following naive code.
Being able to use survreg (or similar) would however be great.
#trucated observations
d$c<-0.1
e<-d[d$y>d$c,]
e$cc<-1
#ss<-survreg(Surv(time=c,time2=y,event=cc,type="right")~x,
dist="gaussian",data=e)
a<-function(x,y,alpha=0,beta=1,sigma=1,s)
{
sum(
log((dnorm(y,mean=alpha+beta*x,sd=sigma^2)/(1-pnorm(q=s,mean=alpha+beta*x,sd=sigma^2))))
)
}
ll<-function(x,y,alpha=0,beta=1,sigma=1,s)
{
aa<-function(z){a(z[1],z[2],alpha=alpha,beta=beta,sigma=sigma,s=s)}
v<-apply(cbind(x,y),1,aa)
sum(v)
}
logl<-function(pars)
{ -a(x=e$x,y=e$y,alpha=pars[1],beta=pars[2],sigma=pars[3],s=0.1)
}
ss<-optim(c(0,0,1),logl)
print(s)
print(ss)
Best Regards
Per Jensen
[[alternative HTML version deleted]]
Dear R-list
I have been trying to make survreg fit a normal regression model with left
truncated data, but unfortunately I am not able to figure out how to do it.
The following survreg-call seems to work just fine when the observations are
right censored:
library(survival)
n<-100000
#censored observations
x<-rnorm(n)
y<-rnorm(n,mean=x)
d<-data.frame(x,y)
d$ym<-pmin(y,0.5)
d$c<-ifelse(d$ym==0.5,0,1)
s<-survreg(Surv(ym,c)~x, dist="gaussian",data=d)
The truncated observations can be fitted with the following naive code.
Being able to use survreg (or similar) would however be great.
#trucated observations
d$c<-0.1
e<-d[d$y>d$c,]
e$cc<-1
#ss<-survreg(Surv(time=c,time2 =y,event=cc,type="right")~x,
dist="gaussian",data=e)
a<-function(x,y,alpha=0,beta=1,sigma=1,s)
{
sum(
log((dnorm(y,mean=alpha+beta*x,sd=sigma^2)/(1-pnorm(q=s,mean=alpha+beta*x,sd=sigma^2))))
)
}
ll<-function(x,y,alpha=0,beta=1,sigma=1,s)
{
aa<-function(z){a(z[1],z[2],alpha=alpha,beta=beta,sigma=sigma,s=s)}
v<-apply(cbind(x,y),1,aa)
sum(v)
}
logl<-function(pars)
{ -a(x=e$x,y=e$y,alpha=pars[1],beta=pars[2],sigma=pars[3],s=0.1)
}
ss<-optim(c(0,0,1),logl)
print(s)
print(ss)
Thx for any help
Per Jensen
[[alternative HTML version deleted]]