# Linear regression with trimmed data

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)  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)  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)  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") }  • 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. – whuber May 2 '16 at 18:34
• @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. – giordano May 2 '16 at 21:44
• 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). – Glen_b May 3 '16 at 0:44
• @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? – giordano May 3 '16 at 6:23
• 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? – Glen_b May 3 '16 at 8:37