Sunday, March 1, 2015

Generalized, linear, and generalized least squares models (LM, GLM, GLS) using R Language

#lm linear regression, normal error, constant variance
     Y = a + bX + E  a Linear Predictor
#glm generalize linear model, non-normal error, non-constant variance
       LogY = a + bX + E  
       Y = e^a*e^bX + E   a Multiplicater, exponential and Logthermic
     in glm, individual slope gives an estimate of multiplitive change
     in the reponse variable for one unit change in corresponding explanatory variable
#gls: generalise least square model, correlated error, spatial, temperal/pattern/trends

airquality

plot(Ozone~Wind,airquality)

model1=lm(Ozone~Wind,airquality)
plot(model1)
coef(model1)

#prediction for Wind speed at 19 and 20 mph

coef(model1)[1]
Ozone1=coef(model1)[1]+coef(model1)[2]*19
Ozone2=coef(model1)[1]+coef(model1)[2]*20

Ozone1
Ozone2

##poisson is generalize linear model

model2=glm(Ozone~Wind,airquality,family=poisson)
glm(model2)

# Coefficients:
#(Intercept)         Wind
  #   96.873       -5.551

Ozone1.glm=exp(coef(model2)[1]+coef(model2)[2])*19
Ozone2.glm=exp(coef(model2)[1]+coef(model2)[2])*20
Ozone1.glm
Ozone2.glm

plot(Ozone~Wind,airquality)

Ozone1.glm/Ozone2.glm

# 0.95

exp(coef(model2))[2]  #exp(-5.551)

###gls

library(nlme)

model3.gls=gls(Ozone~Wind,airquality)

model3=gls(Ozone~Wind,airquality,na.action=na.exclude)
head(airquality)
?airquality

paste(1973,airquality$Month,airquality$Day, sep=",")
as.Date(paste(1973,airquality$Month,airquality$Day, sep=","))

airquality$Date
paste(1973,airquality$Month,airquality$Day, sep=",")

library(lattice)
xyplot(Ozone~Date,airquality)

model4=gls(Ozone~Wind*Date,airquality,na.action=na.exclude)
air2=subset(airquality, complete.cases(Ozone))

model4=gls(Ozone~Wind*Date,air2)
plot(ACF(model5=~Date),alpha=0.5)

model6=(update(model5,correlation=corAR1())

library(MuMIn)

AICc(model5,model6)

summary(model6)

No comments:

Post a Comment