AJPA Submission

This page gives the "R" code for analyses reported in a manuscript submitted to the American Journal of Physical Anthropology.

If you do not already have "R" you can get it from:

http://cran.r-project.org/

You will also need to download two specific "R" libraries, these being mvtnorm and mvord.

Terry Data

First things first, anything shown in fixed-width green text is intended to be copied and pasted into your R "Gui" (graphical user interface). Below is the Terry data in very fine print, because who would want to read this? Just copy and paste to your Rgui.

Terry=structure(list(ID=structure(c(691L,405L,418L,439L,527L,583L,630L,713L,717L,38L,61L,86L,113L,130L,143L,158L,174L,194L,201L,206L,210L,218L,223L,233L,237L,240L,246L,247L,248L,249L,251L,252L,253L,254L,255L,256L,258L,259L,260L,262L,263L,264L,265L,266L,267L,268L,270L,272L,273L,274L,275L,276L,277L,278L,279L,281L,282L,283L,284L,285L,
286L,290L,291L,292L,293L,296L,297L,299L,300L,301L,302L,303L,305L,306L,307L,308L,309L,310L,311L,312L,313L,315L,316L,317L,318L,319L,320L,321L,322L,324L,325L,326L,327L,328L,329L,330L,333L,338L,339L,340L,341L,342L,343L,345L,346L,347L,348L,349L,350L,351L,354L,355L,357L,358L,359L,360L,361L,364L,365L,366L,367L,368L,369L,371L,372L,373L,
374L,377L,381L,382L,384L,385L,386L,387L,388L,389L,390L,391L,392L,393L,394L,396L,397L,398L,400L,402L,403L,406L,407L,410L,411L,412L,413L,414L,419L,420L,421L,422L,424L,425L,426L,427L,428L,429L,430L,431L,433L,434L,435L,436L,437L,441L,442L,443L,444L,445L,447L,449L,450L,451L,452L,453L,454L,455L,457L,458L,459L,460L,461L,463L,464L,465L,
466L,467L,468L,469L,471L,472L,473L,474L,475L,476L,477L,479L,480L,482L,483L,485L,486L,487L,488L,490L,491L,492L,493L,494L,495L,496L,498L,499L,500L,501L,502L,503L,505L,506L,507L,508L,509L,510L,511L,513L,514L,515L,516L,517L,518L,519L,520L,521L,522L,523L,524L,525L,528L,529L,530L,531L,532L,533L,534L,535L,536L,537L,538L,540L,541L,542L,
543L,544L,545L,546L,547L,548L,549L,551L,552L,553L,554L,555L,556L,557L,558L,559L,561L,562L,563L,564L,565L,566L,568L,569L,570L,571L,572L,573L,574L,575L,576L,578L,579L,580L,581L,584L,585L,586L,588L,592L,593L,594L,595L,596L,597L,600L,601L,602L,603L,604L,606L,607L,608L,609L,612L,613L,614L,615L,616L,617L,618L,619L,621L,622L,623L,624L,
626L,631L,632L,633L,636L,637L,638L,639L,641L,642L,643L,644L,645L,646L,647L,648L,649L,650L,654L,656L,657L,658L,659L,660L,662L,663L,665L,666L,667L,670L,671L,672L,673L,674L,676L,677L,679L,682L,683L,686L,687L,688L,689L,693L,694L,696L,697L,698L,699L,700L,702L,703L,704L,705L,706L,707L,709L,710L,714L,715L,716L,718L,719L,720L,721L,722L,
723L,724L,725L,727L,728L,729L,730L,731L,732L,733L,736L,738L,741L,743L,744L,745L,747L,748L,749L,750L,751L,752L,753L,755L,757L,758L,759L,761L,763L,764L,766L,767L,768L,769L,770L,771L,772L,773L,775L,776L,777L,779L,780L,781L,782L,784L,785L,786L,787L,788L,789L,790L,2L,3L,4L,7L,10L,11L,12L,13L,15L,16L,17L,18L,19L,20L,21L,23L,24L,
26L,27L,28L,30L,31L,32L,33L,35L,37L,39L,40L,41L,42L,43L,45L,47L,48L,49L,50L,51L,52L,53L,54L,56L,58L,60L,62L,63L,64L,65L,67L,68L,69L,71L,72L,73L,74L,75L,76L,77L,78L,79L,80L,81L,83L,84L,87L,88L,89L,90L,93L,94L,95L,96L,97L,98L,99L,100L,101L,102L,103L,105L,106L,107L,109L,110L,111L,114L,115L,116L,117L,118L,119L,120L,
121L,122L,123L,127L,128L,131L,132L,133L,135L,137L,139L,140L,141L,144L,145L,148L,149L,150L,151L,152L,153L,154L,155L,156L,157L,159L,160L,161L,162L,163L,164L,167L,168L,169L,172L,175L,176L,177L,178L,180L,181L,182L,183L,185L,186L,187L,188L,189L,190L,192L,193L,195L,196L,197L,198L,199L,203L,204L,207L,208L,209L,211L,212L,213L,214L,215L,
216L,217L,219L,220L,221L,224L,225L,226L,227L,228L,229L,230L,231L,5L,6L,8L,9L,14L,22L,25L,29L,34L,36L,44L,55L,57L,59L,66L,70L,82L,91L,108L,112L,124L,125L,129L,134L,136L,138L,142L,146L,147L,165L,166L,170L,171L,173L,179L,184L,200L,205L,222L,232L,234L,235L,236L,238L,239L,241L,242L,243L,244L,245L,250L,257L,261L,271L,280L,287L,
288L,289L,295L,298L,304L,314L,323L,331L,332L,334L,335L,336L,337L,344L,352L,353L,356L,363L,370L,375L,376L,378L,379L,380L,383L,395L,399L,401L,404L,408L,409L,416L,417L,423L,432L,438L,440L,446L,448L,456L,462L,470L,478L,481L,489L,497L,504L,512L,526L,539L,550L,560L,567L,577L,582L,587L,589L,590L,591L,598L,599L,605L,610L,611L,620L,625L,
628L,629L,634L,635L,640L,651L,652L,653L,655L,664L,668L,675L,678L,680L,681L,684L,685L,690L,695L,701L,708L,711L,712L,726L,734L,735L,737L,739L,740L,742L,746L,754L,756L,760L,762L,765L,774L,778L,783L,791L),.Label=c("1000","1001","1004","1006","1007R","100RR","1010","1016R","101R","1020","1021","1023","1027","102R","1030","1031","1032",
"1033","1034","1035","1039","103R","1047","1048","1052R","1054","1057","1058","105R","1062","1063","1064","1068","106R","1076","1080R","1083","109","1092","1093","1094","1095","1096","1103R","1104","1105","1111","1113","1115","1116","1117","1119","1120","1122","1124R","1127","112R","1135","1138R","1139","114","1141","1143","1145","1148",
"114R","1150","1153","1155","1158RR","1159","1163","1164","1165","1166","1167","1168","1172","1173","1174","1176","117R","1183","1186","1187","119","1195","1196","1198","1199","11R","12","1200","1201","1203","1205","1206","1209","1210","1211","1212","1215","1222","123","1230","1231","1236","123R","1243","1244","1247","124R","125",
"1250","1251","1254","1255","1259","1264","1265","1268","1275","1287","1288R","1289R","129","1290","1298","12R","130","1306","1311","1314","1316R","1319","131R","1322","132R","1331","1332","1333","1337RR","134","1341","1343","1345R","1348R","1349","1350","1351","1354","1356","1362","1363","1368","1377","1379","138","1388","1390","1392",
"1394","1396","1397","139R","13R","1401","1402","1406","1415R","1417R","1419","142R","143","1443","1458","1466","1469","146R","1476","1483","1490","1495","14R","1500","1501","1503","1504","1506","1507","150R","1512","1517","152","1521","1523","1531","1533","1539","153R","154","1544","1551","1559","155R","156","1563","1565","1566",
"157","1572","1573","1574","1575","1579","1580","1587","159","1591","1597","1598","159R","160","1600","1602","1606","1607","1622","1624","1627","1632","164R","167","16R","16RR","171R","172","177R","17R","182","183R","187R","18R","190R","191R","192","194","195","196","19R","200","201","202","203","204","205","206R","207","208","209",
"20R","210","211","212","213","214","215","216","217","219","21R","220","221","222","224","226","227","228","229","22RR","230","231","233","234","235","237","23R","240R","241R","242","244","245","246","247R","24R","250","251","253R","255","256","257","258","259","25R","260","261","262","263","264","265","267","268","269","26R",
"271","272","273","274","275","277","278","279","27R","280","281","282","285","286","287","288","289R","28R","291","293R","294R","29R","302R","304","305","306","307","308","309","30R","310","313","314","315","316","317","318","319R","31R","320","321","321R","322","323","324","326","327","329","32R","330","331","332","335","336",
"337","33R","340","341","343","345","346R","346RR","347","348R","349R","34R","353","355","35R","360","362","365","366","367","368","369","370","371","372","373","374R","375","377","379","37R","381","382R","385","388","38R","39","391","392","393R","393RR","394","395","396","397","398","399","399R","39R","40","400","401","402","403",
"405R","406","409","411","412","414","415","416","419","41R","421","423","424","425","428","42R","43","430RR","431","432","433","434","436","437R","438","43R","440","441","443","444","447","450","451","451R","452","453","456","457","458","45R","460","462","464","465","466","468","469","46R","470","471","473","474","475","477",
"478","47R","480","481","481R","482","483","484","485","486","487","488","48R","490","491","492","494","495","496","498","49R","501","502","504","505","506","509","50R","511","512","515","516","517","518","519","51R","520","522","523","524","525","528","529","530","532","535","537","538","539","53R","54","540","541","544",
"545","546","548","549","551","555","557","559","55R","560","561","562","563","564","565","566","567","568","569","56R","570","571","572","573","574","575","576","578","579","57R","580","581","583","584","586","587","58RR","591","592","593","594","595","596","597","598","599","59R","601","603","606","608","60R","61","610",
"615","617","617R","619","621R","622R","62R","631","632","634","635","636","639","63R","63RR","640","644","645","648","649","64R","653","654","655","656","657R","658R","660","666","672","674","675","676","678","679","67R","680","681","683","684","685R","686","689","689R","68R","69","690","691","696","69R","6RR","703","705","706",
"707","70R","711","712","714","717","718","719","720","723","724","726","728R","729R","72R","731","734RR","735","736","742","743","745","745R","747","748","74R","750","751","755","75R","760","761","764","765","766","767","76R","770","772","774R","775","779R","77R","780","781","78R","792R","793","795","796","798","7R","8","800R",
"802","805","807R","811","812","815","816","817","818R","819","821","822","824","831","832","834R","835","836","837R","83R","84","840","844","847","85","850","851","853","856","859","862","869","871","873R","876","880","881","882","885","886","887","88R","88RR","891","896RR","898","89R","900R","902","903R","905","906","907","90R",
"912","913","915","920","921","926","927","928R","929","92R","932","934","936","93R","940","944R","948","949","94R","952","955","957","959","961","965","966","967","96R","970","978","979","97R","981","983","987","988","98R","992","994","995","996","997","998","999","9R"),class="factor"),Obs=structure(c(3L,3L,3L,3L,3L,1L,3L,1L,3L,3L,
3L,3L,3L,2L,3L,3L,3L,3L,3L,2L,2L,3L,3L,1L,2L,2L,3L,2L,3L,1L,1L,3L,3L,3L,3L,3L,3L,3L,3L,1L,2L,3L,3L,3L,3L,3L,2L,3L,2L,3L,3L,3L,2L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,2L,2L,3L,3L,3L,1L,3L,3L,3L,2L,1L,3L,3L,3L,3L,2L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,2L,2L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,2L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,1L,3L,3L,
3L,2L,1L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,3L,1L,3L,3L,3L,2L,3L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,1L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,1L,3L,1L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,
3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,
3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,2L,3L,3L,2L,3L,3L,2L,3L,2L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,3L,1L,2L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,2L,3L,
3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,2L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,3L,1L,3L,3L,3L,3L,3L,3L,3L,3L,3L),.Label=c("DJW","LWK","NPH"),class="factor"),Sex=c(1,0,0,0,1,1,0,1,0,1,0,1,1,1,0,1,1,1,0,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,1,0,1,1,1,0,1,1,1,0,1,0,1,1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,0,0,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,0,0,1,0,1,1,0,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,0,0,0,0,1,0,1,1,0,0,1,1,1,1,1,1,1,0,1,0,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,0,1,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,1,1,0,0,0,1,1,0,1,1,1,0,0,1,1,1,1,0,1,0,1,1,0,1,1,1,1,0,0,0,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,1,1,1,0,1,1,0,1,1,1,1,1,0,1,1,0,1,1,1,0,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,0,0,
1,1,1,1,1,1,0,1,1,0,1,1,1,0,0,0,1,1,1,0,0,1,0,1,0,0,0,0,0,1,0,1,1,0,0,0,1,0,0,0,1,1,0,0,0,1,1,0,1,0,0,0,0,0,0,1,1,1,0,0,0,1,0,1,0,1,1,0,1,0,1,1,1,1,0,1,0,1,0,0,1,0,0,0,1,0,1,0,0,1,1,1,1,1,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,1,1,1,0,0,1,1,0,1,1,0,0,1,0,0,1,1,0,0,0,0,0,0,1,0,0,0,0,1,1,1,1,1,1,1,0,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,1,0,0,0,0,1,0,1,0,1,1,0,0,1,0,1,0,0,0,0,1,1,1,1,0,0,1,1,0,0,1,0,1,0,0,0,1,0,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,1,0,1,1,1,0,0,0,0,1,
0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,1),VA=c(5,5,1,1,4,5,1,2,1,5,1,5,5,5,1,5,5,5,1,5,5,2,2,5,1,5,5,5,5,2,5,5,5,5,5,2,5,2,5,5,5,5,3,5,5,5,5,5,5,5,2,2,5,1,5,5,5,5,5,5,5,5,5,5,5,
2,3,2,1,5,5,5,5,5,5,5,5,4,2,5,5,5,1,5,5,1,5,5,5,1,5,5,5,1,4,5,5,1,5,3,5,5,5,5,5,5,5,5,5,2,5,1,5,4,5,5,1,1,5,5,4,5,5,5,1,5,5,2,5,5,5,5,5,5,5,5,4,5,5,5,1,5,3,5,5,5,4,5,5,5,5,2,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,2,5,5,5,5,4,5,5,5,5,2,2,2,5,1,5,5,2,5,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,1,1,1,5,2,5,1,1,5,2,5,5,1,5,5,5,2,5,5,5,5,3,1,5,5,5,5,5,1,5,5,5,5,1,1,1,2,1,2,1,5,5,1,1,5,5,5,5,5,5,5,1,5,1,3,5,5,5,5,3,1,5,5,5,5,5,5,5,5,5,5,1,5,1,5,1,1,5,5,5,5,5,5,5,5,5,1,1,5,5,
1,1,1,5,5,1,5,5,5,1,1,5,5,5,5,1,5,1,5,5,1,5,5,5,5,1,1,1,1,5,5,1,5,5,5,5,5,1,5,5,2,5,4,5,5,5,1,5,5,5,3,1,5,5,1,5,5,5,5,5,1,5,5,1,5,5,5,1,2,1,5,5,5,5,1,5,5,5,1,5,5,1,5,3,1,1,5,5,5,1,5,1,5,5,5,5,5,5,1,5,5,1,5,5,5,2,2,1,5,5,5,1,2,5,1,5,2,1,2,1,1,5,1,5,5,1,2,1,5,1,1,2,5,5,1,1,2,2,5,1,5,1,1,1,1,1,1,4,1,5,2,1,1,5,1,5,1,2,5,1,5,1,5,5,5,5,2,5,1,5,1,1,5,1,1,1,5,1,5,1,2,2,5,5,4,5,1,1,1,1,1,2,5,1,1,5,1,1,1,5,1,1,2,1,5,4,5,1,1,5,5,1,5,5,1,1,5,1,1,5,5,1,1,2,1,1,1,2,
1,2,1,1,5,5,5,5,5,5,5,1,2,5,2,1,1,1,1,1,5,5,5,1,1,1,5,5,1,1,1,1,5,1,5,1,5,5,1,1,5,1,5,1,2,1,1,5,5,5,5,2,1,5,5,1,1,5,1,5,1,1,1,5,1,5,5,5,1,1,1,1,1,1,5,1,1,1,1,1,5,5,5,1,5,5,2,1,1,1,2,5,1,1,1,3,1,1,1,1,1,2,1,1,5,1,5,2,1,5,1,1,1,1,1,5,1,1,2,1,1,5,5,1,1,1,3,4,5,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,5,5,1,2,1,1,5,2,1,1,1,5,5,1,1,5,1,1,1,1,1,1,1,2,2,1,1,1,5,1,1,1,1,5,5,1,1,1,1,1,5,5,1,1,5,5,1,1,2,1,1,1,1,1,1,1,1,1,5,1,1,2,5,2,1,1,1,5,1,
1,5,5,1,2,1,1,1,1,2,5,1,5,1,1,1,1,1,1,2,1,1,1,1,5),SPC=c(5,5,1,1,5,5,1,1,1,5,1,5,5,5,1,5,5,5,1,5,5,1,1,5,1,1,5,5,5,1,5,5,5,5,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,5,4,5,5,1,5,5,5,5,5,5,5,5,5,5,5,5,4,1,2,5,5,5,2,5,5,5,5,5,5,5,5,5,1,5,5,1,5,5,5,1,5,5,5,1,5,5,5,1,5,1,5,5,5,5,5,5,5,5,5,5,5,1,5,4,5,5,1,1,5,5,4,5,5,5,1,5,5,5,3,5,5,5,5,5,5,5,5,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,1,5,5,1,2,5,5,5,1,5,5,5,5,5,5,5,5,
5,5,5,1,1,2,5,5,5,1,1,5,1,5,5,1,5,5,5,5,5,5,5,5,1,1,5,5,5,5,5,1,5,5,5,5,1,1,1,2,1,5,1,5,5,1,1,5,5,5,5,5,5,5,1,5,1,1,5,5,5,5,1,1,5,5,5,5,5,5,5,5,5,5,1,5,1,5,1,1,5,5,5,5,5,5,5,5,5,1,1,5,5,1,1,1,5,5,1,5,5,2,1,1,5,5,5,5,1,5,1,5,5,1,5,5,5,5,1,1,1,2,5,5,1,5,5,5,5,5,1,5,5,5,5,5,5,5,5,1,5,5,5,5,1,5,5,1,5,5,4,5,5,1,5,5,1,5,5,5,1,4,1,5,5,5,5,1,5,5,5,1,4,5,1,5,1,1,1,5,5,5,1,1,1,5,5,5,5,5,5,1,5,5,1,5,5,5,1,1,1,5,5,5,1,1,5,1,5,1,1,1,1,1,5,1,5,2,1,5,1,5,1,1,1,5,5,1,
1,2,5,5,1,5,1,1,1,1,1,1,4,1,5,1,1,1,5,1,5,1,5,5,1,5,1,5,5,5,5,2,5,1,5,1,1,5,1,1,1,5,1,5,1,1,5,5,5,5,5,1,1,1,1,1,1,5,1,1,5,1,1,1,5,1,1,5,1,5,5,5,1,1,5,5,1,5,5,1,1,5,1,1,5,5,1,1,2,1,1,1,5,1,1,1,1,5,5,5,5,5,5,5,1,5,5,1,1,1,1,1,1,5,5,5,1,1,1,5,5,1,1,1,1,5,1,5,1,5,5,1,1,5,1,5,1,1,1,1,5,5,5,5,2,1,5,5,1,1,5,1,5,1,1,1,5,1,5,5,5,1,1,1,1,1,1,5,1,1,1,1,1,5,5,5,1,5,5,5,1,1,1,1,5,1,1,1,2,1,1,1,1,1,2,1,1,5,1,5,1,1,5,1,1,1,1,1,5,1,1,1,1,1,5,5,1,1,1,5,4,5,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,5,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,5,1,5,1,1,5,1,1,1,1,5,5,1,1,5,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,5,5,1,1,1,1,1,5,5,1,1,5,5,1,1,1,1,1,1,1,1,1,1,1,1,5,2,1,1,5,1,1,1,1,5,1,1,5,5,1,1,1,1,1,1,1,5,1,5,1,1,1,2,1,1,1,1,1,1,1,5),MA=c(5,5,1,1,5,5,1,1,1,5,1,5,5,5,1,5,5,5,1,5,5,1,1,5,1,5,5,5,5,1,5,5,5,5,5,5,5,1,5,5,5,4,5,5,5,5,2,5,5,5,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,5,4,1,1,5,5,5,2,5,5,5,5,5,5,5,5,5,1,5,5,1,5,5,5,1,5,5,5,1,5,5,5,1,5,3,5,5,4,5,2,5,5,5,
5,5,5,1,5,4,5,5,1,1,5,5,3,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,1,5,5,1,2,5,5,5,1,5,5,5,5,5,5,5,5,5,5,5,1,1,5,5,5,5,1,1,5,1,5,5,1,5,5,5,5,5,5,5,5,1,1,5,5,5,5,5,2,5,5,5,5,5,1,1,1,1,5,1,5,5,1,1,5,5,5,5,5,5,5,3,5,1,3,5,5,5,5,3,1,5,5,5,5,5,5,5,5,5,5,1,5,1,5,1,1,5,5,5,5,5,5,5,5,5,1,1,5,2,1,1,1,5,5,1,5,5,2,1,1,5,5,5,5,2,5,1,5,5,1,5,5,5,5,1,1,1,5,5,5,1,5,5,5,5,5,1,5,5,5,5,5,
5,5,5,1,5,5,5,4,1,5,5,1,5,5,5,5,5,1,5,5,1,5,5,5,1,3,1,5,5,5,5,1,5,5,5,1,5,5,1,5,3,1,1,5,5,5,1,1,1,5,5,5,5,5,5,1,5,5,1,5,5,5,1,1,1,5,5,5,1,1,5,1,5,1,1,1,1,1,5,1,5,5,1,5,1,5,1,1,3,5,5,1,1,1,5,5,1,5,1,1,1,1,1,1,5,1,5,1,1,1,5,1,5,1,5,5,1,5,1,5,5,5,5,2,4,1,5,1,1,5,1,1,1,5,1,5,1,1,5,5,5,5,5,1,1,1,1,1,1,5,1,1,5,1,1,1,5,1,1,5,1,5,5,5,1,1,5,5,1,5,5,1,1,4,1,1,5,5,1,1,3,1,1,1,5,1,1,1,1,5,5,5,5,5,5,5,1,5,5,1,1,1,1,1,1,5,5,5,1,1,1,5,5,1,1,1,1,5,2,5,1,5,5,1,1,5,1,5,
1,1,1,1,5,5,5,5,2,1,5,5,1,1,5,1,5,1,1,1,5,1,5,5,5,1,1,1,1,1,1,5,1,1,1,1,1,5,5,5,1,5,5,5,1,1,1,1,5,1,1,1,2,1,1,1,1,1,1,1,1,5,1,5,1,1,5,1,1,1,1,1,5,1,1,1,1,1,5,5,1,2,1,5,4,5,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,5,5,1,5,1,1,5,1,1,1,1,5,5,1,1,5,1,1,1,1,1,1,1,1,1,1,1,1,5,1,1,1,1,5,5,1,1,1,1,1,5,5,2,1,5,5,1,1,2,1,1,1,1,1,1,1,1,1,5,3,1,1,5,1,1,1,1,5,1,1,5,5,1,1,1,1,1,1,1,5,1,5,1,1,1,1,1,1,1,1,1,1,1,5)),row.names=c(NA,-774L),class="data.frame")

Observers

table(Terry[,2])

Picking the best (lowest BIC) logistic regression for the original coding and for the collapsed coding

The data above is for 774 individuals from the Terry Collection that have complete data on the Phenice characteristics. We will use the Bayesian Information Criterion (BIC) to pick the best model among the saturated model (three predictors), the three models with two predictors, the three models with one predictor, and the null model with no predictors.

all.BICS = function (dat=Terry) 
{
options(warn=-1)

# Fit models

full.model=glm(Sex~VA+SPC+MA,data=dat,family='binomial')

model1=glm(Sex~SPC+MA,data=dat,family='binomial')
model2=glm(Sex~VA+MA,data=dat,family='binomial')
model3=glm(Sex~VA+SPC,data=dat,family='binomial')

model4=glm(Sex~VA,data=dat,family='binomial')
model5=glm(Sex~SPC,data=dat,family='binomial')
model6=glm(Sex~MA,data=dat,family='binomial')

no.model=glm(Sex~1,data=dat,family='binomial')

# Get BIC for each model

BICs=vector()
BICs[1]=round(BIC(full.model),4)
BICs[2]=round(BIC(model1),4)
BICs[3]=round(BIC(model2),4)
BICs[4]=round(BIC(model3),4)
BICs[5]=round(BIC(model4),4)
BICs[6]=round(BIC(model5),4)
BICs[7]=round(BIC(model6),4)
BICs[8]=round(BIC(no.model),4)

# Store the (text) formula for each model

models=vector()
models[1]=noquote(deparse(formula(full.model)))
models[2]=noquote(deparse(formula(model1)))
models[3]=noquote(deparse(formula(model2)))
models[4]=noquote(deparse(formula(model3)))
models[5]=noquote(deparse(formula(model4)))
models[6]=noquote(deparse(formula(model5)))
models[7]=noquote(deparse(formula(model6)))
models[8]=noquote(deparse(formula(no.model)))

# Print the table

cat('\n')
for(i in 1:8) {cat(paste(BICs[i],models[i],'\n'))}

# Print the (best) model with the lowest BIC

cat('\nBest model:\n')
cat(models[which.min(BICs)])
cat('\n\n')
}


all.BICS()


The best logistic regression model is Sex~VA+SPC (to  the exclusion of MA).  Now let's get the confusion matrix for the 5-point scale using the logistic regression model Sex~VA+SPC:

logist = function (dat=Terry) 
{
N=NROW(dat)
cat(paste('\nData object name = ',deparse(substitute(dat)),'\n\n'))
Sexed_As=c('"Female"','"Male"')[as.numeric(predict(glm(Sex~VA+SPC,data=dat,family=binomial),type='response')>0.5)+1]
Sex=c('Female','Male')[dat$Sex+1]
out=table(data.frame(Sex,Sexed_As))
correct=round(sum(diag(out))/N*100,2)
print(out)
cat(paste('\nPercent correct = ',correct,'%\n\n',sep=''))
}

logist()

Logistic regression with "leave one out" cross-validation

Below is the "leave one out" logistic regression. See https://www.r-bloggers.com/calculate-leave-one-out-prediction-for-glm/ for a faster way using the "foreach" library (particularly if you have parallel processors).

LOO = function(dat=Terry[,3:5]){

N = NROW(dat)
Sexed = vector()
cat('\n')
for(i in 1:N){
 model = glm(Sex~VA+SPC,family='binomial',data=dat[-i,])
 Sexed[i] = c('"Female"','"Male"')[as.numeric(predict(model,newdata=dat[i,],type='response')>0.5)+1]
 cat(paste('\r',i,'of',N))
 flush.console()
}

Sex=c('Female','Male')[dat$Sex+1]
cat('\n\nLeave one out cross-validation\n')
table(data.frame(Sex,Sexed))
}

LOO()

Notice that the confusion matrix is unchanged from the "biased" one

sum(apply(Terry[Terry$Sex==0,4:5]==rep(1,2),1,sum)==2) # 291 females scored as "Female" for both the VA and SPC
sum(apply(Terry[Terry$Sex==1,4:5]==rep(5,2),1,sum)==2) # 372 males scored as "Male" for both the VA and SPC

0.632+ Estimator (some code from https://rdrr.io/cran/sortinghat/src/R/errorest-632plus.r)

my.632.plus = function () 
{
x = data.matrix(Terry[,4:5])
y = Terry[,3]

# Find apparent error

train_out=glm(y~x,family='binomial')
newdata=data.frame(x)
N = NROW(newdata)
classify_out=rep(1,N)
p=1/(1+exp(predict(train_out, newdata)))
for(i in 1:N) if(p[i]>0.5) classify_out[i]=0

apparent = (N-sum(diag(table(data.frame(y,classify_out)))))/N

    
LOO.boot=function () 
{
set.seed(42)

loo = matrix(NA,nr=50,nc=N)

for(i in 1:50){
   training = sample(1:774,replace=T)
   test = which(!(1:774 %in% training))
   newdata = x[test,]
   N=NROW(newdata)
   newdata = cbind(rep(1,N),newdata)
   model = as.vector(coef(glm(y[training]~x[training,],family='binomial')))
   model = rep(1,N)%o%model
   clas = apply(model*newdata,1,sum)
   clas=1/(1+exp(clas))
   clas = as.numeric(clas<0.5)
   loo[i,test] = clas!=y[test]
}

loo=colMeans(loo,na.rm=T)
return(mean(loo))
}

# Get the LOO bootstrap

loo_boot=LOO.boot()

n=length(y)

p_k <- as.vector(table(y))/n
q_k <- as.vector(table(classify_out))/n
gamma_hat <- drop(p_k %*% (1 - q_k))
R_hat <- (loo_boot - apparent)/(gamma_hat - apparent)
w_hat <- 0.632/(1 - 0.368 * R_hat)
six.32 = (1 - w_hat) * apparent + w_hat * loo_boot
six.32.plus = round(100*(1-six.32),2)
cat(paste('\n100% - .632+ error% = ',six.32.plus,'%',sep=''))
cat('\n')
}

my.632.plus()

This gives a percent correct of 97.98% as versus the "biased" percent correct of 98.06%.  Given that this is only 0.08% less, we will no longer concern ourselves with the notion that the model might be "over-trained."

Multivariate probit

This requires the library "mvord.".

mv.probit=function (dat=Terry[,3:5])
{
library(mvord)
options(warn=-1)
traits=colnames(dat)[-1]
n.traits=length(traits)
trait.names=vector()
for(i in 1:n.traits) trait.names[i]=traits[i]
if(n.traits==2){
mod=as.formula(paste('MMO2(',trait.names[1],',',trait.names[2],')','~Sex',sep=''))}
else{
mod=as.formula(paste('MMO2(',trait.names[1],',',trait.names[2],',',
    trait.names[3],')~Sex',sep=''))
}
cat(paste('\nModel is:','\n'))
print(mod,showEnv=F)
cat('\n')
flush.console()

model=mvord(mod,data=dat)

thresh1 = model$theta[[1]]-model$beta[1]
thresh2 = model$theta[[2]]-model$beta[2]
if(n.traits==3){
r = error_structure(model,type='corr')[[1]][lower.tri(diag(3))]
thresh3 = model$theta[[3]]-model$beta[3]
thetas = c(thresh1,thresh2,thresh3,model$beta[4:6],r)
last.name = max(which(names(thetas)!=''))
names(thetas)=c(names(thetas)[1:last.name],
     paste('r(',trait.names[1],', ',trait.names[2],')',sep=''),
     paste('r(',trait.names[1],', ',trait.names[3],')',sep=''),
     paste('r(',trait.names[2],', ',trait.names[3],')',sep=''))
}
else{
r = error_structure(model,type='corr')[[1]][1,2]
thetas=c(thresh1,thresh2,model$beta[3:4],r)
last.name = max(which(names(thetas)!=''))
names(thetas)=c(names(thetas)[1:last.name],paste('r(',trait.names[1],', ',trait.names[2],')',sep=''))

}

cat('Note: Intercepts subtracted from all thresholds, so first thresholds are NOT = 0!\n\n')
print(round(thetas,3))
cat('\n')
thetas<<-as.numeric(thetas)
}

mv.probit()

Predictions from the Multivariate Probit for the Training Data

patterns=function (dat=Terry[,3:5])
{
# This function counts data patterns
Sex=rep(c(rep(0,5),rep(1,5)),5)
VA=rep(1:5,10)
SPC=rep(1:5,each=10)
counts=as.numeric(ftable(dat))
out=cbind(Sex,VA,SPC,counts)
out=out[out[,4]>0,]
return(out)
}

mv.probit.predict = function (prior=F)
{
library(mvtnorm)
Patterns=patterns()
N.pats=NROW(Patterns)
lower=rep(NA,2)
upper=rep(NA,2)
thresh=matrix(NA,nr=2,nc=4)

lower=rep(NA,2)
upper=rep(NA,2)
x=thetas
thresh[1,]=x[1:4]
thresh[2,]=x[5:8]
mus=x[9:10]
corrs=x[11]
cor=diag(2)
cor[upper.tri(cor)]=corrs
cor[lower.tri(cor)]=corrs

prob.values = function(pat=1){

if(prior==F) prob.f=0.5
else {prob.f=347/774}

mu=rep(0,2)

mu=Patterns[pat,1]*mus
for(i in 1:2){
if(Patterns[pat,i+1]==1){
    lower[i]=-Inf
    upper[i]=thresh[i,1]
}
if(Patterns[pat,i+1]==2){
    lower[i]=thresh[i,1]
    upper[i]=thresh[i,2]
}
if(Patterns[pat,i+1]==3){
    lower[i]=thresh[i,2]
    upper[i]=thresh[i,3]
}
if(Patterns[pat,i+1]==4){
    lower[i]=thresh[i,3]
    upper[i]=thresh[i,4]
}
if(Patterns[pat,i+1]==5){
    lower[i]=thresh[i,4]
    upper[i]=Inf
}
}
prob1 = as.numeric(pmvnorm(mean=c(0,0),lower=lower,upper=upper,corr=cor,algorithm=Miwa)[1])
p.cond.f = prob.f*prob1
prob2 = as.numeric(pmvnorm(mean=mus,lower=lower,upper=upper,corr=cor,algorithm=Miwa)[1])
p.cond.m = (1-prob.f)*prob2
post.f=p.cond.f/(p.cond.f+p.cond.m)
return(post.f)
}

p=vector()

female=rep(0,N.pats)
male=rep(0,N.pats)

for(i in 1:N.pats){
p = prob.values(pat=i)
if(p>0.5) female[i]=Patterns[i,4]
if(p<0.5) male[i]=Patterns[i,4]
}

out=cbind(Patterns[,c(1,4)],female,male)

sexing=matrix(0,nc=2,nr=2)
for(i in 1:25){
if(out[i,1]==0){
    sexing[1,1]=sexing[1,1]+out[i,3]
    sexing[1,2]=sexing[1,2]+out[i,4]
}
else{
    sexing[2,2]=sexing[2,2]+out[i,4]
    sexing[2,1]=sexing[2,1]+out[i,3]
}
}
rownames(sexing)=c('Female','Male')
colnames(sexing)=c('"Female"','"Male"')
if(prior==F) cat('\nUniform prior\n')
else{cat('\nPrior from training sample\n')}
return(sexing)
}

mv.probit.predict(F)
mv.probit.predict(T)

Applying the Training Data to a Simulated Testing Dataset using Multivariate Probit

simulate = function (nf=25,nm=975)
{
library(mvtnorm)
x=thetas
thresh=rbind(c(-Inf,x[1:4],Inf),c(-Inf,x[5:8],Inf))

r = matrix(c(1,rep(x[11],2),1),nc=2)

set.seed(42)

scores=matrix(NA,nc=2,nr=1000)
latent = rbind(rmvnorm(nf,sigma=r,mean=c(0,0)),
               rmvnorm(nm,sigma=r,mean=x[9:10]))

for(j in 1:2){
    for(i in 1:1000){
       scores[i,j] = findInterval(latent[i,j],thresh[j,])
}
}
scores = cbind(c(rep(0,nf),rep(1,nm)),scores)
colnames(scores)=c('Sex','VA','SPC')
return(data.frame(scores))
}

estimate.pf = function (first=T)
{

post=vector()
if(first==T){ score=patterns(simulate(25,975))}
else{score=patterns(simulate2(680,320))}

freqs=score[,4]
N.pats=NROW(score)
probs=matrix(NA,nr=N.pats,nc=2)
score[,-4]

lower=rep(NA,2)
upper=rep(NA,2)
thresh=matrix(NA,nr=2,nc=6)

x=thetas
thresh[1,]=c(-Inf,x[1:4],Inf)
thresh[2,]=c(-Inf,x[5:8],Inf)
mus=x[9:10]

R=matrix(c(1,rep(x[11],2),1),nc=2)

for(i in 1:N.pats){
probs[i,1]=as.numeric(pmvnorm(mean=c(0,0),
             lower=c(thresh[1,score[i,2]],thresh[2,score[i,3]]),
             upper=c(thresh[1,score[i,2]+1],thresh[2,score[i,3]+1]),
             corr=R,algorithm=Miwa)[1])
probs[i,2]=as.numeric(pmvnorm(mean=mus,
             lower=c(thresh[1,score[i,2]],thresh[2,score[i,3]]),
             upper=c(thresh[1,score[i,2]+1],thresh[2,score[i,3]+1]),
             corr=R,algorithm=Miwa)[1])
}

LnLK = function(x){
zlnlk=0
for(i in 1:N.pats) {
pf = x*probs[i,1]
pm = (1-x)*probs[i,2]
zlnlk=zlnlk+freqs[i]*log(pf+pm)
}
return(zlnlk)
}

mle=optimize(LnLK,lower=0.0001,upper=0.9999,maximum=T)$max

SE = function(x){
Hess=0
for(i in 1:N.pats) {
pf = probs[i,1]
pm = probs[i,2]
Hess = Hess + freqs[i]*((pf-pm)^2/(x*pf+(1-x)*pm)^2)
}
return(sqrt(1/Hess))
}

se=SE(mle)

cat(paste('\nEstimated proportion females = ',(round(mle,4)),'\n',sep=''))
cat(paste('Standard error = ',round(se,4),'\n',sep=''))
cat(paste('95% Confidence interval (SE) = ',round(mle-1.96*se,4),' to ',round(mle+1.96*se,4),'\n',sep=''))

find.CL = function(x){
   return(LnLK(mle)-LnLK(x)-1.920729)
}

low=uniroot(find.CL,lower=0.000001,upper=mle)$root
hi=uniroot(find.CL,lower=mle,upper=1)$root

cat(paste('95% Confidence interval (chi-square) = ',round(low,4),' to ',round(hi,4),'\n',sep=''))

LnLK2 = function(x,max.lnlk=0){
zlnlk=0
for(i in 1:N.pats) {
pf = x*probs[i,1]
pm = (1-x)*probs[i,2]
zlnlk=zlnlk+freqs[i]*log(pf+pm)
}
return(zlnlk-max.lnlk)
}

mle=optimize(LnLK2,lower=0.0,upper=1,maximum=T)
max.lnlk=mle$objective
mle=mle$max

pdf=function(x,denom=1){
   return(exp(LnLK2(x,max.lnlk=max.lnlk))/denom)
}
denom = integrate(pdf,lower=0.0,upper=1)$val

cut.pt = function(x,height=30){
   pdf(x,denom=denom)-height
}

find.95 = function(height){

left = uniroot(cut.pt,lower=0.0,upper=mle,height=height)$root
right = uniroot(cut.pt,lower=mle,upper=1,height=height)$root
left.area = integrate(pdf,lower=0.0,upper=left,denom=denom)$val
right.area = integrate(pdf,lower=right,upper=1,denom=denom)$val
tails = left.area + right.area
return(tails-0.05)
}

mle.height=pdf(mle,denom=denom)

height = uniroot(find.95,lower=mle.height/100000,upper=mle.height*0.999999)$root
left = uniroot(cut.pt,lower=0.0,upper=mle,height=height)$root
right = uniroot(cut.pt,lower=mle,upper=1,height=height)$root

cat(paste('Left and right 95% HPD bounds = ',round(left,4),' ',round(right,4),'\n',sep=''))

for(i in 1:N.pats) {

pf = mle*probs[i,1]
pm = (1-mle)*probs[i,2]

post[i]=pf/(pf+pm)
}

female=rep(0,N.pats)
male=rep(0,N.pats)

for(i in 1:N.pats){
if(post[i]>0.5) female[i]=freqs[i]
else {male[i]=freqs[i]}
}

out = cbind(score,post,female,male)

sexing=matrix(0,nc=2,nr=2)
for(i in 1:N.pats){
if(out[i,1]==0){
    sexing[1,1]=sexing[1,1]+out[i,6]
    sexing[1,2]=sexing[1,2]+out[i,7]
}
else{
    sexing[2,2]=sexing[2,2]+out[i,7]
    sexing[2,1]=sexing[2,1]+out[i,6]
}
}
rownames(sexing)=c('Female','Male')
colnames(sexing)=c('"Female"','"Male"')

cat('\n')
print(sexing)

cat(paste('\nPercent correctly sexed = ',round(sum(diag(sexing))/10,2),'%\n',sep=''))

}

estimate.pf()

Applying the Training Data to a Simulated Testing Dataset using Logistic Regression

test.logist = function ()
{
training=Terry[,3:5]
testing=simulate(25,975)
Sex=c('Female','Male')[testing$Sex+1]
Sexed.As=c('Female','Male')[(as.numeric(predict(glm(Sex~.,data=training),newdata=testing))>0.5)+1]

out=table(data.frame(Sex,Sexed.As))

print(out)
perc=sum(diag(out))/10
cat(paste('\nPercent correct = ',round(perc,2),'\n\n'))
cat('Estimated proportions of the sexes\n')
print(apply(out,2,sum)/1000)
}

test.logist()

Simulating a training set for two traits with less sexual dimorphism than the Phenice characteristics

simulate2 = function (nf=680,nm=320)
{
library(mvtnorm)
thresh=c(-Inf,0,.25,.35,.6,Inf)
thresh=rbind(thresh,thresh)

r = matrix(c(1,rep(.5,2),1),nc=2)

set.seed(42)

scores=matrix(NA,nc=2,nr=1000)
latent = rbind(rmvnorm(nf,sigma=r,mean=c(0,0)),
               rmvnorm(nm,sigma=r,mean=c(.6,.6)))

for(j in 1:2){
    for(i in 1:1000){
       scores[i,j] = findInterval(latent[i,j],thresh[j,])
}
}
scores = cbind(c(rep(0,nf),rep(1,nm)),scores)
colnames(scores)=c('Sex','Trait.1','Trait.2')
return(data.frame(scores))
}

test.logist2 = function ()
{
training=simulate2(500,500)
testing=simulate2(680,320)
Sex=c('Female','Male')[testing$Sex+1]
Sexed.As=c('Female','Male')[(as.numeric(predict(glm(Sex~.,data=training),newdata=testing))>0.5)+1]

out=table(data.frame(Sex,Sexed.As))

print(out)
perc=sum(diag(out))/10
cat(paste('\nPercent correct = ',round(perc,2),'\n\n'))
cat('Estimated proportions of the sexes\n')
print(apply(out,2,sum)/1000)
}

test.logist2()

mv.probit(simulate2(500,500))

estimate.pf(F)

Trivariate multivariate probit for the Phenice characteristics

mv.probit(Terry[,-(1:2)])