Raj at 05:18, 28 October 2016

2016-10-28T05:18:44Z

← Older revision		Revision as of 05:18, 28 October 2016
Line 1:		Line 1:
	[[:Category:R]]		[[Category:R]]

	<source line>		<source line>

Raj at 09:30, 6 August 2013

2013-08-06T09:30:15Z

← Older revision		Revision as of 09:30, 6 August 2013
Line 1:		Line 1:
	[[Category:R]]		[[:Category:R]]

	<source line>		<source line>

Raj at 09:26, 6 August 2013

2013-08-06T09:26:32Z

← Older revision		Revision as of 09:26, 6 August 2013
Line 1:		Line 1:
	[[:Category:R]]		[[Category:R]]

	<source line>		<source line>

Raj: Created page with ":Category:R # oneway.r # code for 1-way ANOVA examples in R diet <- read.csv('diet.csv', as.is=T) diet <- diet[diet$diet !='lopro',] # I use as.is=T to f..."

2013-08-06T09:22:38Z

Created page with "Category:R <source line> # oneway.r # code for 1-way ANOVA examples in R diet <- read.csv('diet.csv', as.is=T) diet <- diet[diet$diet !='lopro',] # I use as.is=T to f..."

New page

[[:Category:R]]

<source line>

# oneway.r
# code for 1-way ANOVA examples in R

diet <- read.csv('diet.csv', as.is=T)
diet <- diet[diet$diet !='lopro',]

# I use as.is=T to force me to explicitly declare
# factors. default conversion does not convert
# factors with numeric levels.

diet$diet.f <- factor(diet$diet)

# anova using lm:
diet.lm <- lm(longevity ~ diet.f, data=diet)

# lm() has lots of helper / reporting functions:
coef(diet.lm)
# coefficients
vcov(diet.lm)
# and their variance-covariance matrix
sqrt(diag(vcov(diet.lm)))
# se's of the coefficients
anova(diet.lm)
# ANOVA table using type I = sequential SS
summary(diet.lm)
# lots of information
plot(predict(diet.lm), resid(diet.lm))
# plot of residuals vs predicted values

# model comparison by hand:
diet.m0 <- lm(longevity ~ +1, data=diet)
# intercept only model
anova(diet.m0, diet.lm)
# change in SS from 1st to 2nd model

drop1(diet.lm)
# drop one term at a time = type II SS

# orthogonal contrasts
diet.means <- tapply(diet$longevity, diet$diet, mean)
# mean for each diet
diet.helm <- contr.helmert(5)
# matrix of Helmert coeff. for 5 groups
diet.c1 <- t(diet.helm) %*% diet.means
# estimate of each contrast
diet.ss1 <- 49*diet.c1^2/apply(diet.helm^2,2,sum)
# SS for each contrast
sum(diet.ss1)

# second set of contrasts were hand entered
# using diet.2nd <- rbind(c(-2,-1,0,1,2), ...)

# third set:
diet.trt <- -contr.treatment(5)
diet.trt[1,] <- 1

diet.c3 <- t(diet.trt) %*% diet.means
# estimate of each contrast
diet.ss3 <- 49*diet.c3^2/apply(diet.trt^2,2,sum)
# SS for each contrast
sum(diet.ss3)

# C beta tests
diet.c1 <- t(diet.helm) %*% diet.means
# estimate of contrast

X1 <- model.matrix(~ -1 + diet.f, data=diet)
# X matrix using cell means parameterization
diet.cm1 <- t(diet.helm) %*% solve(t(X1) %*% X1) %*%
diet.helm

t(diet.c1) %*% solve(diet.cm1) %*% diet.c1

# how do the explicit contrasts of cell means
# compare to R using those contrasts?

contrasts(diet$diet.f) <- contr.helmert
# tell R to use helmert contrasts
# default is contr.treatment,
# which is drop first level of the factor
# contr.SAS is SAS-like (drop last level)

diet.lmh <- lm(longevity ~ diet.f,data=diet)
coef(diet.lmh)
# coefficients are not the same as the
# hand-computed contrasts!

# the R coefficients are related to the
# contrasts among cell means

apply(diet.helm^2,2,sum)
# sum of squared coefficients

coef(diet.lmh)[2:5] * c(2,6,12,20)
# these are the same as diet.c1

# although R calls them contrasts, they really are not.
# they are columns of the X matrix for a regression.
model.matrix(diet.lmh)

# can do C beta from the lm() information

diet.mse <- 4687.7/240 # MSE = est of sigma^2
diet.clm <- coef(diet.lm)[2:5]
# extract coefficients for diet factor
diet.vclm <- vcov(diet.lm)[2:5,2:5]
# and their VC matrix
diet.mse * t(diet.clm) %*% solve(diet.vclm) %*%
diet.clm
# SS for diet
t(diet.clm) %*% solve(diet.vclm) %*% diet.clm/5
# F statistic for diet (5 is # rows in C)

</source>

Anova-Oneway - Revision history

Raj at 05:18, 28 October 2016

Raj at 09:30, 6 August 2013

Raj at 09:26, 6 August 2013

Raj: Created page with ":Category:R # oneway.r # code for 1-way ANOVA examples in R diet <- read.csv('diet.csv', as.is=T) diet <- diet[diet$diet !='lopro',] # I use as.is=T to f..."