ben-domingue

make_grs<-function(nm,df,fam.id="family",id="subjectID",
                   threshold.seq,
                   plink="/hrsshare/analytic_sample/hrs_geno_final",
                   resid.fm=NULL) {
  dir.nm<-paste0("/tmp/",nm)
  system(paste("mkdir",dir.nm))
  getwd()->dir.save
  setwd(dir.nm)
  if (!is.null(resid.fm)) {
    paste(nm,"~",resid.fm)->resid.fm

ben-domingue

#get genetic risk score
read.table("/hrsshare/alz_context/pgs/alz_nc.profile",header=TRUE)->score
score[,c(2,6)]->score
names(score)<-c("subjectID","polygenic_score")

#rand fat files
load("/hrsshare/v2_hrs_linked.Rdata")
x[,c("subjectID","hhidpn","family","raracem","rahispan")]->ids
ids[ids$raracem=="1.white/caucasian" & ids$rahispan=="0. not hispanic",]->ids
merge(ids,score)->x

ben-domingue

pick_lambda<-function(file.name, #file.name needs to be the root of the plink files.
                      #for example, if you have plink_file.bim, etc, then pass file.name="plink_file"
                      work.dir=".",
                      prop.train=10, #this number is inverted and becomes the proportion of the training SNPs used to calculate A-hat, this is the basis for the blackout window.
                      prop.test=500, #this is inverted and multiplied by the number of training SNPs. each of the L draws of test snps is then of this size.
                      lambda.seq=seq(.01,.25,length.out=45), #this is the sequence of lambdas that gets computed
                      L=50, #this is the number of draws from the test snps that are taken.
                      blackout.window=10, #minimum number of indicies that the chosen SNPs will differ by,
                      num.proc=1 #if you specify more, the sampling of L draws of test snps will be done in parallel.
                      ) { 

ben-domingue

herit_helper<-function(twindata,var.nm1,var.nm2) { #data is pair-level data
  require(OpenMx)
  source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")
  library(foreign)
  ######################################################################################################
  #first univariate
  ACE_h2<-function(twindata,var.nm) {
    #twindata <- read.dta("C:/Users/boardman/Docs/Research/Data/MIDUS/wehby/twinsonly.dta")
    #twindata<-read.dta("/home/bd/Dropbox/jason_other/wehby/data/twinsonly.dta")
    #summary(twindata)

ben-domingue

#1
library("statnet")
data("faux.magnolia.high")
fmh <- faux.magnolia.high
fmh #discuss directedness

#
plot(fmh, displayisolates = FALSE)

#2

ben-domingue

#these two functions create a fairly well formatted regression table, completw with stars.


get.stars<-function(x) { #x is a mer object
  summary(x)@coefs->tab
  tab[,3]->tstat
  star<-rep("",length(tstat))
  star[abs(tstat)>1.96]<-"*"
  star[abs(tstat)>2.58]<-"**"
  star[abs(tstat)>3.29]<-"***"

ben-domingue

interplay<-function(x,std.time.in.item=FALSE,nspl=4,plot.den=TRUE,top.plot=TRUE,...) {
    ##x needs to have columns:
    ## item [item id]
    ## id [person id]
    ## diff [item difficulty]
    ## th [person theta]
    ## pv [irt-based p-value]
    ## rt [response time in metric you want to analyze]
    ##resp [item response]
    #####################################################################

ben-domingue

out<-list()
for (s1 in seq(0,.2,by=.01)) {
    print(s1)
    test<-numeric()
    for (i in 1:250) {
        N=1000
        h=sqrt(.6)
        e=sqrt(1-h^2)
        m=.35
        s0=sqrt(1-m^2)

ben-domingue

#this contains the funtion that does the ML estimation
ratio.test<-function(df,
                     hess=FALSE, #if FALSE, compute empirical hessian
                     ctl.list=list(maxit=1000,reltol=sqrt(.Machine$double.eps)/1000)
                     ) {
    ##
    anal_vcov <- function(pars, df){
        neg_hess <- function(pars,df) { #computes the negative hessian
            for (nm in names(pars)) assign(nm,pars[[nm]])
            mu <- m*df$E + pi0*df$g + pi1*df$g*df$E

ben-domingue

simfun<-function(N,a1,beta.inter=.1) {
    ff<-function(i,N,a1,a2,beta.inter) {
        add.error<-function(t,alpha) {
            ve<-var(t)/alpha-var(t)
            e<-rnorm(N,mean=0,sd=sqrt(ve))
            t+e
        }
        g<-rnorm(N)
        e<-rnorm(N)
        y<-.2*g+.5*e+beta.inter*g*e+rnorm(N)