শিক্ষক.কমের R পরিচিতি কোর্সের শেষ লেকচারে আমরা R ব্যবহার করে কীভাবে Regression analysis ও সেটা করার পর Regression modelএর diagnostics করা যায় সেটা নিয়ে আলোচনা করবো।
[কোর্সের সব লেকচার দেখতে হলে কোর্সের মূল পাতায় দেখুন]
রিগ্রেশন এ্যানালাইসিস
এই লেকচারে আমরা যে ড্যাটাটি নিয়ে কাজ করবো সেটা হলো state.x77। এটি R এর একটি ডিফল্ট ড্যাটাসেট। ড্যাটাসেটটি লোড করতে আমরা একটি নতুন data.frame তৈরি করে সেটা attach করে নেব।
data <- data.frame(state.x77) # attaching the dataset attach(data)
ড্যাটাটির ভ্যারিয়েবলগুলোর দিকে একটু নজর বুলোতে আমরা str ও summary কমান্ডটি ব্যবহার করতে পারি।
# summerize the dataset str(data) summary(data)
নিউমেরিক সামারি দেখার সাথে সাথে আমরা গ্রাফিকালি ড্যাটাটির ওপর চোখ বোলাতে পারি scatter plot matrix ও correlation matrix plot এর ওপর।
# Checking the data # scatter plot matrix pairs (data) # correlation matrix library("corrplot") matr <- cor(data) corrplot(matr, method = "circle")
আমরা একটা মডেল অনুমান করে নিতে পারি যেখানে Life.Exp ডিপেনডেন্ট ভ্যারিয়েবল ও Population, Income, Illiteracy, Murder, HS.Grad, Frost এগুলো ইনডিপেনডেন্ট ভ্যারিয়েবল হিসেবে ধরে নিতে পারি। এই মডেলটি রান করতে আমাদের lm কমান্ডটি ব্যবহার করতে হবে।
# regression model regressmod <- lm (Life.Exp ~ Population + Income + Illiteracy + Murder + HS.Grad + Frost)
এখানে মডেলটি regressmod নামে একটি অবজেক্টের ভেতরে রাখা হয়েছে। এর সামারি দেখতে আমরা summary কমান্ডটি ব্যবহার করবো।
summary (regressmod)
রিগ্রেশন ডায়াগনস্টিক্স
রিগ্রেশন মডেলটি রান করার পর আমরা সামারিতে R square, adjusted R square, F value, p-value (F-value) এগুলো পেয়েছি।
এখন এই মডেলটি autocorrelation ও multicollinearity-মুক্ত কিনা সেটা দেখতে আমরা car প্যাকেজটি লোড করে যথাক্রমে dwt ও vif কমান্ডদুটি ব্যবহার করবো।
# Regression diagnostics # Autocorrelation test: Dubin Watson d # The d should be around 2 library(car) dwt(regressmod) # Multicollinearity test # Variance inflation factor # library(car) vif (regressmod)
Heteroscedasticityর প্রভাব রয়েছে কিনা সেটা দেখতে আমরা Breusch-Pagan টেস্ট করে দেখতে পারি। এটা করতে আমরা lmtest প্যাকেজটি ব্যবহার করে bptest কমান্ডটি ব্যবহার করতে পারি।
# Breusch Pagan Heteroscedasticity test # H_0: Homoscedasticity. library(lmtest) bptest(regressmod)
এখন এই মডেলের residual থেকে আমরা সেটা normally distributed কিনা সেটা দেখতে পাবো। এটা করতে আমরা প্রথমে রেসিডুয়ালগুলো বের করে নেব। এটা করতে একটা নতুন অবজেক্ট তৈরি করবো আমরা resid কমান্ডটি ব্যবহার করে।
# Normality test # Generate a residual variable data$residuals <- resid(regressmod)
Residualগুলো normally distributed কিনা সেটা দেখতে আমরা গ্রাফিকাল টুল যেমন Q-Q plot বা histogramএর সাহায্য নিতে পারি। এটা করতে আমরা প্রথমে প্লট এরিয়াটিকে একটি সারি ও দুইটি কলামে বিভক্ত করে নেই।
par(mfrow=c(1,2))
এরপর Q-Q plot ও histogram plot করতে নিন্মোক্ত কমান্ডগুলোর সাহায্য নিতে পারি।
qqnorm(data$residuals) qqline(data$residuals) hist(data$residuals, breaks=100, xlab="Residuals", ylab="Frequency", main="Histogram of Residuals") dev.off()
গ্রাফিকাল পদ্ধতিতে ফলাফল পরিস্কার বোঝা না গেলে আমরা Shapiro Wilk test এর মাধ্যমে Residual সিরিজটা normally distributed কিনা সেটা জানতে পারি।
# H_0: data is normally distributed. p-value more than 0.05 indicates acceptance of null at 5% level shapiro.test(data$residuals)
মিসিং ভ্যালু ও নমিনাল ও ক্যাটেগরিকাল ড্যাটা কোডিং
Rএ কোন ড্যাটা ইমপোর্ট করলে সেটাকে প্রথমে সঠিক কোডিং করে নিতে হয়। প্রথমত, খেয়াল রাখা দরকার যে ড্যাটাতে কোন missing value আছে কিনা। সাধারণত ড্যাটা ইনপুটের সময় missing value চিহ্নিত করতে -99 বা -999 লেখা হয়। অনেক সময় এটা 0 হিসেবেও ড্যাটাতে থাকে। R-এর missing value চিহ্নিত করতে এগুলোকে NA হিসেবে কনভার্ট করে নিতে হবে।
এই উদাহরণটিতে মিসিং ভ্যালু কোড না করে রিগ্রেশন রান করা হয়েছে।
# Coding missing values x <- c(12, 14, 11, 9, 8, 15, 18, 25, 33, 12) y <- c(134, 321, 231, 222, 114, 144, -99, 99, 128, 222) regressmod1 <- lm(y ~ x, na.action=na.omit) summary(regressmod1)
এখানে দেখছি আমাদের R-square ১১% এর মতো।
এখন আমরা মিসিং ভ্যালুকে NA দিয়ে প্রতিস্থাপন করবো।
# Replacing -99 with NA is.na(x) is.na(y) y[y==-99] <- NA is.na(y) regressmod1 <- lm(y ~ x)
এভাবে missing value গুলো চিহ্নিত করার পর regression model রান করার সময় এই ভ্যালুগুলো বাদ দেয়া যায় na.action=na.omit অপশনটি ব্যবহার করে।
এখন আমরা রিগ্রেশন মডেল রান করালে আগের চাইতে ভিন্ন একটা আউটপুট পাবো।
# Dealing with the missing values in regression regressmod1 <- lm(y ~ x, na.action=na.omit)
ক্যাটাগরিকাল/নমিনাল/ডামি ভ্যারিয়েবল কোডিং
Rএ কাজ করার সময় discrete ড্যাটাগুলোও সঠিকভাবে কোডিং নেয়া প্রয়োজনীয়। এখানে অনলাইন থেকে একটি ড্যাটা ডাউনলোড করে আমরা দেখবো কীভাবে একটি ক্যাটাগরিকাল ভ্যারিয়েবলকে কোডিং করা যায়। ড্যাটাটির ভ্যারিয়েবল লিস্ট ও সেগুলোর প্রকার দেখতে প্রথমেই আমরা ডাউনলোডেড অবজেক্টটিতে str কমান্ড প্রয়োগ করে দেখি।
# Coding categorial variables data1 <- read.csv("http://www.ats.ucla.edu/stat/data/hsb2.csv") str(data1)
R-এ ক্যাটাগরিকাল/নমিনাল/ডামি ভ্যারিয়েবলগুলোকে factor ভ্যারিয়েবল বলা হয়। ডাউনলোড করা ড্যাটাটিতে race ভ্যারিয়েবলটিতে চারটি ক্যাটাগরি চিহ্নিত করতে আমাদের factor কমান্ডটির সাহায্য নিতে হবে।
data1$race.fact <- factor(data1$race, labels=c("Hispanic", "Asian", "African-Am", "Caucasian"))
এই কমান্ডটি যথাযথভাবে কাজ করলো কিনা সেটা দেখতে, is.factor কমান্ডটি ব্যবহার করা যেতে পারে। যেমন,
is.factor(race.fact)