# install.packages("psych") # install.packages("QuantPsyc") # Grave = Alt + 126 ex1 <- read.table("lecture2ex1.csv", sep=",", header=TRUE) library(psych) describe(ex1) corr.test(ex1) out1 <- lm(stress ~ homepress + workpress, data=ex1) summary(out1) anova(out1) library(QuantPsyc) lm.beta(out1) # out1s <- lm(scale(stress) ~ scale(homepress) + scale(workpress), data=ex1) # summary(out1s) ex2 <- read.table("lecture2ex2.csv", sep=",", header=TRUE) describe(ex2) corr.test(ex2) out2 <- lm(cute ~ face + shape + habit + close, data=ex2) summary(out2) lm.beta(out2) anova(out2) out2a <- lm(cute ~ I(face - 50) + I(shape - mean(shape)) + I(habit - 50) + I(close - 50), data=ex2) summary(out2a) lm.beta(out2a) anova(out2a) ex3 <- read.table("lecture2ex3.csv", sep=",", header=TRUE) head(ex3, 15) ex3[,"sex"] <- factor(ex3[,"sex"], labels=c("female", "male")) ex3[,"region"] <- factor(ex3[,"region"], labels=c("north", "center", "northeast", "south")) ex3[,"region"] <- relevel(ex3[,"region"], ref="south") head(ex3, 15) describe(ex3) with(ex3, table(sex)) with(ex3, table(region)) corr.test(ex3[,c("chilimg", "age")]) out3 <- lm(chilimg ~ I(age - 40) + sex + region, data=ex3) summary(out3) lm.beta(out3) anova(out3) out3noregion <- lm(chilimg ~ I(age - 40) + sex, data=ex3) summary(out3noregion) anova(out3noregion, out3) out3nosex <- lm(chilimg ~ I(age - 40) + region, data=ex3) anova(out3nosex, out3) out3noage <- lm(chilimg ~ sex + region, data=ex3) anova(out3noage, out3) ex3expanded <- model.matrix(~ chilimg + age + sex + region, data=ex3) head(ex3expanded) ex3expanded <- data.frame(ex3expanded) out3b <- lm(chilimg ~ I(age - 40) + sexmale + regionnorth + regioncenter + regionnortheast, data=ex3expanded) summary(out3b) lm.beta(out3b) ############################### # install.packages("Matching") library(Matching) set.seed(123321) n <- 200 nrep <- 1000 result <- matrix(NA, nrep, 3) for(i in seq_len(nrep)) { age <- runif(n, 18, 60) male <- runif(n, 0, 1) > 0.5 probsmoke <- 0.35 - age/100 + 0.2*male randnum <- runif(n, 0, 1) issmoke <- probsmoke > randnum healthscore <- 20 + 5*issmoke - 0.1*scale(age) + male + rnorm(n, 0, sqrt(30)) dat <- data.frame(age, issmoke, healthscore) out1 <- lm(healthscore ~ age + issmoke + male, data=dat) result[i, 1] <- out1$coef[3] glm1 <- glm(issmoke ~ age + male, family=binomial, data=dat) X <- glm1$fitted Y <- healthscore Tr <- issmoke out2 <- Match(Y=Y, Tr=Tr, X=X, M=1); result[i, 2] <- out2$est result[i, 3] <- mean(issmoke) } apply(result, 2, mean) apply(result, 2, sd) t.test(age ~ issmoke, dat) table(male, issmoke)