library(GenABEL)
load("erf.RData")
tra <- erf@phdata$qtbas
#erf@phdata$qt <- 23 + erf@phdata$qtbas*5 + 1.6*erf@phdata$sex + 0.25*erf@phdata$age - 0.0023*erf@phdata$age*erf@phdata$age + rnorm(695)
age2 <- erf@phdata$age*erf@phdata$age
summary(tra)
summary(glm(tra~sex+age+age2,data=erf@phdata))
summary(glm(tra~sex+age,data=erf@phdata))

# real data analysis
ls()
names(erf@phdata)
gkin <- ibs(erf[,autosomal(erf)],w="freq")
gkin[1:5,1:5]
h2 <- polygenic(tra,kin=gkin,data=erf)
h2$h2an
tras <- qtscore(tra,data=erf)
tras$lam
descriptives.scan(tras,sort="Pc1df")
grs <- trascore(h2.a$pgres,data=erf,clam=FALSE)
grs$lam
descriptives.scan(grs,sort="Pc1df")
mms <- mmscore(h2.a,data=erf)
descriptives.scan(mms,sort="Pc1df")

# small data analysis
load("erfsmall.RData")
tra <- erf
erfs <- erf[1:150,]
gkins <- ibs(erfs[,autosomal(erfs)],w="freq")
h2 <- polygenic(qt,kin=gkins,data=erfs)
h2$h2an
h2.a <- polygenic(qt~sex+age,kin=gkins,data=erfs)
h2.a$h2an
qts <- qtscore(qt~sex+age,data=erfs)
qts$lam
descriptives.scan(qts,sort="Pc1df")
grs <- qtscore(h2.a$pgres,data=erfs,clam=FALSE)
gts$lam
descriptives.scan(grs,sort="Pc1df")
mms <- mmscore(h2.a,data=erfs)
descriptives.scan(mms,sort="Pc1df")