I'm trying to figure out if there is a way in R to get the loglikelihood of
a distribution fit to a set of data where the parameter values are fixed.
For example, I want to simulate data from a given alternate lognormal
distribution and then I will fit it to a lognormal distribution with null
parameter values to see what the likelihood of the null distribution is
given random data from the alternate distribution.
I have been using fitdistrplus for other purposes but I cannot use it to
fix both parameter values.
Here is an example of what I've been working with...
nullmu<-1.66 #set null mu
altmu<-1.58 #set alt mu
sd.log<-0.25 #set common sigma
cens.time<-6 #if simulated times are greater than this turn them into right
censored times
#simulating lognormal data (time) from altnative dist
(sim<-rlnorm(n=samplesize, meanlog=altmu, sdlog=sd.log))
#if the time was > cens.time replace time with cens.time
(sim[which(sim>cens.time)]<-cens.time)
sim
#create a variable indicating censoring
(cens<-sim)
cens[which(sim==cens.time)]<-NA
cens
#create the data frame to be passed to fitdistcens and fitdist
(x<-data.frame(left=sim,right=cens))
#if there is censored data use fitdistcens else use fitdist
ifelse(length(which(is.na(cens)))>0,
simfit<-fitdistcens(censdata=x, distr="lnorm"),
simfit<-fitdist(data=x[,1], distr="lnorm")
)
#Now I can get the loglikelihood of the MLE fitted distribution
simfit$loglik
#I want to get the loglikelihood of the distribution with the null
parameterization
#This is what I can't get to work
#I can't seem to find any function that allows me to set both parameter
values
#so I can get to loglikelihood of the of the parameterization given the data
nulldist<-fitdistcens(censdata=x, distr="lnorm", start=list(meanlog=nullmu,
sdlog=sd.log)
#Then I want to do a likelihood ratio test between the two distributions
pchisq((-2*simfit$loglik--2*nulldist$loglik), df=2, lower.tail=FALSE)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
