Fitting an integral function given a set of data points

I have a set of measures of the resistivity of a given material at different thicknesses and I’m trying to fit them using the Fuchs-Sondheimer model. My code is:

data = {{8.1, 60.166323}, {8.5, 47.01784}, {14, 52.534961}, {15,     50.4681111501753}, {20, 39.0704975714401}, {30,     29.7737879177201}, {45, 22.4406}, {50, 15.2659673601299}, {54,     18.189933218482}, {73, 14.8377093467966}, {100,     15.249523361101}, {137, 15.249523361101}, {170,     10.7190970441753}, {202, 15.249523361101}, {230, 10.9744085456615}}  G[d_, l_, p_] := NIntegrate[(y^(-3) - y^(-5)) (1 - Exp[-yd/l])/(1 - pExp[-yd/l]), {y,0.01, 1000}];  nlm  = NonlinearModelFit[data, 1/(1 - (3 l/(2 d)) G [d, l, p]) , {{l, 200}, {p, 4}}, d, Method -> NMinimize] 

However it returns me these errors:

NIntegrate::inumr: The integrand ((1-E^(-(yd/l))) (-(1/y^5)+1/y^3))/(1-pExp[-(yd/l)]) has evaluated to non-numerical values for all sampling points in the region with boundaries {{0.01,1000}}. 
NonLinearModelFit: the function value is not a real number at {l,p} = {200.,4.} 

I think that the problem is in the way in which I defined the integral function G[d, l, p], because I had to fit a different set of data points with a different function of only one variable which I defined through the NIntegrate function and it gave me no error. Could anyone please help me?