# Loading wages data
wages <- read.csv("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/wages_pp-1.txt",head=TRUE)
head(wages)
## id lnw exper ged postexp black hispanic hgc hgc.9 uerate ue.7
## 1 31 1.491 0.015 1 0.015 0 1 8 -1 3.215 -3.785
## 2 31 1.433 0.715 1 0.715 0 1 8 -1 3.215 -3.785
## 3 31 1.469 1.734 1 1.734 0 1 8 -1 3.215 -3.785
## 4 31 1.749 2.773 1 2.773 0 1 8 -1 3.295 -3.705
## 5 31 1.931 3.927 1 3.927 0 1 8 -1 2.895 -4.105
## 6 31 1.709 4.946 1 4.946 0 1 8 -1 2.495 -4.505
## ue.centert1 ue.mean ue.person.cen ue1
## 1 0.00 3.215 0.00 3.215
## 2 0.00 3.215 0.00 3.215
## 3 0.00 3.215 0.00 3.215
## 4 0.08 3.215 0.08 3.215
## 5 -0.32 3.215 -0.32 3.215
## 6 -0.72 3.215 -0.72 3.215
summary(wages)
## id lnw exper ged
## Min. : 31 Min. :0.708 Min. : 0.001 Min. :0.0000
## 1st Qu.: 3194 1st Qu.:1.591 1st Qu.: 1.609 1st Qu.:0.0000
## Median : 6582 Median :1.842 Median : 3.451 Median :0.0000
## Mean : 6301 Mean :1.897 Mean : 3.957 Mean :0.2719
## 3rd Qu.: 9300 3rd Qu.:2.140 3rd Qu.: 5.949 3rd Qu.:1.0000
## Max. :12543 Max. :4.304 Max. :12.700 Max. :1.0000
## postexp black hispanic hgc
## Min. : 0.0000 Min. :0.0000 Min. :0.000 Min. : 6.000
## 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.: 8.000
## Median : 0.0000 Median :0.0000 Median :0.000 Median : 9.000
## Mean : 0.9076 Mean :0.2529 Mean :0.241 Mean : 8.948
## 3rd Qu.: 0.1168 3rd Qu.:1.0000 3rd Qu.:0.000 3rd Qu.:10.000
## Max. :12.2600 Max. :1.0000 Max. :1.000 Max. :12.000
## hgc.9 uerate ue.7 ue.centert1
## Min. :-3.00000 Min. : 1.795 Min. :-5.2050 Min. :-15.310
## 1st Qu.:-1.00000 1st Qu.: 5.395 1st Qu.:-1.6050 1st Qu.: -2.905
## Median : 0.00000 Median : 7.000 Median : 0.0000 Median : -0.600
## Mean :-0.05248 Mean : 7.733 Mean : 0.7333 Mean : -1.029
## 3rd Qu.: 1.00000 3rd Qu.: 9.400 3rd Qu.: 2.4000 3rd Qu.: 0.600
## Max. : 3.00000 Max. :23.705 Max. :16.7050 Max. : 13.005
## ue.mean ue.person.cen ue1
## Min. : 2.895 Min. :-8.403 Min. : 2.895
## 1st Qu.: 6.066 1st Qu.:-1.350 1st Qu.: 6.195
## Median : 7.154 Median :-0.142 Median : 8.300
## Mean : 7.733 Mean : 0.000 Mean : 8.762
## 3rd Qu.: 8.654 3rd Qu.: 1.133 3rd Qu.:10.200
## Max. :17.128 Max. :13.421 Max. :23.705
# Creating data to work on using id uerate lmw exper variables
# To do Bivariate linear mixed model need to first melt the data
require(dplyr)
## Loading required package: dplyr
## Warning: package 'dplyr' was built under R version 3.5.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
require(reshape2)
## Loading required package: reshape2
## Loading required package: reshape2
swages <- melt(wages[, c("id", "uerate", "lnw", "exper")],
measure.vars = c("lnw","exper"))
head(swages)
## id uerate variable value
## 1 31 3.215 lnw 1.491
## 2 31 3.215 lnw 1.433
## 3 31 3.215 lnw 1.469
## 4 31 3.295 lnw 1.749
## 5 31 2.895 lnw 1.931
## 6 31 2.495 lnw 1.709
# Creating dummy variables for lnw and exper
swages <- within(swages, {
Dl <- as.integer(variable == "lnw")
De <- as.integer(variable == "exper")
})
head(swages)
## id uerate variable value De Dl
## 1 31 3.215 lnw 1.491 0 1
## 2 31 3.215 lnw 1.433 0 1
## 3 31 3.215 lnw 1.469 0 1
## 4 31 3.295 lnw 1.749 0 1
## 5 31 2.895 lnw 1.931 0 1
## 6 31 2.495 lnw 1.709 0 1
## print the data for variables of interest for id 31
swages[swages$id == 31, ]
## id uerate variable value De Dl
## 1 31 3.215 lnw 1.491 0 1
## 2 31 3.215 lnw 1.433 0 1
## 3 31 3.215 lnw 1.469 0 1
## 4 31 3.295 lnw 1.749 0 1
## 5 31 2.895 lnw 1.931 0 1
## 6 31 2.495 lnw 1.709 0 1
## 7 31 2.595 lnw 2.086 0 1
## 8 31 4.795 lnw 2.129 0 1
## 6403 31 3.215 exper 0.015 1 0
## 6404 31 3.215 exper 0.715 1 0
## 6405 31 3.215 exper 1.734 1 0
## 6406 31 3.295 exper 2.773 1 0
## 6407 31 2.895 exper 3.927 1 0
## 6408 31 2.495 exper 4.946 1 0
## 6409 31 2.595 exper 5.965 1 0
## 6410 31 4.795 exper 6.984 1 0
require(nlme)
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## fit the model
# using the 0 notation to allow for a separate intercept for each outcome indicate
# by the effect of the two dummy variables, De and Dl.
m <- lme(value ~ 0 + Dl + Dl:uerate + De + De:uerate, data = swages,
random = ~0 + Dl + Dl:uerate + De + De:uerate | id,
weights = varIdent(form = ~1 |variable),
control = lmeControl(maxIter = 200, msMaxIter = 200, niterEM = 50,msMaxEval = 400))
require(lattice)
## Loading required package: lattice
## Loading required package: lattice
## create data frame of residuals, fitted values, and variable
diagnos <- data.frame(Resid = resid(m, type = "normalized"), Fitted = fitted(m),Variable = swages$variable)
## overal QQ normal plot
qqmath(~Resid, data = diagnos, distribution = qnorm)
## separate by variable
qqmath(~Resid | Variable, data = diagnos, distribution = qnorm)
## we could add a line to indicate 'normal' with a bit more work
qqmath(~Resid | Variable, data = diagnos, distribution = qnorm, prepanel = prepanel.qqmathline,
panel = function(x) {
panel.qqmathline(x)
panel.qqmath(x)
})
## overall plot
xyplot(Resid ~ Fitted, data = diagnos)
## separate by variable
xyplot(Resid ~ Fitted | Variable, data = diagnos)
