## ----echo=FALSE,eval=FALSE---------------------------------------------------- # options(width=80) ## ----eval=FALSE--------------------------------------------------------------- # library(catdata) # data(medcare) # # library(MASS) # library(pscl) ## ----eval=FALSE--------------------------------------------------------------- # medmale <- medcare[medcare$male==1,] # medmale <-medmale[medmale$ofp<=30,] ## ----eval=FALSE--------------------------------------------------------------- # set.seed(5) # # subs<-1:nrow(medmale) # # reps <- 50 # # squerror1<-squerror3<-squerror4<-squerror5<-squerror6<-squerror7<-c() # abserror1<-abserror3<-abserror4<-abserror5<-abserror6<-abserror7<-c() # # # for(i in 1:reps){ # learn <- sample(subs,600) # test <- subs[-learn] # # ########################### # med=glm(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # family=poisson,data=medmale[learn,]) # l <- predict(med, newdata=medmale[test,], type="response") # # a<-rep(0,length(medmale[test,1])) # # for(j in 1:length(medmale[test,1])){ # while(ppois(a[j],lambda=l[j])<0.5){ # a[j] <- a[j] + 1 # } # } # diffs <- a - medmale[test,1] # squerror1[i]<-mean(diffs^2) # abserror1[i]<-mean(abs(diffs)) # # ###################### # med3=glm.nb(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # data=medmale[learn,]) # l <- predict(med3, newdata=medmale[test,], type="response") # a<-rep(0,length(medmale[test,1])) # # for(j in 1:length(medmale[test,1])){ # # while((pnbinom(a[j],mu=l[j],size=med3$theta))<0.5){ # a[j] <- a[j]+1 # } # } # # diffs <- a - medmale[test,1] # squerror3[i]<-mean(diffs^2) # abserror3[i]<-mean(abs(diffs)) # # ################## # med4=zeroinfl(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school|1, # data=medmale[learn,]) # pii <- 1-predict(med4, newdata=medmale[test,], type="zero") # mui <- predict(med4, newdata=medmale[test,], type="count") # a<-rep(0,length(medmale[test,1])) # # for(j in 1:length(medmale[test,1])){ # cdis <- 0 # # while(cdis<0.5){ # # cdis <- cdis + pii[j]*exp(-mui[j])*((mui[j]^a[j])/factorial(a[j])) + (1-pii[j])* # I(a[j]==0) # a[j] <- a[j]+1 # } # a[j]<-a[j] - 1 # } # diffs <- a - medmale[test,1] # squerror4[i]<-mean(diffs^2) # abserror4[i]<-mean(abs(diffs)) # # ################################ # med5=zeroinfl(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # data=medmale[learn,]) # # pii <- 1-predict(med5, newdata=medmale[test,], type="zero") # mui <- predict(med5, newdata=medmale[test,], type="count") # a<-rep(0,length(medmale[test,1])) # # for(j in 1:length(medmale[test,1])){ # cdis <- 0 # # while(cdis<0.5){ # # cdis <- cdis + pii[j]*exp(-mui[j])*((mui[j]^a[j])/factorial(a[j])) + (1-pii[j])* # I(a[j]==0) # a[j] <- a[j]+1 # } # a[j]<-a[j] - 1 # } # # diffs <- a - medmale[test,1] # # squerror5[i]<-mean(diffs^2) # abserror5[i]<-mean(abs(diffs)) # # #################### # med6=hurdle(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school|1, # data=medmale[learn,]) # mui <- predict(med6, newdata=medmale[test,], type="count") # gammai <- predict(med6, newdata=medmale[test,], type="zero") # pii2<-1-gammai*(1-exp(-mui)) # # # fa<-function(z,a){ # ((z^a)/factorial(a))*exp(-z) } # # a<-rep(0,length(medmale[test,1])) # # for(j in 1:length(medmale[test,1])){ # cdis <- pii2[j] # # if(cdis<0.5){ # while(cdis<0.5){ # #while((1-cumdist)>eps){ # cdis <- cdis + fa(z=mui[j],a=(a[j]+1))*gammai[j] # a[j]<-a[j]+1 # #a <- a+1 # } # } # else{a[j]<-0} # } # # diffs <- a - medmale[test,1] # # squerror6[i]<-mean(diffs^2) # abserror6[i]<-mean(abs(diffs)) # # ####################### # med7=hurdle(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # data=medmale[learn,]) # mui <- predict(med7, newdata=medmale[test,], type="count") # gammai <- predict(med7, newdata=medmale[test,], type="zero") # pii2<-1-gammai*(1-exp(-mui)) # # fa<-function(z,a){ # ((z^a)/factorial(a))*exp(-z) } # # a<-rep(0,length(medmale[test,1])) # # for(j in 1:length(medmale[test,1])){ # cdis <- pii2[j] # # if(cdis<0.5){ # while(cdis<0.5){ # #while((1-cumdist)>eps){ # cdis <- cdis + fa(z=mui[j],a=(a[j]+1))*gammai[j] # a[j]<-a[j]+1 # } # } # else{a[j]<-0} # } # # diffs <- a - medmale[test,1] # # squerror7[i]<-mean(diffs^2) # abserror7[i]<-mean(abs(diffs)) # # } ## ----fig.width=12,eval=FALSE-------------------------------------------------- # par(mgp=c(0,3,0)) # boxplot(abserror1,abserror3,abserror4,abserror5,abserror6,abserror7, # names=c("Poisson","Negative\nBinomial","Zero\nInflated 1","Zero\nInflated 2", # "Hurdle 1","Hurdle 2"),cex=1.3,cex.axis=1.6) ## ----eval=FALSE--------------------------------------------------------------- # library(catdata) # data(medcare) # #attach(medcare) # # library(MASS) # library(pscl) # # medmale <- medcare[medcare$male==1,] # medmale <-medmale[medmale$ofp<=30,] # set.seed(5) # # subs<-1:nrow(medmale) # # reps <- 50 # # eps <- 1e-5 # # lrps <- lrps3 <- lrps4 <- lrps5 <- lrps6 <- lrps7 <- rep(0,reps) # # for(i in 1:reps){ # learn <- sample(subs,600) # test <- subs[-learn] # # # med=glm(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # family=poisson,data=medmale[learn,]) # # l <- predict(med, newdata=medmale[test,], type="response") # # # for(j in 1:length(medmale[test,1])){ # for(a in 0:100){ # #while((1-ppois(a,lambda=l[j]))>eps){ # lrps[i] <- lrps[i] + (ppois(a,lambda=l[j]) - I(medmale[test,1][j] <= a))^2 # #a <- a+1 # } # } # # # med3=glm.nb(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # data=medmale[learn,]) # l <- predict(med3, newdata=medmale[test,], type="response") # # for(j in 1:length(medmale[test,1])){ # #a <- 0 # for (a in 0:100){ # #while((1-pnbinom(a,mu=l[j],size=med3$theta))>eps){ # lrps3[i] <- lrps3[i] + (pnbinom(a,mu=l[j],size=med3$theta) - # I(medmale[test,1][j] <= a))^2 # #a <- a+1 # } # } # # # med4=zeroinfl(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school|1, # data=medmale[learn,]) # pii <- 1-predict(med4, newdata=medmale[test,], type="zero") # mui <- predict(med4, newdata=medmale[test,], type="count") # # # for(j in 1:length(medmale[test,1])){ # a <- 0 # cumdist<-(1-pii[j]) # #cumdist <- 0 # for (a in 0:100){ # #while((1-cumdist)>eps){ # cumdist <- cumdist + pii[j]*exp(-mui[j])*(mui[j]^a)*(factorial(a)^(-1)) # lrps4[i] <- lrps4[i] + as.numeric((cumdist - I(medmale[test,1][j] <= a))^2) # #a <- a+1 # } # } # # # med5=zeroinfl(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # data=medmale[learn,]) # pii <- 1-predict(med5, newdata=medmale[test,], type="zero") # mui <- predict(med5, newdata=medmale[test,], type="count") # # for(j in 1:length(medmale[test,1])){ # a <- 0 # cumdist<-(1-pii[j]) # #cumdist <- 0 # for (a in 0:100){ # #while((1-cumdist)>eps){ # cumdist <- cumdist + pii[j]*exp(-mui[j])*(mui[j]^a)*(factorial(a)^(-1)) # lrps5[i] <- lrps5[i] + as.numeric((cumdist - I(medmale[test,1][j] <= a))^2) # #a <- a+1 # } # } # # # med6=hurdle(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school|1, # data=medmale[learn,]) # mui <- predict(med6, newdata=medmale[test,], type="count") # gammai <- predict(med6, newdata=medmale[test,], type="zero") # pii2<-1-gammai*(1-exp(-mui)) # # #rr<-coef(med6,model="zero") # # #pii<-exp(rr)/(1+exp(rr)) # # #ga<-(1-pii)/(1-exp(-mui)) # # fa<-function(z,a){ # ((z^a)/factorial(a))*exp(-z) } # # for(j in 1:length(medmale[test,1])){ # cumdist <- pii2[j] # # for(a in 1:100){ # #while((1-cumdist)>eps){ # cumdist <- cumdist + fa(z=mui[j],a=a)*gammai[j] # lrps6[i] <- lrps6[i] + as.numeric((cumdist - I(medmale[test,1][j] <= a))^2) # # #a <- a+1 # } # } # # # # med7=hurdle(ofp ~ hosp+healthpoor+healthexcellent+numchron+age+married+school, # data=medmale[learn,]) # mui <- predict(med7, newdata=medmale[test,], type="count") # gammai <- predict(med7, newdata=medmale[test,], type="zero") # pii2<-1-gammai*(1-exp(-mui)) # # fa<-function(z,a){ # ((z^a)/factorial(a))*exp(-z) } # # for(j in 1:length(medmale[test,1])){ # cumdist <- pii2[j] # # for(a in 1:100){ # #while((1-cumdist)>eps){ # cumdist <- cumdist + fa(z=mui[j],a=a)*gammai[j] # lrps7[i] <- lrps7[i] + as.numeric((cumdist - I(medmale[test,1][j] <= a))^2) # # #a <- a+1 # } # } # # } # # # lrps <-lrps/length(medmale[test,1]) # lrps3 <-lrps3/length(medmale[test,1]) # lrps4 <-lrps4/length(medmale[test,1]) # lrps5 <-lrps5/length(medmale[test,1]) # lrps6 <-lrps6/length(medmale[test,1]) # lrps7 <-lrps7/length(medmale[test,1]) ## ----fig.width=12,eval=FALSE-------------------------------------------------- # par(mgp=c(0,3,0)) # boxplot(lrps,lrps3,lrps4,lrps5,lrps6,lrps7, # names=c("Poisson","Negative\nBinomial","Zero\nInflated 1","Zero\nInflated 2", # "Hurdle 1","Hurdle 2"),cex=1.3,cex.axis=1.6) # ## ----echo=FALSE,results='hide',eval=FALSE------------------------------------- # detach(package:pscl)