###################################################
### chunk number 1: 
###################################################
library(lattice)
options(width=65,show.signif.stars=FALSE)
options(SweaveHooks=list(fig=function() {
ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
ltheme$strip.background$col <- "transparent" ## change strip bg
ltheme$plot.symbol$col <- "black"
ltheme$axis.text$cex <- 1.2
ltheme$par.xlab.text$cex <- 1.2
ltheme$par.ylab.text$cex <- 1.2
lattice.options(default.theme = ltheme)      ## set as default
}))
source("formats.R")


###################################################
### chunk number 2: oxboys eval=FALSE
###################################################
## data(Oxboys,package="SMIR")


###################################################
### chunk number 3: oxboys
###################################################
data(Oxboys,package="nlme")


###################################################
### chunk number 4: oxboys.plot eval=FALSE
###################################################
## library(lattice)
## print(xyplot(height~I(age+13),xlab="age",group=Subject,data=Oxboys,type="b",lty=3))


###################################################
### chunk number 5: 
###################################################
library(lattice)
print(xyplot(height~I(age+13),xlab="age",group=Subject,data=Oxboys,type="b",lty=3))


###################################################
### chunk number 6: oxboys.gq20
###################################################
library(npmlreg)
data(Oxboys,package="nlme")
oxboys.gq20 <- allvc(height~age,random=~1|Subject,data=Oxboys,
random.distribution="gq",k=20,verbose=FALSE,plot.opt=0)
summary(oxboys.gq20)


###################################################
### chunk number 7: 
###################################################
sa2 <- (coef(oxboys.gq20)[3])^2
s2 <- oxboys.gq20$sdev$sdev^2


###################################################
### chunk number 8: 
###################################################
library(nlme)
#source('misc.R')
oxboys.lme <- lme(height~age,random=~1|Subject,data=Oxboys)


###################################################
### chunk number 9: oxboys.gq500 eval=FALSE
###################################################
## oxboys.gq500 <- update(oxboys.gq20,k=500)


###################################################
### chunk number 10: 
###################################################
summary(oxboys.gq500)


###################################################
### chunk number 11: oxboys.quad eval=FALSE
###################################################
## Oxboys <- transform(Oxboys,a2=age^2)
## oxboys.gq500quad <- update(oxboys.gq500,.~.+I(age^2))


###################################################
### chunk number 12: 
###################################################
summary(oxboys.gq500quad)


###################################################
### chunk number 13: 
###################################################
oxboys.np2 <- update(oxboys.gq500quad,random.distribution="np",k=2)
oxboys.np3 <- update(oxboys.np2,k=3)


###################################################
### chunk number 14: oxboys.np4_10
###################################################
oxboys.np4 <- update(oxboys.np3,k=4)
oxboys.np5 <- update(oxboys.np4,k=5)
oxboys.np6 <- update(oxboys.np5,k=6)
oxboys.np7 <- update(oxboys.np6,k=7)
oxboys.np8 <- update(oxboys.np7,k=8)
oxboys.np9 <- update(oxboys.np8,k=9)
oxboys.np10 <- update(oxboys.np9,k=10)
oxboys.gq20quad <- update(oxboys.np8,k=20,random.distribution="gq")


###################################################
### chunk number 15: oxboys.np8
###################################################
summary(oxboys.np8)


###################################################
### chunk number 16: oxboys.random.coefficients
###################################################
#oxboys.np8lin <- update(oxboys.np8,.~age)
oxboys.np8lin <- update(oxboys.np8,height~age,random=~1|Subject)
oxboys.np8linlin <- update(oxboys.np8lin,random=~age|Subject)
oxboys.np8quad <- update(oxboys.np8,height~age+I(age^2),random=~1|Subject)
oxboys.np8quadlin <- update(oxboys.np8quad,random=~age|Subject)
oxboys.np8quadquad <- update(oxboys.np8,random=~(age+I(age^2))|Subject)
#oxboys.np8rclin <- update(oxboys.np8rc,height~age+I(age^2),random=~age|Subject)
#oxboys.np8quadrcquad <- update(oxboys.np8rc,height~age+I(age^2),random=~(age+I(age^2))|Subject)


###################################################
### chunk number 17: oxboys.quadratic.linear
###################################################
#oxboys.np8quadrclin <- update(oxboys.np8rc,height~age+I(age^2),random=~age|Subject)
summary(oxboys.np8quadlin)


###################################################
### chunk number 18: oxboys.cubic
###################################################
oxboys.np8cubiclin <- update(oxboys.np8quadrclin,.~age+I(age^2)+I(age^3))
#summary(oxboys.np8cubicrclin)


###################################################
### chunk number 19: oxboys.seasonal
###################################################
Oxboys$Season <- Oxboys$Occasion
levels(Oxboys$Season) <- levels(Oxboys$Season)[c(1,2,3,4,1,2,3,4,1)]
oxboys.np8seasonrclin <- update(oxboys.np8quadlin,.~.+Season,tol=0.4)
round(coef(oxboys.np8seasonrclin),3)
round(oxboys.np8seasonrclin$disparity,2)


###################################################
### chunk number 20: 
###################################################
oxboys.post.prob <- post(oxboys.np8quadrclin,level="upper")$prob
boy.cluster <- which(round(oxboys.post.prob)==1,arr.ind=TRUE)[,2]
tabn <- apply(round(oxboys.post.prob),2,sum)
ord <- names(sort(tabn))
oxboys.post.prob <- oxboys.post.prob[names(boy.cluster),]
dimnames(oxboys.post.prob) <- list("boy  "=names(boy.cluster),"      K"=unique(boy.cluster))
print.pp <- function(x) print.noquote(ifelse(x<0.01,"    ",fc0(x)))
print.pp(oxboys.post.prob)


###################################################
### chunk number 21: oxboys.fitted.plot
###################################################
Oxboys$fitted <- predict(oxboys.np8quadrclin)
Oxboys$cluster <- factor(boy.cluster[Oxboys$Subject%in%names(boy.cluster)])
print(xyplot(height~age, group=Subject, data=Oxboys, 
             subscripts=TRUE,
             scales=list(x=list(at=c(-1,-0.5,0,0.5,1.0),labels=c(11,11.5,12,1.5,13))),
             panel=function(x, y, subscripts, ...) {
                 panel.superpose(x, y, subscripts,
                                 type="b",
                                 #pch=unclass(Oxboys$cluster), 
                                 lty=3, ...)#col=unclass(Oxboys$cluster),...)
                 panel.superpose(x, Oxboys$fitted,lty=1,
                                 subscripts, type="l", lwd=1.2,
,...)#                                 col=unclass(Oxboys$cluster),...)
             })) 


###################################################
### chunk number 22: 
###################################################
Oxboys$fitted <- predict(oxboys.np8quadrclin)
Oxboys$cluster <- factor(boy.cluster[Oxboys$Subject%in%names(boy.cluster)])
print(xyplot(height~age, group=Subject, data=Oxboys, 
             subscripts=TRUE,
             scales=list(x=list(at=c(-1,-0.5,0,0.5,1.0),labels=c(11,11.5,12,1.5,13))),
             panel=function(x, y, subscripts, ...) {
                 panel.superpose(x, y, subscripts,
                                 type="b",
                                 #pch=unclass(Oxboys$cluster), 
                                 lty=3, ...)#col=unclass(Oxboys$cluster),...)
                 panel.superpose(x, Oxboys$fitted,lty=1,
                                 subscripts, type="l", lwd=1.2,
,...)#                                 col=unclass(Oxboys$cluster),...)
             })) 


###################################################
### chunk number 23: 
###################################################
betablok <- read.csv("Data/betablock.csv")


###################################################
### chunk number 24: betablok eval=FALSE
###################################################
## data(betablok,package="SMIR")


###################################################
### chunk number 25: betablok.glm
###################################################
library(npmlreg)
(betablok.gq1 <- allvc(cbind(r,(n-r))~treat,data=betablok,random=~1|centre,
family=binomial,random.distribution="gq",k=1,verbose=FALSE,plot.opt=0))


###################################################
### chunk number 26: betablok.gq
###################################################
betablok.gq2 <- update(betablok.gq1,k=2)
betablok.gq3 <- update(betablok.gq1,k=3)


###################################################
### chunk number 27: 
###################################################
betablok.gq4 <- update(betablok.gq1,k=4)
betablok.gq5 <- update(betablok.gq1,k=5)
betablok.gq6 <- update(betablok.gq1,k=6)
betablok.gq7 <- update(betablok.gq1,k=7)
betablok.gq8 <- update(betablok.gq1,k=8)
betablok.gq9 <- update(betablok.gq1,k=9)
betablok.gq10 <- update(betablok.gq1,k=10)


###################################################
### chunk number 28: 
###################################################
betablok.gq20 <- update(betablok.gq1,k=20)


###################################################
### chunk number 29: segq
###################################################
distreatk1 <- update(betablok.gq1,.~.-treat)$disparity
distreatk2 <- update(betablok.gq2,.~.-treat)$disparity
distreatk3 <- update(betablok.gq3,.~.-treat)$disparity
distreatk4 <- update(betablok.gq4,.~.-treat)$disparity
distreatk5 <- update(betablok.gq5,.~.-treat)$disparity
distreatk6 <- update(betablok.gq6,.~.-treat)$disparity
distreatk7 <- update(betablok.gq7,.~.-treat)$disparity
distreatk8 <- update(betablok.gq8,.~.-treat)$disparity
distreatk9 <- update(betablok.gq9,.~.-treat)$disparity
distreatk10 <- update(betablok.gq10,.~.-treat)$disparity
distreatk20 <- update(betablok.gq20,.~.-treat)$disparity
sebeta1 <- abs(coef(betablok.gq1)[2]/sqrt(distreatk1-betablok.gq1$disparity))
sebeta2 <- abs(coef(betablok.gq2)[2]/sqrt(distreatk2-betablok.gq2$disparity))
sebeta3 <- abs(coef(betablok.gq3)[2]/sqrt(distreatk2-betablok.gq3$disparity))
sebeta4 <- abs(coef(betablok.gq4)[2]/sqrt(distreatk2-betablok.gq4$disparity))
sebeta5 <- abs(coef(betablok.gq5)[2]/sqrt(distreatk2-betablok.gq5$disparity))
sebeta6 <- abs(coef(betablok.gq6)[2]/sqrt(distreatk2-betablok.gq6$disparity))
sebeta7 <- abs(coef(betablok.gq7)[2]/sqrt(distreatk2-betablok.gq7$disparity))
sebeta8 <- abs(coef(betablok.gq8)[2]/sqrt(distreatk2-betablok.gq8$disparity))
sebeta9 <- abs(coef(betablok.gq9)[2]/sqrt(distreatk2-betablok.gq9$disparity))
sebeta10 <- abs(coef(betablok.gq10)[2]/sqrt(distreatk2-betablok.gq10$disparity))
sebeta20 <- abs(coef(betablok.gq20)[2]/sqrt(distreatk2-betablok.gq20$disparity))


###################################################
### chunk number 30: 
###################################################
(betablok.gq20.1 <- update(betablok.gq20,.~1))


###################################################
### chunk number 31: betablok.gq100 eval=FALSE
###################################################
## betablok.gq100 <- update(betablok.gq20,k=100)
## summary(betablok.gq100)


###################################################
### chunk number 32: betablok.np
###################################################
betablok.np1 <- update(betablok.gq1,random.distribution="np")
betablok.np2 <- update(betablok.gq2,random.distribution="np")
betablok.np3 <- update(betablok.np2,k=3)
betablok.np4 <- update(betablok.np2,k=4)


###################################################
### chunk number 33: betablok.np5-8
###################################################
betablok.np5 <- update(betablok.np2,k=5)
betablok.np6 <- update(betablok.np2,k=6)
betablok.np7 <- update(betablok.np2,k=7)


###################################################
### chunk number 34: betablok.npses
###################################################
distreatk2 <- update(betablok.np2,.~.-treat)$disparity
distreatk3 <- update(betablok.np3,.~.-treat)$disparity
distreatk4 <- update(betablok.np4,.~.-treat)$disparity
distreatk5 <- update(betablok.np5,.~.-treat)$disparity
distreatk6 <- update(betablok.np6,.~.-treat)$disparity
distreatk7 <- update(betablok.np7,.~.-treat)$disparity
sebeta1q <- abs(coef(betablok.np1)[2]/sqrt(distreatk1-betablok.np1$disparity))
sebeta2q <- abs(coef(betablok.np2)[1]/sqrt(distreatk2-betablok.np2$disparity))
sebeta3q <- abs(coef(betablok.np3)[1]/sqrt(distreatk3-betablok.np3$disparity))
sebeta4q <- abs(coef(betablok.np4)[1]/sqrt(distreatk4-betablok.np4$disparity))
sebeta5q <- abs(coef(betablok.np5)[1]/sqrt(distreatk5-betablok.np5$disparity))
sebeta6q <- abs(coef(betablok.np6)[1]/sqrt(distreatk6-betablok.np6$disparity))
sebeta7q <- abs(coef(betablok.np7)[1]/sqrt(distreatk7-betablok.np7$disparity))


###################################################
### chunk number 35: betablok.np5.summary
###################################################
summary(betablok.np5)


###################################################
### chunk number 36: betablok.np5rc
###################################################
betablok.np5rc <- update(betablok.np5,random=~treat|centre)
summary(betablok.np5rc)
round(betablok.np5rc$disparity,2)


###################################################
### chunk number 37: 
###################################################
print.pp2 <- function(x) print.noquote(ifelse(x>0.99,round(x),ifelse(x<0.01,"     ",round(x,2))))


###################################################
### chunk number 38: betablok.postprob
###################################################
betablok.post.prob <- apply(betablok.np5$post.prob,2,tapply,betablok$centre,mean)
betablok.post.prob <- betablok.post.prob[,rev(order(betablok.np5$masses))]
center.cluster <- sort(apply(betablok.post.prob,1,which.max))
betablok.post.prob <- betablok.post.prob[names(center.cluster),]
dimnames(betablok.post.prob) <- list("center  "=names(center.cluster),
"        K"=rev(order(betablok.np5$masses)))
print.pp2(betablok.post.prob)


###################################################
### chunk number 39: 
###################################################
treateff <- c(betablok.np8rc$coef[1]+betablok.np8rc$coef[10:16],betablok.np8rc$coef[1])
effmean <- sum(treateff*betablok.np8rc$masses)
treatsd <- sum(treateff^2*betablok.np8rc$masses)-effmean^2


###################################################
### chunk number 40: obesity
###################################################
Obese <- read.table("Data/obese.dat",skip=7,header=T)[,c(2,6,7,8,9)]
Obese$pattern <- apply(Obese[,2:4],1,paste,collapse="")
Obese[,2:4] <- data.frame(lapply(Obese[,2:4],function(x) ifelse(x==9,NA,x)))
Obese <- transform(Obese, 
o77 = o77 - 1,
o79 = o79 - 1,
o81 = o81 - 1)
expanded.Obese <- data.frame(lapply(Obese,
    function(x) rep(x, Obese$a8)))
Obesity <- reshape(expanded.Obese,varying=list(names(expanded.Obese)[2:4]),
                   times=c(8,10,12),timevar="age",
                   v.names = "status",direction="long")
Obesity$id2 <- rep(1:1014,3)
Obesity$sex <- Obesity$sex-1
Obesity$age1 <- (Obesity$age-10)/2
Obesity$age2 <- (Obesity$age1^2)*3-2


###################################################
### chunk number 41: obesity
###################################################
Obesity <- na.omit(Obesity)
library(npmlreg)


###################################################
### chunk number 42:  eval=FALSE
###################################################
## data(Obesity,package="SMIR")


###################################################
### chunk number 43: obesity.np1
###################################################
library(npmlreg)
obesity.np1 <- allvc(cbind(status,(1-status))~sex*(age1+age2),
                     data=Obesity,
                     tol=0.2,
family=binomial,random=~1|id2,k=1,
plot.opt=0,verbose=FALSE)


###################################################
### chunk number 44: 
###################################################
summary(obesity.np1)
fc2(obesity.np1$disparity)


###################################################
### chunk number 45: obesity.np2 eval=FALSE
###################################################
## obesity.np2 <- update(obesity.np1,k=2)


###################################################
### chunk number 46: 
###################################################
summary(obesity.np2)
fc2(obesity.np3$disparity)


###################################################
### chunk number 47: obesity.np3 eval=FALSE
###################################################
## obesity.np3 <- update(obesity.np1,k=3)


###################################################
### chunk number 48: 
###################################################
summary(obesity.np3)
fc2(obesity.np3$disparity)


###################################################
### chunk number 49: obesity.np4 eval=FALSE
###################################################
## obesity.np4 <- update(obesity.np1,k=4)


###################################################
### chunk number 50: 
###################################################
summary(obesity.np4)
fc2(obesity.np4$disparity)


###################################################
### chunk number 51: obesity.np2models eval=FALSE
###################################################
## obesity.np2a <- update(obesity.np2,.~.-sex:age2)
## obesity.np2b <- update(obesity.np2a,.~.-sex:age1)
## obesity.np2c <- update(obesity.np2b,.~.-age2)
## obesity.np2d <- update(obesity.np2c,.~.-sex)
## obesity.np2e <- update(obesity.np2d,.~.-age1)


###################################################
### chunk number 52: 
###################################################
print(obesity.np2d,digits=3)


###################################################
### chunk number 53: 
###################################################
## invlogit <- function(x) 1/(1+exp(-x))
## Obesity$fval <- invlogit(obesity.np2d$ebp)
## tab1 <- with(Obesity,tapply(fval,list(id,age1,status),mean,na.rm=T))
## obese.pp$id <- as.numeric(sapply(strsplit(rownames(obese.pp),"\\."),"[",1))
## obese.pp$age <- as.numeric(sapply(strsplit(rownames(obese.pp),"\\."),"[",2))
## obese.pp$sex <- obese.pp$id%%2
## obese.pp$patt <- floor((obese.pp$id+1)/2)
## Obesepp <- matrix(NA,nrow=52,ncol=2)
## Obesepp[obese.pp$id[1:36],1] <- round(obese.pp[1:36,1],3)
## Obesepp[obese.pp$id[37:72],2] <- round(obese.pp[37:72,2],3)
## #Obesepp[obese.pp$id[73:108],3] <- round(obese.pp[73:108,3],3)
## Obesepp <- as.data.frame(Obesepp)
## Obesepp$patt <- apply(Obese[,2:4],1,paste,collapse=" ")
## Obesepp$patt <- gsub("9","-",Obesepp$patt)
## Obesepp$patt <- gsub("1","0",Obesepp$patt)
## Obesepp$patt <- gsub("2","1",Obesepp$patt)
## Obesepp <- Obesepp[,c(4,1,2,3)]
## names(Obesepp) <- c("Pattern","1","2","3")
## Obesepp$sex <- rep(c(0,1),26)
## Obese.pp.male <- subset(Obesepp, sex==0)


###################################################
### chunk number 54:  eval=FALSE
###################################################
## ## These fitted probabilities do not match GLIMs
## Obesity$fval <- fitted(obesity.np2d)
## Obesity$Pattern <- gsub("9"," -",as.character(Obesity$pattern))
## Obesity$Pattern <- gsub("1"," 0",as.character(Obesity$Pattern))
## Obesity$Pattern <- gsub("2"," 1",as.character(Obesity$Pattern))
## ptab <- sort(with(Obesity,tapply(fval,list("Pattern   "=Pattern),mean,na.rm=T)))
## cbind(round(ptab[1:13],3))
## cbind(round(ptab[14:26],3))


###################################################
### chunk number 55: obesity_fitted eval=FALSE
###################################################
## inv.logit <- function(eta) exp(eta)/(1+exp(eta))
## fval <- inv.logit(as.numeric(names(table(obesity.np2d$linear.predictors))))
## age <- rep(c(8,10,12),2)
## group <- gl(2,3)
## prob <- rep(obesity.np2d$masses,c(3,3))
## paop <- as.vector(apply(matrix(prob*fval,nrow=3,ncol=2),1,sum))
## print(xyplot(fval~age,groups=group,ylab="proportion",
##              scales=list(y=list(at=seq(0,1,by=0.1),labels=seq(0,1,by=0.1))),
##              panel=function(x,y,...){
##                panel.superpose.2(x,y,type="l",lty=2,lwd=1.3,...)
##               panel.lines(c(8,10,12),paop,lwd=1.3,...)
##              }))


###################################################
### chunk number 56: 
###################################################
inv.logit <- function(eta) exp(eta)/(1+exp(eta))
fval <- inv.logit(as.numeric(names(table(obesity.np2d$linear.predictors))))
age <- rep(c(8,10,12),2)
group <- gl(2,3)
prob <- rep(obesity.np2d$masses,c(3,3))
paop <- as.vector(apply(matrix(prob*fval,nrow=3,ncol=2),1,sum))
print(xyplot(fval~age,groups=group,ylab="proportion",
             scales=list(y=list(at=seq(0,1,by=0.1),labels=seq(0,1,by=0.1))),
             panel=function(x,y,...){
               panel.superpose.2(x,y,type="l",lty=2,lwd=1.3,...)
              panel.lines(c(8,10,12),paop,lwd=1.3,...)
             }))


###################################################
### chunk number 57:  eval=FALSE
###################################################
## data(woolson,package="SMIR")


###################################################
### chunk number 58: Woolson
###################################################
woolson <- read.table("Data/lairdAR.dat",skip=6)
names(woolson) <- c("x","y","sex","age")


###################################################
### chunk number 59: 
###################################################
woolson <- transform(woolson,
                     a1=(age-10)/2,
                     a2=3*((age-10)/2)^2 - 2)
woolson <- transform(woolson,
                     sa1=sex*a1,
                     sa2=sex*s2,
                     t = ifelse(age>8,1,0))
woolson <- transform(woolson,
                     tsex = t * sex,
                     ta1= t*a1,
                     tsa1=t*sex*a1)
woolson <- woolson[order(woolson$age,1:48),]
wt <- c(150,15,8,8,8,9,7,20,154,14,13,19,2,6,6,21)
Wt <- rep(wt,3)
##
Woolson <- data.frame(lapply(woolson,
    function(x) rep(x, Wt)))
Woolson$Child <- rep(1:460,3)


###################################################
### chunk number 60: obesity.allvc1 eval=FALSE
###################################################
## library(npmlreg)
## obesity.allvc1 <- allvc(y~sex+t+x+tsex+ta1 +tsa1,random=~t|Child,
## random.distribution='np',k=1,data=Woolson,family=binomial,verbose=FALSE,plot.opt=0)


###################################################
### chunk number 61: 
###################################################
summary(obesity.allvc1)


###################################################
### chunk number 62: obesity.allvc2 eval=FALSE
###################################################
## obesity.allvc2 <- update(obesity.allvc1,k=2)#,tol=0.1)


###################################################
### chunk number 63: 
###################################################
summary(obesity.allvc2)


###################################################
### chunk number 64: obesity.allvc3 eval=FALSE
###################################################
## obesity.allvc3 <- update(obesity.allvc1,k=3)


###################################################
### chunk number 65: 
###################################################
summary(obesity.allvc3)


###################################################
### chunk number 66: obesity.allvc4
###################################################
obesity.allvc3a <- update(obesity.allvc3,random=~1|Child)


###################################################
### chunk number 67: 
###################################################
summary(obesity.allvc3a)


###################################################
### chunk number 68: obesity.allvc3b
###################################################
obesity.allvc3b <- update(obesity.allvc3a,.~.-t)


###################################################
### chunk number 69: 
###################################################
summary(obesity.allvc3b)


###################################################
### chunk number 70: obesity.allvc3c
###################################################
obesity.allvc3c <- update(obesity.allvc3b,.~.-x)


###################################################
### chunk number 71: 
###################################################
summary(obesity.allvc3c)


###################################################
### chunk number 72: 
###################################################
lsat <- read.table('Data/lsat.txt')
names(lsat) <- c("y1","y2","y3","y4","y5")
lsat$wt6 <- c(3,6,2,11,1,1,3,4,1,8,0,16,
0,3,2,15,10,29,14,81,3,28,15,80,16,
56,21,173,11,61,28,298)
lsat$wt7 <- c(12,19,1,7,3,19,3,17,10,5,3,7,
7,23,8,28,7,39,11,34,14,51,15,90,6,
25,7,35,18,136,32,308)


###################################################
### chunk number 73: lsat eval=FALSE
###################################################
## data(lsat,package="SMIR")


###################################################
### chunk number 74: lsat.total.scores
###################################################
tapply(lsat$wt7,list("Number of items correct "=apply(lsat[,1:5],1,sum)),sum)


###################################################
### chunk number 75: 
###################################################
sapply(lsat[,1:5],function(x) sum(x*lsat$wt7)/1000)


###################################################
### chunk number 76: 
###################################################
y <- NULL
for (i in 1:dim(lsat)[1]){
y <- c(y,unlist(rep(lsat[i,1:5],lsat[i,7])))
}
Lsat <- expand.grid(Item=1:5,Person=1:sum(lsat[,7]))
Lsat$y <- y
Lsat <- transform(Lsat,Item=factor(Item),Person=factor(Person))


###################################################
### chunk number 77: lsat.gq20.rasch eval=FALSE
###################################################
## library(npmlreg)
## lsat.gq20.rasch <- allvc(y~Item-1,data=Lsat,family=binomial,
## random=~1|Person,random.distribution="gq",k=20,plot.opt=0,verbose=FALSE)


###################################################
### chunk number 78: 
###################################################
summary(lsat.gq20.rasch)


###################################################
### chunk number 79: lsat.np2.rasch eval=FALSE
###################################################
## lsat.np2.rasch <- allvc(y~Item,data=Lsat,family=binomial,random=~1|Person,
## random.distribution="np",k=2,plot.opt=0,verbose=FALSE)


###################################################
### chunk number 80: 
###################################################
summary(lsat.np2.rasch)


###################################################
### chunk number 81: lsat.np3.rasch eval=FALSE
###################################################
## lsat.np3.rasch <- update(lsat.np2.rasch,k=3)


###################################################
### chunk number 82: 
###################################################
summary(lsat.np3.rasch)


###################################################
### chunk number 83: lsat.np4.rasch eval=FALSE
###################################################
## lsat.np4.rasch <- update(lsat.np2.rasch,k=4)


###################################################
### chunk number 84: 
###################################################
summary(lsat.np4.rasch)


###################################################
### chunk number 85: lsat.np5-6.rasch eval=FALSE
###################################################
## lsat.np5.rasch <- update(lsat.np2.rasch,k=5)
## lsat.np6.rasch <- update(lsat.np2.rasch,k=6)


###################################################
### chunk number 86: 
###################################################
lsat2m <- sum(lsat.np2.rasch$masses*lsat.np2.rasch$mass.points)
lsat2.sigma <- sqrt(sum(lsat.np2.rasch$masses*lsat.np2.rasch$mass.points^2) -
                    lsat2m^2)
lsat2.coef <- c("Item1"=1.868,1.868+coef(lsat.np2.rasch)[1:4])
lsat3m <- sum(lsat.np3.rasch$masses*lsat.np3.rasch$mass.points)
lsat3.sigma <- sqrt(sum(lsat.np3.rasch$masses*lsat.np3.rasch$mass.points^2) -
                    lsat3m^2)
lsat3.coef <- c("Item1"=1.868,1.868+coef(lsat.np3.rasch)[1:4])
lsat4m <- sum(lsat.np4.rasch$masses*lsat.np4.rasch$mass.points)
lsat4.sigma <- sqrt(sum(lsat.np4.rasch$masses*lsat.np4.rasch$mass.points^2) -
                    lsat4m^2)
lsat4.coef <- c("Item1"=1.868,1.868+coef(lsat.np4.rasch)[1:4])
lsat5m <- sum(lsat.np5.rasch$masses*lsat.np5.rasch$mass.points)
lsat5.sigma <- sqrt(sum(lsat.np5.rasch$masses*lsat.np5.rasch$mass.points^2) -
                    lsat5m^2)
lsat5.coef <- c("Item1"=1.868,1.868+coef(lsat.np5.rasch)[1:4])
lsat6m <- sum(lsat.np6.rasch$masses*lsat.np6.rasch$mass.points)
lsat6.sigma <- sqrt(sum(lsat.np6.rasch$masses*lsat.np6.rasch$mass.points^2) -
                    lsat6m^2)
lsat6.coef <- c("Item1"=1.868,1.868+coef(lsat.np6.rasch)[1:4])