hwk3-8output - > library(MASS) > > visc <‐

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: > library(MASS) > > visc <‐ read.table("visc.dat",header=T) > visc p f v 1 0 0 26 2 0 12 28 3 0 24 30 4 0 36 32 5 0 48 34 6 0 60 37 7 10 0 18 8 10 12 19 9 10 24 20 10 10 36 21 11 10 48 24 12 10 60 24 13 20 0 12 14 20 12 14 15 20 24 14 16 20 36 16 17 20 48 17 18 20 60 17 19 30 12 12 20 30 24 12 21 30 36 13 22 30 48 14 23 30 60 14 > > # Part a > > visc$P <‐ factor(visc$p) > visc$F <‐ factor(visc$f) > > bc1 <‐ boxcox(v ~ P + F, data = visc, + lambda = seq(‐1, 1, length = 100),plotit=T) > > bc2 <‐ boxcox(v ~ P + F, data = visc, + lambda = seq(‐.4, .25, length = 500),plotit=T) > > bc2$x[bc2$y==max(bc2$y)] [1] ‐0.06392786 > max(bc2$y) [1] ‐15.68125 > lim <‐ max(bc2$y)‐0.5*qchisq(0.90,1) > cbind(bc2$x,bc2$y,lim) lim [1,] ‐0.4000000000 ‐18.89343 ‐17.03402 [2,] ‐0.3986973948 ‐18.87204 ‐17.03402 <snip> [97,] ‐0.2749498998 ‐17.07117 ‐17.03402 [98,] ‐0.2736472946 ‐17.05520 ‐17.03402 [99,] ‐0.2723446894 ‐17.03930 ‐17.03402 [100,] ‐0.2710420842 ‐17.02348 ‐17.03402 [101,] ‐0.2697394790 ‐17.00773 ‐17.03402 <snip> [257,] ‐0.0665330661 ‐15.68139 ‐17.03402 [258,] ‐0.0652304609 ‐15.68126 ‐17.03402 [259,] ‐0.0639278557 ‐15.68125 ‐17.03402 [260,] ‐0.0626252505 ‐15.68135 ‐17.03402 [261,] ‐0.0613226453 ‐15.68157 ‐17.03402 <snip> [412,] 0.1353707415 ‐17.01294 ‐17.03402 [413,] 0.1366733467 ‐17.02966 ‐17.03402 [414,] 0.1379759519 ‐17.04648 ‐17.03402 [415,] 0.1392785571 ‐17.06339 ‐17.03402 <snip> [499,] 0.2486973948 ‐18.76751 ‐17.03402 [500,] 0.2500000000 ‐18.79074 ‐17.03402 > > lm1 <‐ lm(I((v^(‐.064)‐1)/(‐.064))~F+P,data=visc) > summary(lm1) Call: lm(formula = I((v^(‐0.064) ‐ 1)/(‐0.064)) ~ F + P, data = visc) Residuals: Min 1Q Median 3Q Max ‐0.0338447 ‐0.0152421 0.0004949 0.0123468 0.0297721 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.93298 0.01527 192.017 < 2e‐16 *** F12 0.08691 0.01775 4.897 0.000236 *** F24 0.11142 0.01775 6.278 2.03e‐05 *** F36 0.17952 0.01775 10.115 8.10e‐08 *** F48 0.24735 0.01775 13.937 1.34e‐09 *** F60 0.26417 0.01775 14.885 5.64e‐10 *** P10 ‐0.32018 0.01323 ‐24.205 7.98e‐13 *** P20 ‐0.60175 0.01323 ‐45.490 < 2e‐16 *** P30 ‐0.74751 0.01408 ‐53.081 < 2e‐16 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.02291 on 14 degrees of freedom Multiple R‐squared: 0.9962, Adjusted R‐squared: 0.994 F‐statistic: 458.7 on 8 and 14 DF, p‐value: 1.361e‐15 > anova(lm1) Analysis of Variance Table Response: I((v^(‐0.064) ‐ 1)/(‐0.064)) Df Sum Sq Mean Sq F value Pr(>F) F 5 0.11187 0.02237 42.621 5.543e‐08 *** P 3 1.81447 0.60482 1152.149 < 2.2e‐16 *** Residuals 14 0.00735 0.00052 ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plot(lm1) Waiting to confirm page change... Waiting to confirm page change... Error in object$coefficients : $ operator is invalid for atomic vectors > > lm2 <‐ lm(I((v^(‐.064)‐1)/(‐.064))~f+p,data=visc) > summary(lm2) Call: lm(formula = I((v^(‐0.064) ‐ 1)/(‐0.064)) ~ f + p, data = visc) Residuals: Min 1Q Median 3Q Max ‐0.100021 ‐0.030925 0.003916 0.039741 0.099562 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.9064884 0.0250576 115.992 < 2e‐16 *** f 0.0045798 0.0005746 7.971 1.23e‐07 *** p ‐0.0254541 0.0010447 ‐24.366 2.42e‐16 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.05469 on 20 degrees of freedom Multiple R‐squared: 0.9691, Adjusted R‐squared: 0.966 F‐statistic: 313.2 on 2 and 20 DF, p‐value: 8.039e‐16 > anova(lm2,lm1) Analysis of Variance Table Model 1: I((v^(‐0.064) ‐ 1)/(‐0.064)) ~ f + p Model 2: I((v^(‐0.064) ‐ 1)/(‐0.064)) ~ F + P Res.Df RSS Df Sum of Sq F Pr(>F) 1 20 0.059828 2 14 0.007349 6 0.052479 16.662 1.211e‐05 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > logLik(lm1) 'log Lik.' 59.92375 (df=10) > logLik(lm2) 'log Lik.' 35.80976 (df=4) > logLR <‐ 2*(logLik(lm1)‐logLik(lm2)) > logLR [1] 48.22799 attr(,"nall") [1] 23 attr(,"nobs") [1] 23 attr(,"df") [1] 10 attr(,"class") [1] "logLik" > 1‐pchisq(logLR,6) [1] 1.063935e‐08 attr(,"nall") [1] 23 attr(,"nobs") [1] 23 attr(,"df") [1] 10 attr(,"class") [1] "logLik" > > bc3 <‐ boxcox(v ~ p + f, data = visc, + lambda = seq(‐2, 1, length = 500),plotit=T) > bc4 <‐ boxcox(v ~ p + f, data = visc, + lambda = seq(‐1.1, ‐.1, length = 500),plotit=T) > > bc4$x[bc4$y==max(bc4$y)] [1] ‐0.6410822 > max(bc4$y) [1] ‐36.46484 > > lm3 <‐ lm(I((v^(‐.64)‐1)/(‐.64))~f+p,data=visc) > summary(lm3) Call: lm(formula = I((v^(‐0.64) ‐ 1)/(‐0.64)) ~ f + p, data = visc) Residuals: Min 1Q Median 3Q Max ‐0.024160 ‐0.003169 0.001036 0.004740 0.012133 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.360e+00 3.965e‐03 343.089 < 2e‐16 *** f 8.266e‐04 9.091e‐05 9.092 1.53e‐08 *** p ‐4.604e‐03 1.653e‐04 ‐27.857 < 2e‐16 *** ‐‐‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.008654 on 20 degrees of freedom Multiple R‐squared: 0.9761, Adjusted R‐squared: 0.9738 F‐statistic: 409.3 on 2 and 20 DF, p‐value: < 2.2e‐16 > > logLik(lm1) 'log Lik.' 59.92375 (df=10) > logLik(lm3) 'log Lik.' 78.21632 (df=4) > #these loglikelihoods are based ont the assumptions that lambda is known, > # so they do not form the basis of an appropriate LRT to compare models > # from parts (a) and (c). That is, the following LRT is NOT appropriate. > # Instead, see hwk3‐7c.sas > #logLR <‐ 2*(logLik(lm1)‐logLik(lm3)) > #logLR > #1‐pchisq(logLR,6) > > ...
View Full Document

This note was uploaded on 11/13/2011 for the course STAT 8620 taught by Professor Hall during the Fall '11 term at University of Florida.

Ask a homework question - tutors are online