##################################################################################### #Nonlinear mixed effects model to analyze the purchasing data (Application Section) #nlme function in R is used. # #Dependent variable is mjpt (marijuana purchase). Independent variable is price. ##################################################################################### #mjptdata is a data set with grouping by each participant (id). mjptdata<-groupedData(mjpt~price|id, data=original.data) #inivals is initial estimates of l, a, and b from function 'nls' (nonlienar model # without random effects). prenlmeout <- nls(mjpt~((l*price^b)/exp(a*price)), data=mjptdata, start=c(l=1,b=0,a=0)) lini<-summary(prenlmeout)$parameters[1] bini<-summary(prenlmeout)$parameters[2] aini<-summary(prenlmeout)$parameters[3] inivals<-c(lini,bini,aini) #Nonlinear mixed effects model fitting nlmeout <- nlme(mjpt~((l*price^b)/exp(a*price)),data=mjptdata, fixed=l+b+a~1, random=l~1, start=inivals, weights=varPower(form=~price), method="ML") #Parameter estimates lvals<-nlmeout$coefficients$fixed[[1]] bvals<-nlmeout$coefficients$fixed[[2]] avals<-nlmeout$coefficients$fixed[[3]] #Covariance matrix estimates for the estimates of l, a and b scaleval<-(summary(nlmeout)$tTable[[4]])^2/summary(nlmeout)$varFix[[1]] v1s<-summary(nlmeout)$varFix[[1]]*scaleval #Error variance of l estimate v3s<-summary(nlmeout)$varFix[[5]]*scaleval #Error variance of b estimate v2s<-summary(nlmeout)$varFix[[9]]*scaleval #Error variance of a estimate #And, corresponding covariances c12s<-summary(nlmeout)$varFix[[7]]*scaleval c13s<-summary(nlmeout)$varFix[[4]]*scaleval c23s<-summary(nlmeout)$varFix[[6]]*scaleval #Derived values, Pmax, Omax and average elasticity (Note: pbar=mean of prices) pmaxs<-(1+bvals)/avals omaxs<-(lvals*pmaxs^bvals)/exp(avals*pmaxs)*pmaxs elasts<-bvals-avals*pbar #Variance estimates for derived values varOmax<-((1+bvals)^2*((1+bvals)/avals)^(2*bvals)*exp(-2*(1+bvals))* (-2*avals*(1+bvals)*c12s*lvals+avals^2*v1s+(1+bvals)^2*lvals^2*v2s+2*avals*lvals* (avals*c13s-(1+bvals)*c23s*lvals)*log((1+bvals)/avals)+avals^2*lvals^2*v3s* log((1+bvals)/avals)^2))/avals^4 varPmax<-(-2*avals*(1+bvals)*c23s+(1+bvals)^2*v2s+avals^2*v3s)/avals^4 varElast<--2*c23s*pbar+pbar^2*v2s+v3s