2
$\begingroup$

I would like to know how experts deal with real data. Even if statistical text books uses real data I'm always surprised how good the real data are and at the end of the exercises the residuals are normal distributed and have constant variances.

I produced here artificial data based on a simple linear model $y\sim x$ which I treated to get fake real data. One treatment is trimming. Below I wrote the code of the treatment. Let's assume we have the data set and we want to see if variable $x$ has an effect on $y$. $x$ is age between 20 and 80 years and $y$ is a continuous variable. Here are some information about this dataset data1:

> head(data1)
   x         y
1 20  2983.026
2 21  6026.804
3 22 48976.006
4 23 10948.566
5 24 12637.740
6 25 62337.722


> dim(data1)
[1] 541   2

> summary(data1)
       x               y         
 Min.   :20.00   Min.   :  1561  
 1st Qu.:35.00   1st Qu.: 17578  
 Median :53.00   Median : 44783  
 Mean   :51.11   Mean   : 63214  
 3rd Qu.:67.00   3rd Qu.: 89139  
 Max.   :80.00   Max.   :410210


par(mfrow=c(1,2))
plot(y ~ x, data=data1)
hist(data1$y)

figure 1

Next step is to model this data. First I tried a simple regression:

> fit1   <- lm(y ~ x, data=data1)
> summary(fit1)

Call:
lm(formula = y ~ x, data = data1)

Residuals:
   Min     1Q Median     3Q    Max 
-99053 -30495  -7793  15710 303476 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -44813       6560  -6.831 2.28e-11 ***
x               2114        121  17.469  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Residual standard error: 50950 on 539 degrees of freedom
Multiple R-squared:  0.3615,    Adjusted R-squared:  0.3603 
F-statistic: 305.2 on 1 and 539 DF,  p-value: < 2.2e-16

Residual plots:

f_ra(fit1)  # function to display Tukey-Anscomb-plot and qqnorm-plot (see below)

figure 2

As already figure 1 showed and definitly the residuals analyses in figure 2 confirms the data has to be transformed. I wouldn't do that if I would need only the estimates but since I need p-value respectively confidence intervals the residuals should satisfy equal variance at each x and normality. I tried to take the log and the cube root of $y$ and decided for the cube root. Here are the results:

> fit2 <- lm(y^(1/3) ~ x, data=data1)
> summary(fit2)

Call:
lm(formula = y^(1/3) ~ x, data = data1)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.033  -6.070  -0.694   5.843  32.475 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.08906    1.18977    9.32   <2e-16 ***
x            0.48108    0.02194   21.93   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.24 on 539 degrees of freedom
Multiple R-squared:  0.4714,    Adjusted R-squared:  0.4705 
F-statistic: 480.7 on 1 and 539 DF,  p-value: < 2.2e-16

And here again the residual plots:

f_ra(fit2)

figure3

I'm very happy with this solution. Nevertheless, the Tukey-Anscomb-plot (residuals against fitted) shows some structure (V-notch at x=35) and also the trim is visible. Also the residuals seems to be still right tailed.

Is there a way to cope with this V-notch or I don't have to care?

Is the qqnorm-plot good enough or should I try to find a better fit?

I could use bootstrap but if the data is very big and the model is much more complex (e.g. mixed models) so that one computation (one regression) takes severel minitues bootstrap wouldn't feasible.

Here is the code of the simulated data:

# Number of measurements for age x
n  <- 10
# error
sd <- 10
# coefficients
beta0 <- 10
beta1 <- 0.5

# data
x <- rep(20:80,n)
y <- beta0 + beta1*x + rnorm(length(x),0,sd)

# ideal data
data0 <- data.frame(x=x,y=y)

# make real data
data1 <- data0
# - make skewed (due to 
data1$y <- data1$y^3
# - trim (due to threshold)
data1 <- subset(data1, y > 1500)
# - a weird lacking of data
data1 <- subset(data1, !(x>40 & x <= 50 & y > 40000))

# function to plot Tukey-Anscomb-plot and qqnorm-plot
f_ra <- function(xfit) {
    plot(fitted(xfit), resid(xfit))
    abline(h=0, lty=2, col="red")
    qqnorm(resid(xfit))
    qqline(resid(xfit), lty=2, col="red")
}
$\endgroup$
  • 2
    $\begingroup$ NB: your "trimming" is usually known as (lower) truncation. The code shows you did more than that: you selectively changed the response in the middle of the $x$ range. This is an unusual (and strong) violation of many of the implicit assumptions in your models: linearity, homoscedasticity, and common distributional shape of the response. If you made that change to study something you suspect may be happening in real data, perhaps you want to look into changepoint analysis. $\endgroup$ – whuber May 2 '16 at 18:34
  • $\begingroup$ @whuber Thanks for the correction. Indeed, trimming was only one of the violation. This is the crude reality with, for example, epidemiological data: there is may be a threshold which produce truncated data, there are inclompete parts etc. I would like to know if there are somekind like a policy or technicques how to cope with this. Thanks for the link to changepoint analysis. This does not answer my posted question but another one I have in my rucksack. $\endgroup$ – giordano May 2 '16 at 21:44
  • 1
    $\begingroup$ You can try to model what you actually think is going on with the data rather than bash a square peg into a round hole (i.e. force it into an ordinary regression). $\endgroup$ – Glen_b May 3 '16 at 0:44
  • $\begingroup$ @Glen_b I'm with you. Let's assume I thought about that but couldn't find nothing and the only I have are the data as they are. How deal with a square peg? How to make inference (confidence intervals) based only on the data without pre-knowledge? $\endgroup$ – giordano May 3 '16 at 6:23
  • $\begingroup$ If the data look like that, and the data are all I have to go on (presumably I'm doing some form of sample splitting to pick my model and again to evaluate, fitting on the remainder), then why would I choose a regression model rather than something else? The conditional distribution looks skewed with nonconstant variance and the relationship doesn't look straight. Why wouldn't I look to a GLM or GAM type model rather than an obviously inadequate one? $\endgroup$ – Glen_b May 3 '16 at 8:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.