R and SAS code for Math 485 #R code HW2: 7 mean(sample(1:100,5)) mean(sample(1:100,5)) mean(sample(1:100,5)) mean(sample(1:100,5)) mean(sample(1:100,5)) mean(sample(1:100,45)) mean(sample(1:100,45)) mean(sample(1:100,45)) mean(sample(1:100,45)) mean(sample(1:100,45)) #R code HW3: 5) source("http://lagrange.math.siu.edu/Olive/sipack.txt") accisimf(n=50,p=0.01) accisimf(n=100,p=0.01) accisimf(n=150,p=0.01) accisimf(n=200,p=0.01) accisimf(n=250,p=0.01) accisimf(n=300,p=0.01) accisimf(n=350,p=0.01) accisimf(n=400,p=0.01) accisimf(n=450,p=0.01) *The program below is for Math 485 HW4.1; options ls = 70; data aspirin; * Agresti (1996, p. 268); input group mi count @@; cards; 1 1 189 1 2 10845 2 1 104 2 2 10933 ; proc print; proc freq; weight count; tables group*mi/chisq expected measures; run; *The program below is for Math 485 HW5.2; options ls = 70; data crabs; * Agresti (1996, p. 272); input width cases satell; cards; 22.69 14 5 23.84 14 4 24.77 28 17 25.84 39 21 26.79 22 15 27.74 24 20 28.67 18 15 30.41 14 14 ; proc logistic; model satell/cases = width; output out = predict p = pi_hat; proc print data = predict run; #HW6.5 R #a) copy Gladstone data into R source("http://lagrange.math.siu.edu/Olive/regdata.txt") y <- cbrainx[,10] x1 <- cbrainx[,1] x2 <- cbrainx[,3] x3 <- cbrainx[,4] x4 <- cbrainx[,5] x5 <- cbrainx[,6] x6 <- cbrainx[,7] x7 <- cbrainx[,8] x8 <- cbrainx[,9] x9 <- cbrainx[,11] outf <- glm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9,family=binomial) summary(outf) drop1(outf,test="Chi") #b) Test whether LR model is better than the null model #could get G^2_o and G^2_full from anova(outf) outn<-glm(y~1,family=binomial) anova(outn,outf,test="Chi") #c) reduced model with breadth, cause, height, size outr<-glm(y~x2+x3+x7+x9,family=binomial) anova(outr,outf,test="Chi") #d) ESP <- predict(outf) ESPR <- predict(outr) plot(ESPR,ESP) abline(0,1) #HW7.2 R code Copy Gladstone data into R with the source command. #HW7.2 a) source("http://lagrange.math.siu.edu/Olive/regdata.txt") y <- cbrainx[,10] age <- cbrainx[,1] lnag <- log(age) brd <- cbrainx[,3] fcau <- as.factor(cbrainx[,4]) ceph <- cbrainx[,5] circ <- cbrainx[,6] hdht <- cbrainx[,7] ht <- cbrainx[,8] len <- cbrainx[,9] siz <- cbrainx[,11] lnsz <- log(siz) outf <- glm(y~age+lnag+brd+fcau+ceph+circ+hdht+ht+len+siz+lnsz,family=binomial) summary(outf) #full model #HW 7.2b back <- step(outf) #backward elimination is the default outn <- glm(y~1,family=binomial) #null model with just a constant #HW 7.2 c) AIC(outn) AIC(outf) ##hit Enter #HW7.2 d) forward selection library(MASS) outfs<-stepAIC(outn,direction="forward",scope=list(upper=outf,lower=outn)) outfs #hit enter #outf = glm(y ~ x1+x2+x3+x4+x5+x6+x7+x8,family=binomial) #outn <- glm(y ~ 1,family=binomial) #forwards = step(outn,scope=list(lower=formula(outn),upper=formula(outf)), #direction="forward") #R code for R help handout #Poisson regression n <- 100 q <- 5 y <- 0 * 1:n beta <- 0 * 1:q beta[1:3] <- 1 alpha <- -2.5 x <- matrix(rnorm(n * q), nrow = n, ncol = q) x <- 0.5*x + 1 SP <- alpha + x%*%beta y <- rpois(n,lambda=exp(SP)) #fit the Poisson regression out <- glm(y ~ x[, 1] + x[, 2] + x[, 3] + x[, 4] + x[,5], family = poisson) out summary(out) #make a response plot ESP <- x %*% out$coef[-1] + out$coef[1] Y <- y plot(ESP,Y) lines(lowess(ESP,Y)) llrplot(x,y) #GAM analog of Poisson regression library(mgcv) x1 <- x[,1] x2 <- x[,2] x3 <- x[,3] out <- gam(y ~ s(x1) + s(x2) + s(x3),family = poisson) out summary(out) plot(out) #make a response plot EAP<-predict.gam(out) plot(EAP,Y) lines(lowess(EAP,Y)) #make an EE plot plot(EAP,ESP) abline(0,1) #binary logistic regression mv <- 1 + 0 * y pv <- exp(SP)/(1 + exp(SP)) for(i in 1:n) y[i] <- rbinom(1, size = 1, prob = pv[i]) out <- glm(y ~ x[, 1] + x[, 2] + x[, 3] + x[, 4] + x[, 5], family = binomial, weights = mv) out summary(out) lressp(x,y) #GAM analog of binary logistic regression x1 <- x[,1] x2 <- x[,2] x3 <- x[,3] out <- gam(y ~ s(x1) + s(x2) + s(x3),family = binomial, weights = mv) out summary(out) plot(out) #loading data into R cbrn <- matrix(scan(),nrow=267,ncol=13,byrow=T) dim(cbrn) y <- cbrn[,12] out <- glm(y~cbrn[,3]+cbrn[,8],family=binomial) summary(out) *The program below is for Math 485 HW10.2 = Math 583 Applied Robust Statistics 13.21; options ls = 70; data cbrain; infile "J:cbrain2.txt"; input obs x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 sex x12; proc rsquare; model sex = x2 x3 x4 x5 x6 x7 x9 x10 / cp; proc reg; model sex = x2 x3 x4 x5 x6 x7 x9 x10; output out = a p = pred r = resid; run; *This program is for HW 11.1 and comes from Agresti's webpage. data cmh; input center smoke cancer count @@; cards; 1 1 1 126 1 1 2 100 1 2 1 35 1 2 2 61 2 1 1 908 2 1 2 688 2 2 1 497 2 2 2 807 3 1 1 913 3 1 2 747 3 2 1 336 3 2 2 598 4 1 1 235 4 1 2 172 4 2 1 58 4 2 2 121 5 1 1 402 5 1 2 308 5 2 1 121 5 2 2 215 6 1 1 182 6 1 2 156 6 2 1 72 6 2 2 98 7 1 1 60 7 1 2 99 7 2 1 11 7 2 2 43 8 1 1 104 8 1 2 89 8 2 1 21 8 2 2 36 ; proc freq; weight count; tables center*smoke*cancer / cmh; run;