pdf("g11.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class553")
setwd("~/Taught1/960-553/Data")
<pre>
#<a name=block1> Block 1</a>
#**********************************************************/
# Stuart (1953 Biometrika 40, 105-110), quoted by Stuart */
# (1955 Biometrika 42, 412-416), and Agresti (1983 Bio-  */
# metrics 39(2) 505-510), present an ordinal categoriza- */
# tion of right and left eye vision of Royal Ordinance   */
# female factory workers 1943-46.                        */
#*********************************************************/
eyes<-read.table("eyes.dat",col.names=c("r","l","count"))
#*****************************************************/
# Cohen's Kappa                                      */
#*****************************************************/
library(psych)#For cohen.kappa; Calculates Cohen's kappa to measure agreement.
cohen.kappa(xtabs(count~r+l,data=eyes))
library(TheorStat)# For fun.testpoly and fun.testplotpoly; Explores the relationship between multivariate normality and the polychoric correlation via simulation and plots the results.
out<-fun.testpoly(nmc=1000)
fun.plottestpoly(out)
#<a name=block2> Block 2</a>
#**********************************************************/
#Cronbach's alpha                                         */
#**********************************************************/
#*********************************************************/
# Data from SAS Usere's Guide, documentation for CORR pr-*/
# ocedure.  Data are from 35 fish of a certain species   */
# from a lake in Finland.  Variables are weight, length, */
# max height as a percentage of Length3, and max width as*/
# a percentage of Length3.                               */
#*********************************************************/
fish <- read.table("fish.dat",na.strings=".",
   col.names=c("Weight","Length3","HtPct","WidthPct"))
#How reliable are the four as measures of a fish's linear 
#size?  
#Rescale weight to put on dimension scale.
fish$weightrt<-fish$Weight^(1/3)
# Recover height from percentages.
fish$Height<-fish$HtPct*fish$Length3*0.01
fish$Width<-fish$WidthPct*fish$Length3*0.01
library(psych)#For alpha; Calculate Cronbach's alpha to measure reliability.
alpha(fish[,c("Length3","weightrt","Height","Width")])
xx<-fish[-14,c(2,5:7)]
yy<-xx
for(j in 1:4) yy[[j]]<-(yy[[j]]-mean(yy[[j]]))/sd(yy[[j]])
alpha(xx)
alpha(yy)
#<a name=block3> Block 3</a>
#************************************************************/
#              Cow Mastitis Example                         */
# Data from Andrews and Herzberg (1985), Data: A Collection */
# of Problems from Many Fields for the Student and Research */
# Worker, Table 61.  Observations represent infections in   */
# udders of cows given a placebo or one of eight anti-      */
# biotics.  First two variables identify the table.  Third  */
# and fourth identify cow, fifth identifies herdd, sixth    */
# treatment, and next eight represent level of infection at */
# each of eight sites.  Treatment assignment is non-random. */
# http://lib.stat.cmu.edu/datasets/Andrews/T61.1            */
#************************************************************/
cows<-as.data.frame(scan("T61.1",flush=T,what=list(ex=0,   
  ta=0,c1=0,c2=0,herd=0,treaetment=0,q11=0,q12=0,q21=0)))
#*************************************************************/
# Explore relation between the first four treatments and     */
# infection severity in udder 2 location 1 in two ways: both */
# variables treated as ordered, and both as unordered.       */
#*************************************************************/
cat('\n Exact Tests for Assoc. betw. treatment and infection  \n')
#*************************************************************/
#Do a "toy" example of logistic regression.                  */
#*************************************************************/
#If we knew the intercept term, we would use it as an offset.*/
#The p-value depends on the value of this offset.            */
#Investigate actual size of unconditional test.              */
#*************************************************************/
cat('\n asy  \n')
glm(c(1,1,0,0)~c(-1,1,-1,1),weight=c(7,3,3,7),family=binomial)
library(TheorStat)# For waldsize; Calculate true size of Wald statistic for the comparison of two proportions.
out<-waldsize(10,10,201)
plot(out[[1]],out[[2]],xlab="Common Probability",
  ylab="Size of nominal .05 test",type="l",
  main="Size of Wald Test from Logistic Regression")
#<a name=block4> Block 4</a>
#*************************************************************/
#*  Prostate data                                            */
#* From Andrews and Herzberg (1985), Data: A Collection of   */
#* Problems from Many Fields for the Student and Research    */
#* Worker, Table 46.  Observations represent subjects in a   */
#* prostate cohort study, randomized to one of four dose     */
#* levels of diethylstilbestrol.  Rx records dose in four    */
#* ordered categories, with 1 being placebo.  Disease stage  */
#* is 3 or 4.  monfol is months of followup.  surv is 0 if   */
#* alive after 50 mo, and codes cause of death otherwise.    */
#* http://lib.stat.cmu.edu/datasets/Andrews/T46.1            */
#* The on-line version of the data set adds 3 fields before  */
#* the first field in the book.  Variables of interest are   */
#* stage, rx, monfol, and surv in fields 5, 6, 10, 11 of the */
#* online version, resp.  Causes of death are given by var-  */
#* ious positive integers in surv; I recode these to 1.  The */
#* data file has more columns than we are interested in.  Put*/
#* place-holding variables in for variables not of interest  */
#* between variables of interest.  Data were previously pub- */
#* lished by Byar and Corle (1977 Chronic Disease) and Byar  */
#* and Green (1980 Bull. Cancer).  Lower value of dichoto-   */
#* mized dose begins with blank to make it alphabetized      */
#* before high.                                              */
#*************************************************************/
# If we omit the ones past where the data of interest stops
# in the "what" list, and use flush=T, R will ignore them. 
prostate<-read.table("T46.1")[,c(5,6,10,11)]
dimnames(prostate)[[2]]<-c("stage","rx","monfol","surv")
prostate$alive<-(prostate$surv==0)+0                      
prostate$dose<-c(" low","high")[1+(prostate$rx>1)]         
stage4<-prostate[prostate$stage==4,]
#*************************************************************/
#* Causes of death are given by various positive integers in */
#* surv; I recode these to 1.                                */
#*************************************************************/
prostate$alive<-(prostate$surv==0)+0
#*******************************************************/
#* Test whether survival is the same across treatments.*/
#*******************************************************/
chisq.test(stage4$rx,stage4$alive,simulate.p.value=TRUE)
# One can score tables by inverse probability.
fisher.test(stage4$rx,stage4$alive)
</pre>