# 2014.11.21 'The Heritability of Otitis Media, A Twin and Triplet Study' (Casselbrant et al), JAMA. 1999;282:2125-2130 http://jama.ama-assn.org/ set.no a b c z x=scan() 1 0.00000 0.04646 NA 2 2 0.05162 0.00000 NA 2 3 0.00000 0.09130 NA 2 4 0.09594 0.00000 NA 2 5 0.04968 0.04968 NA 2 6 0.00000 0.10791 NA 2 7 0.05399 0.05399 NA 2 8 0.08533 0.04192 NA 2 9 0.04136 0.08999 NA 2 10 0.04136 0.10635 NA 2 11 0.05847 0.09850 NA 2 12 0.10581 0.05663 0.07824 2 13 0.08881 0.08881 NA 2 14 0.17428 0.04566 NA 2 15 0.14453 0.08047 NA 2 16 0.06761 0.16549 NA 2 17 0.11491 0.12578 NA 2 18 0.00000 0.24394 NA 2 19 0.16949 0.09605 NA 2 20 0.05344 0.23664 NA 2 21 0.14898 0.14317 NA 2 22 0.00000 0.29493 NA 2 23 0.12482 0.17134 NA 2 24 0.23180 0.09955 NA 2 25 0.08598 0.25723 NA 2 26 0.10146 0.24432 0.03003 2 27 0.11662 0.24162 NA 2 28 0.30645 0.05914 NA 2 29 0.30729 0.08333 NA 2 30 0.10000 0.31511 NA 2 31 0.09802 0.32931 NA 2 32 0.26603 0.16837 NA 2 33 0.18287 0.32221 NA 2 34 0.45220 0.05311 NA 2 35 0.41667 0.09568 NA 2 36 0.09836 0.42012 NA 2 37 0.06414 0.48568 NA 2 38 0.22798 0.32623 NA 2 39 0.30296 0.25481 NA 2 40 0.39074 0.17004 NA 2 41 0.30248 0.26718 NA 2 42 0.32093 0.26634 NA 2 43 0.29124 0.29978 NA 2 44 0.42562 0.18839 NA 2 45 0.21677 0.41258 NA 2 46 0.22757 0.40203 NA 2 47 0.50673 0.13677 NA 2 48 0.24482 0.45118 NA 2 49 0.20313 0.50744 NA 2 50 0.23400 0.50041 NA 2 51 0.24090 0.53110 NA 2 52 0.40377 0.40028 NA 2 53 0.40254 0.50710 NA 2 54 0.26294 0.65562 NA 2 55 0.39004 0.55572 NA 2 56 0.59776 0.38134 NA 2 57 0.43266 0.55673 NA 2 58 0.32767 0.70868 0.25000 2 59 0.69925 0.41269 NA 2 60 0.00000 0.04117 NA 1 61 0.03228 0.03228 NA 1 62 0.03717 0.03717 NA 1 63 0.09752 0.00000 NA 1 64 0.09838 0.00000 NA 1 65 0.07907 0.02108 NA 1 66 0.06975 0.03084 NA 1 67 0.10703 0.00000 NA 1 68 0.06287 0.06287 NA 1 69 0.06444 0.06444 NA 1 70 0.02562 0.12445 NA 1 71 0.07740 0.07740 NA 1 72 0.13907 0.03477 NA 1 73 0.17406 0.00000 NA 1 74 0.04831 0.13420 NA 1 75 0.13450 0.05623 NA 1 76 0.03343 0.17425 NA 1 77 0.00000 0.20950 0.69704 1 78 0.10958 0.10958 NA 1 79 0.14985 0.08086 NA 1 80 0.05415 0.19120 NA 1 81 0.12329 0.12329 NA 1 82 0.17656 0.08457 NA 1 83 0.13803 0.12317 NA 1 84 0.08595 0.19656 NA 1 85 0.14574 0.14574 NA 1 86 0.15556 0.14583 NA 1 87 0.08226 0.22662 NA 1 88 0.15455 0.15455 NA 1 89 0.06063 0.25929 NA 1 90 0.09010 0.25702 NA 1 91 0.12958 0.21962 NA 1 92 0.19747 0.16203 NA 1 93 0.21547 0.15240 NA 1 94 0.09020 0.28196 NA 1 95 0.19175 0.22287 NA 1 96 0.17554 0.25755 NA 1 97 0.28947 0.14812 NA 1 98 0.31811 0.12224 NA 1 99 0.21977 0.22467 NA 1 100 0.06292 0.45618 NA 1 101 0.27197 0.27652 NA 1 102 0.23021 0.34164 NA 1 103 0.24155 0.35532 NA 1 104 0.22527 0.37731 NA 1 105 0.22548 0.37921 NA 1 106 0.17035 0.44525 NA 1 107 0.31745 0.38563 NA 1 108 0.34685 0.40480 NA 1 109 0.35294 0.42502 NA 1 110 0.46055 0.36495 NA 1 111 0.62998 0.21311 NA 1 112 0.33485 0.50917 NA 1 113 0.51975 0.40637 NA 1 114 0.52364 0.41404 NA 1 115 0.45992 0.49782 NA 1 116 0.53920 0.44400 NA 1 117 0.42180 0.58195 NA 1 118 0.61628 0.39535 NA 1 119 0.67966 0.40367 0.61315 1 120 0.62105 0.51296 NA 1 121 0.55990 0.65226 NA 1 122 0.80210 0.61965 NA 1 123 0.67050 0.76305 NA 1 124 0.75076 0.75982 NA 1 125 0.77882 0.78349 NA 1 126 0.82672 0.74788 NA 1 # # WIDE file ds = matrix(x,ncol=5,byrow=T) colnames(ds)=c("set.no", "a", "b", "c", "z") ; ds[1:20,]; tail(ds) sum(is.na(ds[,4])) ; sum(!is.na(ds[,4])) set.no a b c z # TALL file #proportion of time with MEE family=c(); proportion=c(); zygosity =c(); for ( i in 1:length(ds[,1]) ) { if( is.na(ds[i,4]) ) { family =c( family, rep(ds[i,1],2) ) ; proportion =c( proportion, ds[i,2:3] ) ; zygosity =c( zygosity, rep(ds[i,5],2) ) ; } if( !is.na(ds[i,4]) ) { family =c( family, rep(ds[i,1],3) ) ; proportion =c( proportion, ds[i,2:4] ) ; zygosity =c( zygosity, rep(ds[i,5],3) ) ; } } DS=cbind(family,proportion,zygosity) length(DS[,1]) DS[1:20,]; tail(DS) setwd("/Users/jameshanley/Dropbox/Courses/bios602/MultilevelModels/Otitis") write.csv(DS,"otitisDataTall.txt",row.names=F) ################################################## Variance components via R JAGS D=read.csv("otitisDataTall.txt") MODEL=readLines(,22) model { #data for( i in 1:n) { proportion[i] ~ dnorm(mu[i],inv.var.within) mu[i] <- mu.overall + alpha[ family[i] ] } # random.effects for( f in 1:n.families ) { alpha[f] ~ dnorm(0, inv.var.family) } #priors mu.overall ~ dnorm(0, 0.001) sigma.family ~ dunif(0,1000) sigma.within ~ dunif(0,1000) #structural var.family <- sigma.family * sigma.family var.within <- sigma.within * sigma.within icc <- var.family / (var.family + var.within) inv.var.family <- 1/var.family inv.var.within <- 1/var.within } writeLines(MODEL,'otitisModel.txt') library('rjags') par(mfrow=c(3,1),mar = c(4.5,4,3,1)) for (z in 1:2) { Z = z # 1 egg or 2 n=sum(D[,3]==Z) ; # 136 # 121 who=D[,3]==Z n.families = length( unique(D[who,1]) ) if(Z==1) FAMILY=D[who,1]-59 ; if(Z==2) FAMILY=D[who,1] ; FAMILY INITS=list( mu.overall = 0.5, alpha = rnorm(n.families), sigma.family = .1, sigma.within = .1) jags <- jags.model('otitisModel.txt', data = list(n=n, n.families=n.families, family=FAMILY , proportion=D[who,2]) , inits=INITS,n.chains = 1,n.adapt = 4000) update(jags, 20000) results=jags.samples(jags, c('mu.overall', 'var.family', 'var.within', 'icc'), 1000) str(results) print(summary(results$mu[1,,1]) ) print(summary(results$var.family[1,,1]) ) print(summary(results$var.within[1,,1]) ) print( summary( results$icc[1,,1] ) ) hist( results$icc[1,,1], main= paste(toString(Z),"egg(s)"),xlim=c(0,1), xlab="ICC" ) if(Z==1) ICC1=results$icc[1,,1] if(Z==2) ICC2=results$icc[1,,1] } # end z # posterior for difference in icc's icc.diff = ICC1 - ICC2 hist( icc.diff, main= "ICC Difference",xlim=c(0,1) ) round( quantile(icc.diff,c(0.025, 0.50, 0.975)) , 2)