R Under development (unstable) (2024-09-14 r87148 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(ISwR) > ps.options(height=3.5, width=4.4, pointsize=8, horiz=F) > par(mar=c(4,4,3,2)+.1) > set.seed(310367) > x <- y <- c(7, 9, NA, NA, 13) > all(is.na(x) == is.na(y)) & all((x == y)[!is.na(x)]) [1] TRUE > x <- factor(c("Huey", "Dewey", "Louie", "Huey")) > y <- c("blue", "red", "green") > x [1] Huey Dewey Louie Huey Levels: Dewey Huey Louie > y[x] [1] "red" "blue" "green" "red" > juul.girl <- juul[juul$age >=7 & juul$age < 14 & juul$sex == 2,] > summary(juul.girl) age menarche sex igf1 tanner Min. : 7.000 Min. :1.000 Min. :2 Min. : 95.0 Min. :1.000 1st Qu.: 8.885 1st Qu.:1.000 1st Qu.:2 1st Qu.:214.2 1st Qu.:1.000 Median :10.560 Median :1.000 Median :2 Median :302.5 Median :1.000 Mean :10.557 Mean :1.132 Mean :2 Mean :351.8 Mean :1.951 3rd Qu.:12.207 3rd Qu.:1.000 3rd Qu.:2 3rd Qu.:473.2 3rd Qu.:3.000 Max. :13.940 Max. :2.000 Max. :2 Max. :914.0 Max. :5.000 NA's :5 NA's :11 NA's :5 NA's :119 NA's :66 testvol Min. : NA 1st Qu.: NA Median : NA Mean :NaN 3rd Qu.: NA Max. : NA NA's :351 > sapply(1:10, function(i) mean(rexp(20))) [1] 0.7969888 0.8648102 0.7280609 1.3195396 1.1868598 1.2246834 0.8794734 [8] 0.9722758 1.0604842 0.8782972 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > x <- 1:10 > z <- append(x, 1.23, after=7) > z [1] 1.00 2.00 3.00 4.00 5.00 6.00 7.00 1.23 8.00 9.00 10.00 > z <- c(x[1:7],1.23,x[8:10]) > z [1] 1.00 2.00 3.00 4.00 5.00 6.00 7.00 1.23 8.00 9.00 10.00 > v <- 1.23; k <- 7 > i <- seq(along=x) > z <- c(x[i <= k], v, x[i > k]) > z [1] 1.00 2.00 3.00 4.00 5.00 6.00 7.00 1.23 8.00 9.00 10.00 > write.table(thuesen, file="foo.txt") > # edit the file > read.table("foo.txt", na.strings=".") blood.glucose short.velocity 1 15.3 1.76 2 10.8 1.34 3 8.1 1.27 4 19.5 1.47 5 7.2 1.27 6 5.3 1.49 7 9.3 1.31 8 11.1 1.09 9 7.5 1.18 10 12.2 1.22 11 6.7 1.25 12 5.2 1.19 13 19.0 1.95 14 15.1 1.28 15 6.7 1.52 16 8.6 NA 17 4.2 1.12 18 10.3 1.37 19 12.5 1.19 20 16.1 1.05 21 13.3 1.32 22 4.9 1.03 23 8.8 1.12 24 9.5 1.7 > write.table(thuesen, file="foo.txt", na=".") > read.table("foo.txt", na.strings=".") blood.glucose short.velocity 1 15.3 1.76 2 10.8 1.34 3 8.1 1.27 4 19.5 1.47 5 7.2 1.27 6 5.3 1.49 7 9.3 1.31 8 11.1 1.09 9 7.5 1.18 10 12.2 1.22 11 6.7 1.25 12 5.2 1.19 13 19.0 1.95 14 15.1 1.28 15 6.7 1.52 16 8.6 NA 17 4.2 1.12 18 10.3 1.37 19 12.5 1.19 20 16.1 1.05 21 13.3 1.32 22 4.9 1.03 23 8.8 1.12 24 9.5 1.70 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > 1 - pnorm(3) [1] 0.001349898 > 1 - pnorm(42, mean=35, sd=6) [1] 0.1216725 > dbinom(10, size=10, prob=0.8) [1] 0.1073742 > punif(0.9) # this one is obvious... [1] 0.9 > 1 - pchisq(6.5, df=2) [1] 0.03877421 > pnorm(-2) * 2 [1] 0.04550026 > qnorm(1-.01/2) [1] 2.575829 > qnorm(1-.005/2) [1] 2.807034 > qnorm(1-.001/2) [1] 3.290527 > qnorm(.25) [1] -0.6744898 > qnorm(.75) [1] 0.6744898 > rbinom(10, 1, .5) [1] 0 1 0 0 0 0 1 0 0 1 > ifelse(rbinom(10, 1, .5) == 1, "H", "T") [1] "T" "T" "T" "T" "H" "T" "T" "T" "H" "T" > c("H", "T")[1 + rbinom(10, 1, .5)] [1] "T" "H" "T" "T" "H" "T" "T" "H" "T" "T" > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > x <- 1:5 ; y <- rexp(5,1) ; opar <- par(mfrow=c(2,2)) > plot(x, y, pch=15) # filled square > plot(x, y, type="b", lty="dotted") > plot(x, y, type="b", lwd=3) > plot(x, y, type="o", col="blue") > par(opar) > plot(rnorm(10),type="o", pch=21, bg="white") > x1 <- rnorm(20) > x2 <- rnorm(10)+1 > q1 <- qqnorm(x1, plot.it=F) > q2 <- qqnorm(x2, plot.it=F) > xr <- range(q1$x, q2$x) > yr <- range(q1$y, q2$y) > qqnorm(x1, xlim=xr, ylim=yr) > points(q2, col="red") > hist(react) > library(MASS) > truehist(react,h=1,x0=.5) > z <- runif(5) > curve(quantile(z,x), from=0, to=1) > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > qqnorm(react) > t.test(react) One Sample t-test data: react t = -7.7512, df = 333, p-value = 1.115e-13 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: -0.9985214 -0.5942930 sample estimates: mean of x -0.7964072 > wilcox.test(react) Wilcoxon signed rank test with continuity correction data: react V = 9283.5, p-value = 2.075e-13 alternative hypothesis: true location is not equal to 0 > wilcox.test(vital.capacity~group, data=vitcap) Wilcoxon rank sum test with continuity correction data: vital.capacity by group W = 30.5, p-value = 0.01783 alternative hypothesis: true location shift is not equal to 0 Warning message: In wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...) : cannot compute exact p-value with ties > attach(intake) ; opar <- par(mfrow=c(2,2)) > plot(post ~ pre) ; abline(0,1) > plot((post+pre)/2, post - pre, + ylim=range(0,post-pre)); abline(h=0) > hist(post-pre) > qqnorm(post-pre) > detach(intake) > par(opar) > shapiro.test(react) Shapiro-Wilk normality test data: react W = 0.95701, p-value = 2.512e-08 > shapiro.test(react[-c(1,334)]) Shapiro-Wilk normality test data: react[-c(1, 334)] W = 0.96869, p-value = 1.377e-06 > qqnorm(react[-c(1,334)]) > attach(ashina) > t.test(vas.active, vas.plac, paired=TRUE) Paired t-test data: vas.active and vas.plac t = -3.2269, df = 15, p-value = 0.005644 alternative hypothesis: true mean difference is not equal to 0 95 percent confidence interval: -71.1946 -14.5554 sample estimates: mean difference -42.875 > t.test((vas.active-vas.plac)[grp==1], + (vas.plac-vas.active)[grp==2]) Welch Two Sample t-test data: (vas.active - vas.plac)[grp == 1] and (vas.plac - vas.active)[grp == 2] t = -3.2517, df = 13.97, p-value = 0.005807 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -130.56481 -26.76853 sample estimates: mean of x mean of y -53.50000 25.16667 > t.test(rnorm(25))$p.value #repeat 10x [1] 0.08718793 > t.test(rt(25,df=2))$p.value #repeat 10x [1] 0.8077125 > t.test(rexp(25), mu=1)$p.value #repeat 10x [1] 0.4056805 > x <- replicate(5000, t.test(rexp(25), mu=1)$p.value) > qqplot(sort(x),ppoints(5000),type="l",log="xy") > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > fit <- lm(metabolic.rate ~ body.weight, data=rmr) > summary(fit) Call: lm(formula = metabolic.rate ~ body.weight, data = rmr) Residuals: Min 1Q Median 3Q Max -245.74 -113.99 -32.05 104.96 484.81 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 811.2267 76.9755 10.539 2.29e-13 *** body.weight 7.0595 0.9776 7.221 7.03e-09 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 157.9 on 42 degrees of freedom Multiple R-squared: 0.5539, Adjusted R-squared: 0.5433 F-statistic: 52.15 on 1 and 42 DF, p-value: 7.025e-09 > 811.2267 + 7.0595 * 70 # , or: [1] 1305.392 > predict(fit, newdata=data.frame(body.weight=70)) 1 1305.394 > qt(.975,42) [1] 2.018082 > 7.0595 + c(-1,1) * 2.018 * 0.9776 # , or: [1] 5.086703 9.032297 > confint(fit) 2.5 % 97.5 % (Intercept) 655.883819 966.5695 body.weight 5.086656 9.0324 > summary(lm(log(ab)~age, data=malaria)) Call: lm(formula = log(ab) ~ age, data = malaria) Residuals: Min 1Q Median 3Q Max -4.0753 -1.0622 0.1181 1.1012 2.7335 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.83697 0.38021 10.092 <2e-16 *** age 0.10350 0.03954 2.618 0.0103 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.478 on 98 degrees of freedom Multiple R-squared: 0.06536, Adjusted R-squared: 0.05582 F-statistic: 6.853 on 1 and 98 DF, p-value: 0.01025 > plot(log(ab)~age, data=malaria) > rho <- .90 ; n <- 100 > x <- rnorm(n) > y <- rnorm(n, rho * x, sqrt(1 - rho^2)) > plot(x, y) > cor.test(x, y) Pearson's product-moment correlation data: x and y t = 17.673, df = 98, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.8158420 0.9124999 sample estimates: cor 0.8724529 > cor.test(x, y, method="spearman") Spearman's rank correlation rho data: x and y S = 23676, p-value < 2.2e-16 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.8579298 > cor.test(x, y, method="kendall") Kendall's rank correlation tau data: x and y z = 9.8457, p-value < 2.2e-16 alternative hypothesis: true tau is not equal to 0 sample estimates: tau 0.6678788 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > walk <- unlist(zelazo) # or c(..,recursive=TRUE) > group <- factor(rep(1:4,c(6,6,6,5)), labels=names(zelazo)) > summary(lm(walk ~ group)) Call: lm(formula = walk ~ group) Residuals: Min 1Q Median 3Q Max -2.7083 -0.8500 -0.3500 0.6375 3.6250 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.1250 0.6191 16.355 1.19e-12 *** grouppassive 1.2500 0.8755 1.428 0.1696 groupnone 1.5833 0.8755 1.809 0.0864 . groupctr.8w 2.2250 0.9182 2.423 0.0255 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.516 on 19 degrees of freedom Multiple R-squared: 0.2528, Adjusted R-squared: 0.1348 F-statistic: 2.142 on 3 and 19 DF, p-value: 0.1285 > t.test(zelazo$active,zelazo$ctr.8w) # first vs. last Welch Two Sample t-test data: zelazo$active and zelazo$ctr.8w t = -3.0449, df = 8.6632, p-value = 0.01453 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.8878619 -0.5621381 sample estimates: mean of x mean of y 10.125 12.350 > t.test(zelazo$active,unlist(zelazo[-1])) # first vs. rest Welch Two Sample t-test data: zelazo$active and unlist(zelazo[-1]) t = -2.386, df = 9.0865, p-value = 0.04058 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.22070952 -0.08811401 sample estimates: mean of x mean of y 10.12500 11.77941 > fit <- lm(volume~method+subject, data=lung) > anova(fit) Analysis of Variance Table Response: volume Df Sum Sq Mean Sq F value Pr(>F) method 2 1.08111 0.54056 6.4953 0.01557 * subject 5 2.18278 0.43656 5.2457 0.01271 * Residuals 10 0.83222 0.08322 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(fit) Call: lm(formula = volume ~ method + subject, data = lung) Residuals: Min 1Q Median 3Q Max -0.35556 -0.16389 -0.03889 0.17361 0.32778 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.17222 0.19232 16.494 1.4e-08 *** methodB 0.28333 0.16656 1.701 0.11975 methodC 0.60000 0.16656 3.602 0.00483 ** subject2 -0.83333 0.23555 -3.538 0.00538 ** subject3 0.10000 0.23555 0.425 0.68016 subject4 -0.06667 0.23555 -0.283 0.78293 subject5 -0.03333 0.23555 -0.142 0.89027 subject6 -0.60000 0.23555 -2.547 0.02900 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.2885 on 10 degrees of freedom Multiple R-squared: 0.7968, Adjusted R-squared: 0.6546 F-statistic: 5.603 on 7 and 10 DF, p-value: 0.00768 > kruskal.test(walk ~ group) Kruskal-Wallis rank sum test data: walk by group Kruskal-Wallis chi-squared = 6.8805, df = 3, p-value = 0.0758 > wilcox.test(zelazo$active,zelazo$ctr.8w) # first vs. last Wilcoxon rank sum test with continuity correction data: zelazo$active and zelazo$ctr.8w W = 3, p-value = 0.03493 alternative hypothesis: true location shift is not equal to 0 Warning message: In wilcox.test.default(zelazo$active, zelazo$ctr.8w) : cannot compute exact p-value with ties > wilcox.test(zelazo$active,unlist(zelazo[-1])) # first vs. rest Wilcoxon rank sum test with continuity correction data: zelazo$active and unlist(zelazo[-1]) W = 18, p-value = 0.02224 alternative hypothesis: true location shift is not equal to 0 Warning message: In wilcox.test.default(zelazo$active, unlist(zelazo[-1])) : cannot compute exact p-value with ties > friedman.test(volume ~ method | subject, data=lung) Friedman rank sum test data: volume and method and subject Friedman chi-squared = 5.8261, df = 2, p-value = 0.05431 > wilcox.test(lung$volume[lung$method=="A"], + lung$volume[lung$method=="C"], paired=TRUE) # etc. Wilcoxon signed rank test with continuity correction data: lung$volume[lung$method == "A"] and lung$volume[lung$method == "C"] V = 0, p-value = 0.05906 alternative hypothesis: true location shift is not equal to 0 Warning message: In wilcox.test.default(lung$volume[lung$method == "A"], lung$volume[lung$method == : cannot compute exact p-value with zeroes > attach(juul) > tapply(sqrt(igf1),tanner, sd, na.rm=TRUE) 1 2 3 4 5 3.203440 3.242068 3.479512 2.566712 3.164611 > plot(sqrt(igf1)~jitter(tanner)) > oneway.test(sqrt(igf1)~tanner) One-way analysis of means (not assuming equal variances) data: sqrt(igf1) and tanner F = 254.17, num df = 4.00, denom df = 162.43, p-value < 2.2e-16 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > binom.test(0, 10, p=.20, alt="less") Exact binomial test data: 0 and 10 number of successes = 0, number of trials = 10, p-value = 0.1074 alternative hypothesis: true probability of success is less than 0.2 95 percent confidence interval: 0.0000000 0.2588656 sample estimates: probability of success 0 > binom.test(0, 13, p=.20, alt="less") Exact binomial test data: 0 and 13 number of successes = 0, number of trials = 13, p-value = 0.05498 alternative hypothesis: true probability of success is less than 0.2 95 percent confidence interval: 0.0000000 0.2058167 sample estimates: probability of success 0 > binom.test(0, 14, p=.20, alt="less") Exact binomial test data: 0 and 14 number of successes = 0, number of trials = 14, p-value = 0.04398 alternative hypothesis: true probability of success is less than 0.2 95 percent confidence interval: 0.0000000 0.1926362 sample estimates: probability of success 0 > prop.test(c(210,122),c(747,661)) 2-sample test for equality of proportions with continuity correction data: c(210, 122) out of c(747, 661) X-squared = 17.612, df = 1, p-value = 2.709e-05 alternative hypothesis: two.sided 95 percent confidence interval: 0.05138139 0.14172994 sample estimates: prop 1 prop 2 0.2811245 0.1845688 > M <- matrix(c(23,7,18,13),2,2) > chisq.test(M) Pearson's Chi-squared test with Yates' continuity correction data: M X-squared = 1.6243, df = 1, p-value = 0.2025 > fisher.test(M) Fisher's Exact Test for Count Data data: M p-value = 0.1737 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.6936416 8.4948588 sample estimates: odds ratio 2.339104 > prop.test(M) 2-sample test for equality of proportions with continuity correction data: M X-squared = 1.6243, df = 1, p-value = 0.2025 alternative hypothesis: two.sided 95 percent confidence interval: -0.08462185 0.50657307 sample estimates: prop 1 prop 2 0.5609756 0.3500000 > tbl <- c(42, 157, 47, 62, 4, 15, 4, 1, 8, 28, 9, 7) > dim(tbl) <- c(2,2,3) > dimnames(tbl) <- list(c("A","B"), + c("not pierced","pierced"), + c("ok","broken","cracked")) > ftable(tbl) ok broken cracked A not pierced 42 4 8 pierced 47 4 9 B not pierced 157 15 28 pierced 62 1 7 > fisher.test(tbl["B",,]) # slice analysis Fisher's Exact Test for Count Data data: tbl["B", , ] p-value = 0.1116 alternative hypothesis: two.sided > fisher.test(tbl["A",,]) Fisher's Exact Test for Count Data data: tbl["A", , ] p-value = 1 alternative hypothesis: two.sided > fisher.test(margin.table(tbl,2:3)) # marginal Fisher's Exact Test for Count Data data: margin.table(tbl, 2:3) p-value = 0.3211 alternative hypothesis: two.sided > p <- seq(0,1,0.001) > pval <- sapply(p,function(p)binom.test(3,15,p=p)$p.value) > plot(p,pval,type="l") > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > power.t.test(power=.8,delta=.30,sd=.20) Two-sample t test power calculation n = 8.06031 delta = 0.3 sd = 0.2 sig.level = 0.05 power = 0.8 alternative = two.sided NOTE: n is number in *each* group > power.t.test(power=.8,delta=.30,sd=.20,alt="one.sided") Two-sample t test power calculation n = 6.298691 delta = 0.3 sd = 0.2 sig.level = 0.05 power = 0.8 alternative = one.sided NOTE: n is number in *each* group > (qnorm(.975)+qnorm(.8))^2*2*(.2/.3)^2 # approx. formula [1] 6.976782 > power.t.test(n=8, delta=.30, sd=.20) # power with eq.size Two-sample t test power calculation n = 8 delta = 0.3 sd = 0.2 sig.level = 0.05 power = 0.7965441 alternative = two.sided NOTE: n is number in *each* group > d2 <- .30 * sqrt(2/8) / sqrt(1/6+1/10) # corr.f.uneq. size > power.t.test(n=8, delta=d2, sd=.20) Two-sample t test power calculation n = 8 delta = 0.2904738 sd = 0.2 sig.level = 0.05 power = 0.7707066 alternative = two.sided NOTE: n is number in *each* group > power.prop.test(power=.9, p1=.6, p2=.75) Two-sample comparison of proportions power calculation n = 202.8095 p1 = 0.6 p2 = 0.75 sig.level = 0.05 power = 0.9 alternative = two.sided NOTE: n is number in *each* group > power.prop.test(power=.8, p1=.6, p2=.75) Two-sample comparison of proportions power calculation n = 151.8689 p1 = 0.6 p2 = 0.75 sig.level = 0.05 power = 0.8 alternative = two.sided NOTE: n is number in *each* group > curve(dt(x-3, 25), from=0, to=5) > curve(dt(x, 25, 3), add=TRUE) > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > attach(thuesen) > f <- cut(blood.glucose, c(4, 7, 9, 12, 20)) > levels(f) <- c("low", "intermediate", "high", "very high") > bcmort2 <- within(bcmort,{ + period <- area <- cohort + levels(period) <- rep(c("1991-2001","1981-1991"), each=2) + levels(area) <- rep(c("Cph+Frb","Nat"),2) + }) > summary(bcmort2) age cohort bc.deaths p.yr period 50-54:4 Study gr. :6 Min. : 9.0 Min. : 25600 1991-2001:12 55-59:4 Nat.ctr. :6 1st Qu.: 51.0 1st Qu.: 86434 1981-1991:12 60-64:4 Hist.ctr. :6 Median :104.0 Median : 169261 65-69:4 Hist.nat.ctr.:6 Mean :213.2 Mean : 396520 70-74:4 3rd Qu.:436.2 3rd Qu.: 781897 75-79:4 Max. :545.0 Max. :1067778 area Cph+Frb:12 Nat :12 > ashina.long <- reshape(ashina, direction="long", + varying=1:2, timevar="treat") > ashina.long <- within(ashina.long, { + m <- matrix(c(2,1,1,2),2) + id <- factor(id) + treat <- factor(treat) + grp <- factor(grp) + period <- factor(m[cbind(grp,treat)]) + rm(m) + }) > within(ashina.long, + period2 <- ifelse(treat != "active", + as.numeric(grp), 3 - as.numeric(grp)) + ) grp treat vas id period period2 1.active 1 active -167 1 2 2 2.active 1 active -127 2 2 2 3.active 1 active -58 3 2 2 4.active 1 active -103 4 2 2 5.active 1 active -35 5 2 2 6.active 1 active -164 6 2 2 7.active 1 active -3 7 2 2 8.active 1 active 25 8 2 2 9.active 1 active -61 9 2 2 10.active 1 active -45 10 2 2 11.active 2 active -38 11 1 1 12.active 2 active 29 12 1 1 13.active 2 active 2 13 1 1 14.active 2 active -18 14 1 1 15.active 2 active -74 15 1 1 16.active 2 active -72 16 1 1 1.plac 1 plac -102 1 1 1 2.plac 1 plac -39 2 1 1 3.plac 1 plac 32 3 1 1 4.plac 1 plac 28 4 1 1 5.plac 1 plac 16 5 1 1 6.plac 1 plac -42 6 1 1 7.plac 1 plac -27 7 1 1 8.plac 1 plac -30 8 1 1 9.plac 1 plac -47 9 1 1 10.plac 1 plac 8 10 1 1 11.plac 2 plac 12 11 2 2 12.plac 2 plac 11 12 2 2 13.plac 2 plac -9 13 2 2 14.plac 2 plac -1 14 2 2 15.plac 2 plac 3 15 2 2 16.plac 2 plac -36 16 2 2 > stroke.trim <- function(t1, t2) + subset(transform(stroke, + entry=t1, exit=pmin(t2, obsmonths), + dead=dead & obsmonths <= t2), + entry < exit) > stroke2 <- do.call(rbind, mapply(stroke.trim, + c(0,0.5,2,12), c(0.5,2,12,Inf), SIMPLIFY=F)) > table(stroke$dead) FALSE TRUE 344 485 > table(stroke2$dead) FALSE TRUE 1956 485 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > summary(lm(log(bwt) ~ log(bpd) + log(ad), data=secher)) Call: lm(formula = log(bwt) ~ log(bpd) + log(ad), data = secher) Residuals: Min 1Q Median 3Q Max -0.35074 -0.06741 -0.00792 0.05750 0.36360 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -5.8615 0.6617 -8.859 2.36e-14 *** log(bpd) 1.5519 0.2294 6.764 8.09e-10 *** log(ad) 1.4667 0.1467 9.998 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1068 on 104 degrees of freedom Multiple R-squared: 0.8583, Adjusted R-squared: 0.8556 F-statistic: 314.9 on 2 and 104 DF, p-value: < 2.2e-16 > summary(lm(log(bwt) ~ log(ad), data=secher)) Call: lm(formula = log(bwt) ~ log(ad), data = secher) Residuals: Min 1Q Median 3Q Max -0.58560 -0.06609 0.00184 0.07479 0.48435 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.4446 0.5103 -4.791 5.49e-06 *** log(ad) 2.2365 0.1105 20.238 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1275 on 105 degrees of freedom Multiple R-squared: 0.7959, Adjusted R-squared: 0.794 F-statistic: 409.6 on 1 and 105 DF, p-value: < 2.2e-16 > pairs(tlc) > summary(lm(log(tlc) ~ ., data=tlc)) Call: lm(formula = log(tlc) ~ ., data = tlc) Residuals: Min 1Q Median 3Q Max -0.32648 -0.12520 0.02708 0.12424 0.42314 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.157899 0.565698 -2.047 0.050163 . age -0.003647 0.003864 -0.944 0.353272 sex 0.086184 0.082015 1.051 0.302324 height 0.017328 0.004032 4.298 0.000188 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1898 on 28 degrees of freedom Multiple R-squared: 0.5927, Adjusted R-squared: 0.549 F-statistic: 13.58 on 3 and 28 DF, p-value: 1.181e-05 > opar <- par(mfrow=c(2,2)) > plot(lm(log(tlc) ~ ., data=tlc), which=1:4) > > drop1(lm(log(tlc) ~ ., data=tlc)) Single term deletions Model: log(tlc) ~ age + sex + height Df Sum of Sq RSS AIC 1.0087 -102.626 age 1 0.03210 1.0408 -103.623 sex 1 0.03978 1.0485 -103.388 height 1 0.66546 1.6742 -88.413 > drop1(lm(log(tlc) ~ . - age, data=tlc)) Single term deletions Model: log(tlc) ~ (age + sex + height) - age Df Sum of Sq RSS AIC 1.0408 -103.623 sex 1 0.05132 1.0921 -104.083 height 1 0.71317 1.7540 -88.923 > > par(mfrow=c(1,1)) > plot(log(tlc) ~ height, data=tlc) > par(mfrow=c(2,2)) > plot(lm(tlc ~ ., data=tlc), which=1:4) # slightly worse > par(opar) > summary(lm(sqrt(igf1) ~ age, data=juul2, subset=(age >= 25))) Call: lm(formula = sqrt(igf1) ~ age, data = juul2, subset = (age >= 25)) Residuals: Min 1Q Median 3Q Max -4.936 -1.184 0.113 1.036 4.021 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 18.89742 0.46904 40.29 <2e-16 *** age -0.10871 0.01034 -10.52 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.742 on 126 degrees of freedom (9 observations deleted due to missingness) Multiple R-squared: 0.4676, Adjusted R-squared: 0.4633 F-statistic: 110.6 on 1 and 126 DF, p-value: < 2.2e-16 > anova(lm(sqrt(igf1) ~ age + weight + height, + data=juul2, subset=(age >= 25))) Analysis of Variance Table Response: sqrt(igf1) Df Sum Sq Mean Sq F value Pr(>F) age 1 16.386 16.3861 3.7610 0.06132 . weight 1 2.454 2.4536 0.5632 0.45848 height 1 0.004 0.0036 0.0008 0.97715 Residuals 32 139.419 4.3568 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(lm(dl.milk ~ . - no, data=kfm)) Call: lm(formula = dl.milk ~ . - no, data = kfm) Residuals: Min 1Q Median 3Q Max -1.74201 -0.81173 -0.00926 0.78326 2.52646 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -11.681839 4.361561 -2.678 0.010363 * sexgirl -0.499532 0.312672 -1.598 0.117284 weight 1.349124 0.322450 4.184 0.000135 *** ml.suppl -0.002233 0.001241 -1.799 0.078829 . mat.weight 0.006212 0.023708 0.262 0.794535 mat.height 0.072278 0.030169 2.396 0.020906 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.075 on 44 degrees of freedom Multiple R-squared: 0.5459, Adjusted R-squared: 0.4943 F-statistic: 10.58 on 5 and 44 DF, p-value: 1.03e-06 > summary(lm(dl.milk ~ . - no - mat.weight, data=kfm)) Call: lm(formula = dl.milk ~ . - no - mat.weight, data = kfm) Residuals: Min 1Q Median 3Q Max -1.77312 -0.81196 -0.00683 0.76988 2.52240 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -12.112571 3.997860 -3.030 0.00405 ** sexgirl -0.494675 0.308875 -1.602 0.11626 weight 1.372524 0.306612 4.476 5.14e-05 *** ml.suppl -0.002313 0.001190 -1.943 0.05824 . mat.height 0.076363 0.025560 2.988 0.00454 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.064 on 45 degrees of freedom Multiple R-squared: 0.5452, Adjusted R-squared: 0.5047 F-statistic: 13.48 on 4 and 45 DF, p-value: 2.658e-07 > summary(lm(dl.milk ~ . - no - mat.weight - sex, data=kfm)) Call: lm(formula = dl.milk ~ . - no - mat.weight - sex, data = kfm) Residuals: Min 1Q Median 3Q Max -2.06540 -0.74758 -0.02408 0.67488 2.79882 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -13.064926 4.020073 -3.250 0.00216 ** weight 1.464781 0.306231 4.783 1.81e-05 *** ml.suppl -0.002237 0.001209 -1.850 0.07074 . mat.height 0.077600 0.025979 2.987 0.00451 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.082 on 46 degrees of freedom Multiple R-squared: 0.5192, Adjusted R-squared: 0.4879 F-statistic: 16.56 on 3 and 46 DF, p-value: 1.953e-07 > summary(lm(dl.milk ~ weight + mat.height, data=kfm)) Call: lm(formula = dl.milk ~ weight + mat.height, data = kfm) Residuals: Min 1Q Median 3Q Max -2.19598 -0.82149 0.01822 0.75582 2.83375 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -11.92014 4.07325 -2.926 0.00527 ** weight 1.42862 0.31338 4.559 3.67e-05 *** mat.height 0.07063 0.02636 2.680 0.01013 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.109 on 47 degrees of freedom Multiple R-squared: 0.4835, Adjusted R-squared: 0.4615 F-statistic: 22 on 2 and 47 DF, p-value: 1.811e-07 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > ashina.long <- reshape(ashina, direction="long", + varying=1:2, timevar="treat") > ashina.long <- within(ashina.long, { + m <- matrix(c(2,1,1,2),2) + id <- factor(id) + treat <- factor(treat) + grp <- factor(grp) + period <- factor(m[cbind(grp,treat)]) + rm(m) + }) > fit.ashina <- lm(vas ~ id + period + treat, data=ashina.long) > drop1(fit.ashina, test="F") Single term deletions Model: vas ~ id + period + treat Df Sum of Sq RSS AIC F value Pr(>F) 19679 241.49 id 15 51137 70816 252.47 2.4254 0.05287 . period 1 1505 21184 241.85 1.0709 0.31830 treat 1 11603 31282 254.32 8.2550 0.01228 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(fit.ashina) Analysis of Variance Table Response: vas Df Sum Sq Mean Sq F value Pr(>F) id 15 51137 3409.2 2.4254 0.05287 . period 1 4608 4608.0 3.2783 0.09171 . treat 1 11603 11603.3 8.2550 0.01228 * Residuals 14 19679 1405.6 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > attach(ashina) > dd <- vas.active - vas.plac > t.test(dd[grp==1], -dd[grp==2], var.eq=T) Two Sample t-test data: dd[grp == 1] and -dd[grp == 2] t = -2.8731, df = 14, p-value = 0.01228 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -137.39089 -19.94244 sample estimates: mean of x mean of y -53.50000 25.16667 > t.test(dd[grp==1], dd[grp==2], var.eq=T) Two Sample t-test data: dd[grp == 1] and dd[grp == 2] t = -1.0348, df = 14, p-value = 0.3183 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -87.05756 30.39089 sample estimates: mean of x mean of y -53.50000 -25.16667 > attach(tb.dilute) > anova(lm(reaction ~ animal + logdose)) Analysis of Variance Table Response: reaction Df Sum Sq Mean Sq F value Pr(>F) animal 5 35.208 7.042 8.2641 0.002527 ** logdose 2 72.396 36.198 42.4817 1.295e-05 *** Residuals 10 8.521 0.852 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > ld <- c(0.5, 0, -0.5)[logdose] > anova(lm(reaction ~ animal + ld)) Analysis of Variance Table Response: reaction Df Sum Sq Mean Sq F value Pr(>F) animal 5 35.208 7.042 6.4353 0.00496 ** ld 1 68.880 68.880 62.9489 7.067e-06 *** Residuals 11 12.036 1.094 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(lm(reaction ~ animal + ld)) Call: lm(formula = reaction ~ animal + ld) Residuals: Min 1Q Median 3Q Max -1.72917 -0.22917 -0.08333 0.23437 2.41667 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 11.0000 0.6039 18.214 1.45e-09 *** animal2 1.0000 0.8541 1.171 0.266408 animal3 2.1667 0.8541 2.537 0.027637 * animal4 2.1667 0.8541 2.537 0.027637 * animal5 3.5833 0.8541 4.195 0.001497 ** animal6 4.0833 0.8541 4.781 0.000571 *** ld 4.7917 0.6039 7.934 7.07e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.046 on 11 degrees of freedom Multiple R-squared: 0.8963, Adjusted R-squared: 0.8398 F-statistic: 15.85 on 6 and 11 DF, p-value: 7.816e-05 > 4.7917 + 0.6039 * qt(c(.025,.975), 11) [1] 3.462525 6.120875 > # or: > confint(lm(reaction ~ animal + ld))["ld",] 2.5 % 97.5 % 3.462408 6.120925 > > slopes <- reaction[logdose==0.5] - reaction[logdose==-0.5] > t.test(slopes) One Sample t-test data: slopes t = 16.429, df = 5, p-value = 1.525e-05 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 4.041914 5.541420 sample estimates: mean of x 4.791667 > > anova(lm(reaction ~ animal*ld)) Analysis of Variance Table Response: reaction Df Sum Sq Mean Sq F value Pr(>F) animal 5 35.208 7.042 3.9264 0.0630805 . ld 1 68.880 68.880 38.4076 0.0008133 *** animal:ld 5 1.276 0.255 0.1423 0.9753664 Residuals 6 10.760 1.793 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > a <- gl(2, 2, 8) > b <- gl(2, 4, 8) > x <- 1:8 > y <- c(1:4,8:5) > z <- rnorm(8) > model.matrix(~ a:b) ; lm(z ~ a:b) (Intercept) a1:b1 a2:b1 a1:b2 a2:b2 1 1 1 0 0 0 2 1 1 0 0 0 3 1 0 1 0 0 4 1 0 1 0 0 5 1 0 0 1 0 6 1 0 0 1 0 7 1 0 0 0 1 8 1 0 0 0 1 attr(,"assign") [1] 0 1 1 1 1 attr(,"contrasts") attr(,"contrasts")$a [1] "contr.treatment" attr(,"contrasts")$b [1] "contr.treatment" Call: lm(formula = z ~ a:b) Coefficients: (Intercept) a1:b1 a2:b1 a1:b2 a2:b2 0.2537 -0.3222 -0.1409 -0.6558 NA > model.matrix(~ a * b) ; lm(z ~ a * b) (Intercept) a2 b2 a2:b2 1 1 0 0 0 2 1 0 0 0 3 1 1 0 0 4 1 1 0 0 5 1 0 1 0 6 1 0 1 0 7 1 1 1 1 8 1 1 1 1 attr(,"assign") [1] 0 1 2 3 attr(,"contrasts") attr(,"contrasts")$a [1] "contr.treatment" attr(,"contrasts")$b [1] "contr.treatment" Call: lm(formula = z ~ a * b) Coefficients: (Intercept) a2 b2 a2:b2 -0.06848 0.18134 -0.33360 0.47449 > model.matrix(~ a:x) ; lm(z ~ a:x) (Intercept) a1:x a2:x 1 1 1 0 2 1 2 0 3 1 0 3 4 1 0 4 5 1 5 0 6 1 6 0 7 1 0 7 8 1 0 8 attr(,"assign") [1] 0 1 1 attr(,"contrasts") attr(,"contrasts")$a [1] "contr.treatment" Call: lm(formula = z ~ a:x) Coefficients: (Intercept) a1:x a2:x 0.12999 -0.09671 0.00482 > model.matrix(~ a * x) ; lm(z ~ a * x) (Intercept) a2 x a2:x 1 1 0 1 0 2 1 0 2 0 3 1 1 3 3 4 1 1 4 4 5 1 0 5 0 6 1 0 6 0 7 1 1 7 7 8 1 1 8 8 attr(,"assign") [1] 0 1 2 3 attr(,"contrasts") attr(,"contrasts")$a [1] "contr.treatment" Call: lm(formula = z ~ a * x) Coefficients: (Intercept) a2 x a2:x 0.02594 0.32161 -0.07464 0.04477 > model.matrix(~ b * (x + y)) ; lm(z ~ b * (x + y)) (Intercept) b2 x y b2:x b2:y 1 1 0 1 1 0 0 2 1 0 2 2 0 0 3 1 0 3 3 0 0 4 1 0 4 4 0 0 5 1 1 5 8 5 8 6 1 1 6 7 6 7 7 1 1 7 6 7 6 8 1 1 8 5 8 5 attr(,"assign") [1] 0 1 2 3 4 5 attr(,"contrasts") attr(,"contrasts")$b [1] "contr.treatment" Call: lm(formula = z ~ b * (x + y)) Coefficients: (Intercept) b2 x y b2:x b2:y 0.002029 -0.128603 0.066866 -0.058803 NA NA > attach(secretin) > model1 <- lm(gluc ~ person * time) > model2 <- lm(gluc ~ person + time) > model3 <- lm(gluc ~ person * time20plus + time) > model4 <- lm(gluc ~ person * time20plus + time.comb) > tt <- c(20,30,60,90,0)[time] > plot(fitted(model4)~tt,pch=as.character(person)) > bp.obese <- transform(bp.obese,sex=factor(sex, labels=c("M","F"))) > plot(log(bp) ~ log(obese), pch=c(20,21)[sex], data=bp.obese) > summary(lm(log(bp) ~ sex, data=bp.obese)) Call: lm(formula = log(bp) ~ sex, data = bp.obese) Residuals: Min 1Q Median 3Q Max -0.30067 -0.07468 -0.01583 0.09436 0.50927 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.84396 0.02037 237.851 <2e-16 *** sexF -0.01569 0.02701 -0.581 0.563 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1351 on 100 degrees of freedom Multiple R-squared: 0.003365, Adjusted R-squared: -0.006601 F-statistic: 0.3376 on 1 and 100 DF, p-value: 0.5625 > summary(lm(log(bp) ~ sex + log(obese), data=bp.obese)) Call: lm(formula = log(bp) ~ sex + log(obese), data = bp.obese) Residuals: Min 1Q Median 3Q Max -0.21260 -0.09068 -0.00472 0.05675 0.43846 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.79180 0.02249 213.05 < 2e-16 *** sexF -0.06393 0.02744 -2.33 0.0218 * log(obese) 0.31235 0.07366 4.24 5.02e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1249 on 99 degrees of freedom Multiple R-squared: 0.1565, Adjusted R-squared: 0.1395 F-statistic: 9.187 on 2 and 99 DF, p-value: 0.0002188 > summary(lm(log(bp) ~ sex*log(obese), data=bp.obese)) Call: lm(formula = log(bp) ~ sex * log(obese), data = bp.obese) Residuals: Min 1Q Median 3Q Max -0.21455 -0.09020 -0.00741 0.05710 0.42968 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.80589 0.02896 165.930 <2e-16 *** sexF -0.09047 0.04393 -2.059 0.0421 * log(obese) 0.22798 0.13159 1.733 0.0863 . sexF:log(obese) 0.12310 0.15895 0.774 0.4405 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1252 on 98 degrees of freedom Multiple R-squared: 0.1617, Adjusted R-squared: 0.136 F-statistic: 6.3 on 3 and 98 DF, p-value: 0.0005933 > vitcap2 <- transform(vitcap2,group=factor(group, + labels=c("exp>10", + "exp<10", "unexp"))) > attach(vitcap2) > plot(vital.capacity~age, pch=(20:22)[group]) > vit.fit <- lm(vital.capacity ~ age*group) > summary(vit.fit) Call: lm(formula = vital.capacity ~ age * group) Residuals: Min 1Q Median 3Q Max -1.24497 -0.36929 0.01977 0.43681 1.13953 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.18344 0.99358 8.236 3.28e-12 *** age -0.08511 0.01967 -4.327 4.44e-05 *** groupexp<10 -1.95341 1.10481 -1.768 0.0810 . groupunexp -2.50315 1.04184 -2.403 0.0187 * age:groupexp<10 0.03858 0.02327 1.658 0.1014 age:groupunexp 0.05450 0.02107 2.587 0.0116 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5942 on 78 degrees of freedom Multiple R-squared: 0.422, Adjusted R-squared: 0.385 F-statistic: 11.39 on 5 and 78 DF, p-value: 2.871e-08 > drop1(vit.fit, test="F") Single term deletions Model: vital.capacity ~ age * group Df Sum of Sq RSS AIC F value Pr(>F) 27.535 -81.689 age:group 2 2.4995 30.035 -78.391 3.5402 0.03376 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > for (i in 1:3) abline(lm(vital.capacity ~ age, + subset=as.numeric(group)==i), lty=i) > legend(20, 3.5 ,legend=levels(group), pch=20:22, lty=1:3) > juul.prepub <- subset(juul, tanner==1) > > summary(lm(sqrt(igf1)~age, data=juul.prepub, subset= sex==1)) Call: lm(formula = sqrt(igf1) ~ age, data = juul.prepub, subset = sex == 1) Residuals: Min 1Q Median 3Q Max -7.8319 -1.7113 -0.1784 1.4845 11.1109 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.50382 0.46198 18.41 <2e-16 *** age 0.64564 0.05357 12.05 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.678 on 190 degrees of freedom (99 observations deleted due to missingness) Multiple R-squared: 0.4333, Adjusted R-squared: 0.4303 F-statistic: 145.3 on 1 and 190 DF, p-value: < 2.2e-16 > summary(lm(sqrt(igf1)~age, data=juul.prepub, subset= sex==2)) Call: lm(formula = sqrt(igf1) ~ age, data = juul.prepub, subset = sex == 2) Residuals: Min 1Q Median 3Q Max -6.2532 -1.3766 -0.0526 1.1144 6.1815 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.7055 0.9047 11.833 < 2e-16 *** age 0.4761 0.1020 4.667 8.18e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.18 on 117 degrees of freedom (105 observations deleted due to missingness) Multiple R-squared: 0.157, Adjusted R-squared: 0.1498 F-statistic: 21.79 on 1 and 117 DF, p-value: 8.184e-06 > > summary(lm(sqrt(igf1)~age*factor(sex), data=juul.prepub)) Call: lm(formula = sqrt(igf1) ~ age * factor(sex), data = juul.prepub) Residuals: Min 1Q Median 3Q Max -7.832 -1.586 -0.085 1.337 11.111 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.50382 0.43127 19.718 <2e-16 *** age 0.64564 0.05001 12.911 <2e-16 *** factor(sex)2 2.20169 1.12350 1.960 0.0509 . age:factor(sex)2 -0.16952 0.12721 -1.333 0.1836 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.5 on 307 degrees of freedom (204 observations deleted due to missingness) Multiple R-squared: 0.3969, Adjusted R-squared: 0.391 F-statistic: 67.33 on 3 and 307 DF, p-value: < 2.2e-16 > summary(lm(sqrt(igf1)~age+factor(sex), data=juul.prepub)) Call: lm(formula = sqrt(igf1) ~ age + factor(sex), data = juul.prepub) Residuals: Min 1Q Median 3Q Max -7.7224 -1.5711 -0.0889 1.3310 11.1234 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.70903 0.40335 21.59 <2e-16 *** age 0.61945 0.04604 13.45 <2e-16 *** factor(sex)2 0.75669 0.29445 2.57 0.0106 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.503 on 308 degrees of freedom (204 observations deleted due to missingness) Multiple R-squared: 0.3934, Adjusted R-squared: 0.3894 F-statistic: 99.86 on 2 and 308 DF, p-value: < 2.2e-16 > summary(fit.aicopt <- step(lm(dl.milk ~ . - no, data=kfm))) Start: AIC=12.83 dl.milk ~ (no + sex + weight + ml.suppl + mat.weight + mat.height) - no Df Sum of Sq RSS AIC - mat.weight 1 0.0793 50.919 10.911 50.840 12.833 - sex 1 2.9492 53.789 13.652 - ml.suppl 1 3.7408 54.580 14.383 - mat.height 1 6.6318 57.471 16.963 - weight 1 20.2269 71.067 27.580 Step: AIC=10.91 dl.milk ~ sex + weight + ml.suppl + mat.height Df Sum of Sq RSS AIC 50.919 10.911 - sex 1 2.9023 53.821 11.682 - ml.suppl 1 4.2735 55.193 12.940 - mat.height 1 10.0999 61.019 17.958 - weight 1 22.6741 73.593 27.326 Call: lm(formula = dl.milk ~ sex + weight + ml.suppl + mat.height, data = kfm) Residuals: Min 1Q Median 3Q Max -1.77312 -0.81196 -0.00683 0.76988 2.52240 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -12.112571 3.997860 -3.030 0.00405 ** sexgirl -0.494675 0.308875 -1.602 0.11626 weight 1.372524 0.306612 4.476 5.14e-05 *** ml.suppl -0.002313 0.001190 -1.943 0.05824 . mat.height 0.076363 0.025560 2.988 0.00454 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.064 on 45 degrees of freedom Multiple R-squared: 0.5452, Adjusted R-squared: 0.5047 F-statistic: 13.48 on 4 and 45 DF, p-value: 2.658e-07 > opar <- par(mfrow=c(2,2)) > plot(fit.aicopt, which=1:4) > kfm[32,] no dl.milk sex weight ml.suppl mat.weight mat.height 32 37 7.22 girl 5.336 590 58 160 > summary(kfm) no dl.milk sex weight ml.suppl Min. : 1.00 Min. : 4.440 boy :25 Min. :4.120 Min. : 0.00 1st Qu.: 28.75 1st Qu.: 6.555 girl:25 1st Qu.:4.977 1st Qu.: 16.25 Median : 50.50 Median : 7.660 Median :5.352 Median : 57.50 Mean : 52.36 Mean : 7.504 Mean :5.319 Mean : 96.00 3rd Qu.: 79.75 3rd Qu.: 8.428 3rd Qu.:5.637 3rd Qu.:103.75 Max. :105.00 Max. :10.430 Max. :6.578 Max. :590.00 mat.weight mat.height Min. :47.00 Min. :153.0 1st Qu.:55.00 1st Qu.:162.0 Median :58.00 Median :167.0 Mean :59.96 Mean :167.4 3rd Qu.:64.50 3rd Qu.:172.0 Max. :80.00 Max. :185.0 > summary(update(fit.aicopt, ~ . - sex)) Call: lm(formula = dl.milk ~ weight + ml.suppl + mat.height, data = kfm) Residuals: Min 1Q Median 3Q Max -2.06540 -0.74758 -0.02408 0.67488 2.79882 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -13.064926 4.020073 -3.250 0.00216 ** weight 1.464781 0.306231 4.783 1.81e-05 *** ml.suppl -0.002237 0.001209 -1.850 0.07074 . mat.height 0.077600 0.025979 2.987 0.00451 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.082 on 46 degrees of freedom Multiple R-squared: 0.5192, Adjusted R-squared: 0.4879 F-statistic: 16.56 on 3 and 46 DF, p-value: 1.953e-07 > plot(update(fit.aicopt, ~ . - sex - ml.suppl), which=1:4) > par(opar) > juulyoung <- subset(juul, age < 25) > juulyoung <- transform(juulyoung, + sex=factor(sex), tanner=factor(tanner)) > fit.untf <- lm(igf1 ~ age * sex * tanner, data=juulyoung, + na.action=na.exclude) > plot(fitted(fit.untf) ~ age, data=juulyoung, + col=c("red","green")[sex]) > fit.log <- update(fit.untf, log(igf1) ~ .) > fit.sqrt <- update(fit.untf, sqrt(igf1) ~ .) > opar <- par(mfrow=c(2,2)) > plot(fit.untf, which=1:4) > plot(fit.log, which=1:4) > plot(fit.sqrt, which=1:4) > par(opar) > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > summary(glm(mal~age+log(ab), binomial, data=malaria)) Call: glm(formula = mal ~ age + log(ab), family = binomial, data = malaria) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.57234 0.95184 2.702 0.006883 ** age -0.06546 0.06772 -0.967 0.333703 log(ab) -0.68235 0.19552 -3.490 0.000483 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 116.652 on 99 degrees of freedom Residual deviance: 98.017 on 97 degrees of freedom AIC: 104.02 Number of Fisher Scoring iterations: 5 > attach(graft.vs.host) > type <- factor(type,labels=c("AML", "ALL", "CML")) > m1 <- glm(gvhd~rcpage+donage+type+preg+log(index), binomial) > m1a <- glm(gvhd~rcpage+donage+type+preg+index, binomial) > summary(m1) Call: glm(formula = gvhd ~ rcpage + donage + type + preg + log(index), family = binomial) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -5.47420 2.59647 -2.108 0.0350 * rcpage 0.01597 0.08205 0.195 0.8457 donage 0.10724 0.08262 1.298 0.1943 typeALL -0.42824 1.25389 -0.342 0.7327 typeCML 1.61517 1.30743 1.235 0.2167 preg 1.66236 1.18103 1.408 0.1593 log(index) 1.81874 0.88919 2.045 0.0408 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 51.049 on 36 degrees of freedom Residual deviance: 26.252 on 30 degrees of freedom AIC: 40.252 Number of Fisher Scoring iterations: 6 > summary(m1a) Call: glm(formula = gvhd ~ rcpage + donage + type + preg + index, family = binomial) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -6.08898 2.67548 -2.276 0.0229 * rcpage 0.03273 0.08224 0.398 0.6906 donage 0.09267 0.07974 1.162 0.2452 typeALL -0.37184 1.26596 -0.294 0.7690 typeCML 1.80883 1.34681 1.343 0.1793 preg 1.92847 1.18157 1.632 0.1027 index 0.68422 0.38065 1.797 0.0723 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 51.049 on 36 degrees of freedom Residual deviance: 25.964 on 30 degrees of freedom AIC: 39.964 Number of Fisher Scoring iterations: 6 > drop1(m1, test="Chisq") Single term deletions Model: gvhd ~ rcpage + donage + type + preg + log(index) Df Deviance AIC LRT Pr(>Chi) 26.252 40.252 rcpage 1 26.291 38.291 0.0382 0.84495 donage 1 28.082 40.082 1.8292 0.17622 type 2 28.959 38.959 2.7064 0.25841 preg 1 28.447 40.447 2.1947 0.13849 log(index) 1 32.343 44.343 6.0904 0.01359 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > drop1(update(m1, ~ . - rcpage), test="Chisq") Single term deletions Model: gvhd ~ donage + type + preg + log(index) Df Deviance AIC LRT Pr(>Chi) 26.291 38.291 donage 1 28.832 38.832 2.5411 0.11092 type 2 29.403 37.403 3.1126 0.21092 preg 1 28.812 38.812 2.5218 0.11228 log(index) 1 32.630 42.630 6.3397 0.01181 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > drop1(update(m1, ~ . - rcpage - type), test="Chisq") Single term deletions Model: gvhd ~ donage + preg + log(index) Df Deviance AIC LRT Pr(>Chi) 29.403 37.403 donage 1 33.654 39.654 4.2506 0.0392375 * preg 1 31.068 37.068 1.6653 0.1968922 log(index) 1 41.723 47.723 12.3198 0.0004482 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > drop1(update(m1, ~ . - rcpage - type - preg), test="Chisq") Single term deletions Model: gvhd ~ donage + log(index) Df Deviance AIC LRT Pr(>Chi) 31.068 37.068 donage 1 37.740 41.740 6.671 0.0097993 ** log(index) 1 45.476 49.476 14.407 0.0001473 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(m2 <- glm(gvhd~donage + log(index), binomial)) Call: glm(formula = gvhd ~ donage + log(index), family = binomial) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -5.45399 2.08147 -2.620 0.00879 ** donage 0.14594 0.06465 2.257 0.02399 * log(index) 2.17773 0.78986 2.757 0.00583 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 51.049 on 36 degrees of freedom Residual deviance: 31.068 on 34 degrees of freedom AIC: 37.068 Number of Fisher Scoring iterations: 5 > confint(m2) Waiting for profiling to be done... 2.5 % 97.5 % (Intercept) -10.38849705 -1.9550010 donage 0.03262764 0.2946854 log(index) 0.89146245 4.0706703 > ## normal approximation: > est <- coefficients(summary(m2))[,1] > se <- coefficients(summary(m2))[,2] > est + cbind(qnorm(.025)*se, qnorm(.975)*se) [,1] [,2] (Intercept) -9.53358695 -1.3743858 donage 0.01922664 0.2726562 log(index) 0.62964464 3.7258196 > confint.default(m2) 2.5 % 97.5 % (Intercept) -9.53358695 -1.3743858 donage 0.01922664 0.2726562 log(index) 0.62964464 3.7258196 > counts <- c(13,40,157,40,21,61) > total <- c(108,264,375,310,181,162) > age <- gl(3,1,6) > type <- gl(2,3,6) > anova(glm(counts/total~age+type,weights=total, binomial), + test="Chisq") Analysis of Deviance Table Model: binomial, link: logit Response: counts/total Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 5 136.910 age 2 134.857 3 2.054 <2e-16 *** type 1 1.273 2 0.780 0.2591 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > juul.girl <- transform(subset(juul,age>8 & age<20 & + complete.cases(menarche)), + menarche=factor(menarche)) > logit.menarche <- glm(menarche~age+I(age^2)+I(age^3), + binomial, data=juul.girl) Warning message: glm.fit: fitted probabilities numerically 0 or 1 occurred > probit.menarche <- glm(menarche~age+I(age^2)+I(age^3), + binomial(probit), data=juul.girl) Warning message: glm.fit: fitted probabilities numerically 0 or 1 occurred > summary(logit.menarche) Call: glm(formula = menarche ~ age + I(age^2) + I(age^3), family = binomial, data = juul.girl) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -766.0804 377.7631 -2.028 0.0426 * age 175.2400 87.1896 2.010 0.0444 * I(age^2) -13.4140 6.6902 -2.005 0.0450 * I(age^3) 0.3435 0.1707 2.013 0.0442 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 719.39 on 518 degrees of freedom Residual deviance: 191.80 on 515 degrees of freedom AIC: 199.8 Number of Fisher Scoring iterations: 12 > summary(probit.menarche) Call: glm(formula = menarche ~ age + I(age^2) + I(age^3), family = binomial(probit), data = juul.girl) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -328.53729 184.42237 -1.781 0.0748 . age 74.94081 42.66126 1.757 0.0790 . I(age^2) -5.73565 3.27880 -1.749 0.0802 . I(age^3) 0.14722 0.08373 1.758 0.0787 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 719.39 on 518 degrees of freedom Residual deviance: 191.98 on 515 degrees of freedom AIC: 199.98 Number of Fisher Scoring iterations: 13 > Age=seq(8,20,.1) > newages <- data.frame(age=Age) > p.logit <- predict(logit.menarche,newdata=newages,type="resp") > p.probit <- predict(probit.menarche,newdata=newages,type="resp") > matplot(Age,cbind(p.probit,p.logit),type="l") > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > library(survival) Attaching package: 'survival' The following object is masked from 'package:ISwR': lung > attach(graft.vs.host) > plot(survfit(Surv(time,dead)~gvhd)) > survdiff(Surv(time,dead)~gvhd) Call: survdiff(formula = Surv(time, dead) ~ gvhd) N Observed Expected (O-E)^2/E (O-E)^2/V gvhd=0 20 7 11.99 2.07 6.55 gvhd=1 17 11 6.01 4.14 6.55 Chisq= 6.6 on 1 degrees of freedom, p= 0.01 > summary(coxph(Surv(time,dead) ~ gvhd)) # for comparison Call: coxph(formula = Surv(time, dead) ~ gvhd) n= 37, number of events= 18 coef exp(coef) se(coef) z Pr(>|z|) gvhd 1.2419 3.4620 0.5144 2.414 0.0158 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 gvhd 3.462 0.2888 1.263 9.489 Concordance= 0.671 (se = 0.052 ) Likelihood ratio test= 6.19 on 1 df, p=0.01 Wald test = 5.83 on 1 df, p=0.02 Score (logrank) test = 6.55 on 1 df, p=0.01 > summary(coxph(Surv(time,dead) ~ + gvhd + log(index) + donage + rcpage + preg)) Call: coxph(formula = Surv(time, dead) ~ gvhd + log(index) + donage + rcpage + preg) n= 37, number of events= 18 coef exp(coef) se(coef) z Pr(>|z|) gvhd 1.29074 3.63548 0.69131 1.867 0.0619 . log(index) -0.49623 0.60882 0.35579 -1.395 0.1631 donage 0.03440 1.03500 0.04000 0.860 0.3899 rcpage -0.01230 0.98778 0.04085 -0.301 0.7634 preg 1.13641 3.11555 0.57058 1.992 0.0464 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 gvhd 3.6355 0.2751 0.9378 14.093 log(index) 0.6088 1.6425 0.3031 1.223 donage 1.0350 0.9662 0.9569 1.119 rcpage 0.9878 1.0124 0.9118 1.070 preg 3.1155 0.3210 1.0183 9.533 Concordance= 0.748 (se = 0.065 ) Likelihood ratio test= 14.33 on 5 df, p=0.01 Wald test = 14.86 on 5 df, p=0.01 Score (logrank) test = 17.48 on 5 df, p=0.004 > attach(melanom) > cox1 <- coxph(Surv(days, status==1) ~ + log(thick) + sex + strata(ulc)) > new <- data.frame(sex=2, thick=c(0.1, 0.2, 0.5)) > svfit <- survfit(cox1,newdata=new) > plot(svfit[2], ylim=c(.985, 1)) > summary(coxph(Surv(obsmonths, dead)~age+sex, data=stroke)) Call: coxph(formula = Surv(obsmonths, dead) ~ age + sex, data = stroke) n= 829, number of events= 485 coef exp(coef) se(coef) z Pr(>|z|) age 0.049289 1.050524 0.004374 11.268 <2e-16 *** sexMale 0.022332 1.022584 0.100717 0.222 0.825 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.051 0.9519 1.0416 1.060 sexMale 1.023 0.9779 0.8394 1.246 Concordance= 0.66 (se = 0.013 ) Likelihood ratio test= 158.6 on 2 df, p=<2e-16 Wald test = 137.7 on 2 df, p=<2e-16 Score (logrank) test = 136.3 on 2 df, p=<2e-16 > summary(coxph(Surv(obsmonths, dead)~sex, data=stroke)) Call: coxph(formula = Surv(obsmonths, dead) ~ sex, data = stroke) n= 829, number of events= 485 coef exp(coef) se(coef) z Pr(>|z|) sexMale -0.36571 0.69371 0.09615 -3.803 0.000143 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 sexMale 0.6937 1.442 0.5746 0.8376 Concordance= 0.547 (se = 0.011 ) Likelihood ratio test= 14.97 on 1 df, p=1e-04 Wald test = 14.47 on 1 df, p=1e-04 Score (logrank) test = 14.63 on 1 df, p=1e-04 > with(stroke, tapply(age,sex,mean)) Female Male 73.32549 64.38558 > stroke.trim <- function(t1, t2) + subset(transform(stroke, + entry=t1, exit=pmin(t2, obsmonths), + dead=dead & obsmonths <= t2), + entry < exit) > stroke2 <- do.call(rbind, mapply(stroke.trim, + c(0,0.5,2,12), c(0.5,2,12,Inf), SIMPLIFY=F)) > summary(coxph(Surv(entry, exit, dead)~age+sex, data=stroke2)) Call: coxph(formula = Surv(entry, exit, dead) ~ age + sex, data = stroke2) n= 2441, number of events= 485 coef exp(coef) se(coef) z Pr(>|z|) age 0.049289 1.050524 0.004374 11.268 <2e-16 *** sexMale 0.022332 1.022584 0.100717 0.222 0.825 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.051 0.9519 1.0416 1.060 sexMale 1.023 0.9779 0.8394 1.246 Concordance= 0.66 (se = 0.013 ) Likelihood ratio test= 158.6 on 2 df, p=<2e-16 Wald test = 137.7 on 2 df, p=<2e-16 Score (logrank) test = 136.3 on 2 df, p=<2e-16 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > bcmort2 <- within(bcmort,{ + period <- area <- cohort + levels(period) <- rep(c("1991-2001","1981-1991"), each=2) + levels(area) <- rep(c("Cph+Frb","Nat"),2) + }) > bcfit <- glm(bc.deaths ~ (age + period + area)^2, poisson, + offset=log(p.yr), data=bcmort2) > summary(bcfit) Call: glm(formula = bc.deaths ~ (age + period + area)^2, family = poisson, data = bcmort2, offset = log(p.yr)) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -8.71065 0.20674 -42.133 < 2e-16 *** age55-59 0.59227 0.23373 2.534 0.011277 * age60-64 1.08881 0.22321 4.878 1.07e-06 *** age65-69 1.24619 0.21945 5.679 1.36e-08 *** age70-74 1.62200 0.21709 7.472 7.93e-14 *** age75-79 1.93431 0.23583 8.202 2.36e-16 *** period1981-1991 0.76507 0.15539 4.924 8.50e-07 *** areaNat -0.36892 0.20220 -1.825 0.068070 . age55-59:period1981-1991 -0.34415 0.14996 -2.295 0.021735 * age60-64:period1981-1991 -0.55380 0.14786 -3.746 0.000180 *** age65-69:period1981-1991 -0.58432 0.14731 -3.967 7.29e-05 *** age70-74:period1981-1991 -0.65739 0.14776 -4.449 8.62e-06 *** age75-79:period1981-1991 -0.56839 0.16360 -3.474 0.000512 *** age55-59:areaNat 0.69382 0.22708 3.055 0.002248 ** age60-64:areaNat 0.52311 0.21633 2.418 0.015600 * age65-69:areaNat 0.49641 0.21251 2.336 0.019492 * age70-74:areaNat 0.38044 0.21047 1.808 0.070675 . age75-79:areaNat 0.33547 0.22865 1.467 0.142325 period1981-1991:areaNat -0.29262 0.08796 -3.327 0.000879 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1039.855 on 23 degrees of freedom Residual deviance: 5.458 on 5 degrees of freedom AIC: 202.4 Number of Fisher Scoring iterations: 4 > drop1(bcfit, test="Chisq") Single term deletions Model: bc.deaths ~ (age + period + area)^2 Df Deviance AIC LRT Pr(>Chi) 5.458 202.41 age:period 5 34.792 221.74 29.334 1.994e-05 *** age:area 5 17.279 204.23 11.822 0.0373179 * period:area 1 16.775 211.72 11.317 0.0007678 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > confint(bcfit, parm="period1981-1991:areaNat") Waiting for profiling to be done... 2.5 % 97.5 % -0.4663402 -0.1213920 > stroke.trim <- function(t1, t2) + subset(transform(stroke, + entry=t1, exit=pmin(t2, obsmonths), + dead=dead & obsmonths <= t2), + entry < exit) > stroke2 <- do.call(rbind, mapply(stroke.trim, + c(0,0.5,2,12), c(0.5,2,12,Inf), SIMPLIFY=F)) > summary(glm(dead~sex+age+factor(entry), poisson, + offset=log(exit-entry), data=stroke2)) Call: glm(formula = dead ~ sex + age + factor(entry), family = poisson, data = stroke2, offset = log(exit - entry)) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.252978 0.353585 -12.028 <2e-16 *** sexMale 0.021233 0.100668 0.211 0.833 age 0.050231 0.004386 11.452 <2e-16 *** factor(entry)0.5 -1.588098 0.127052 -12.500 <2e-16 *** factor(entry)2 -3.220317 0.126616 -25.434 <2e-16 *** factor(entry)12 -3.878493 0.123296 -31.457 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 3745.4 on 2440 degrees of freedom Residual deviance: 2443.2 on 2435 degrees of freedom AIC: 3425.2 Number of Fisher Scoring iterations: 7 > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > ## IGNORE_RDIFF_BEGIN > girls <- subset(juul2, age<20 & age>5 & sex==2) > boys <- subset(juul2, age<20 & age>5 & sex==1) > young <- subset(juul2, age<20 & age>5) > stval <- c(alpha=exp(5.3),beta=exp(0.42),gamma=0.15) > fit.boys <- nls(height~alpha*exp(-beta*exp(-gamma*age)), + start=stval, data=boys) > fit.girls <- nls(height~alpha*exp(-beta*exp(-gamma*age)), + start=stval, data=girls) > fit.young <- nls(height~alpha*exp(-beta*exp(-gamma*age)), + start=stval, data=young) > ms.pooled <- (deviance(fit.boys) + deviance(fit.girls))/(499+625) > ms.diff <- (deviance(fit.young) - + deviance(fit.boys) - deviance(fit.girls))/3 > ms.diff/ms.pooled [1] 90.58118 > fit.young2 <- nls(height~(alpha+da*(sex==1))* + exp(-(beta+db*(sex==1))* + exp(-(gamma+dg*(sex==1))*age)), + start=c(alpha=exp(5.3),beta=exp(0.42),gamma=0.15, + da=0, db=0, dg=0), data=young) > summary(fit.young2) Formula: height ~ (alpha + da * (sex == 1)) * exp(-(beta + db * (sex == 1)) * exp(-(gamma + dg * (sex == 1)) * age)) Parameters: Estimate Std. Error t value Pr(>|t|) alpha 180.77121 1.93392 93.474 < 2e-16 *** beta 1.18704 0.06218 19.091 < 2e-16 *** gamma 0.16391 0.01015 16.144 < 2e-16 *** da 62.03692 11.26004 5.509 4.46e-08 *** db -0.01107 0.06477 -0.171 0.864 dg -0.08488 0.01306 -6.500 1.20e-10 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 6.528 on 1124 degrees of freedom Number of iterations to convergence: 6 Achieved convergence tolerance: 1.981e-06 (11 observations deleted due to missingness) > anova(fit.young, fit.young2) Analysis of Variance Table Model 1: height ~ alpha * exp(-beta * exp(-gamma * age)) Model 2: height ~ (alpha + da * (sex == 1)) * exp(-(beta + db * (sex == 1)) * exp(-(gamma + dg * (sex == 1)) * age)) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 1127 59475 2 1124 47896 3 11580 90.581 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > e1 <- subset(philion, experiment==1) > fit <- nls(sqrt(response) ~ sqrt(ymax / (1 + (dose/ec50)^exp(la))), + start=list(ymax=28, ec50=.3, la=0), data=e1, + lower=c(.1,.0001,-Inf), algorithm="port") > summary(fit) Formula: sqrt(response) ~ sqrt(ymax/(1 + (dose/ec50)^exp(la))) Parameters: Estimate Std. Error t value Pr(>|t|) ymax 25.4331 4.7681 5.334 0.00177 ** ec50 0.2087 0.1439 1.450 0.19730 la -0.1732 0.2399 -0.722 0.49747 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5311 on 6 degrees of freedom Algorithm "port", convergence message: relative convergence (4) > confint(fit) Waiting for profiling to be done... 2.5% 97.5% ymax 15.16409523 41.0317406 ec50 0.03655499 1.3770101 la -0.72605274 0.6841776 > # p <- profile(fit, alphamax=.2) > # alphamax semantics got changed in R 2.8 > p <- profile(fit) > par(mfrow=c(3,1)) > plot(p) > confint(p) 2.5% 97.5% ymax 15.16043812 41.031201 ec50 0.05332678 1.377884 la -0.72610588 0.684268 > e1 <- subset(philion, experiment==1) > fit1 <- nls(sqrt(response) ~ sqrt(ymax / (1 + dose/b)^exp(la)), + start=list(ymax=28, b=.3, la=0), data=e1, + lower=c(.1,.0001,-Inf), algorithm="port") > summary(fit1) Formula: sqrt(response) ~ sqrt(ymax/(1 + dose/b)^exp(la)) Parameters: Estimate Std. Error t value Pr(>|t|) ymax 23.38316 4.00491 5.839 0.00111 ** b 0.29938 0.32259 0.928 0.38919 la -0.02689 0.42805 -0.063 0.95196 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5474 on 6 degrees of freedom Algorithm "port", convergence message: relative convergence (4) > fit2 <- nls(sqrt(response) ~ sqrt(ymax / (1 + + dose/(d50/(2^(1/exp(la))-1)))^exp(la)), + start=list(ymax=28, d50=.3, la=0), data=e1, + lower=c(.1,.0001,-Inf), algorithm="port") > summary(fit2) Formula: sqrt(response) ~ sqrt(ymax/(1 + dose/(d50/(2^(1/exp(la)) - 1)))^exp(la)) Parameters: Estimate Std. Error t value Pr(>|t|) ymax 23.38316 4.00491 5.839 0.00111 ** d50 0.31080 0.18567 1.674 0.14517 la -0.02689 0.42805 -0.063 0.95196 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5474 on 6 degrees of freedom Algorithm "port", convergence message: relative convergence (4) > ## IGNORE_RDIFF_END > dd <- seq(0,1,,200) > yy <- predict(fit, newdata=data.frame(dose=dd)) > y1 <- predict(fit2, newdata=data.frame(dose=dd)) > matplot(dd,cbind(yy,y1)^2, type="l") > rm(list=ls()) > while(search()[2] != "package:ISwR") detach() > > proc.time() user system elapsed 3.10 0.12 3.23