Full_Name: Richard Cai
Version: R 2.3.1
OS: Windows XP SP2
Submission from: (NULL) (220.233.184.74)
My R scrpit code was interrupted when I was trying to load it into R's built-in
script editor.
#################### R code for question 4 ########################
dat <- read.table("Drosophila.txt",header=T)
library(lattice)
trellis.device(color=F)
attach(dat)
########## scatter plot of the data ################
data_plot <- xyplot(Weight~density)
print(data_plot)
dev.copy2eps(file="D:/CaiYu/UNE/Statistics/Stat
356/Assignment1/q4_data_plot.eps",horizontal=F)
######### plot of group means against density ##########
group_mean <- sapply(split(Weight,density),mean)
density_level <- as.numeric(levels(as.factor(density)))
group_plot <- xyplot(group_mean ~ density_level)
windows()
print(group_plot)
dev.copy2eps(file="D:/CaiYu/UNE/Statistics/Stat
356/Assignment1/q4_group_plot.eps",horizontal=F)
######### fitted model #######################
fitted_model <- lm(Weight ~density, data=dat)
print(anova(fitted_model))
print(summary(fitted_model)$coefficients)
residual_SS <- anova(fitted_model)$Sum[2]
df_residual_SS <- anova(fitted_model)$Df[2]
########### One-way analysis of variance #########
aov <- aov(Weight ~ as.factor(density), data=dat)
print(anova(fitted_model))
print(summary(aov))
pure_error_SS <- anova(aov)$Sum[2]
df_pure_error_SS <- anova(aov)$Df[2]
########## lack of fittness ###############
lack_of_fittness_SS <- residual_SS - pure_error_SS
df_lack_of_fittness_SS <- df_residual_SS - df_pure_error_SS
F_lack_of_fittness <-
(lack_of_fittness_SS/df_lack_of_fittness_SS)/(pure_error_SS/df_pure_error_SS)
P_value_lack_of_fittness <- pf(F_lack_of_fittness, df_lack_of_fittness_SS,
df_pure_error_SS, lower.tail=F)
print(P_value_lack_of_fittness)
########## testing hypothesis beta1=0 ###############
print(summary(fitted_model)$coefficients)
########### prediction for density=15 #############
new1 <- data.frame(density=15)
print(predict(lm(Weight~density), new1, se.fit = TRUE))
print(predict(lm(Weight~density), new1, interval="prediction"))
print(predict(lm(Weight~density), new1, interval="confidence"))
########### prediction for density=40 #############
new2 <- data.frame(density=40)
print(predict(lm(Weight~density), new2, se.fit = TRUE))
print(predict(lm(Weight~density), new2, interval="prediction"))
print(predict(lm(Weight~density), new2, interval="confidence"))
############# prediciton plot #################
windows()
pred.w.plim <- predict(fitted_model, data.frame(density=seq(1,40,1)),
interval="prediction")
pred.w.clim <- predict(fitted_model, data.frame(density=seq(1,40,1)),
interval="confidence")
larval_density <- seq(1,40,1)
matplot(larval_density,cbind(pred.w.clim, pred.w.plim[,-1]),
lty=c(1,2,2,3,3), col=rep(1,5), type="l", ylab="predicted Weight")
dev.copy2eps(file="D:/CaiYu/UNE/Statistics/Stat
356/Assignment1/q4_prediction_plot.eps",horizontal=F)
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
