5.4
(a)
library(GAD)
## Loading required package: matrixStats
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
cut <- as.fixed(cut)
FR <- as.fixed(FR)
model<-aov(finish~cut+FR+cut*FR)
GAD::gad(model)
## Analysis of Variance Table
##
## Response: finish
## Df Sum Sq Mean Sq F value Pr(>F)
## cut 3 2125.11 708.37 24.6628 1.652e-07 ***
## FR 2 3160.50 1580.25 55.0184 1.086e-09 ***
## cut:FR 6 557.06 92.84 3.2324 0.01797 *
## Residual 24 689.33 28.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\alpha\beta_{ij}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\). all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i
Since P-value of Cut:FR=0.01797<0.05, reject null hypothesis, significant interaction between Cut and Feed Rate. So herein we stop exploration on the main factors, even if we got P-value of Cut=1.652e-07<0.05 and P-value of Feed Rate= 1.086e-09<0.05, indicating both significant impacts on means. We need to plot the interaction.
Interaction plot:
f1 <- finish[1:3]
f2 <- finish[4:6]
f3 <- finish[7:9]
f4 <- finish[10:12]
f5 <- finish[13:15]
f6 <- finish[16:18]
f7 <- finish[19:21]
f8 <- finish[22:24]
f9 <- finish[25:27]
f10 <- finish[28:30]
f11 <- finish[31:33]
f12 <- finish[34:36]
FR1 <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
avg1 <- ave(f1)
avg2 <- ave(f2)
avg3 <- ave(f3)
avg4 <- ave(f4)
avg5 <- ave(f5)
avg6 <- ave(f6)
avg7 <- ave(f7)
avg8 <- ave(f8)
avg9 <- ave(f9)
avg10 <- ave(f10)
avg11 <- ave(f11)
avg12 <- ave(f12)
?plo
## No documentation for 'plo' in specified packages and libraries:
## you could try '??plo'
plot(FR1, c(f1,f2,f3,f4),main="Scatterplot with average", xlab="Feed Rate",ylab="Finish", pch=19, col="red", ylim=range(c(min(finish), max(finish))))
lines(FR1, c(avg1,avg2,avg3,avg4), type = "b", col = "red", lwd = 3, pch = 15)
points(FR1, c(f5,f6,f7,f8),pch = 19,col="blue")
lines(FR1, c(avg5,avg6,avg7,avg8), type = "b", col = "blue", lwd = 3, pch = 15)
points(FR1, c(f9,f10,f11,f12),pch = 19,col="green")
lines(FR1, c(avg9,avg10,avg11,avg12), type = "b", col = "green", lwd = 3, pch = 15)
legend(x=0.21,y=80,c("Feed rate=0.2 in/min","Feed rate=0.25 in/min","Feed rate=0.30 in/min"),pch=c(19,19),col = c("red","blue","green"))

From interaction plot, we verify that do have significant interaction.
(c)
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
f1 <- finish[1:12]
f2 <- finish[13:24]
f3 <- finish[25:36]
sdev <-c(rep(sd(f1),12),rep(sd(f2),12),rep(sd(f3),12))
avg <- c(ave(f1), ave(f2), ave(f3))
plot(FR, finish, ylim=range(c(avg-sdev, avg+sdev)), main="Scatterplot with average and error bar", xlab="Feed Rate ", ylab="Surface Finish", pch=19)
points(FR,avg,pch = 15,col="red")
arrows(FR, avg-sdev, FR, avg+sdev, length=0.05, angle=90, code=3, col="red")
lines(FR, avg, type = "b", col = "red", lwd = 3, pch = 15)

##?lines()
avg
## [1] 81.58333 81.58333 81.58333 81.58333 81.58333 81.58333 81.58333
## [8] 81.58333 81.58333 81.58333 81.58333 81.58333 97.58333 97.58333
## [15] 97.58333 97.58333 97.58333 97.58333 97.58333 97.58333 97.58333
## [22] 97.58333 97.58333 97.58333 103.83333 103.83333 103.83333 103.83333
## [29] 103.83333 103.83333 103.83333 103.83333 103.83333 103.83333 103.83333
## [36] 103.83333
For feed rate 0.20, average is 81.58333; For feed rate 0.25, average is 97.58333; For feed rate 0.30, average is 103.83333.
(d)
Please see part (a)
————————————————————————————————
5.34
All fixed effect, we can use ANOVA or GAD
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
block1 <- c(rep(seq(1,3),12))
dat1 <- cbind(finish,cut,FR,block1)
dat1 <- as.data.frame(dat1)
dat1$cut <- as.factor(dat1$cut)
dat1$FR <- as.factor(dat1$FR)
dat1$block1 <- as.factor(dat1$block1)
str(dat1)
## 'data.frame': 36 obs. of 4 variables:
## $ finish: num 74 64 60 79 68 73 82 88 92 99 ...
## $ cut : Factor w/ 4 levels "0.15","0.18",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ FR : Factor w/ 3 levels "0.2","0.25","0.3": 1 1 1 1 1 1 1 1 1 1 ...
## $ block1: Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
model<-aov(finish~cut+FR+block1+cut*FR,data=dat1)
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## cut 3 2125.1 708.4 30.637 4.89e-08 ***
## FR 2 3160.5 1580.2 68.346 3.64e-10 ***
## block1 2 180.7 90.3 3.907 0.03532 *
## cut:FR 6 557.1 92.8 4.015 0.00726 **
## Residuals 22 508.7 23.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\gamma_{k}+\alpha\beta_{ij}+\epsilon_{ijkl}\)
\(\gamma_{k}\) is block effect.
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\). all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i
Here is an important concept: block do not interact with other factors, so we have the linear equation and hypothesis above.
Since P-value of Cut: FR=0.00726<0.05, reject null hypothesis, significant interaction between Cut and Feed Rate. We stop exploration of main effect, even if P-value of Cut=4.89e-08<0.05 and P-value of Feed Rate= 3.64e-10<0.05, indicating significant impact of both on means. We need to plot interaction.
Interaction plot is the same as I plotted in Question 5.4.
Since we introduce the block effect, nuisance variability will introduce the sum square block (SSB) and decrease the sum square error(SSE). MSE is likely to drop, f increases and p-value decreases, all the p-value for Cut:FR, Cut and FR have been reduced. Blocking was useful in this experiment.
variance component: Since I=4,J=3,K=3,L=1
We can use these following equations:
\(\sigma_{\alpha\beta}^2=\frac{MSAB-MSE}{KL}=\frac{92.8-23.1}{3*1}=23.2333\)
\(\sigma_{\alpha}^2=\frac{MSA-MSAB}{JKL}=\frac{708.4-92.8}{3*3X1}=68.4\).
\(\sigma_{\beta}^2=\frac{MSB-MSAB}{IKL}=\frac{1580.2-92.8}{4*3X1}=123.95\).
\(\sigma_{\gamma}^2=\frac{MSC-MSE}{IJL}=\frac{90.3-23.1}{4*3X1}=5.6\).
E[MSE]=\(\sigma^2=23.1\).
We further confirmed the above mentioned P-value related conclusions by these variances’ estimations, variance=0 means no interaction or no impact on means, vice versa.
————————————————————————————————
Problem 13.5
Question assumes no interaction. If we count interaction in:
density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position+temp*position)
GAD::gad(model)
## Analysis of Variance Table
##
## Response: density
## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 945342 472671 1155.518 0.0008647 ***
## position 1 7160 7160 15.998 0.0017624 **
## temp:position 2 818 409 0.914 0.4271101
## Residual 12 5371 448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (model)




\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\alpha\beta_{ij}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\), all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i
Question assumes no interaction. Since P-value of temp:position=0.4271101>0.05, we confirm no interaction.
P-value of Temp=0.0008647<0.05, P-value of Position=0.0017624<0.05, we reject both null hypothesis, both of them significantly affect the means.
Model’s adequacy: From norm plot and residual plot, we know that it is normally distributed and variances are roughly equal.
Analysis of variance: Since I=3,J=2,K=3, since no interaction,\(\sigma_{\alpha\beta}^2=0\)
We can use these following equations:
\(\sigma_{\alpha\beta}^2=\frac{MSAB-MSE}{K}=\frac{409-448}{3}<0\), so \(\sigma_{\alpha\beta}^2=0\).
\(\sigma_{\alpha}^2=\frac{MSA-MSAB}{JK}=\frac{472671-409}{2*3}=78710.3333\).
\(\sigma_{\beta}^2=\frac{MSB-MSAB}{IK}=\frac{7160-409}{3*3}=750.1111\).
E[MSE]=\(\sigma^2=448\).
We further confirmed the above mentioned P-value related conclusions by these variances’ estimations, variance=0 means no interaction or no impact on means, vice versa.
Question assumes no interaction.If we do not count interaction in:
density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position)
GAD::gad(model)
## Analysis of Variance Table
##
## Response: density
## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 945342 472671 1069.257 4.924e-16 ***
## position 1 7160 7160 16.197 0.001254 **
## Residual 14 6189 442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (model)




\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\), all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
Question assumes no interaction.
P-value of Temp=4.924e-16<0.05, P-value of Position=0.001254<0.05, we reject both null hypothesis, both of them significantly affect the means.
Model’s adequacy: From norm plot and residual plot, we know that it is normally distributed and variances are roughly equal.
Analysis of variance: Since I=3,J=2,K=3 and no interaction per se.
We can use these following equations:
\(\sigma_{\alpha}^2=\frac{MSA-MSE}{JK}=\frac{472671-442}{2*3}=78704.8333\).
\(\sigma_{\beta}^2=\frac{MSB-MSE}{IK}=\frac{7160-442}{3*3}=746.4444\).
E[MSE]=\(\sigma^2=442\).
We further confirmed the above mentioned P-value related conclusions by these variances’ estimations, variance=0 means no interaction or no impact on means, vice versa.
————————————————————————————————
Problem 13.6
dat1 <- c(50,49,50,50,48,51,52,52,51,51,51,51,53,50,50,54,52,51,49,51,50,48,50,51,48,49,48,48,49,48,52,50,50,52,50,50,51,51,51,51,50,50,52,50,49,53,48,50,50,51,50,51,48,49,47,46,49,46,47,48)
operator <- c(rep(c(rep(1,3),rep(2,3)),10))
operator
## [1] 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1
## [39] 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2
part <- c(rep(1,6),rep(2,6),rep(3,6),rep(4,6),rep(5,6),rep(6,6),rep(7,6),rep(8,6),rep(9,6),rep(10,6))
part
## [1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5
## [26] 5 5 5 5 5 6 6 6 6 6 6 7 7 7 7 7 7 8 8 8 8 8 8 9 9
## [51] 9 9 9 9 10 10 10 10 10 10
operator <- as.fixed(operator)
part <- as.random(part)
model<-aov(dat1~operator+part+operator*part)
GAD::gad(model)
## Analysis of Variance Table
##
## Response: dat1
## Df Sum Sq Mean Sq F value Pr(>F)
## operator 1 0.417 0.4167 0.6923 0.4269
## part 9 99.017 11.0019 7.3346 3.216e-06 ***
## operator:part 9 5.417 0.6019 0.4012 0.9270
## Residual 40 60.000 1.5000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (model)




Since test were made in random order indicated by question, there is not block effect, herein the measurements 1 2 3 do not mean blocking. This a mixed effect.
\(\overline{y_{ijk}}=\mu+\alpha_{i}+\beta_{j}+\alpha\beta_{ij}+\epsilon_{ijk}\)
\(H_{0}\): \(\alpha_{i}=0\), all i
\(H_{1}\): \(\alpha_{i}\neq0\), some i
\(H_{0}\): \(\beta_{j}=0\), all i
\(H_{1}\): \(\beta_{j}\neq0\), some i
\(H_{0}\): \(\alpha\beta_{ij}=0\), all i
\(H_{1}\): \(\alpha\beta_{ij}\neq0\), some i
P-value for operator:part=0.9270>0.05, fail to reject null hypothesis, no significant interaction. So analyze main effect. P-value for operator=0.4269>0.05, no impact by operator on means. P-value for part=3.216e-06<0.05, there is impact by part on means.
Model’s adequacy: From norm plot and residual plot, we can draw that it is roughly noramlly distributed, and variances are roughly equal. Assumptions are verified.
Analysis of variance: Since I=2,J=10,K=3, E[MSE]=\(\sigma^2\)=1.5000, E[MSAB]=\(\sigma^2+K\sigma_{\alpha\beta}^2\), \(0.6019=0.1500+3*\sigma_{\alpha\beta}^2\), so \(\sigma_{\alpha\beta}^2\)<0, assume \(\sigma_{\alpha\beta}^2\)=0.
We can ALSO use these following equations:
\(\sigma_{\alpha\beta}^2=\frac{MSAB-MSE}{K}=\frac{0.6019-1.5000}{3}<0\), so \(\sigma_{\alpha\beta}^2=0\).
\(\sigma_{\alpha}^2=\frac{MSA-MSAB}{JK}=\frac{0.4167-0.6019}{10*3}<0\), so \(\sigma_{\alpha}^2\)=0.
\(\sigma_{\beta}^2=\frac{MSB-MSAB}{IK}=\frac{11.0019-0.6019}{2*3}=1.7333\),
We further confirmed the above mentioned P-value related conclusions by these variances’ estimations, variance=0 means no interaction or no impact on means, vice versa.
————————————————————————————————
Raw Code:
library(GAD)
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
cut <- as.fixed(cut)
FR <- as.fixed(FR)
model<-aov(finish~cut+FR+cut*FR)
GAD::gad(model)
f1 <- finish[1:3]
f2 <- finish[4:6]
f3 <- finish[7:9]
f4 <- finish[10:12]
f5 <- finish[13:15]
f6 <- finish[16:18]
f7 <- finish[19:21]
f8 <- finish[22:24]
f9 <- finish[25:27]
f10 <- finish[28:30]
f11 <- finish[31:33]
f12 <- finish[34:36]
FR1 <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
avg1 <- ave(f1)
avg2 <- ave(f2)
avg3 <- ave(f3)
avg4 <- ave(f4)
avg5 <- ave(f5)
avg6 <- ave(f6)
avg7 <- ave(f7)
avg8 <- ave(f8)
avg9 <- ave(f9)
avg10 <- ave(f10)
avg11 <- ave(f11)
avg12 <- ave(f12)
?plo
plot(FR1, c(f1,f2,f3,f4),main="Scatterplot with average", xlab="Feed Rate",ylab="Finish", pch=19, col="red", ylim=range(c(min(finish), max(finish))))
lines(FR1, c(avg1,avg2,avg3,avg4), type = "b", col = "red", lwd = 3, pch = 15)
points(FR1, c(f5,f6,f7,f8),pch = 19,col="blue")
lines(FR1, c(avg5,avg6,avg7,avg8), type = "b", col = "blue", lwd = 3, pch = 15)
points(FR1, c(f9,f10,f11,f12),pch = 19,col="green")
lines(FR1, c(avg9,avg10,avg11,avg12), type = "b", col = "green", lwd = 3, pch = 15)
legend(x=0.21,y=80,c("Feed rate=0.2 in/min","Feed rate=0.25 in/min","Feed rate=0.30 in/min"),pch=c(19,19),col = c("red","blue","green"))
plot (model)
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
f1 <- finish[1:12]
f2 <- finish[13:24]
f3 <- finish[25:36]
sdev <-c(rep(sd(f1),12),rep(sd(f2),12),rep(sd(f3),12))
avg <- c(ave(f1), ave(f2), ave(f3))
plot(FR, finish, ylim=range(c(avg-sdev, avg+sdev)), main="Scatterplot with average and error bar", xlab="Feed Rate ", ylab="Surface Finish", pch=19)
points(FR,avg,pch = 15,col="red")
arrows(FR, avg-sdev, FR, avg+sdev, length=0.05, angle=90, code=3, col="red")
lines(FR, avg, type = "b", col = "red", lwd = 3, pch = 15)
##?lines()
avg
finish <- c(74,64,60,79,68,73,82,88,92,99,104,96,92,86,88,98,104,88,99,108,95,104,110,99,99,98,102,104,99,95,108,110,99,114,111,107)
cut <- c(rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3),rep(0.15,3),rep(0.18,3),rep(0.2,3),rep(0.25,3))
FR <- c(rep(0.2,12),rep(0.25,12),rep(0.3,12))
block1 <- c(rep(seq(1,3),12))
dat1 <- cbind(finish,cut,FR,block1)
dat1 <- as.data.frame(dat1)
dat1$cut <- as.factor(dat1$cut)
dat1$FR <- as.factor(dat1$FR)
dat1$block1 <- as.factor(dat1$block1)
str(dat1)
model<-aov(finish~cut+FR+block1+cut*FR,data=dat1)
summary(model)
density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position+temp*position)
GAD::gad(model)
plot
density <- c(570,565,583,1063,1080,1043,565,510,590,528,547,521,988,1026,1004,526,538,532)
temp <- c(rep(800,3),rep(825,3),rep(850,3),rep(800,3),rep(825,3),rep(850,3))
position <- c(rep (1,9),rep(2,9))
temp <- as.fixed(temp)
position <- as.random(position)
model<-aov(density~temp+position)
GAD::gad(model)
plot (model)
dat1 <- c(50,49,50,50,48,51,52,52,51,51,51,51,53,50,50,54,52,51,49,51,50,48,50,51,48,49,48,48,49,48,52,50,50,52,50,50,51,51,51,51,50,50,52,50,49,53,48,50,50,51,50,51,48,49,47,46,49,46,47,48)
operator <- c(rep(c(rep(1,3),rep(2,3)),10))
operator
part <- c(rep(1,6),rep(2,6),rep(3,6),rep(4,6),rep(5,6),rep(6,6),rep(7,6),rep(8,6),rep(9,6),rep(10,6))
part
operator <- as.fixed(operator)
part <- as.random(part)
model<-aov(dat1~operator+part+operator*part)
GAD::gad(model)
plot (model)