title: Is the observed allele frequency significantly less than the predicted?

text: **Question** : How can I construct a test to determine if the observed "mountain"-allele frequency (Fig 1) is significantly lower in the central to southern mountains than predicted (Fig 2) by the ecological selection model ( _see below for details_ )? **Problem** : My initial thought was to regress the model residuals against latitude:longitude and altitude (which results in only the interaction between latitude and longitude as significant). The problem is that the residuals (Fig 3) may reflect variation not explained by the model and/or that their is something biological happening, e.g. allele has not had time to spread south to its potential or there is some barrier to gene flow. If you compare the observed (Fig1) versus expected (Fig 2) mountain-allele frequencies there is an obvious difference especially in the central to southern mountains of Sweden and Norway. I accept that the model may not explain all of the variation, but can I come up with a reasonable test to explore the idea that mountain-allele has not reached its potential in the central to southern mountains? **Background** : I have an bi-allelic AFLP marker whose frequency distribution appears to be associated with mountain (and latitude:longitude) versus lowland habitats on the Scandinavian peninsula (Fig 1). The "mountain"-allele is nearly fixed in the north, which is mountainous. It is nearly absent or fixed for the "lowland"-allele in the south, which lacks mountains. As one moves north to south in the mountains the "mountain"-allele occurs in lower frequency. This difference in the "mountain"-allele frequency from north to south may simply be due to phylogeography or historical processes, as the region was colonized from both the north and the south. For example, if the mountain-allele originates in the northern population, maybe it has not had time to expand fully into the southern population, assuming there are no barriers to gene flow (the northern and southern populations are panmictic for mtDNA and 12 microsatellites). My working hypothesis is that the "mountain"-allele frequency is a result of ecological selection (null hypothesis is neutral selection). For my ecological selection model, I have used a generalized additive model (GAM) with the binomial allele _frequency_ as the response variable (129 sites sampled across Fennoscandinavia with typically 10 to 20 individuals sampled at each site) and several climatic and growing season variables as the predictor variables. The model results are as follows (TMAX04-06 = max temperature in April-June, Phen_NPPMN = mean growing season vegetation productivity, PET_HE_YR = annual potential evapotranspiration, Dist_Coast = distance to coast):               Family: binomial      Link function: logit           Formula: Binomial_WW1 ~ s(TMAX_04) + s(TMAX_05) + s(TMAX_06) + s(Phen_NPPMN) +      s(PET_HE_YR) + s(Dist_Coast)          Parametric coefficients:                  Estimate Std. Error z value Pr(>|z|)         (Intercept) -0.74372    0.04736   -15.7   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Approximate significance of smooth terms:                  edf Ref.df Chi.sq  p-value         s(TMAX_04)    3.8100  4.812 25.729 9.43e-05 ***     s(TMAX_05)    0.8601  1.000  5.887  0.01526 *       s(TMAX_06)    0.8862  1.000  7.644  0.00569 **      s(Phen_NPPMN) 6.2177  7.375 39.028 3.16e-06 ***     s(PET_HE_YR)  3.1882  4.147 18.039  0.00145 **      s(Dist_Coast) 2.2882  2.857  9.725  0.01906 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           R-sq.(adj) =  0.909   Deviance explained = 89.7%     REML score = 326.73  Scale est. = 1         n = 129      ![enter image description here](http://i.stack.imgur.com/uJO5H.jpg) Figure 1. Observed "mountain"-allele frequency for bi-allelic AFLP marker. Contour lines 0.1 frequency intervals, color-shading is altitude with blues for lowest and reds highest. ![enter image description here](http://i.stack.imgur.com/sDN5x.jpg) Figure 2. Predicted"mountain"-allele frequency for bi-allelic AFLP marker. Contour lines 0.1 frequency intervals, color-shading is altitude with blues for lowest and reds highest. ![enter image description here](http://i.stack.imgur.com/nOv1P.png) Figure 3. Ecological selection model (using GAM) residuals broken down by the entire study area (Fennoscandinavia) and separately for Norway, Sweden and Finland. The red dashed lines represent a secondary contact zone between the northern and southern populations inferred from other AFLP markers and stable isotope analysis of feathers grown on their separate wintering grounds in Africa. The thin black dotted line is the center of the zone.


title: What method is used in Google's correlate?

text: Here is a recent Google correlate query:   http://www.google.com/trends/correlate/search?e=internet+usage&t=weekly# As you can see in the search box at that link, I entered "internet usage" and Google did the rest. It shows a value of 0.9298 as the "correlation" with the query "data mining". However, when I read **page 2 of the Google white paper [PDF]** , it says: > >     The objective of Google Correlate is to surface the queries in >     the database whose spatial or temporal pattern is most highly correlated >     with a target pattern. Google Correlate employs a novel approximate > nearest >     neighbor (ANN) algorithm over millions of candidate queries in an online >     search tree to produce results similar to the batch-based approach > employed >     by Google Flu Trends but in a fraction of a second. For additional > details, >     please see the Methods section below.... >   So, my question is:   Is Google using a normal Pearson or Spearman correlation to find this stuff or are they using something else? If so, can you explain the general technique? ================== Also, notice in the plot that the search for "internet usage" (and "data mining") drops during the summer months and really dives around Christmas. I would guess that kids and their homework have something to do with this.


title: Estimating a confidence interval around the average of an interacted treatment effect

text: I've got a causal inference situation with heterogeneous treatment effects. In addition to simply estimating coefficients, I'd like to get an average effect of the treatment at the covariate values that I observed in my dataset. For example: $$ y = \alpha+\beta_1 T + \beta_2Z + \beta_3 (T\times Z) + X'\gamma + \epsilon $$ I want to know the effect of the treatment $T$ (binary) across the gradient of $Z$ (continuous) that I observe. The expectation is the average of linear combinations -- that's easy. But the spread of the estimate seems a bit tricky. I think I've come up with a reasonable method -- basically take the parameter covariance matrix only for the treatment and the interacted variables, and use that to get the variances of each linear combination of data-parameter groups that I observe, sum those variances, get an SE, and make CI's like normal. it should work because $$ var\left(\displaystyle\sum X'\beta\right) = \displaystyle\sum var\left(X'\beta\right) = \displaystyle\sum X'V_pX $$ where $\beta$, $X$, and $V_p$ are only the subsets of the $\beta$'s, X's and the parameter covariance matrix that correspond to either the treatment or the interaction(s). Here's is silly little example in R:               #Population model     set.seed(413)     d = data.frame(         T = rbinom(1000,1,.5)         ,Z = runif(1000)         ,x1 = rnorm(1000)         ,x2 = rnorm(1000)         )     y = with(d,10+5*T+5*Z+-10*T*Z + 10*x1+10*x2+rnorm(1000,sd=10))          #Fit linear model     library(mgcv) #using mgcv because the predict lpmatrix functionality is useful, though this is parametric     m = gam(y~T*Z+x1+x2,data=d)     summary(m)     #Caluclate Average treatment effect on the treated     c1 = m$model$T*m$coef[names(m$coef)=="T"]     c2 = m$model$Z*m$model$Z*m$coef[names(m$coef)=="T:Z"]     ATET = (c1+c2)[m$model$T==1]     hist(ATET); abline(v = mean(ATET),lty=2,lwd=2)     lpm = predict(m,type='lpmatrix')     lpm = lpm[,grepl("T",colnames(lpm))]     vcv = vcov(m)     vcv = vcv[grepl("T",colnames(vcv)),grepl("T",colnames(vcv))]     variances = diag((lpm)%*%vcv%*%t(lpm))     variances = variances[m$model$T==1]     se = sqrt(sum(vars))     abline(v=mean(ATET)+se*1.96,col="red",lty=2,lwd=2)     abline(v=mean(ATET)-se*1.96,col="red",lty=2,lwd=2)      My basic question: is this right? Does my 95% CI have the advertised coverage properties? I guess that I am more of an applied coder than a proper statistician, so I'd appreciate reassurance from people who are more confident in such matters. EDIT: One more question about this -- the parameter covariance matrix that I'm using (cluster-robust) comes with a degree of freedom correction already in it. Does that justify the assumption that this interval is normally distributed? Or should I use $t$ with some number of degrees of freedom? Or would that be double-counting somehow?


title: Logistic Regression using dummy variables in MATLAB

text: I have a model where categorical (mutually exclusive) variables predict bankruptcy. Chi-square is significant. How can I code a logistic regression model in `MATLAB` that proves that some of these variables explain the bankruptcies better. There are 4 variables: which implies 3 dummy variables.               X = [  0 0 0           0 1 0            0 0 0            0 0 1            0 0 1            1 0 0             1 0 0 ......... ] ;     Y = [1 0 0 1 1 0 0 .....]; % 1 means bankruptcy      I think I should use `glmfit` (from http://matlabdatamining.blogspot.com/2009/03/logistic-regression.html) but I wasn't sure if using dummy indicators would require any additional interpretation/inputs. I use latest version and have Statistics and Optimization Toolbox.


title: Dependent variable is a function of independent variables; can I sensibly include them in a regression?

text: We've created a survey asking students, among other things, their GPA (=weighted average of grades) and their marks in some specific courses (which count towards GPA). We wanted to see which regressors influence the GPA using a simple OLS model. **Is it sensible to use a formula like this?**               GPA ~ grade_maths + grade_statistics + grade_privatelaw + ... + {other regressors, like study habits or origin}      Of course, the grade regressors turn out to be highly significant (some more than others, and _not directly related to the weight they have in GPA_ ), while few of the other ones are... Is this a case of endogeneity, i.e. does a regression like this violate strict exogeneity? With this regression we want to get an quick overview of which variables will likely be useful in the regressions that follow, which for example try to find out whether being good in quantitative courses tends to help get good grades in law and other ideas like this...


title: Difference between point-biserial and rank-biserial correlations

text: According to the wikipedia article the point-biserial correlation is just Pearson correlation where one variable is continuous but the other is dichotomous (e.g. Yes/No, Male/Female). However the article later introduces rank-biserial correlation, which is a correlation measure between a dichotomous variable and a ordinal/ranked variable. What is the difference? Is rank-biserial correlation related to Pearson correlation?


title: Inverse of false discovery rate (FDR)

text: Suppose we have a set $S$ consisting of $p$ features, and a subset $S_+$ of the features are positive. If $Q$ is any subset of $S$, define the false positive rate as the proportion of features in $Q$ which are not positive: $$FPR[Q] = 1 - \frac{|Q \cap S_+|}{|Q|}$$ where $|\cdot|$ denotes cardinality. If $Q$ is a function of the data, $d$, then we can define the false discovery rate as the expected false positive rate: $$FDR[Q(d)] = \mathbb{E}_d[FPR[Q(d)]].$$ Now suppose that I have a method for ranking the features in $S$ by likelihood of significance. I will report the top $r$ features most likely to be significant, based on my data, $Q_r(d)$. Formally, I have a family of set- valued functions $$Q_1(d) \subset Q_2(d) \subset Q_3(d) \subset \cdots \subset Q_p(d)$$ where $|Q_r(d)| = r.$ What I want to know is the maximum $r$ such that the set $Q_r$ has a false discovery rate less than a certain critical value, $q_{crit}$. That is, I want to know what is the value $$IFDR_{q_{crit}} = \max_r \{r \in \{1,...,d\}: FDR[Q_r] \leq q_{crit}\}$$ Is there a name for this 'inverse false discovery rate' function? If not, can you suggest a name better than 'inverse false discovery rate'?


title: Reducing false positive rate

text: I am running a model for a problem in insurance domain. The final results show some false positive x and some false negative y. I am using SAS Enterprise Miner for this. Can somebody suggest me how to reduce false positive? I know for this i have to increase the false negative. I want to know two things:   1. Is there any option in e-miner where I can give more weight to false negative and less to false positive?   2. Is there any general approach in modeling which tells us any ways to reduce false negatives or is it just a hit and trial approach?


title: Is it possible to specify a lmer model without any fixed effects?

text: In R, how do I specify lmer model without global fixed effect? For example, if I say something like               lmer(y ~ (1 | group) + (0 + x | group), data = my_df)      the fitted model will be $y_{ij} = a + \alpha_{i} + \beta_{i} x_{ij}$ How do I fit model $y_{ij} = \alpha_{i} + \beta_{i} x_{ij}$?


title: James-Stein estimator: How to calculate $\sigma^2$ in shrinkage factor?

text: I had a question on calcuating James-Stein Shrinkage factor in the 1977 Scientific American paper by Bradley Efron and Carl Morris. I gathered the data for the baseball players and it is given below:               Name, avg45, avgSeason         Clemente, 0.400, 0.346         Robinson, 0.378, 0.298         Howard, 0.356, 0.276         Johnstone, 0.333, 0.222         Berry, 0.311, 0.273         Spencer, 0.311, 0.270         Kessinger, 0.289, 0.263         Alvarado, 0.267, 0.210         Santo, 0.244, 0.269         Swoboda, 0.244, 0.230         Unser, 0.222, 0.264         Williams, 0.222, 0.256         Scott, 0.222, 0.303         Petrocelli, 0.222, 0.264         Rodriguez, 0.222, 0.226         Campaneris, 0.200, 0.285         Munson, 0.178, 0.316         Alvis, 0.156, 0.200      avg45 is the average after 45 at bats and is denoted as y in the article. avgSeason is the end of the season average. The James-Stein estimator for the average ($z$) is given by $ z = \bar{y} + c (y-\bar{y})$ and the the shrinkage factor $c$ is given by (page 5 of Scientific American 1977 article) $ c = 1 - \frac{(k-3) \sigma^2} {\sum (y - \bar{y})^2} $ $k$ = number of unknown means. Here there are 18 players so $k = 18$. I can calculate $\sum (y - \bar{y})^2$ using avg45 values. But I don't know how to calculate $\sigma^2$. The authors say $c = 0.212$ for the given data set. I tried using both $\sigma_{x}^2$ and $\sigma_{y}^2$ for $\sigma^2$ but they don't give the correct answer of $c = 0.212$ Can anybody be kind enough to let me know how to calculate $\sigma^2$ for this data set?


title: Multiplicative unobservable component in state space model

text: I'm new here and wondering if anyone could give me some hints on how to estimate the time varying coefficient and state variable together. Here is my model: observation equation: $Y(t)= A(t)X(t)+ w(t)$,   state equation: $X(t)=\phi X(t-1)+v(t)$, here I have time varying coefficient $A(t)$, it doesn't depend on any predetermined parameter $\theta$, for example. If I treat $A(t)$ as another state variable, then it is nonlinear state space, I have no idea how to estimate multiplicative state variables. Any hints? Thank you.


title: When forecasting sequential data is it best to use auto-regressive models or build a more traditional n x p dataset with features?

text: I'm familiar with the use of auto-regressive models when it comes to forecasting a single vector of time-series data. Is anybody familiar with a more traditional modeling approach, i.e. - creating sets of features such as indicators for day of the week, time of the day, day of the month, holiday, and then running a model such as a regression or random forest on this? Are there pros/cons to each of these methods? I'm basically tasked with forecasting hourly requests based on A LOT of historical data. There are strong intraday trends as well as pretty strong weekly trends. So far we have been using the average of the past 4 data points of the same hour in the same day (so for this Friday at 4 PM we would average the last 4 Friday's counts at 4 PM) and this works surprisingly well. Is it even worth building a more sophisticated model? Would I have to continually retrain it between every hour or would a few months worth of this hourly data be enough to forecast for several days before retraining? I'm sure I left some questions out so any suggestions and literature you could point me to would be much appreciated.


title: Learning an interpretable model

text: I am working on problems in the field of medical imaging where the need for a simple and interpretable model is important from a clinical perspective. This means that I have to explain the algorithm's prediction to non-experts (well non-experts in Mathematics). My questions is twofold:   1. As far as I know, an interpretable model is a model that assigns weights to each feature or weak classifier to form a strong classifier. Are there other methods that make a model interpretable?   2. Which classification methods are interpretable? I know that the linear SVM and AdaBoost are so. Are there other methods?


title: Periodogram vs. spectral density diagram of a time series

text: Could someone explain to me the difference between a periodogram and spectral density diagram? ![enter image description here](http://i.stack.imgur.com/0Z3SB.png) The first diagram is produced with this block of code:               FF = abs(fft(datalist)/sqrt(128))^2     f = (0:63)/128     plot(f,FF[2:65],type="l",xlab="Frequenz",ylab="Spektrum")      and the second one with this code:               x.spec<-spectrum(datalist,log=c("no"))      In both diagrams I have two peaks, which are clearly higher than the others, but in the first they are almost the same, but in the second (periodogram) the second peak is obviously higher than the first peak. Why are these two diagrams different? In a periodogram, the first peak shows that we have something periodic in our time series, or second one or both of them? How can I interpret these two peaks? **EDIT** in the first diagram for the second peak:               x=0.109375  y=36657.41193        x=0.1171875 y=36731.11184      in the second diagram for the second peak:               x=0.128  y=88176.01878       and My Data: 471 379 484 479 527 548 576 534 443 375 475 514 516 527 445 403 487 382 510 451 562 569 575 528 450 351 467 504 505 520 441 407 460 421 504 475 569 555 575 516 460 359 460 496 492 529 490 465 410 460 475 509 549 564 571 515 375 392 458 450 488 510 459 475 410 514 496 563 550 572 539 396 388 491 447 471 515 454 467 363 499 486 522 583 567 543 459 370 437 500 483 529 451 408 463 395 511 468 559 551 576 527 418 380 471 496 494 515 453 415 454 431 503 490 575 564 571 512 400 371 484


title: Probability of overlap of areas in 2D space

text: Assuming I observe, in a unit square, $n_1$ circles of area $A_1$ (non- overlapping amongst themselves) and $n_2$ circles of area $A_2$ (again, non- overlapping) and that each of the centres is uniformly distributed across the square (well, I guess that can't be precisely true given the non-overlapping constraint), is there an expression for the proportion of the first circles that overlap with at least one of the second type of circles. The actual situation is observed cancer cells and pores on a microscope slide, to which the above seems a good approximation. In the 1D case, I'd like to think that for each line segment (1D equivalent of a circle) of the first group, there's an interval of size $A_1+ A_2$ before the end of each segment of the second type to place the start of the interval so that they would overlap, and so $n_2(A_1+A_2)$ of the unit interval would result in an overlap. That seems to be true for the first segment of the first group, but it won't necessarily hold for further segments as the non- overlapping criterion breaks the independence, so I'm probably looking at this (and the 2D case) from the wrong angle. ### Updated thoughts Maybe this is the way to look at it? There's a total area $n_1A_1$ covered by type 1, and $n_2A_2$ covered by type 2, out of a total area of 1. So an expected area of $n_1A_1n_2A_2$ is covered by both, and divide this by $A_1$ to give the number of type 1's that would be needed to cover that area, and divide by $n_1$ to get the proportion of type 1's that are in the double- covered area. So if 7% of a slide's area is covered by spores, then 7% of cells are hit by spores?


title: "Exotic" text representations methods?

text: I'm looking to the different methods of representing a text into a machine- readable format. However, until now, I only found "Bag of Words" approachs with a lot of variations (boolean BoW, weighted BoW, Bag of Topics...). What is troubling me is that these methods don't take context into account (some approach this goal by counting n-grams, wich seems a bit light to me). I would like to know if other approaches are also used (I found no mentions of graph representation, for instance). Papers/books/websites presenting these methods would be a huge bonus.


title: design and analysis of experiments

text: Kindly recommend some teaching material for "design and analysis of experiments" for the students of biostatistics and epidemiolody. Pls note that the students have only little knowledge on statistics. It will be useful if the book contains many applications rather than theory.


title: One class classification with libsvm. Accuracy results in 0%

text: A quick recap for what I want to do, I want to determine if a text is written by the same author or not. Thus I use one-class classification.   In my training set (18 samples), it looks like this (for simplifying, I used x as data value):               1 1:x 2:x "until" 200:x     1 1:x 2:x "until" 200:x      In my testing set (3 samples), it looks like this (for simplifying, I used y as data value):               1 1:y 2:y "until" 200:y      For data preparation (training and testing set), I set upper and lower scaling limit to +-1               -l -1 -u 1      For training, I use svm_type is one class svm, kernel type is sigmoid. Yet the accuracy is 0%               optimization finished, #iter = 13     obj = 22.901769047004553, rho = 5.476401914859387     nSV = 11, nBSV = 6     Accuracy = 0.0% (0/21) (classification)      Can someone show me what I did wrong here?   Thank you very much.


title: Finding Repeating Combinations

text: I have a table where               Field1 = Name {John, Mark, Brian} and      Field2 = Ball Color {Red, Orange, Yellow}      Each person has 6 balls of different colors so for example here are all the records for John and Mark: **Name** **Ball Color** John Red John Orange John Yellow John Green John Blue John Indigo Mark Orange Mark Yellow Mark Green Mark Blue Mark Indigo Mark Red How can I find every possible combination of colors chosen (regardless of order - order doesn't matter) and the number of times each combination repeats. Note: I don't mean the theoretical combination, I actually have a dataset that I want to analyze.


title: Abnormally high accuracy with repeated 10-fold CV and ordinal regression

text: I am using repeated 10-fold CV to calculate the accuracy of my ordinal regression model. I have 6 predictors, 10 ordered response categories, and a total of 1166 data points. For the ordinal model, I have defined accuracy as _1 - loss_ , with _loss_ being a simple linear function of the distance between observed and predicted classes, assuming classes are equidistant. I chose the number of repetitions for the 10-fold CV by examining the stability of the results as advised here, and decided to use 5 repetitions. Finally, I am using this measure of accuracy to compare the quality of the predictions given different predictor values. More precisely, I have a range of possible storm wave conditions and associated parameters, 10 levels of increasing damage for ships caught in the storm, and I am looking to find the wave conditions which have more likely triggered the damage. So I run the model for each one and compare the accuracy. My problem is, for all conditions tested my accuracy is suspiciously high i.e. > 80%. I am not expecting this kind of accuracy given the quality of some of my predictors (for example, the wave velocity has been estimated numerically on a coarse grid with a number of simplifying assumptions). Also, there is not much variation in accuracy between the model outputs for different wave conditions, maximum +/- 1%. Again I would have expected more than this. I am wondering if there is something I am missing here, regarding my estimation of accuracy. What could be the cause of this issue?


title: PCA of hyperspectral image data

text: I have some hyperspectral image data similar to this, and I want to do a PCA on it. The problem: I've never done a PCA, and its specially difficult for me to do it on 3D data. How can I do it in MATLAB or R? [See also Analyzing 3D data: What can be done? ]


title: Predicting zeros with count data model

text: I'm working with count data (ticket sales, to be specific) and I'm having trouble fitting a model to it. I've tried a linear one and a transformed linear one, but the residuals end up being non-normal and having a non- constant variance. I've been trying to play around with other models (negative binomial, poisson, etc.) but I can't seem to get anything to fit very well. Here's the histogram for my response variable and what the residual plots look like for my linear model: ![enter image description here](http://i.stack.imgur.com/WKmyT.png) ![enter image description here](http://i.stack.imgur.com/1wLRS.png) ![enter image description here](http://i.stack.imgur.com/OAlpU.png)


title: Combining labeled and unlabeled data for training

text: When training my algorithm, if I can get some i.e. data my future test data that has no labels can it improve my algorithm's efficiency, is there any mathematical proof for it? PS: I think semi-supervised methods and transductive learning are something similar to what I look for.


title: Calculation of posterior distribution of a Gaussian process

text: Suppose I have a set of points $Y$, $X$ and I model them by using a Gaussian process $\mathcal{GP}(m,K)$. Let the noisy observations be given by $$y^i = f^i + \sigma^2$$ $p(f)$ gives the prior on the Gaussian process. Then what will be the posterior $p(f|Y,Z)$ which is again a Gaussian process with a mean $(m',K')$. How can $(m',K')$ be derived?


title: Assessing the ability of two quantitative methods measuring the same thing

text: Suppose I have two scientific approaches to measuring the same physical quantity taking values between $[0,1]$. One (say method A) is generally accepted to be more accurate under certain conditions, but it is more of a gross measure (low spatial resolution). The other method is newer in the field, has the a higher resolution, but may have some accuracy problems at lower levels. What is the best way to determine how well the two methods agree? I've visualized it with a scatter plot, and fit the data to a straight line and can get a slope close to 1.0, reasonably small intercept, and $R^2$ quite near 1.0 -- but isn't that more of a correlation measure? I would appreciate any suggestions on what types of statistical methods to explore.


title: Sample size vs. number of samples

text: This question came up in an archaeology class, and it seems like it should be something very basic. Suppose I want to count the number of marbles in a large bucket of sediment. I only have the manpower to process 20 liters of sediment in total. Is it better to take:   20 samples of 1 liter each   5 samples of 4 liters each   1 sample of 20 liters   ...? Does it even matter, if the same total amount of sediment ends up being counted? I am assuming that that this is with replacement, so the sediment in a sample is tossed back in after its marbles are counted.


title: Ordinal regression vs. summated score

text: It seems quite common to see scientific researches using a summated score to represent an overall scale that consists of many individual scale items (ex. Likert scale items). I assume that regular statistical tools (ex. ANOVA, linear regression, etc.) can be used on these data, which are treated as quantitative data. However I've came upon another tool when dealing with ordinal data - ordinal regression (or ordered logistic regression). My initial question was: would it be better if I just used ordinal regression rather than calculating summated score first _then_ use multiple regression? But I found out that ordinal regression seemed applicable only when all individual scales were included to fit a dependent variable. In cases where scales are made up of subscales, then generating summated scores is a must (or perhaps not?), meaning that data are treated as quantitative already. I'm not sure if what I understand here is correct or not. Verification is greatly appreciated.


title: Does high log-likelihood imply high R^2

text: Does a high LL value imply that the model has a high R^2? I'm a very beginner to statistics so please excuse my naivete.


title: What does "Conditioning on the margins of ____" mean?

text: What does it mean to "condition on the margin of ____"? I lack statistical/mathematical training and phrases like these leave me bewildered and unsure where to look. In this post (Which are differences between the hypergeometric distribution and chi-square distribution), there is an example where @Glen_b uses the phrase: > That is, your example data are like so: >   >   >                 Drawn    Not drawn    Total >   >        Red      160        6520        6680 >        Black    222       11938       12160 >   >        Total    382       18458       18840 >   > > Looking at the number of reds drawn as a random variable, that has a > hypergeometric distribution (though there formulated in terms of white and > black balls drawn from an urn rather than red and black balls drawn from a > universe). > > [ **Conditioning on the margins gives the hypergeometric** - this is also > the situation used for Fisher's exact test based on the hypergeometric, and > one of the situations for which the usual 2x2 chi-square test of > association/test of independence applies. If you don't condition on both > margins, you don't have a hypergeometric; but that's what you normally do in > the specific balls-in-urns model you describe.]


title: Pairwise comparisons after significant interaction results: parametric or non?

text: I've run a 2-way ANOVA on growth rate data (grams/day) with the factors of year type (good and poor) and site (A and B). Though the data themselves are non-normal and do not have homogeneous variances, the residuals fall pretty nicely along the qq plot, and they are not heteroskedastic. Running the test shows that there is an interaction between year-type and site. It's been suggested to me that I now must run a series of pairwise comparisons to look for differences because of this interaction effect, which I assumed I'd need to do anyway. My stats program (Sigmaplot11, which includes the SigmaStat package) automatically does post-hoc tests for any significant results (in this case is reports: "All Pairwise Multiple Comparison Procedures (Holm-Sidak method): Overall significance level = 0.05). I assume that this is the correct way to conduct the pairwise comparisons rather than making them separately? Here's why I ask: Using these built-in post-hoc tests, I get significant results in both good years and poor (p<0.001 in each case). However, doing these tests separately (good vs. good and poor vs. poor) as Mann-Whitney rank sum tests (the raw data are non-normal with heterogeneous variances), I get no significant difference between sites within poor years. t-tests show the same as the ANOVA. I assume this is because I'm not comparing means in the case of the Mann-Whitney, but which post-hoc method should I be using in this case?! **EDIT** : following @whuber's informative comment, I've taken a look at my other 2-way ANOVAs run in a similar manner. I've done another similar test today comparing adult weight. The t-test from the multiple pairwise comparisons after the 2-way ANOVA shows no difference (t=1.547, p=0.122), but a t-test run outside of the ANOVA shows highly significant difference (t=-4.739, p<0.001). As @chl pointed out, I would expect these 2 t-tests to have the same t-values; note that these 2 tests use exactly the same original data. Any idea why this might be or how I can interpret this? Thanks for any suggestions you can provide!! **EDIT #2** : Just to update anyone who's interested, I've taken a look more closely at the numbers behind the test, and it looks like somehow the software is not doing what I'm asking. It lists a table called the Least Square Means, listing each site, its mean and its standard error of the mean. The overall site means listed in this table are not correct. However, for each site within year type, it does have the correct means in the Least Square Means table. I'm giving up on the main factor of site comparison within this test, and sticking to the t-value given in the t-test for this main factor comparison. I'm still not sure what's going on exactly, but from comments and answers provided (thanks again!), I feel that this is a safe move, and I must move on with other things. Thanks for all of your help!


title: Autocorrelation for regression

text: I am attempting to use a significant autocorrelation (where it lies outside a 95% interval around 0) indicating periodicity of a signal and use it as predictive variable in a regression. If, for example, there's a significant autocorrelation at frequency 3, does this mean that I can use the time series as predictor of itself by placing together the last signal with a 3 period lag? Or as I have read autocorrelation does not indicate where in time the significant periodicity occurs just the how frequently it occurs. However if I have 3 periods lag, I have 3 choices - is this correct?


title: Accounting for spatial autocorrelation in neural network analysis in R

text: I have been exploring the use of neural networks to identify relationships between large numbers of differing species and the presence of a specific bacteria in the environment. Previously I had been using mixed effects models, as I could account for spatial autocorrelation in my data. I was wondering however if anybody knew whether it was possible to incorporate location data to account for spatial autocorrelation in neural network analysis? Many thanks for any help.


title: Failure Rate of Contraception over n years

text: According to multiple sources, the probability that, after one year of typical use, a woman who uses male condoms will get pregnant is $18\%$. Presumably, this assumes that the woman will have sex at an average frequency, but I found it hard to find information about this. What's surprising is that I calculated the probability of unintended pregnancy over four years using male condoms to be very close to $55\%$. I calculated this by noting that, in order not to be pregnant after four years, contraception has to be successful $4$ times. In other words, the probability of successful contraception over four years is $(1 - .18)^4 = .4521$. Therefore, the probability of contraception failing and producing a pregnancy sometime during those four years is $54.79\%$. Is this technique sound?


title: Simulated single value based on multiple chains in RJAGS

text: I am using RJAGS to simulate the posterior distribution of event that a certain candidate will win the presidential election. I need to find the actual percentage that one of the candidates will have. What I do have is the values of each iteration for each of my three chains for both the proportion of the first candidate winning and the second candidate winning. I don't know how to obtain a single result since I have all of these values for each iteration in the the chains. Should I calculate just an average of all those values for each candidate? Thanks for the help


title: Identification of peer/neighborhood effects in a multilevel framework

text: My question concerns estimation of “peer effects“ or “neighborhood effects” in a multilevel framework. The idea of such an effect is that the behavior of a household (on level-1) is influenced by the behavior of others in the same cluster/neighborhood (level-2) who are perceived as peers. Manski (1995:127) termed this an “endogenous effect”, i.e. that “the propensity of an individual to behave in some way varies with the prevalence of the behavior in the group’’, and pointed out that identifying such an effect is difficult if not impossible. When going through literature and searching the net, I have found many different examples of how people try to estimate this, but also criticism of some of the same strategies. _The reference for Charles F. Manski (1995) is “Identification Problems in the Social Sciences”, Harvard University Press, but the main arguments are also included in the following article: Charles F. Manski (1993) “Identification of Endogenous Social Effects: The Reflection Problem, The Review of Economic Studies, Vol. 60, No. 3 (Jul., 1993), pp. 531-542”._ In a multilevel framework where the outcome is a binary choice, say to buy product A or not, it would concern how the outcome on the household level is influenced by the outcome of others on level-2. The data I am using is a representative sample of 32.000 Indian households, from 1450 villages or neighborhoods if urban. It is cross-sectional observational data. The villages/neighborhoods are primary sampling units used by the Census of India. People in the same village or neighborhood are close enough to be potential influencers. Manski (1995, 1993) discussed three hypotheses to explain how groups (or societies, or environments) may affect individuals: endogenous, contextual, and correlated effects. _Endogenous effects_ : The individual behavior tends to vary with the prevalence of that behavior in the group. In my example, this effect arises when households’ choices of purchasing product A tends to vary with the average choice of all households in the group (i.e. village or urban neighborhood). _Contextual effects_ : The household behavior tends to vary with the distribution of background characteristics in the village/neighborhood. There is a contextual effect when the propensity of an individual to behave in some ways varies with the mean of exogenous variables. _Correlated effects_ : Individuals tend to behave similarly because they have similar individual characteristics in the group (i.e., household education or wealth) or become part of similar institutional environments (i.e., households are part of the same social community – “caste” – or the same social clubs, etc.). My main hypothesis concerns an endogenous effect, as I described in the original posting, but I also hypothesize contextual as well as correlated effects. Manski referred to the so-called reflection problem that arises when “a researcher observes the distribution of behavior in a population and wishes to infer whether the average behavior in some group influences the behavior of the individuals that compose the group”. This ‘‘reflection problem’’ arises out of the presence of the mean outcome of the group as a regressor in an equation that also includes the effects of context (e.g. the mean wealth of households in the village/neighborhood) as well as individual characteristics or other group-level variables (for example, institutional environments). Manski showed that the endogenous effects (in Manksi’s terminology) are unidentifiable (in particular it is not possible to distinguish endogenous and contextual effects, although in some situations it may be possible to determine whether an overall social effect (defined as endogenous plus contextual effects), is present. One way I can think of dealing with this would be to calculate group means but subtract the values for each household and include these as level-2 variables, but I have a feeling that there will still be endogeneity problems that come along with doing so that I am not seeing yet. My question is: how I can best identify whether the propensity of a household to buy product A varies with the prevalence of buying product A in the village/neighborhood? KML


title: Bayesian Hierarchical Model - Exact Conjugate Solution?

text: I was hoping to get some help. In understand how to compute an exact numerical solution (http://www.cs.berkeley.edu/~jordan/courses/260-spring10/lectures/lecture5.pdf) for the following Bayesian model: $$ \tau \sim Ga(\alpha, \beta)$$ $$\mu \sim N(m ,p)$$ $$Y_i \sim N(\mu, \tau)$$ where the data is $Y_i$. I am performing a meta-analysis where my data are $Y_i$ from each study, and $\tau_i$ from each study. The problem is described by a hierarchical model: $$ \tau \sim Ga(\alpha, \beta)$$ $$\mu \sim N(m ,p)$$ $$\theta_i \sim N(\mu, \tau)$$ $$Y_i \sim N(\theta_i, \tau_i)$$ Is there any way to solve this model analytically? Any help/suggestions would be greatly appreciated.


title: Simulating state space model with AR(1) dynamics

text: I asked a question similar to this previously: Check that state space model implemented correctly However I think I have a better handle on it now and want to re-ask it: I simply want to simulate data from a state space model where the state variables follow an AR(1) process (see the code in the first link above). Given the burn in issues (see link below) I assume it's better to determine empirically how many observations are required until the system reaches its theoretical unconditional variance then ensure that I simulate at least 2 times that amount before using the x(t) generated by the AR(1) process in my state space model. http://www.mathworks.co.uk/help/econ/simulate-stationary-arma-processes.html Q1.) If I want to simulate data from my state space model is it necessary to ensure that the AR(1) process is in it's equilibrium state first? Q2.) From the simulated data I will be able to estimate the observation and state error variances as well as the AR(1) parameters and the unconditional mean and variance of the process (averaged over many sample runs). Assuming these empirical values all match their corresponding theoretical ones, can I then be fully satisfied that the state space model which is based on the AR(1) process has been implemented correctly? Q3.) How can I estimate what the likely error bounds should be on the parameters that I propose to estimate in Q2? Baz


title: Standard errors for lasso prediction using R

text: I'm trying to use a LASSO model for prediction, and I need to estimate standard errors. Surely someone has already written a package to do this. But as far as I can see, none of the packages on CRAN that do predictions using a LASSO will return standard errors for those predictions. So my question is: Is there a package or some R code available to compute standard errors for LASSO predictions?


title: Why does non-normally distributed errors compromise the validity of our significance statements?

text: There is a normality assumption when it comes to consider OLS models and that is that the errors be normally distributed. I have been browsing through Cross Validated and it sounds like Y and X don't have to be normal in-order for errors to be normal. My question is why when we have non-normally distributed errors is the validity of our significance statements compromised? Why will confidence intervals be too wide or narrow?


title: Predicting the growth of a social network

text: I am building a predictive model for the growth of the amount of users of a new p2p protocol inspired by bitcoin and I would like to use historical data collected from the growth of major social networks in my model. My specific questions are the following:   1. Does a public database exist containing growth data for major social networks?   2. Does there exist a generally accepted method for estimating contagion given specific network characteristics/assumptions?   3. I need to account for the loss of users over time and thus I need a proper way of measuring node death and its impact upon connected nodes in the model. Is there a generally accepted method to estimate this given some set of parameters? As a side note, I am uncertain, which stack exchange to post this question in or the best way to format the question for this channel. If there is a forum one could recommend in the comments that is better suited, then I would appreciate the advice.


title: Is statistical insignificance fatal

text: I apologize if this question has been done to death, but as a non-statistician I really don't know what is the bottom line takeaway. I am looking at a sample of 30,000 individuals who were the subject of an economic intervention. This intervention produces a mean increase in the annual income of the subjects of $2000. However the result is not statistically significant at usual levels. From a statistical viewpoint, has the intervention failed, and if not what further can be done?


title: Constrained Optimization library for equality and inequality constraints

text: Any recommendations for choice of a constrained optimization library suitable for my optimization function? I am minimizing a i) non-linear function with linear equality and inequality constraints, and ii) have available the gradient and the hessian of the function. If it helps, the function I am minimizing is the Kullback-Liebler divergence. constrOptim only deals with inequality constraints. Quadprog handles quadratics. Trust does not support constraints. So the KL divergence does not fit into these solutions. There are quite a few solutions on the R Cran Task page for Optimization. Iam able to perform the optimization in MATLAB using the fmincon() function which seems to use an interior-point or a trust-region-reflective. Ideally there is a library that is well-suited to the problem defined.


title: Can I use all pairwise comparisons at a factor level to perform correlation analysis?

text: I have 4 data sets, each with thousands of records/instances. There are 2 factors. One factor has 3 instances at each factor level (there are thousands of levels in the data set). The other factor also has thousands of levels, but with varying numbers of instances at each level. For each factor, I would like to measure how similar the values are. The variable whose similarity I want to measure is non-Normally distributed in at least 1 of the data sets. One idea was to create a table with all pairwise comparisons at each factor level and do correlation analysis. i.e., for the first factor, I would get a table that looks like this:               factor_level col1 col2     factor1 val1  val2     factor1 val1  val3     factor1 val2  val3     factor2 val1  val2     factor2 val1  val3     factor2 val2  val3     ...       and that would be used to calculate correlation coefficients (col1 vs. col2). For the first factor, I would like to get one correlation coefficient for the entire data set (to compare to another published data set for which I would perform the same analysis). The table for the second factor looks like this:               factor_level col1 col2     factor1 val1 val2     factor1 val1 val3     factor1 val1 val4     factor1 val2 val3     factor1 val3 val4     factor2 val1  val2     factor2 val1  val3     factor2 val2  val3     ...      For the second factor, I would also like to calculate the correlation at each factor level (in the end, I would like to rank the factor levels in my data set by the similarity of their values and design some experiments to test the highest-ranking ones). After doing this, I got results that seemed to make sense (i.e. the correlations, Spearman's rho, were high in cases where I would expect that based on published data and low where I would expect them to be low). Out of curiosity, I decided to try another method to calculate correlation coefficients. I randomly split the data at each factor level into 2 groups (if there were 3 rows at a factor level, I left one of them out) and paired up 1 row from each group. The table I used to do the correlation analysis looks like this:               factor_level group1_val group2_val     factor1 val1  val1     factor2 val1  val1     factor2 val2  val2     factor2 val3  val3     ...      I repeated this analysis 10 times. That is, I did the random splitting into 2 groups 10 times and calculated the mean Spearman's rho across the 10 runs. The correlation coefficients using this method were generally lower than what I obtained using the all-pairwise comparisons. I'm sure there is an obvious explanation for this. In a way, I am not disturbed by the differences, because the trends are still the same (i.e. those 'control' cases where I expect high or low correlations are still ranked as expected). And the rankings based on correlation coefficients at the second factor level are similar. My question, though, is whether it is considered appropriate to calculate correlation coefficients as I have done (i.e. using all pairwise comparisons). Or whether it's something other people have done. I don't have 2 variables, so it seems wrong to do it as I have. I'm guessing it might be better to use some form of ANOVA, but I would be happy to use my correlation coefficients so long as their use is not too egregious. Thank you in advance for your help.


title: Likelihood function of a Linear probability model

text: What is the Likelihood function of a linear probability model? I know the likelihood function is the joint probability density, but how to construct the likelihood function when we only have the probability $P(Y_i=1|X_i)$ and $P(Y_i=0|X_i)$?


title: Assumptions for intraclass correlation

text: I am having difficulty finding information on the assumptions for the intraclass correlation. Can someone please tell me what they are?


title: What to do when ovtest and linktest in Stata suggest model misspecification?

text: I have a sample that consists of 50 observations. The base model of the OLS- Regression with three control variables, two of them significant, has a $R^2=0.50$ and its F-Value is 7. Both `–ovtest` and `–linktest` (Stata 10.0) indicate no model misspecification. I am testing four Hypotheses. However, after having included the four variables – two of them significant - in the regression, the –ovtest and –linktest (Stata Version 10.0) indicate a model misspecification. I did regression diagnostics, which lead to the removal of three outliers because of their influence and a violation of the assumption that the residuals are normally distributed. With all seven variables in the model, $R^2$ is 0.65 (Huber-White robust standard errors) and the F-Value is 7. Even if I add only the two significant independent variables, both tests are failing, i.e. the ovtest, for example, rejects the H0 that the model has no omitted variables. So my question is: what can I possibly do to improve matters? I read that the `ovtest` is critical for proceeding with the analysis, but I cannot detect any further violation of regression assumptions. Dropping the variables is no alternative, since they are the ones of theoretical interest. Edit: I added three observed vs. fitted plots 1) 50 Observations, 7 variables 2) 47 Observations (without 10,12,48), 3 variables 3) 47 Observations (without 10,12,48), 7 variables ![50 Observations , 7 variables](http://i.stack.imgur.com/f1YXt.jpg) ![47 Observations \(without 10,12,48\), 3 variables](http://i.stack.imgur.com/Kkzq7.jpg) ![47 Observations \(without 10,12,48\), 7 variables](http://i.stack.imgur.com/r6bRT.jpg)


title: What is the effect of dichotomising variables?

text: * When dichotomising variables, what information is lost in the process?    * How does a dichotomisation help in the analyses?


title: Hidden Markov model (forward algorithm) in R

text: I have studied some of these resources and I know that there is an R package called HMM. Could anybody explain the usefulness of the 'forward algorithm' with a simple example in R? Finally, I find a good book which is Hidden Markov Models for Time Series: An Introduction Using R (Chapman & Hall CRC Monographs on Statistics & Applied Probability) Walter Zucchini, Iain L. MacDonald


title: What are alternatives to scatterplots when studying the relationship between two numeric variables?

text: I have questions about a study I was reading. In this study about children and money, the relationship between the age of the children and their weekly pocket money was assessed. The researchers studied the relationship between age and pocket money on a scatter plot and added a linear regression model.   * Why did they use a scatter plot?    * Could they have used something else?


I want to calculate the percentage of workers getting wages between $22$ and $58.$ I don't know what to use. Please help me. I've learned central tendency and interpolation. $$\begin{cases}Wages&&No~of~Workers\\0-10&&20\\10-20&&45\\20-30&&85\\30-40&&160\\40-50&&70\\50-60&&55\\60-70&&35\\70-80&&30\end{cases}$$


title: lme vs lmer vs aov for three factor repeated measures design

text: I'm trying to run a tree-way repeated measures ANOVA on the following data: I have a completely balanced design with three within-subject factors (type, form and ch - channel) and one dependent variable amp (amplitude). I'm inclined to believe the results I got using `aov` function:               res <- aov(amp ~ type*form*ch + Error(sbj/(type*form*ch)), data = p3vals)      Here is the anova table I have:               Error: Within                    Df Sum Sq Mean Sq F value   Pr(>F)         type            1   25.0  24.950  12.315 0.000462 ***     form            1   12.9  12.910   6.372 0.011693 *       ch              1    0.9   0.875   0.432 0.511113         type:form       1    3.1   3.123   1.542 0.214581         type:ch         1    0.9   0.938   0.463 0.496404         form:ch         1    1.3   1.256   0.620 0.431239         type:form:ch    1    3.0   2.960   1.461 0.226974         Residuals    1514 3067.3   2.026                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      However, after reading several pages of similar examples(including Why do lme and aov return different results for repeated measures ANOVA in R? ) I decided to try `lme` and `lmer` functions from name and lme4 packages for further pairwise multiple comparisons using `glht` from multcomp package. In the example mentioned above F values are at least closet those obtained using `aov`, but I cannot figure out how to get any meaningful results.               lme_p3amp = lmer(amp ~ type*form*ch + (1|sbj) + (1|type:sbj) + (1|form:sbj) + (1|ch:sbj), data = p3vals)     > anova(lme_p3amp)     Analysis of Variance Table                  Df  Sum Sq Mean Sq F value     type          1  0.0266  0.0266  0.0433     form          1  2.0782  2.0782  3.3863     ch            1 28.5513 28.5513 46.5227     type:form     1  2.1980  2.1980  3.5815     type:ch       1  2.9789  2.9789  4.8539     form:ch       1  0.9278  0.9278  1.5118     type:form:ch  1  6.6072  6.6072 10.7661      and lme produces the following result:                anova(lme(amp ~ type*form*ch, random=list(sbj=pdBlocked(list(~1, pdIdent(~type-1), pdIdent(~form-1), pdIdent(~ch-1)))), data=p3vals))                  numDF denDF  F-value p-value     (Intercept)      1  1508 32.56485  <.0001     type             1  1508  0.02920  0.8643     form             1  1508  3.54422  0.0599     ch               1  1508  8.05747  0.0046     type:form        1  1508  2.79623  0.0947     type:ch          1  1508  3.78969  0.0518     form:ch          1  1508  1.18037  0.2775     type:form:ch     1  1508  8.40553  0.0038      I'd appreciate if you tell me what is wrong with my code and how can I perform a valid post hoc analysis.


title: How does the size of a dataset affect confounding in both randomised trials and observational studies?

text: I have heard that larger sample sizes in randomised trials lead to a smaller possibility of confounding.   * Why is this true in the case of **randomised trials**?   * Also, how does sample size of an **observational study** affect confounding?


title: What does the Reward function depend on in a Markov Decision Processes (MDPs), in the context of Reinforcement Learning?

text: I was trying to understand MDP's in the context of reinforcement learning, specifically I was trying to understand what the reward function explicitly depends on. I have seen a formulation of the reward function as defined by Andrew Ng in his lecture notes as: $$R: S \times A \mapsto \mathbb{R}$$ Which means that the the reward function depends on the current state and the action take at that state and maps to some real number (the reward). To get a different perspective, I read the interpretation wikipedia had: The process responds at the next time step by randomly moving into a new state s', and giving the decision maker a corresponding reward $R_a(s,s')$. Which seems to be a different interpretation in my opinion since this would make the reward function more of a function of the form: $$R: S \times A \times S\mapsto \mathbb{R}$$ Which in my opinion, seems to be a completely different thing. I was trying to understand if the two formulations were actually the same (and if it was possible to prove their equivalence) in the context of MDP's applied to reinforcement learning.


title: chisq.test & post-hoc analysis?

text: I'm working on a R data frame, such as:               - column 1 : ID       - column 2 : count of event 1       - column 3 : count of event 2        (1 unique ID per row) I try to compare the proportion of events 2 among events (1+2) across all IDs. For now, I am just plotting the proportion estimates and their confidence interval, defined by applying separately a a two sample proportions test (via the function `prop.test`) to each ID. Is there a more rigorous approach to do this? I'm thinking of doing a chisquare test, followed by a post-hoc analysis (testing for significant differences among all pairs of ID). However, I'm not sure whether this is the proper approach and I don't know how to do this in R... Can anyone help?


title: Importance Sampling MC - a couple of questions regarding PDF

text: I implore the good people to quickly glance through this thread over at StackOverflow to get a better idea of my question if the following isn't clear. I have these integrals to evaluate using Importance Sampling-- 1) `x^(-0.5)` ; `x in [0.01,1]` 2) `[1+sinh(2x)ln(x)]^-1` ; `x in [0.8,3]` **Question 1** How does one chose the importance PDF really? For the above functions I just plotted exponential and normal functions (since it just occurred to me that they are similar to the integrals to be evaluated. Otherwise, is there a procedure to arrive at _which_ PDF is to be used? **Question 2** Is it necessary that the chosen importance PDF has to have lesser variance than the one to be evaluated? My understanding is that the PDF chosen should follow and cover the plot of the given integral. **Question 3** So, I have a function `x^2 + y^2`, where `x,y in[-1,1]` . For cases like this, where the integral to be evaluated is a double integral, what PDF should i think about? 2-D normal? If it's a random function, there's no way to generate samples from it (for the lack of built in commands in R- only standard PDFs exist.)


title: Prediction with categorical and numeric variables

text: I used R to estimate a regression with both numeric and categorical variables, and obtained coefficient estimates. However, when I try to make predictions using new data, there appear to be some problems with dimensions. Is there anything that must be done to the code so that R handles both types of variables? Is there a command other than `predict(...)` to combine my coefficients with the new data?


title: Sampling from low probability regions

text: I have a pdf that I need to sample from, but I'm having numerical issues and was hoping someone could provide some insight. My pdf is a product of gaussians $$ P(\theta |d ) = \prod_{i=1}^{N} C_{i} e^{-\frac{1}{2 \sigma_{i}^{2}}(d_{i} - y_{i}(\theta) )^{2} } \,. $$ Lets say for instance that $\theta$ is only one dimensional, making $P(\theta | d)$ a one dimensional pdf. I wish to plot P(\theta | d) as a function of $\theta$, where $\theta$ is on the x-axis of the plot and $P(\theta | d)$ is on the y-axis. In this case I can simply plug in a bunch of various values for $\theta$ and plot the result. The problem is, the value inside the exponent which I'll call $K$, $$ K = - \sum_{i = 1}^{N} \frac{1}{2 \sigma^{2}_{i}}(d_{i} - y_{i}(\theta) )^{2} $$ can vary very much, and never really approaches zero. For example, the range of values for K I got varies from about -20000 to about -500. This is because my data, $d_{i}$, is pretty fairly well corrupted by noise. Is there anything I can do to plot this pdf as a function of $\theta$? Obviously everything I do on my computer (taking $e^{K}$ and getting mad) is pretty much giving me all zeros.


title: probability question

text: I would like to calculate the following conditional probability: I know that the probability of a BLUE ball being drawn is 0.3. I receive a message from A or B who saw the ball that has been drawn. This person tells me that the ball drawn is BLUE. I know the probabilities: A observing BLUE and sending BLUE- 0.16 B observing RED and sending BLUE - 0.09 A observing BLUE and sending BLUE- 0.05 B observing RED and sending BLUE- 0.17 I also know the probabilities: A observing BLUE and sending RED - 0.019 B observing RED and sending RED- 0.323 A observing BLUE and sending RED-0.073 B observing RED and sending RED=0.106 I would like to calculate the conditional probability of the ball drawn actually being BLUE given that I received the message BLUE and the information I know about each player's probability of sending a given message.


title: Question on scaling data

text: I have been reading Elements of Statistical Learning. And I saw some example R codes written online. I realise that very often the explanatory variables $X$ (the data of the predictors) are standardised to ' _zero mean and unit variance_ ', while the data of the response variable $y$ is scaled to have zero mean only. My question is that why don't we also standardised $y$ (response variable) to have not only zero mean, **but also _unit variance (of response variable)_** ?


title: How useful are per capita calculations in level of education of a population?

text: I had a discussion recently in which it was stated that Scotland having more universities in the top 200 per capita than any other country is a good indication of the intelligence of the population of Scotland: **NOTE: I am not questioning the intelligence of Scots, just this calculation's validity** I don't see any way to compare the measurement of a population's intelligence by this method. For example: if Dundee, Scotland's lowest ranking (in the 200) university, fell just 4 places to 201, the entire country's intelligence measured in this manner would drop drastically. Then if it rose one place it would rise again by the same level as the fall of 4 places. If the five universities in the top 200 were at places 1,2,3, 4 and 5 there would be no change. As well as that, funding can be external as in the UK's case the Scottish universities receive **some** EU and UK Research Council funding which, even if Scotland is at a tax/public spending surplus could still cause fluctuations which do not affect the intelligence level of Scots but will most likely affect their university positions in the top 200. Many students, staff and teachers are also not Scottish so really the measure of how intelligent a population is (using this method) comes down to those who make the decisions in the universities. Am I missing something here, is this a useful calculation at all and does it relate to a country's intelligence? Are there any useful indicators of the population of a country's intelligence?


title: How to combine 2 different observations to improve state estimates?

text: ### Context Let $\mathbf{x}_i \in \mathbb{R}^{100}$ and $\mathbf{z}_i \in \mathbb{R}^{20}$ be input vectors with the same corresponding target $\mathbf{y}_i \in \mathbb{R}^{25}$. Using ridge regression we obtain two separate mappings, $\mathbf{A}_\mathbf{x}$ and $\mathbf{A}_\mathbf{z}$, such that $$\mathbf{A}_\mathbf{x}\mathbf{x}_i = \mathbf{y}_i + \mathbf{\epsilon}_\mathbf{x},$$ $$\mathbf{A}_\mathbf{z}\mathbf{z}_i = \mathbf{y}_i + \mathbf{\epsilon}_\mathbf{z},$$ with $\mathbf{\epsilon}$ denoting noise. I wish to minimise the noise (error) by considering both observations $\mathbf{x}_i$ and $\mathbf{z}_i$. ### Questions How can I combine the observations ($\mathbf{x}$ and $\mathbf{z}$) to obtain a lower error on average? I have two approaches, neither of which give satisfactory results on average:   1. Append the two vectors [$\mathbf{x}$ ; $\mathbf{z}$] and learn a new mapping.   2. Use a weighted average of the two estimates. For (1), do I need to do any preprocessing before appending the vectors, such as spatial alignment (PCA?) or some form of value normalisation? For (2), is there a way to automatically determine a suitable weighting? What other ways are likely to yield better results? How can I analyse the complementary nature of the two observations in relation to the state being estimated?


title: How to Interpret GAM P-Values?

text: My name is Hugh, and I'm a PhD student using generalised additive models to do some exploratory analysis. I'm not sure how to interpret the p-values that come from the MGCV package and wanted to check my understanding (I'm using version 1.7-29, and have consulted some of Simon Wood's documentation). I looked for other CV questions first, but the most relevant ones seem to be about general regressions, not GAM p-values in particular. I know there are lots of different arguments to GAM, and the p-values are only approximate. But I'm just starting simple to see if there is any "signal" whatsoever for my covariates. E.g.: **Y ~ s(a, k = 3) + s(b, k = 3) + s(c, k = 3) + s(d, k = 3) + s(e, k = 3)** Approximate p-values of smooth terms: s(a) = 0.000473   s(b) = 1.13e-05   s(c) = 0.000736   **s(d) = 0.887579   s(e) = 0.234017** _R_ ² (adjusted) = 0.62$\quad$ **Deviance explained = 63.7%**   GCV score = 411.17$\quad$ Scale est. = 390.1$\quad$ _n_ = 120 I cut the df columns, etc., due to formatting. I'm interpreting the p-values for each covariate as a test of whether the corresponding smooth function significantly reduces model deviance, where p is the probability of obtaining data at least as 'relatively implausible' as that observed under a null model of 0. This would mean that (e.g. with alpha = 0.05) the smoothed functions did not reduce the deviance for "d" & "e" vs. a null model, whereas they did for the other terms. Hence (d) and (e) do not add significant information to the regression, and the deviance explained is down to (a) (b) (c)? Any advice would be greatly appreciated, and best of luck with your research.


title: How to run a prediction model

text: I have a dataset from a cohort of 500 children under 15 years of age that were screened for measles antibody titre. Each row in the dataset represents an individual child. Covariates included in the dataset for each individual are: age, sex, whether a child ever suffered from measles before (a binary outcome, yes/no), whether a child ever received measles vaccination before (a binary outcome, yes/no), observed apparent sign of malnutrition (a binary outcome, yes/no), and family's socio-economic status (a binary outcome, low/middle). The level of antibody titre per individual is our dependent variable and was categorised as 10, 20, 40, 80, 160, and 320. As antibody level greater than or equals to 40 is considered "protective" and anything less as "unprotective", another dichotomous dependent variable (protected/unprotected) for level of protection was generated therefrom, thereby having 2 dependent variables. The antibody titre is an indication of the level of protection due to either vaccination or natural infection. My questions are   1. Should I employ ordinal and logistic regressions in modelling the categorical antibody titre level (10, 20, 40, 80, 160, and 320) and the dichotomous level of protection (protected/unprotected), respectively, and thereby produce 2 models from a data with 2 outcomes?   2. As one of my hypotheses is that "measles may likely increase to 'X' fold in 'Y' year period under this present condition", how could I run a model on this kind of data to make such predictions?


title: Notion of orthogonality in clustering

text: I'm reading in a paper " _the goal is to cluster data in a way that is meaningful and at the same time as orthogonal as possible to the given classification_ " Could someone have any idea how the notion of orthogonality is used in this case?


title: Data APIs/feeds available as packages in R

text: EDIT: The Web Technologies and Services CRAN task view contains a much more comprehensive list of data sources and APIs available in R. You can submit a pull request on github if you wish to add a package to the task view. * * * I'm making a list of the various data feeds that are already hooked into R or that are easy to setup. Here's my initial list of packages, and I was wondering what else I'm missing. I'm trying to limit this list to "real time" or "close to real time" data feeds/APIs, where the underlying data might change between downloads. There's plenty of lists out there for static datasets, which only require one download. This list is currently biased towards financial/time series data, and I could use some help expanding it into other domains. Free Data:   Data Source - Package   Google Finance historical data - quantmod   Google Finance balance sheets - quantmod   Yahoo Finance historical data - quantmod   Yahoo Finance historical data - tseries   Yahoo Finance current options chain - quantmod   Yahoo Finance historical analyst estimates - fImport   Yahoo Finance current key stats - fImport - seems to be broken   OANDA historic exchange rates/metal prices - quantmod   FRED historic macroeconomic indicators - quantmod   World Bank historic macroeconomic indicators - WDI   Google Trends historic search volume data - RGoogleTrends   Google Docs - RGoogleDocs   Google Storage - RGoogleStorage   Twitter - twitteR   Zillow - Zillow   New York Times - RNYTimes   US Census 2000 - UScensus2000   infochimps - infochimps   datamarket - rdatamarket - requires free account   Factual.com - factualR   Geocode addresses - RDSTK   Map coordinates to political boundaries - RDSTK   Weather Underground - Roll your own   Google News - Roll your own   Earth Sciences netCDF Data - Roll your own   Climate Data - Roll your own   Public health data - Roll your own   OAI Harvester - Open Archives Initiative harvester   RAmazonS3 - S3 Amazon storage server   Rflikr - Flikr api Requires a subscription:   Bloomberg - RBloomberg   LIM - LIM   Trades and Quotes from NYSE - RTAQ   Interactive Brokers - IBrokers


title: Formula for weighted simple linear regression

text: This wiki page Simple linear regression has formulas to calculate $\alpha$ and $\beta$. Could anyone tell me how to derive the formulas in weighted case?


title: Examples of derived data quality features?

text: I have a predictive classification task and the end users want to get **initial** understanding of:   * the quality (not specified further) of each data point in a data set   * whether there are systematic patterns between data quality problems (e.g. rows with missings values have also high number of outliers etc.) What kind of derived data quality features can be used/ experimented with? No metadata (who did the measurements, when etc.) or other information external to the dataset is available. So far I have computed for **each row**   * number of features with missing value/ total number of features   * number of features with univariate outlier (+1.5 IQR)/ total number of features   * multivariate outlier score (LOF algorithm)   * minimum distance to another data point (close to being a duplicate)  I would appreciate any ideas I can experiment with or references that deal with data point quality computations. As 'quality' is not well-defined here, any ideas can help.


title: How to get the expected counts when computing a chi-squared test?

text: I just have 1500 data points in excel. However, I want to do a chi squared test to prove if the data is normaly distributed. My question is: How to calculate the expected range for such a large data set? **UPDATE** How to calculate the relative frequencies for such a large data set?


title: Possible to evaluate GLM in Python/scikit-learn using the Poisson, Gamma, or Tweedie distributions as the family for the error distribution?

text: Trying to learn some Python and SKLearn, but for my work I need to run regressions that use error distributions from the Poisson, Gamma, and especially Tweedie families. I don't see anything in the documentation about them, but they are in several parts of the R distribution, so I was wondering if anyone has seen implementations anywhere for Python. It would be extra cool if you can point me towards SGD implementations of the Tweedie distribution! Thanks.


title: What can't be expressed as a linear model?

text: Say I have outcome variable $Y_i$ and predictors $X_{i1}$ and $X_{i2}$ for some data point $i$. Wikipedia says that a model is linear when: > the mean of the response variable is a linear combination of the parameters > (regression coefficients) and the predictor variables. I _thought_ this meant that a model can be no more complicated than: $Y_i = \beta_1 X_{i1} + \beta_2 X_{i2}$. However, upon further reading, I found out you could handle non-linear "interactions" of the predictors like in $Y_i = \beta_1 X_{i1} + \beta_2 X_{i2} + \beta_3 X_{i1}\ X_{i2}$ by viewing $X_{i1}\ X_{i2}$ as just another predictor (which happens to be dependent on $X_{i1}$ and $X_{i2}$). This seems to mean that you can use any (linear or non-linear) function of the predictors like $\log(X_{i1} / X_{i2}^2)$ or whatever. Conceptually, this type of "recoding" seems like it should work for the coefficients as well. So: What exactly are the limits of linear regression, given you can do this kind of manipulation?


title: Why must one use ANOVA to interpret the impact of interaction effects and factor variables?

text: What is the purpose of the ANOVA table? I once learned that you can only interpret the significance (p-value) of a multi-level discrete variable, or an interaction effect using the ANOVA table. Why? Why can't you use the p-value outputs of the regression? Why do people look at the ANOVA table in practice? GLM               Coefficients:                                         Estimate Std. Error z value Pr(>|z|)     (Intercept)                          -1.9800     1.3697  -1.446    0.148     ConnectivityHIGH                      1.9214     1.6361   1.174    0.240     SusceptibilityHIGH                    0.8636     1.7183   0.503    0.615     ConnectivityHIGH:SusceptibilityHIGH  -0.6555     2.1348  -0.307    0.759      ANOVA                                           Df Deviance Resid. Df Resid. Dev Pr(>Chi)     NULL                                           19     3.6712              Connectivity                 1  2.43379        18     1.2374   0.1187     Susceptibility               1  0.19710        17     1.0403   0.6571     Connectivity:Susceptibility  1  0.09598        16     0.9443   0.7567


title: Skewness in robust regression

text: I am doing a robust regression to predicted a left-skewed outcome. Do I still need to do a transform on my data before doing robust regression, or will that issue be taken care of by the robust method? If so, what transform would you recommend? When I plot my residuals after doing robust regression, they still look left-skewed. I tried doing a transform ($\text{predicted.value}^4$) but even then my residuals are still left-skewed. I'm worried I'm not meeting certain assumptions of robust regression, but I'm not sure since I'm new to this method. Any thoughts? Thanks!


title: Efficent global minimum search of costly to evaluate low dimensional error function

text: Trying to minimize the error from a probabilistic regression model which is composed of hierarchical KDE estimated PDFs. The top level is the result of 2-4 separate meta-PDFs from different data fields the combination of many of their own KDE estimated PDFs. I'm trying to train weights affecting how much each top-level PDF contributes to the final distribution but evaluating the error is costly (a minute over a representative training set of approx 3% of the data) and a brute force sweep over 4 weights at sufficient resolution would be on the order of a week to a month. I've heard about a bunch of different techniques but am at a loss as to which one should be relevant or effective to minimize error here. Anyone know which techniques would be effective in this situation?


title: Binomial GLM and different sample sizes

text: I have a data set which consists of binomial proportions, let's say the success rate of converting a customer depending on the advertisement, the customer age, and various other factors. For some common combinations of covariates, I have a lot of data, and therefore the binomial proportion of successes has low variance. For rare combinations of covariates, however, I have very little data, and therefore the variance of the proportion is high. The magnitude of differences is very large, for example I might have 1 million trials for some combinations of covariates, and only 50 for others. However, I want to include ALL data in my model and weight it appropriately to get the best model fit. I've tried to use R to do binomial (logistic) regression using a generalized linear model:               lrfit <- glm ( cbind(converted,not_converted) ~ advertisement + age, family = binomial)      This is a good start because it automatically weights the observations by the number of trials. However, it's not good enough. Here's why: Let's say you have some observations with 100,000 trials and others with 1,000,000 trials. If you weight by number of trials the latter group is going to receive 10 times the weight. This seems nonsensical, however, because both observations are easily precise enough to receive equal treatment in the model. Clearly you want to penalize groups with only 10 or 100 trials, but as the number of trials gets larger, the weight should stop increasing. Since in weighted least squares the reciprocal of the error variance is used as the weight, my idea would be to use calculate the posterior variance of the proportion (using Jeffrey's prior), then add some constant term to it (this will make sure the variance stops increasing at a certain number of trials) and then use the reciprocal of the sum as the weight. Is this approach reasonable? Am I missing something? Can someone give me more information about this method?


title: Expectation values of functions

text: This question is more to do with interpretation than calculation. I have a model which predicts the probability of a detector 'firing' under a certain intensity of signal, or actually in this case lack of a signal. The detector is made of a number of sub-pixels so the probability of none of the sub-pixels firing in a certain window is $$ P(N=0) = \exp^{-\lambda(x, \gamma)} $$ The detector is made of a number of elements so the probability of none of them firing is given by $$ P(tot = 0) = \Pi P_{(pix)}(N=0) $$ There is some uncertainty on gamma which I assume to be distributed according to some PDF so the expectation value for the total probability is now $$ \int (\Pi P_{(pix)}(\gamma)) P(\gamma)d\gamma $$ My problem is this - I might, and often do, want to look at individual pixels. So I might want to move the integral inside the product. Or alternatively I might wish to calculate the expected mean value in the first equation. $$ \bar \lambda = \int \lambda(x,\gamma)P(\gamma)d\gamma $$ I know what happens numerically if I compare F[E(X)] and E[F(X)] but I struggling with the conceptual difference between the cases. What is the difference between the product of the expectation values for the individual pixels and the expectation value of the product. Thanks


title: How to find patterns and identify changes in them in time series with R

text: This is my first question on stats, just trying to learn the basics of time series analysis with `R`. So any good suggestions about learning resources will be highly appreciated as well as the answer to the question. For the data below, let's say it represents number of website visits per day, I would like to find out:   1. What the weekly pattern is (e.g. the highest number of visits occurs on Thursdays, the lowest on Fridays etc.)   2. Automatically detect changes in that pattern (e.g. in 2008 most of the visits occur on Thursdays, but the in 2009-01-04 the pattern changes to something else) Please let me know if I can provide more details.               > str(daily)     An ‘xts’ object on 2007-02-19 23:32:16/2013-05-05 15:09:17 containing:       Data: num [1:2268, 1] 55 32 70 48 75 50 48 46 36 55 ...      - attr(*, "dimnames")=List of 2       ..$ : NULL           ..$ : chr "cnt"       Indexed by objects of class: [POSIXct,POSIXt] TZ:        xts Attributes:        NULL          > head(daily)                         cnt     2007-02-19 23:32:16  55     2007-02-20 23:58:58  32     2007-02-21 23:40:41  70     2007-02-22 23:01:41  48     2007-02-23 23:53:06  75     2007-02-24 23:47:07  50          plot(daily)      ![PlotDaily](http://i.stack.imgur.com/OzUw5.png) Full dataset: https://dl.dropboxusercontent.com/u/65347419/daily.csv


title: Correlation between two nodes of a single layer MLP for joint-Gaussian input

text: Let's say you have a jointly gaussian vector random variable $\mathbf{x}$, with mean $\mathbf{M}$ and covariance $\mathbf{S}$. I now transform each scalar element of $\mathbf{x}$ , say $x_j$, with a sigmoid: $$y_j = 1/(1+\exp(-x_j))$$ I am interested in the expectation between two variables $y_j$ and $y_j'$ of the resulting distribution, that is $$E\{y_j y_j'\}$$ where $E\{\}$ is the expectation operator. Note:   * The whole PDF of $y$ can be computed just by applying a change of variable. Unfortunately the integrals leading to the expectations are intractable (to the best of my knowledge).   * I'm looking for closed-form approximations, no Markov Chain-Monte Carlo, no variational stuff. That is, approximations to the variable change, to the expectation integrals, to the sigmoid, to the resulting PDF, that allow computing $E\{y_j y_j'\}$.    * Some dead ends: Taylor, often used in papers on the topic, is inaccurate by the mile. Gradshteyn and Ryzhik does not seem to contain the integrals.


title: What is the difference between lmer(.. (1|groupingfactor1) +(1|1|groupingfactor2)) and lmer(..(groupingfactor1|groupingfactor2)?

text: I am quite new to `lmer` and quite confused about how to select the random effects part of the model. The data I am working on is divided into dependent and independent variables. Each variable has two grouping factors. 1). sampling time (sample 1:7) and sampling column (1:24, where 1:3, 4:6 and so are replicates). I have explored these different possibilities, but I am not sure on which is the correct one to use.               fm1 <- lmer(dependent ~ independent +(1 + independent|Sample)+(1+ independent|Column), data=Data)          fm2 <- lmer(dependent ~ independent +(1|Sample)+(1|Column), data=Data)          fm3 <- lmer(dependent ~ independent + (Sample|Column), data=Data)      `fm3` is the model with the lowest AIC and BIC values, but I am not sure that this is the right one to use?


title: How scientifically handle outliers fairly?

text: I conducted a survey, in which participants ranked themselves depending on how often they did something and how good they were at doing that thing. The input were integer values on the survey and _deliberately_ did not have a maximum (or minimum) range. Now I'm running some analysis and calculations so that each participant is given a calculated score. This score takes into account a variety of factors and I'm quite happy with that calculation. The graph below plots the participant scores on the left hand side/Y axis, with each data point representing a participant. The results shown in this graph have been unscientifically butchered to give this nice distribution... (more on this below) ![Distribution of scores](http://i.stack.imgur.com/f8zaM.png) My problem is, that a couple of participants scores were ludicrously high, where most people were between 0-160, some scores were coming out at 2500, 3000. The calculation that produced these scores are just fine, it's just based off survey data input, and some individuals entered unrealistic values, hence a unrealistic score that is an outlier (ie, outside the normal looking distribution of between 0-160). I need to eventually normalize all scores down to a scale between 0 - 10 for a report, but if I do that now the participant that entered 2500 gets a score of 10, and everybody else will get a score of 0 or 1! I don't know how to scientifically deal with these odd results of 2500, 3000, etc. I can manually cut 2500 down to 160, which would still be 10 in the final analysis, but would mean that everybody else gets fair final scores too. However, if I manually cut these scores down on a case by case basis it seems really unscientific and unfair. I could leave those results out of the report, but again, that seems rather unscientific. So, how do I handle these outliers, so that I can re-plot all scores from 0-160 (a reasonable distribution)? * * * Note: I expect that there are really quite well defined solutions to this problem, but I'm afraid my Math experience (and education) is negligible and I'm worried about Googling a solution and implementing something that is quite wrong. Therefore, I've taken the time to state my problem/question in full at risk of asking a duplicate question, just so I've got a higher degree of confidence that you're advising me on the most appropriate way to normalize these results. Additionally, because my Math is quite bad, can I just say that I'm very keen to learn, but would really appreciate it if you can build up your answer in steps as if you were explaining it to a toddler!


title: Which statistical tests can be used to compare two paired/related rankings?

text: To understand my question, can anyone to evaluate if the statistical analysis in the below research (from the 2011 year) was performed correctly? I must to analyze a similar research and I am interested in the correct method. Initially I supposed that in this kind of research must be used the Wilcoxon test, but now I am inclined to the Spearman correlation test. Is this correct? Thank you! Patient-Perceived Changes in the System of Values After Cancer Diagnosis.


title: How to measure the accuracy of multinomial target?

text: I have trained a model that predicts a multinomial output. Would anyone know how the accuracy can be measured. The target takes on the following values: "Yellow","Red,"Green","Blue". I know that for binomial targets, ROC/AUC provide a good solution. But I haven't found anything for my problem. Would appreciate any help. Thanks


title: How to determine trend strength from linear regression slope?

text: I am new in statistics so sorry about my elementary question. I have a data set and simple linear regression equation calculated from this data:               f(x) = ax + b      I would like to know if I can (and if so how) use slope (`a`) from this linear regression to determine strength of data set trend? I need to know is trend is growing or falling and how strong (steep) is. I am not sure how to calculate it because there is not max value of tangent value so I'm not able how to find out if it's a strong or weak trend. I hope I described my problem correctly. **Edit:** I am working on my school project and I would like to accomplish this: there is a dataset, in my situation these datasets contain a business data (monthly sales, number of orders, an average order size etc.). I would like to have a simple indicator for each of these quantatities that is able to visually describe a trend (a strength) of this dataset in this manner:   * Growth (green color) - if growth is strong, e.g. between 10 - 100% of maximum growth   * Low Growth (yellow color) - if growth is weak, e.g. between 0-9% of maximum growth   * No move (gray color) - no move (when regression has no slope or very small)   * Low drop (yellow color) - drop is between 0-9% of maximum drop   * Drop (red color) - drop is strong, 10-100% of maximum drop P.S. In my app I have another visual indicator but it describes a dataset in different way - in particular time, how it works - I compute a median value of whole dataset and in a particular point I compare a value in my dataset with median value and determine a color that describes a status in this particular value. Some simple solution I need for the **whole** dataset. **Additional info** : to describe what I am trying to solve I created these two plots that plot subset of my dataset. In header there is a slope value of this data subset. Both slopes are negative, so trend in this data is decreasing. This information I can determine myself but both trends are not equal. The first trend is stronger than the second (slope value is higher too). I need to know if is there some scale or dimension according to I can normalize slope value and say in my application, e.g. dataset that describes this slope value is not so much important, because slope is in max 5% of its maximum value but this slope is in 75% of its maximum negative value, so be careful and take care of this trend! All what I need to know is: **what is minimum and maximum value of a slope parameter in linear regression equation?**. If I know the scale of the slope value I can compute everything I am asking. I can preprocess data and normalize them or do whatever else what is needed. ![enter image description here](http://i.stack.imgur.com/llG6t.png)


title: Cumulative function to probability density function

text: I have some data points for which I would like to show the density using Python, but the package SciPy is not working for Python 3.3. I am looking for a solution to be able to plot some data in the shape of a normal density. I do not know how to get the nice drawing. When I sort the data points, it would be a cumulative function right? So, given that there is no `density()` function available, what is the best way to plot something that will make a nice "normal-ish" graph. In R I would simply use:               plot(density(data))      EDIT: So I am looking for a statistical solution rather than a software solution. What is the easiest way to make a density (plot) without having a density function available.


title: Prediction with intervals as the independent variable

text: I have sample data that maps intervals to a number:               [3,7] => 1     [6,8] => 2     [6,13] => 3     [7,10] => 3     [10,13] => 4      The dependent variable's values will generally increase as the intervals' values increase in a nonlinear fashion (and they're not part of a finite set, so classification won't work here). What method(s) could be used to predict the mapped value for any given new interval?


title: ezAnova vs. lme for factorial repeated-measures design: results differ, why?

text: I have used `lme` and `ezAnova` to analyse data from a 2$\times$3 repeated- measures experiment. Theoretically those are two different ways to perform the same analysis. However, the resulting $F$-statisics and DF differ and I am lost in why. Here is the exact data and output: I have a data set with 2 independent variables (`marker_lang` and `congruency`) and the dependent variable `RT`: Both IV are repeated and completely crossed (thus 6 conditions overall). The data are not collapsed to cell means, meaning that per condition and subject I have several data points. Here is what I did with ezAnova:               ezANOVA(subset(data.mark.afc, !is.na(afc.RT)), dv=afc.RT, wid=subjectID,      within=.(marker_lang,congruency), within_full=.(marker_lang,congruency), detailed=1, type=3)      And the output:               $Anova                   Effect DFn DFd          SSn        SSd            F            p   p<.05          ges     1            (Intercept)   1  24 84879098.819 1892110.06 1076.6278430 1.881814e-21     * 0.9762134497     2            marker_lang   1  24    36392.804   80595.30   10.8371986 3.071336e-03     * 0.0172922873     3             congruency   2  48    25426.393   47319.45   12.8960382 3.292730e-05     * 0.0121448066     4 marker_lang:congruency   2  48     1160.152   48150.91    0.5782581 5.647333e-01       0.0005606399      Here is what I did with lme:               basemodel <- lme(data=subset(cdata, !is.na(afc.RT)), afc.RT~1,     random=~1|subjectID/congruency/marker_lang, method="ML")          langmodel <- update(basemodel, .~. + marker_lang)          angcongmodel <- update(langmodel, .~. + congruency)          fulmodel <- update(langcongmodel, .~. +marker_lang:congruency)      And the anova-tables for the lme analysis:               anova(fulmodel)                        numDF denDF   F-value p-value     (Intercept)                1  4715 1112.3468  <.0001     marker_lang                1    72   24.8917  <.0001     congruency                 2    48    8.3902  0.0008     marker_lang:congruency     2    72    0.4475  0.6410          anova(basemodel, langmodel, langcongmodel, fulmodel)               Model df      AIC      BIC    logLik   Test   L.Ratio p-value     basemodel         1  5 64203.43 64235.88 -32096.72                              langmodel         2  6 64185.06 64224.00 -32086.53 1 vs 2 20.366313  <.0001     langcongmodel     3  8 64173.27 64225.19 -32078.64 2 vs 3 15.790082  0.0004     fulmodel          4 10 64176.38 64241.28 -32078.19 3 vs 4  0.892535  0.6400      I would expect the $F$, DF, and $p$-values for corresponding effects to be the same, which is not the case. This seems not an issue of different anova-types, as I tried out different types for `ezAnova`. None yield the same result as anova of `fulmodel`. Any help/ideas will be greatly appreciated!


title: What do the "coefficients" in R's HoltWinters function represent?

text: I'm using the HoltWinters function in R and I'm trying to understand what the "coefficients" represent in the object that is returned by that function. They don't seem to match in any obvious way the values returned when you look at the $


title: Using R for dice probabilities

text: Using R, what function(s) would I use to obtain the following probabilities?   * Roll at least one 1 when rolling 2 six-sided dice (2d6) = 11/36   * Roll at least one 1 when rolling 3 six-sided dice (3d6) = 91/216   * Roll at least one 1 when rolling 1d4, 1d6, 1d8, and 1d8 = 801/1536 First I hope my answers above are correct! I did these pretty much manually. I think I need to use binomial distributions and/or probability-generating functions, but not sure if I'm over-complicating things. I've tried using R's *binom() functions but can't seem to arrive at the answers I need.


title: Test to calculate the probability of winning

text: I'm running a predictive model to predict the probability of winning a certain item based on the price that I bid (other factors also). After running the model (ols) in R, I wanted to account for all the variables in my model and develop a graph highlighting the 'predicted probabilities' regarding the primary variables I'm concerned about. So want to have a line graph showing the probability of winning on the y axis, and the bid on the x axis. The following data would result in a graph which shows that the probability of winning decreases as the bid increases.               Bid                  8  6      4     Probability Winning 30% 22%    18%        1. Are predicted probabilities only relevant for logistic regression models or can be equally relevant for linear regression models?   2. What is the reasoning and logic behind going from a model to a probability curve which would show the 'trend' in one variable as predicted by another, while accounting for all other factors. Sorry for the elementary question, I'm just a little clueless. Thanks for the help!


title: Reliability measure with no single item scores

text: Let's say, I have an intelligence test made of 60 questions (items) which produces results for two intelligence constructs on a 0-100 scale (for instance crystallised and fluid intelligence). I only have respondent's final scores on the two constructs. So, I don't have answers to single items. In such a scenario, I cannot calculate Cronbach's Alpha as a measure of internal consistency, because I don't have single items' scores. Is there any reliable alternative measure of reliability that I could use, only having final test scores?


title: Providing variance measures for speedup ratios

text: I'm doing timings of an application which speeds up as we add threads. However, the times vary on repeated runs. What I'd like to have is a speedup chart which shows how many times faster than the baseline the application is, for each thread count I test on. The data I was planning on gathering was 5-10 runs for each number of threads, from which I can get standard deviation for the times at each thread count. But I'm at a loss as to how to represent this deviation in terms of 1 std.dev error bars on a speedup-curve. I'm thinking deviation I want to show is the std. of the speedup calcs, not the times intially right? To make this more concrete:   * Suppose as a baseline, we have the times 10s,10s,10.1s,9.9s (sd = 0.081);    * at 4 threads, we have 2.5,2.5,8,2.8 (sd = 2.72);    * Speedups are 4,4,1.2625,3.5357 (sd = 1.3) So the if I'm plotting the point for 4 threads, does it make sense to have the 1 std. error bar be 1.3 (2.6 in both dir)? ( ignoring skew of course. )