1

I have a dataset consisting of a single variable which happens to be left censored (censoring point is 0). I believe that the latent variable (i.e. the variable before censoring takes place) more or less follows a normal distribution. How can I - using R - find the parameters of this distribution?

Given the wealth of R-packages I'm surprised I have not been able to find any that easily solves the problem at hand. Judging by the name the fitdistcens-function from the fitdistrplus-package might be useful in this context. But if I read the documentation correctly - which I doubt - the function requires two columns, one of which should contain the uncensored data:

censdata A dataframe of two columns respectively named left and right, describing each observed value as an interval. The left column contains either NA for left censored observations, the left bound of the interval for interval censored observations, or the observed value for non-censored observations. The right column contains either NA for right censored observations, the right bound of the interval for interval censored observations, or the observed value for non-censored observations.

Does this mean that the function can't be used for my purpose? If so what are the alternatives?

Help (possibly involving an example) is much appreciated.

6
  • 1
    please add an example dataset, the code you tried and the expected output
    – Ujjwal
    Oct 15, 2014 at 9:02
  • My data consists of a single vector with 200 observation, of which only 40 entries non-zero values.
    – Matt
    Oct 15, 2014 at 11:49
  • Eg y<-c(0,0,0,1.43,0,2.27,0,0,0,0,.84,0)
    – Matt
    Oct 15, 2014 at 11:55
  • Dunno how to use censdata, but if I remember correctly, the likelihood function for censored data is not too complicated. Each censored datum (assuming independence) contributes a factor integral(p(x), x, -infinity, a) where a is the cut-off point (0 as you described it) and p(x) is the density function. So that will be a term which you can express in terms of erf and, I guess, you can at least minimize the log-likelihood numerically even if you can't find an exact solution. I haven't worked out the details but I hope this is enough to get going. Oct 15, 2014 at 18:43
  • @RobertDodier I'm just surprised no one has implemented the MLE procedure already. This really isn't that exotic. But you are absolutely right; if no one comes up with an easier solution that is the way to go.
    – Matt
    Oct 16, 2014 at 7:30

2 Answers 2

0

I ran into the same thing and I believe you just have to specify the minimum value on the right side and NA on the left side when an observation has been censored.

For example, here's a simulation showing that much:

library(fitdistrplus)

# These are the actual parameters of the right-censored, normal distribution to be fit
mu <- 10; sd <- 10; 

# Number of samples to use
n <- 100000; 

# Minimum value to censor at
x.min <- 0

# Create the sample to fit
x <- rnorm(n, mu, sd)


d <- data.frame(
  left= ifelse(x <= x.min, NA, x),   # Assign left side as NA when censored
  right=ifelse(x <= x.min, x.min, x) # Assign right side as x.min when censored
)
fitdistcens(d, 'norm')
# Fitting of the distribution ' norm ' on censored data by maximum likelihood 
# Parameters:
#      estimate
# mean 10.01450
# sd   10.00456
0

I wrote up chunk of R-code that seems to do the job:

#Set censoring point
a<-0 

#Generate random censored data
x<-rnorm(1000,-1,2)
x[x<a]<-a
length(which(x>a))

#Log-likelihood function
ll<-function(u=0,s=1){
-sum(pnorm(a,mean=u,sd=s,log=TRUE)*(x==a)+dnorm(x,mean=u,sd=s,log=TRUE)*(x>a))
}

#Run MLE - change initiation values (u and s)
require("stats4")
est<-mle(ll,start=list(u=0,s=1),,method="L-BFGS-B",lower = c(-Inf, 0))

#Show fit
plot.ecdf(x)
xplot<-seq(from=a,to=max(x),length=100)
lines(xplot,pnorm(xplot,mean=est@coef[1],sd=est@coef[2]),col="green")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.