To part 9a and 9b

From Wiki1

Jump to: navigation, search

Question 9

Part a and b

Take a simple random sample of 1/5 of the radon data using sample(). Fit the varying-intercept model with floor as an individual-level predictor and log uranium as a county level predcitor, and compare your inferences to the model using the entire dataset. (Method was based on: http://www.unc.edu/courses/2007spring/enst/562/001/docs/solutions/assign10.htm)

Repeat step a with a different random sample each time, summarize how the estimates vary.

Level 1:

  logradon = β0i + β1i Floorij + εij

Level 2:

  β0i = β0 + β2log(Uppm)i + u0i 
  
β1i = β1

Composite Equation:

   logradon = β0 + β2log(Uppm)i + u0i + β1 Floorij + εij

Where:

εij ~ N(0, σ2) and u0i ~ N(0, τ2)

R-Code:

  Call up the nmle package: library("nlme", lib.loc="C:/Program Files/R/R-3.0.1/library")

We use "lme" to represent our model using the entire dataset:


> lme(fixed=logradon~1+floor+log(Uppm), random=~1|county, data=merge3, method='ML')->out3
> summary(out3)
Linear mixed-effects model fit by maximum likelihood
 Data: merge3 
       AIC      BIC    logLik
  2157.913 2182.029 -1073.956

Random effects:
 Formula: ~1 | county
        (Intercept)  Residual
StdDev:   0.1457421 0.7689218

Fixed effects: logradon ~ 1 + floor + log(Uppm) 
                 Value  Std.Error  DF  t-value p-value
(Intercept)  1.4627406 0.03756813 833 38.93568       0
floor       -0.6788426 0.06969627 833 -9.74001       0
log(Uppm)    0.7180305 0.09064038  83  7.92175       0
 Correlation: 
          (Intr) floor 
floor     -0.362       
log(Uppm)  0.156 -0.011

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-5.12860182 -0.60172541  0.03181403  0.64952352  3.32903089 

Number of Observations: 919
Number of Groups: 85 

 

We use the sample() function to take 1/5 of our 919 observations randomly, and then use lme to model our data.


> mysample1<-merge3[sample(1:nrow(merge3),919/5,replace=FALSE),]

> lme(fixed=logradon~1+floor+log(Uppm), random=~1|county, data=mysample, method='ML')->out4
> summary(out4)
Linear mixed-effects model fit by maximum likelihood
 Data: mysample 
       AIC      BIC    logLik
  408.4154 424.4629 -199.2077

Random effects:
 Formula: ~1 | county
         (Intercept)  Residual
StdDev: 4.232921e-05 0.7186565

Fixed effects: logradon ~ 1 + floor + log(Uppm) 
                 Value  Std.Error  DF  t-value p-value
(Intercept)  1.4354261 0.06095193 123 23.55013  0.0000
floor       -0.4986576 0.14886944 123 -3.34963  0.0011
log(Uppm)    0.9814369 0.15357413  57  6.39064  0.0000
 Correlation: 
          (Intr) floor 
floor     -0.383       
log(Uppm)  0.297 -0.031

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.56246393 -0.68657626  0.06167267  0.73418155  2.48527625 

Number of Observations: 183
Number of Groups: 59 

> mysample2<-merge3[sample(1:nrow(merge3),919/5,replace=FALSE),]
> lme(fixed=logradon~1+floor+log(Uppm), random=~1|county, data=mysample2, method='ML')->out5
> summary(out5)
Linear mixed-effects model fit by maximum likelihood
 Data: mysample2 
       AIC      BIC    logLik
  425.8058 441.8532 -207.9029

Random effects:
 Formula: ~1 | county
        (Intercept)  Residual
StdDev:   0.2595507 0.7173485

Fixed effects: logradon ~ 1 + floor + log(Uppm) 
                 Value  Std.Error  DF   t-value p-value
(Intercept)  1.4527070 0.07735835 122 18.778929   0e+00
floor       -0.4510375 0.15380505 122 -2.932527   4e-03
log(Uppm)    0.6644033 0.18017837  58  3.687475   5e-04
 Correlation: 
          (Intr) floor 
floor     -0.397       
log(Uppm)  0.259 -0.101

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-2.884649994 -0.557167047 -0.006331528  0.693333489  3.147828643 

Number of Observations: 183
Number of Groups: 60 

> mysample3<-merge3[sample(1:nrow(merge3),919/5,replace=FALSE),]
> lme(fixed=logradon~1+floor+log(Uppm), random=~1|county, data=mysample2, method='ML')->out6
> summary(out6)
Linear mixed-effects model fit by maximum likelihood
 Data: mysample2 
       AIC      BIC    logLik
  425.8058 441.8532 -207.9029

Random effects:
 Formula: ~1 | county
        (Intercept)  Residual
StdDev:   0.2595507 0.7173485

Fixed effects: logradon ~ 1 + floor + log(Uppm) 
                 Value  Std.Error  DF   t-value p-value
(Intercept)  1.4527070 0.07735835 122 18.778929   0e+00
floor       -0.4510375 0.15380505 122 -2.932527   4e-03
log(Uppm)    0.6644033 0.18017837  58  3.687475   5e-04
 Correlation: 
          (Intr) floor 
floor     -0.397       
log(Uppm)  0.259 -0.101

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-2.884649994 -0.557167047 -0.006331528  0.693333489  3.147828643 

Number of Observations: 183
Number of Groups: 60 

<div>
The confidence intervals for the estimates overlap between the sampled models and the model with the entire dataset (or are very very close).
Personal tools