[Rd] R's built-in Script Editor Bug (PR#9137)

[yu_cai99 at yahoo.com]

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)