這篇是以晶圓曝光時間(能量)和線寬(cd)資料分析, 討論資料相關性分析。

蝕刻製程

載入必要套件, 只須執行一次

library(data.table) #載入套件 data.table 才可使用fread()函數

讀取資料

filename<-"H:/Etch分析講義/stat_etch.csv" #完整資料路徑與名稱, 學生要根據自己的資料放的位置與名稱修改這條指令
x<-fread(filename,header="auto")

查看資料

str(x)
## Classes 'data.table' and 'data.frame':   165 obs. of  2 variables:
##  $ 能量: num  28.9 29.3 28.9 28.9 28.7 ...
##  $ 線寬: num  2.13 2.13 2.14 2.13 2.11 ...
##  - attr(*, ".internal.selfref")=<externalptr>

練習

  • 請問能量與線寬為“屬質”或“屬量”變數?
  • 請問能量與線寬的衡量尺度, 名目,順序, 區間或比例?

顯示資料

head(x, 6)
##     能量     線寬
## 1: 28.94 2.133501
## 2: 29.32 2.133610
## 3: 28.85 2.137660
## 4: 28.95 2.128999
## 5: 28.68 2.105443
## 6: 29.82 2.079087
tail(x, 6)
##     能量     線寬
## 1: 24.17 2.047422
## 2: 22.72 2.017259
## 3: 24.93 2.063129
## 4: 25.23 2.086112
## 5: 23.81 2.043845
## 6: 23.57 2.051075

基本統計量

mean(x$能量,na.rm=TRUE)
## [1] 27.146
mean(x$線寬,na.rm=TRUE)
## [1] 2.082079
#使用apply() 求算每個變數的統計量
apply(x,2, mean)
##      能量      線寬 
## 27.146000  2.082079
median(x$能量,na.rm=TRUE)
## [1] 26.2
median(x$線寬,na.rm=TRUE)
## [1] 2.078988
sd(x$能量,na.rm=TRUE)
## [1] 3.528171
sd(x$線寬,na.rm=TRUE)
## [1] 0.04130594
#第25, 50,75個百分位數和最大最小值
quantile(x$能量,type=2,na.rm=TRUE)
##    0%   25%   50%   75%  100% 
## 20.37 24.48 26.20 29.34 36.66
quantile(x$線寬,type=2,na.rm=TRUE)
##       0%      25%      50%      75%     100% 
## 2.001263 2.052519 2.078988 2.104579 2.207242
#其他百分位數
quantile(x$能量,probs=0.45,type=2,na.rm=TRUE)
##   45% 
## 25.93
quantile(x$線寬,probs=0.65,type=2,na.rm=TRUE)
##      65% 
## 2.093503
min(x$能量,na.rm=TRUE);min(x$線寬,na.rm=TRUE)
## [1] 20.37
## [1] 2.001263
max(x$能量,na.rm=TRUE);max(x$線寬,na.rm=TRUE)
## [1] 36.66
## [1] 2.207242

練習

  • 使用apply()函數計算以下統計量
    • min
    • max
    • median

練習

  • 使用quantile函數計算以下百分位數 (type=2)
    • 能量的第20個百分位數
    • 能量的第90個百分位數
    • 線寬的第95個百分位數
  • summary(資料物件): 摘要資訊
summary(x)
##       能量            線寬      
##  Min.   :20.37   Min.   :2.001  
##  1st Qu.:24.48   1st Qu.:2.053  
##  Median :26.20   Median :2.079  
##  Mean   :27.15   Mean   :2.082  
##  3rd Qu.:29.34   3rd Qu.:2.105  
##  Max.   :36.66   Max.   :2.207

  • summary()是一個簡便的函數, 可以顯示資料中每一個變數幾項統計量

  • 這些統計量有最小值(Min), 最大值(Max), 三個四分位數(1st Qu, Median, 3rt Qu)與平均值(Mean)

  • sd(資料向量,na.rm=TRUE): 標準差

sd(x$能量,na.rm=TRUE)
## [1] 3.528171
sd(x$線寬,na.rm=TRUE)
## [1] 0.04130594

問題

  • 根據summary()的輸出回答以下問題
    • 請問能量的全距
    • 請問線寬的四分位距
#使用apply()
apply(x,2,sd)
##       能量       線寬 
## 3.52817109 0.04130594

不良率

lsl<-2.0;usl<-2.2 #規格界線
idx<-which(x$線寬<lsl |x$線寬>usl) # 找出超出規格界線的線寬資料
loss<-length(idx)/length(x$線寬)*100
paste("不良率=", loss ,"%",sep="")
## [1] "不良率=0.606060606060606%"

練習

  • 若規格界限為[2.05,2.15], 請問
    • 超出規格上界的比率?
    • 低於規格下界的比率?
    • 不良率為多少?

常態假設下的不良率估計

mx<-mean(x$線寬);sdx<-sd(x$線寬)
pp<-pnorm(lsl,mx,sdx)+1-pnorm(usl,mx,sdx)
loss<-pp*100
paste("不良率=", loss ,"%",sep="")
## [1] "不良率=2.56083732531527%"

練習

  • 若規格界限為[2.05,2.15], 在常態假設下, 請問
    • 超出規格上界的比率?
    • 低於規格下界的比率?
    • 不良率為多少?

視覺化資料呈現

直方圖

par(mfrow=c(1,2))
brs1<-seq(min(x$能量)-1,max(x$能量)+1,length.out=25) # 能量變數之X軸切割點
brs2<-seq(min(x$線寬)-0.01,max(x$線寬)+0.01,length.out=25) # 線寬變數之X軸切割點
ht1<-hist(x$能量,breaks=brs1,main="能量的直方圖")
lw1<-lowess(ht1$mids, ht1$counts,f=0.3) #平滑函數
lines(lw1$x,lw1$y,lwd=2,col="blue") #繪製平滑曲線
ht2<-hist(x$線寬,breaks=brs2,main="線寬的直方圖")
lines(c(lsl,lsl),c(0,15),col="green4",lwd=2,lty=2)#規格下界
lines(c(usl,usl),c(0,15),col="green4",lwd=2,lty=2)#規格上界
lw2<-lowess(ht2$mids, ht2$counts,f=0.3)
lines(lw2$x,lw2$y,lwd=2,col="blue")

  • seq(a, b, length.out=c): 將 a 到 b的區間, 切隔成 c-1 段, 包含a, b兩點, 共有c+1個等差序列值。

  • lowess() 是一個資料平滑函數, 將波動大的數據取臨邊資料加權平均, 可以使波動較為平順, 目的是要將主要結構呈現出來

  • lowess()中的f值是平滑取樣數, 指取多少比例的資料來平滑, f值越大取樣筆數越多, 結果越平滑, 但可能失真, 不具意義。

  • 由平滑線看出這個變數的資料都不像常態分配

練習

  • 產生以下序列值
    • 從0到200,等分成20段
    • 從0.1 到 0.8, 產生間隔0.1的序列值

練習

  • 將組距改為0.0025,規格界線改為2.05, 2.15重繪線寬的直方圖。

常態平滑曲線

par(mfrow=c(1,2))
brs1<-seq(min(x$能量)-1,max(x$能量)+1,length.out=25)
brs2<-seq(min(x$線寬)-0.01,max(x$線寬)+0.01,length.out=25)
ht1<-hist(x$能量,breaks=brs1,main="能量的直方圖")
lw1<-lowess(ht1$mids, ht1$counts,f=0.3) #平滑函數
lines(lw1$x,lw1$y,lwd=2,col="blue") #繪製平滑曲線
## 常態曲線
mx<-mean(x$能量)
sdx<-sd(x$能量)
stepx<-ht1$mids[2]-ht1$mids[1] #組距
cut1<-seq(min(x$能量)-1,max(x$能量)+1,length.out=250)
norm.curve<-dnorm(cut1,mx,sdx)
norm.curve.adj<-norm.curve*length(x$能量)*stepx
lines(cut1,norm.curve.adj,col="red",lwd=2)

ht2<-hist(x$線寬,breaks=brs2,main="線寬的直方圖")
lines(c(lsl,lsl),c(0,15),col="green4",lwd=2,lty=2)#規格下界
lines(c(usl,usl),c(0,15),col="green4",lwd=2,lty=2)#規格上界
lw2<-lowess(ht2$mids, ht2$counts,f=0.3)
lines(lw2$x,lw2$y,lwd=2,col="blue")

## 常態曲線
my<-mean(x$線寬)
sdy<-sd(x$線寬)
stepy<-ht2$mids[2]-ht2$mids[1]
cut2<-seq(min(x$線寬)-0.01,max(x$線寬)+0.01,length.out=250)
norm.curve<-dnorm(cut2,my,sdy)
norm.curve.adj<-norm.curve*length(x$線寬)*stepy
lines(cut2,norm.curve.adj,col="red",lwd=2)

  • 紅色線為假設資料為常態分配下, 以資料的平均值和標準差推估它的常態曲線

  • 若紅色與藍色貼合表示資料像常態分配, 視覺上, 這兩個變數資料都不像常態分配

盒形圖

par(mfrow=c(2,1),mar=c(3,3,3,3))
boxplot(x$能量,horizontal=TRUE,col="blue",main="能量的盒形圖")
boxplot(x$線寬,horizontal=TRUE,col="purple",main="線寬的盒形圖")

  • 由盒形圖可以看出能量分配為右偏, 有一個異常值,而線寬似乎是對稱的, 但有兩個異常值。
  • 視覺上, 線寬比能量更像常態分配
  • 線寬是製程產出測量值, 而能量為製程輸入變數, 通常有人為操控的因素。

趨勢圖(假設資料是按時間先後排列)

par(mfrow=c(1,2))
plot(x$能量,type="b",pch=19,xlab="run",ylab="能量")
plot(x$線寬,type="b",pch=19,xlab="run",ylab="線寬")

- 這兩張圖發現, 能量趨勢圖與線寬趨勢圖似乎非常相像, 都有下降的趨勢。 - 可以使用lm() 函數驗證是否有線性下降或上升趨勢, lm 是linear model的簡寫

par(mfrow=c(1,2))
runs<-1:length(x$能量)
plot(x$能量,type="b",pch=19,xlab="run",ylab="能量")
out1<-lm(x$能量~runs) #估算線性模式
abline(out1,lwd=2,col="blue")#加入線性趨勢線
plot(x$線寬,type="b",pch=19,xlab="run",ylab="線寬")
out2<-lm(x$線寬~runs)#估算線性模式
abline(out2,lwd=2,col="blue")#加入線性趨勢線

  • 查看lm()函數傳出的數值與意義
out1
## 
## Call:
## lm(formula = x$能量 ~ runs)
## 
## Coefficients:
## (Intercept)         runs  
##    31.66423     -0.05444

  • 這表示能量與批次的線性模型為 \[能量=31.66423-0.0544\times runs\]

  • 截距項為31.66423可以解釋為起始能量值, 斜率為 -0.0544表示能量隨批次的變化量

  • 使用summary()函數可以查看模型的合適性

summary(out1)
## 
## Call:
## lm(formula = x$能量 ~ runs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4491 -1.4952 -0.1273  1.4357  6.0367 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 31.664226   0.374048   84.65   <2e-16 ***
## runs        -0.054436   0.003909  -13.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.391 on 163 degrees of freedom
## Multiple R-squared:  0.5434, Adjusted R-squared:  0.5406 
## F-statistic:   194 on 1 and 163 DF,  p-value: < 2.2e-16
  • 查看模型的幾項指標
    • F-statistic 的p-value
    • Adjusted R-square
    • 係數表(Coefficients)斜率項 (intercept以外)的p-value 值
  • Residuals 是殘差值: 實際值與預測值的差距
  • 任何估計值的殘差值的總和必定為零, 殘差值的變異越小, 表示估計值與實際值越貼近。
  • 先看F-statistic 的p-value是否小於0.05, 若是表示這個模型是顯著的, 其意是模型的殘差變異明顯比實際值的變異小, 統計上是指用這個模型來解釋資料是有意義的, 但不是說是最好的解釋, 或是最好模型。
  • F-statistic: 194 on 1 and 163 DF, p-value: < 2.2e-16中,顯示 p-value小於0.05, 故說明模型有解釋資料的意義
  • 再看Coefficients的Pr(>|t|) 的機率值 (這也是p-value)是否小於0.05, 若是表示接受係數不為零。
  • Adjusted R-squared=0.5406, 模型可以解釋約54%的線寬變異。
  • Residual standard error: 殘差的標準誤
sd(out1$residuals);sd(out1$residuals)*sqrt(164/163)
## [1] 2.384155
## [1] 2.391457
  • 殘差標準誤的計算公式為 \[standard\_error=\sqrt{\frac{ 殘差平方總和}{d.f}}=\sqrt{\frac{\sum residuals^2}{d.f}}\]

-上述報告中有``Residual standard error: 2.391 on 163 degrees of freedom", 自由度(degrees of freedom,d.f.)為163, 等於\(n-2\) (減2是因為有兩個參數截距和斜率需要估計)

  • 標準差的計算公式 \[ s=\sqrt{\frac{ 殘差平方總和}{n-1}}=\sqrt{\frac{\sum residuals^2}{n-1}}\]

所以由標準差要修正成標準誤, 必須乘上係數\(\sqrt{164/163}\)修正。當 \(n\) 很大, \(\sqrt{164/163}\)接近1, 是否修正差異不大。

殘差的分配

直方圖, 盒形圖, 常態機率繪圖

par(mfrow=c(1,3))
ht3<-hist(out1$residuals,nclass=25)
## 常態曲線
mr<-mean(out1$residuals)
sr<-sd(out1$residuals)
stepr<-ht3$mids[2]-ht3$mids[1]
cut3<-seq(min(out1$residuals)-0.01,max(out1$residuals)+0.01,length.out=250)
norm.curve<-dnorm(cut3,mr,sr)
norm.curve.adj<-norm.curve*length(out1$residuals)*stepr
lines(cut3,norm.curve.adj,col="red",lwd=2)

boxplot(out1$residuals) #盒形圖
qqnorm(out1$residuals);qqline(out1$residuals,col="blue",lwd=2) #常態機率繪圖

\[H_0: 係數為0, H_1: 係數不為0\] - 截距(Intercept)和斜率(runs)的p-value都小於0.05, 故說模型中截距=31.664226 , 斜率=-0.05444都是顯著的

out2
## 
## Call:
## lm(formula = x$線寬 ~ runs)
## 
## Coefficients:
## (Intercept)         runs  
##   2.1272766   -0.0005445
  • 能量逐漸下降, 線寬也是
  • 截距項為2.1272766 可以解釋為起始線寬估計值, 斜率為 -0.0005445表示線寬隨批次的變化量 \[能量=2.1272766-0.0005445\times runs\]
summary(out2)
## 
## Call:
## lm(formula = x$線寬 ~ runs)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.093436 -0.019111 -0.000157  0.020606  0.096302 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.1272766  0.0050336  422.62   <2e-16 ***
## runs        -0.0005445  0.0000526  -10.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03218 on 163 degrees of freedom
## Multiple R-squared:  0.3967, Adjusted R-squared:  0.393 
## F-statistic: 107.2 on 1 and 163 DF,  p-value: < 2.2e-16

練習

  • 根據summary(out2)的結果, 請回答:
    • 這個模型的p-value是否小於0.05, 模型顯著嗎?
    • 這個模型對能量變異的解釋比例有多少?
    • 請問估計能量逐批減少量? 是否顯著?
    • 殘差的標準誤(standard error) 為多少?
plot(x$能量,x$線寬,pch=19,type="p")

plot(x$能量,x$線寬,pch=19,type="p")
out3<-lm(x$線寬~x$能量) #線寬與能量間的線性模型估計
abline(out3,lwd=2,col="blue")

summary(out3)
## 
## Call:
## lm(formula = x$線寬 ~ x$能量)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.051184 -0.013388  0.000204  0.012988  0.058739 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.8082108  0.0127352  141.99   <2e-16 ***
## x$能量      0.0100887  0.0004652   21.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02102 on 163 degrees of freedom
## Multiple R-squared:  0.7426, Adjusted R-squared:  0.741 
## F-statistic: 470.2 on 1 and 163 DF,  p-value: < 2.2e-16
  • 查看模型的幾項指標
    • F-statistic 的p-value
    • Adjusted R-square
    • 係數表(Coefficients)斜率項 (intercept以外)的p-value 值
  • 先看F-statistic 的p-value是否小於0.05, 若是表示這個模型是顯著的, 其意是模型的殘差變異明顯比實際值的變異小, 統計上是指用這個模型來解釋資料是有意義的, 但不是說是最好的解釋, 或是最好模型。
  • F-statistic: 470.2 on 1 and 163 DF, p-value: < 2.2e-16中,顯示 p-value小於0.05, 故說明模型有解釋資料的意義
  • 再看Coefficients的Pr(>|t|) 的機率值 (這也是p-value)是否小於0.05, 若是表示接受係數不為零。 \[線寬=1.8082108+0.0100887\times 能量\] 增加一單位能量可以使先寬增加約0.01。
  • Adjusted R-squared=0.741, 模型中能量變化可以解釋約74.1%的線寬變異。

殘差分析

殘差的趨勢圖

  • 主要查看模型殘差(residuals)是否跟時間有關, 有無特殊型態
  • 可以使用一些統計檢定方法驗證, 但超出課程範圍略過
  • 用視覺化方式, 看看是否有無特殊形態?
plot(out3$residuals,xlab="runs",pch=19,col="blue",type="b")

  • 一些特殊形態
par(mfrow=c(2,2),mar=c(3,3,3,3))
runs<-1:100
y<-runs+rnorm(100,0,10)-50.5;
#線性
plot(runs, y,ylab="residuals",main="線性") 
lines(runs,runs-50.5,lty=2,col="blue",lwd=2)
#週期波動
plot(runs, sin(runs/10)+rnorm(100,0,1),ylab="residuals",main="週期波動")
lines(runs,sin(runs/10),lty=2,col="blue",lwd=2)
#分層
plot(runs,c(rnorm(50,2,1),c(rnorm(50,-2,1))),ylab="residuals",main="分層")
lines(runs,c(rep(2,50),rep(-2,50)),lty=2,col="blue",lwd=2)
y<-rnorm(100,1,abs(runs-50.5))
#變異不均
plot(runs, y,ylab="residuals",main="變異不均")
lines(runs,rep(0,100),lty=2,col="blue",lwd=2)

par(mfrow=c(2,2))
ht5<-hist(out3$residuals,nclass=25)
## 常態曲線
mr<-mean(out3$residuals)
sr<-sd(out3$residuals)
stepr<-ht5$mids[2]-ht5$mids[1]
cut5<-seq(min(out3$residuals)-0.005,max(out3$residuals)+0.005,length.out=250)
norm.curve<-dnorm(cut5,mr,sr)
norm.curve.adj<-norm.curve*length(out3$residuals)*stepr
lines(cut5,norm.curve.adj,col="red",lwd=2)

boxplot(out3$residuals,horizontal=TRUE)
qqnorm(out3$residuals);qqline(out3$residuals,col="blue",lwd=2)
plot(out3$fitted, out3$residuals,pch=19,main="fitted vs residuals")

  • 根據殘差分析,我們接受殘差服從常態分配, 且變異無特殊型態(iid)

  • 我們建立線寬與能量的模型為 \[線寬=1.8082108+0.0100887\times 能量\] 增加一單位能量可以使先寬增加約0.01。

  • 討論:

    • 如果欲使線寬為2.1, 請問能量應控制為多少?
    • 假設目前能量設定為\(x_c\), 線寬為2.2, 為了使線寬減少0.1變為2.1, 請問能量應該增加或減少, 變化量為多少?

練習

  • 美化殘差分析的四個圖形
    • 改變直方圖的組距和組數, 將常態曲線改為藍色
    • 盒形圖的方向改為直立, 並使用紫色外框
    • 常態機率繪圖的點改用pch=19
    • 配適值與殘差值圖(fitted vs residuals)的顏色改為綠色