## Settings for RMarkdown http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE,
echo = TRUE, fig.width = 10, fig.height = 8)
options(width = 116, scipen = 10)
setwd("~/statistics/bio509/")
Data loading
## Load gamair package
library(gamair)
## Load hubble dataset
data(hubble)
## Check data
str(hubble)
'data.frame': 24 obs. of 3 variables:
$ Galaxy: Factor w/ 24 levels "IC4182","NGC0300",..: 2 3 4 5 6 8 9 7 10 11 ...
$ y : int 133 664 1794 1594 1473 278 714 882 80 772 ...
$ x : num 2 9.16 16.14 17.95 21.88 ...
head(hubble)
Galaxy y x
1 NGC0300 133 2.00
2 NGC0925 664 9.16
3 NGC1326A 1794 16.14
4 NGC1365 1594 17.95
5 NGC1425 1473 21.88
6 NGC2403 278 3.22
Plotting points
## Plot
plot(y ~ x, data = hubble, xlim = c(0, 25), ylim = c(0, 2000),
ylab = "galaxy's relative velocity in KM/sec",
xlab = "Galaxy's distance in Mega parsecs")
Age calculation
## Age of universe in
kmTOmegaparsec <- 3.09*10^(19)
fit.task2 <- lm(I(y / kmTOmegaparsec) ~ -1 + x, data=hubble)
summary(fit.task2)
Call:
lm(formula = I(y/kmTOmegaparsec) ~ -1 + x, data = hubble)
Residuals:
Min 1Q Median 3Q Max
-2.38e-17 -4.29e-18 -6.15e-19 5.57e-18 1.81e-17
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 2.48e-18 1.28e-19 19.3 1e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.38e-18 on 23 degrees of freedom
Multiple R-squared: 0.942, Adjusted R-squared: 0.939
F-statistic: 373 on 1 and 23 DF, p-value: 1.03e-15
(1 / coef(fit.task2)) / 3600*24*365
x
981834020134355200
Plotting fitted line
## Graphing
plot(y ~ x, data = hubble, xlim = c(0, 25), ylim = c(0, 2000),
ylab = "galaxy's relative velocity in KM/sec",
xlab = "Galaxy's distance in Mega parsecs")
fit.task2.raw <- lm(y ~ -1 + x, data = hubble)
abline(reg = fit.task2.raw)
##
fit.task3 <- lm(log(y) ~ 1 + offset(log(x)), data=hubble)
summary(fit.task3)
Call:
lm(formula = log(y) ~ 1 + offset(log(x)), data = hubble)
Residuals:
Min 1Q Median 3Q Max
-1.1786 -0.0830 0.0236 0.2262 0.4395
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2714 0.0704 60.7 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.345 on 23 degrees of freedom
## plot
plot(log(y) ~ log(x), data = hubble, xlim = c(0,4), ylim = c(0, 8),
ylab = "galaxy's relative velocity in KM/sec (log)",
xlab = "Galaxy's distance in Mega parsecs (log)")
## abline intercept and beta = 1
abline(a = coef(fit.task3), b = 1)
